Qualitative dynamics of wavepackets in turbulent jets
Abstract
It has long been established that turbulent jets comprise largescale coherent structures, now more commonly referred to as “wavepackets” ((1)). These structures exhibit a remarkable spatiotemporal organisation, despite turbulence.
In this work we analyse, from a qualitative point of view, the temporal dynamics of axisymmetric wavepackets educed, experimentally, from subsonic isothermal jets. We use the data presented by (2), where timeseries of the wavepackets are extracted at different streamwise locations. A thorough analysis is performed; statistical tools are used for estimating the embedding and correlation dimensions characterising the dynamical system. System identification is used for computing nonlinear surrogate models. Finally, controloriented linear models are computed.
The goal of the contribution is to assess the extent to which nonlinear models are necessary, or appropriate, for description of the temporal wavepacket dynamics and to provide a complementary perspective to the current modelling.
keywords:
1 Introduction
Following the early studies of MollöChristensen ((3); (4); (5)), a considerable body of work has been devoted to exploring the nature of organised motions that are observed in turbulent jets ((6); (7); (8); (9); (10)). It was hypothesised early on that this component of the flow might be understood in terms of an instability of the turbulent mean ((6); (11); (12); (13); (14); (15)); and the importance of such flow structures for sound radiation has been suggested in numerous studies ((16); (17); (18); (19); (20); (8); (21)).
It is only more recently, however, thanks largely to progress in theory, numerical simulation and experimental diagnostics, that it has been possible to explore these hypotheses in a comprehensive manner. Exhaustive comparison of the results of theory with experimental measurements has confirmed that the average characteristics of coherent structures in turbulent jets are remarkably well described by solutions of the NavierStokes equations linearised about the turbulent mean ((22); (23); (24)); these solutions are synonymous with globally stable modal solutions ((25); (26); (27)), which, physically, amount to hydrodynamic waves that are convectively amplified in the upstream region of the flow, but become neutrally stable and then decay farther downstream. It is for this reason that they have been given the denomination ‘wavepacket’ ((1)). An example of a modal solution, from (27), is shown in figure 1, where it is compared with a realisation taken directly from the Large Eddy Simulation ((28)) that provided the mean flow. The wavepackets can be seen to comprise a hydrodynamic component, with an amplitude envelope as described above, and an acoustic component that takes the form of a directive beam radiating to shallow angles.
Considering the Reynolds number of the jet, , and the fact that it issues from a nozzle with fully turbulent boundary layers, the agreement between the LES realisation and the linear global mode is striking. One is compelled to ask how, despite the nonlinear, orderless character of the turbulence that dominates the fluctuation energy of this flow, such organisation can exist and how a linear model can capture so many of its features.
However, something less clear from the images is the considerably lower acoustic efficiency of the linear wavepacket. As shown in a number of recent studies ((29); (24); (30)), the intensity of the soundfield of the linear solution, once the hydrodynamic fluctuation levels of LES and model have been matched, can be as much as 40dB below that of the LES. Such a large discrepancy, in the face of such compelling organisational similarity, is intriguing and the subject of a number of ongoing research efforts, which we summarise briefly in what follows.
The present understanding is that the linear dynamics do not contain the acoustically important degrees of freedom: the average wavepacket does not generate the average sound ((31); (32); (33); (24); (29)). The essential soundproducing motions are those associated with higherorder statistics of the wavepacket dynamics. These motions have been denoted ‘jitter’ by ((31)), and would appear to be underpinned by nonlinearity ((24); (29); (34); (35)).
The nature of this nonlinearity remains to be clarified, and it motivates the present study. Do wavepackets jitter on account of nonlinear wavewave interactions ((36); (37)), i.e. is the nonlinearity of an intrinsic kind? Or is it rather an extrinsic, stochastic forcing of linear waves by background turbulence that leads to the activation of these higherorder degrees of freedom ((34); (35))? Or might both extrinsic and intrinsic mechanisms be at work? In this work, in light of the foregoing considerations and questions, we are interested in lowdimensional projections of the complete wavepacket dynamics.
A further issue, that motivates both the modelling efforts evoked above and the work we undertake in this paper, is that of control. Given the success of linear models in predicting both the average wavepacket structure, in frequency space, and the realtime evolution in an experimental context (38), is it legitimate to ask if the tools provided by linear control theory might be sufficient? Such possibilities are presently being explored by (39). But if it were to be necessary to directly manipulate the more subtle dynamics associated with the energeticallyunimportantbutacousticallyessential degreesoffreedom, then a nonlinear control framework would be required (machine learning etc. (40)). In both cases the question of the dimension of the space spanned by the dynamics must again be asked: is this small enough for control to be realistic in an experimental context?
And, finally, there is the issue of the physical interpretation of POD modes, which some recent results show may correspond to the highest gain structures of the linear NavierStokes operator subject to stochastic external forcing by turbulence. With this in mind, the various analyses of the dynamics will be performed using both raw and PODfiltered data.
Consideration of the above issues via qualitative dynamical analysis has, to the best of our knowledge, not previously been attempted. We therefore consider the same turbulent jets studied by ((41); (23); (2)), using a variety of nonlinear signalprocessing and system identification techniques, in order to explore dynamics of lowfrequency wavepackets. Particular attention is given to pressure signals dominated by fluctuations at Strouhal numbers , where linear models are found to perform most poorly.
Our analysis starts with a minimal description of the experimental setup in Sec. 2; we refer to (2) for a detailed description. The remainder of the paper is organised in two parts. The behaviour of the timeseries obtained from the experimental campaign is discussed in Sec. 3, where statistical tools and spectral analysis are introduced for a preliminary characterisation of the temporal dynamics.
In the second part, we pursue an alternative approach by analysing reducedorder models of the timeseries, based on system identification; the approach is introduced in Sec. 4, while the dynamical analysis is reported in Sec. 5 and Sec. 6 . More details on the techniques are given in A. The manuscript finalizes with a discussion of the results and concluding remarks in Sec. 7.
2 Experimental data
The study presented in this paper is based on the analysis of the temporal behaviour of axisymmetric wavepackets educed from isothermal turbulent jets. In particular, we consider the nearfield pressure fluctuations measured by (2). For a detailed description of the measurement setup, data acquisition and postprocessing the reader can refer to that paper. We here limit ourselves to a brief presentation of the data, in order to highlight those features, which are relevant to the analysis we perform.
2.1 Experimental setup
We consider an unforced, isothermal, subsonic jet issuing from a round nozzle with a fully turbulent boundary layers. The flow is thus very different to the transitional jet studied by (42). The setup is shown in Fig. 2, and a schematic of the nearfield microphone array used for eduction of the axisymmetric wavepacket signature is sketched in Fig. 2. The experiments were carried out in the anechoïc free jet facility at the Centre d’Etudes Aérodynamiques et Thermiques (CEAT), Institut Pprime, Poitiers (FR). The jet Mach number ranged from to . In what follows, the diameter of the nozzle, , is taken as the reference lengthscale, while the exit velocity, , is the reference speed. Time is scaled as . The corresponding Reynolds number ranges between , where is the density and the air viscosity. The potential core length of the jet ranges between .
Azimuthal rings of six microphones recorded pressure fluctuations from the nearfield of the jet, on a conical surface. An azimuthal ring of three microphones recorded the farfield fluctuations. In this work, we focus on the nearfield fluctuations; two different datasets are analysed, obtained from two measurement campaigns.

The first set of data is obtained by simultaneous measurements at different streamwise locations in the nearfield. At each location, an azimuthal ring, each with microphones, is positioned. The array thus comprises a total of microphones, placed on a conical surface in the nearfield of the jet, as shown in Fig. 2. The spacing between the rings is . At each location a Fourierseries decomposition is performed in the azimuthal direction.

The measurements of the second campaign are performed using a rig of equispaced, azimuthal rings, which are displaced in the streamwise direction over the range with a streamwise resolution of . This set of measurements allows a more detailed computation of the twopoint flow statistics, providing the crossspectraldensity matrix, an eigendecomposition of which provides POD modes.
This POD basis is used for the projection of the data from the first campaign. In this way, a larger domain along is spanned. The temporal dynamics of the axisymmetric mode is shown in Fig. 3, as a function of time and the streamwise direction. This spacetime picture is complementary to the frequencyspace realisation shown in Fig. 1. Both show clearly the highly organised nature of wavepackets in these highReynoldsnumber, fully turbulent jets.
Our study focuses on the temporal dynamics of the axisymmetric mode . The analysis is carried out by considering the timeseries extracted at given streamwise locations, for each of the datasets. By definition, a timeseries is a sequence of equispaced data points, function of time. For our case, the sampling rate is , and the total number of points for each series is .
In the following section, preliminary analysis of the temporal behaviour is performed by means of spectral analysis and estimation – using nonlinear tools – of the embedding and correlation dimensions associated with these wavepackets.
3 Part I: Timeseries analysis
In this section, we consider both the original data and the timeseries postprocessed by means of properorthogonal decomposition (POD). In particular, we study the runs at and Reynolds number , at different locations along the streamwise direction (, ).
The main goal in this section is to characterise the time–series from a dynamicalsystem point of view. Although we are considering a fully turbulent case, we study a projection of the system on a set of linear, axisymmetric wavepackets. Thus, the working hypothesis is that the resulting temporal dynamics is lowdimensional. In principle, a dynamical system can be qualitatively analysed by reconstructing the – unknown – underlying attractor from a sequence of observables of its state, for instance the timeseries. This idea is at the heart of the Takens’ theorem (43), a delay embedding theorem providing the conditions under which the reconstructed dynamics preserves the properties of the original system. In particular, the theorem ensures a diffeomorphism between the original – unknown – phase space and the reconstructed phase space. In that sense, there are two crucial parameters: i) the minimal embedding dimension , i.e the dimension of the projection phase space of the reconstructed dynamics; ii) the correlation dimension , a measure of the fractal dimension of an attractor (see (44)). According to these values one can evaluate the dimensions of the attractor of the system and – if possible – study its projection.
In general, the available time series are univariate; thus, the embedding is typically the first step for the qualitative analysis of a dynamical system. Different procedures of embedding are available. For instance, a singularvalue decomposition of the Hankel matrix based on univariate timeseries can be performed, where the dimensions of the rectangular matrix are the length of the time series and the maximum embedding. A second alternative consists of computing higherorder derivatives of the timeseries, such that each of the derivatives corresponds to one of the columns of the embedding vector. In this work, we rely on the most classical approach: once a maximal embedding dimension is chosen, each coordinate is obtained by a time shift ; thus, introducing the time series , the embedding vector is obtained as
(1) 
where indicates the timespan; note that the last channel will have a total shift of with respect to the first one. Once an embedding vector is obtained, we need to identify the minimal embedding dimension that allows determination of the dimension of the (hypothesised) attractor. A small indicates the possibility of unfolding the attractor. A popular method for estimating the minimum embedding dimension is the false nearest algorithm, (45); (46).
The second dimension of interest is the correlation dimension , providing an estimate of the fractal dimension of an attractor. This quantity is related to the embedding dimension by the Takens’ theorem, stating that ((43)).
In what follows, we adopt statistical tools described from the theoretical point of view in the book by Kantz & Schreiber (45). A thorough discussion of these techniques is beyond the scope of this paper: we will focus mainly on the validation of the results. Unfortunately, these statistical methods are capable of producing results also in the presence of purely random timeseries. Thus, the validation has a twofold goal: assessing the robustness of the results, and confirming the determinism of the dataset.
3.1 Filtering and spectral analysis
In order to analyse the timeseries, we first perform nonlinear noise reduction. Indeed, noise is a dominant, limiting factor for the embedding procedure (45). The denoising algorithm is a movingaverage filter, applied along the trajectory identified in the embedding space, on the assumption that the dynamics is continuous. In Fig. 4, the timeseries is shown at . In particular, we are interested in assessing to what extent the applied filters modify the essential dynamics of the timeseries. In the following, for ease of discussion, the original timeseries is denoted TSRaw while data filtered using POD as the projection basis is denoted TSPOD.
In Fig. 4, the TSRaw data are compared with the TSPOD set in the timedomain. A quantitative assessment is reported in Fig. 4 using Welch’s power spectral density (PSD) estimate. At the significant Strouhal numbers, , it can be see that the POD filter does not strongly modify the dynamics of the timeseries. Differences can be observed only at . At lower , we observe the cutoff due to highpass filtering of the original data applied to remove energy below the anechoïc cutoff frequency of the windtunnel. The two data sets are filtered using the nonlinear filters in Fig. 4 and Fig. 4. By definition, the nonlinear filter does not act on specific bandwidths, but on the whole spectrum, while preserving the foliation of the attractor in the embedded space. Welch analysis confirms this behaviour: in Fig. 4 we can observe that the nonlinear filter acts on the whole range of numbers. The standard deviation of the removed noise is and , for TSRaw and TSPOD, respectively.
The spectral content of the TSPOD data at is compared in Fig. 5. We observe that the maximum fluctuation levels move progressively to lower frequencies as the observation position moves downstream, from at to at . A change in the maximum amplitude can also be observed. An alternative analysis is provided by the spectrogram, shown in Fig. 6 for the TSPOD at . The graph is in two dimensions: the horizontal axis represents the time scale, while the number is reported along the vertical axis. The amplitudes are shown as a colour map: at each instant, the intensity of the fluctuation is shown as a function of . The Welch estimate can be regarded as an average along the timespan of the results shown in the spectrogram. In particular, in Fig. 6, two instants are shown at and , respectively. The instantaneous amplitudes of certain frequencies in the bandwidth can be observed; this behaviour is a timefrequency view of the organised behaviour manifest in spacetime in Fig. 3. These behaviours are the signature of the amplifier nature of the flow: wavepackets are hydrodynamic instability waves that exist on the turbulent mean flow. These waves acquire initial amplitudes and phase at upstream stations – either from disturbances issuing from the turbulent motions within the nozzle, or from distributed turbulence that acts as a volume force – and evolve in the downstream direction according to the stability properties of the linear operator, possibly subject also to distributed forcing from background turbulence; this last issue is one that is presently being studied.
3.2 Embedding dimension and correlation dimension
Assessment of the correlation dimension, , and embedding dimension, , is performed for both TSRaw and TSPOD datasets, filtered nonlinearly. The results are collected in Fig. 8, as a function of the streamwise direction .
Minimal embedding dimension
The minimal embedding dimension is estimated by using the algorithm described by Cao in (46). The algorithm does not strongly depend on the length of the time series and it is capable of estimating the embedding dimension also for time series describing highdimensional attractors. More importantly, it allows to clearly determine whether the signal is deterministic or stochastic, as shown in Fig. 7. The dimension is estimated by analysing the behaviour of the parameter : the saturation indicates the minimal embedding dimension. We chose it using a threshold value , giving . The quantity explicitly accounts for the correlation: random series will be characterised by , for any . For deterministic data, cannot be a constant. In our case, we find that the values of are not constantly unitary, as shown in Fig. 7: this is a first clue that we are analysing deterministic time series.
A critical aspect is the choice of the embedding delay . In principle, the embedding dimension is independent of the time delay; in practice, the minimum embedding dimension estimate may depend on this choice and needs to be verified case by case. In Fig. 8, we identify a confidence interval for the estimates; the upper bound corresponds to a timedelay , while the lower bound is obtained for . The choice of the embedding delay is done by analysing the autocorrelations. points have been used for the computations, although we verified that is scarcely affected by this choice.
The confidence intervals are shown using shadowed areas. Each point corresponds to the value where the parameter saturates, as discussed in Fig. 7. Dotted curves identify the minimal embedding dimension computed by averaging the results obtained with points and . It can be observed that the TSRaw series is characterised by a rather constant value of , oscillating between a minimum value at and maximum values further upstream, at . This estimate follows closely the upper bound of our confidence interval. The lowest value is found at , for the runs at .
The behaviour of the TSPOD series is coherent with the original data, although we found higher values. This apparently counterintuitive behaviour is discussed later (see Sec. 3.2.3). For values , the minimal embedding dimension is practically constant until ; however, note that in this region the data are extrapolated, so relevant dynamics carried by higher order POD modes is missing.
Correlation dimension
The correlation dimension is estimated using the routines included in (47), based on the GrassbergerProcaccia algorithm (45). Convergence tests were performed over the input parameters of the algorithm: the embedding quantities , the time delay, and the Theiler window. The Theiler window accounts for the possible oversampling. In fact, the computation of is based on a correlation sum: in the timeseries, successive elements are generally not independent and can be highly correlated. To limit this effect, that could lead to inconsistent results, a timeshift – i.e. the Theiler window – is introduced to reduce the correlation between points during the pair counting.
A limit characterises the computation of the correlation dimension; in particular, the value of the correlation dimension over a decade cannot exceed , where is the number of points in the time series (EckmannRuelle limit, (48)). The amount of data available allows us to get a value of . This means that when we approach this upper bound, it is not legitimate to conclude that the value corresponds to the dimension of the inner dynamics.
In Fig. 8, the results are shown with squareddotted curves. For each positions of the available dataseries, we analysed two subsequent blocks of data containing a total of points; thus, only values of correlation dimension are considered relevant (a green shadowed area indicates the limit). The final value of is obtained by averaging the results of the blocks. We observe that the timedelay imposed for the embedding procedure does not influence the final result. A Theiler windows of was chosen; we do not observe relevant changes when increasing this parameter.
For the series TSRaw, we observe a behaviour similar to , with in the span . For the TSPOD, consistent with the previous estimate of , we found that the value of are slightly higher. In the upstream region , the value close to the upper bound does not allow us to conclude that this is the effective value of the correlation dimension. However, the dimension with respect to the embedding dimension is bounded within the limits imposed by the Takens theorem.
A decrease is observed at the location extrapolated in the vicinity of the nozzle and far downstream, for locations . Also in this case, low values observed in the TSPOD data are due to the POD filtering suppressing important dynamics in these regions of the flow (see (2)).
A brief discussion on the estimated dimensions
The preliminary analysis of the minimal embedding dimension and the correlation dimension indicates deterministic behaviour, despite the limitations that both algorithms pose. Interestingly, the values of and are higher when considering the PODfiltered data. This result has been found also by computing the embedding dimension with a different algorithm based on the false nearest strategy (results are not reported). We believe that the behaviour is due to the physics of the system described in the different datasets and it is not a numerical artefact. As mentioned in Sec. 2, the simultaneous pressure measurements of TSRaw are projected on a basis of PODs obtained using twopoint flow statistics on a larger number of points in the streamwise direction leading to a finer resolution. Although from the spectral point of view the resulting timeseries are equivalent, except at higher (see Fig. 4), from the dynamical point of view the projection on the PODbasis might lead to a richer dynamics due to the increased accuracy of the measurements. When considering the dimensions outside the range , i.e. where the simultaneous measurements are not available, we observe a rather low correlation dimension, while the embedding dimension is constant until . This effect is due to the extrapolation of the data in this region of the flow. For this reason, in the following sections we focus only on the points .
The embedding space that should be used according to the Takensâ criterion is rather large. In Fig. 9 and Fig. 10, we show three phaseportraits reconstructed in a three dimensional embedding space obtained by imposing a delay of , between each of the three coordinates. In all cases, at a first glance, the reconstructed phase space suggests a toroidal nature for the dynamical system; unfortunately, the projection is poor due to the highdimensionality and cannot be further investigated. Moreover, the first return map in Fig. suggests contamination of the dynamics due to external forcings. Similar results were obtained for all the positions along regardless of which embedding strategy was used (SVD embedding and derivative embedding).
In conclusion, the results are consistent with a determinism of the timeseries under investigation and suggest that the dynamics is relatively lowdimensional. However, the embedding and correlation dimensions are still too large for useful characterisation of the attractor geometry. For this reason, alternative strategies are necessary for exploration of the dynamics of the system. In the next section, system identification is applied: with the same tools, we attempt to replace the original data with surrogate models that capture a part of the dynamics.
4 Part II: A systemidentification modelling approach
In the previous section, the preliminary study of the timeseries – based on statistical tools and Fourier analysis – provided clues that the underlying temporal dynamics of the axisymmetric coherent structures is deterministic. However, the dynamics educed locally is rather rich and cannot be analysed by simply using the standard tools from dynamical system theory. From a global view point, the turbulent jet is characterised by a strong convective behaviour: the intrinsic dynamics is sensitive to external forcing (49). Thus, one might wonder to what extent the behaviour educed from the timeseries analysis is due to intrinsic nonlinear interactions, external activation of degrees of freedom driven by external (stochastic) forcings or a combination of intrinsic and extrinsic mechanisms.
In an effort to understand the underlying dynamics of the system, surrogate models can be introduced by applying system identification (45); (50). By definition, system identification aims at building mathematical representations of dynamical systems from set of observables, i.e. measured data. In other words, given a dynamical system, we can analyse it from its observables and infer a model reproducing the dynamics (or part of it). Starting from the hypothesis that the chosen observables properly represent the system behaviour, this approach can be summarised in the following steps:

Model validation – among all the possible models, only a few allow a proper characterisation of the dynamical system. Validation is necessary in order to discard the models unable to reproduce a dynamics close to the real one.

Prediction or dynamical analysis – according to the goal of the identification process, we can aim at minimising the prediction error given a time horizon or analysing the dynamical properties of the system.
The last aspect is particularly relevant: according to the goal of the identification process, the choices of the inputs/outputs as well as the model structure are of great importance.
In the following, we will often refer to noise indicating with this term the external forcings driving the system dynamics. As already said, turbulent jets are convectively unstable, so strongly sensitive to external forcings; measurement noise will be generally indicated as error or residual in the context of the identification problem.
4.1 Choice of the model, inputs and outputs
For our application, the chosen model consists of nonlinear Volterra series of polynomials. Hereafter, the outputs will be denoted by and will be represented by the measurements taken at different locations along the streamwise direction. By including also the inputs, indicated by , the model reads
(2) 
This relation is called the equation error model. On the righthand side, the first term relating the past outputs at the present is referred to as the autoregressive term; the second term reproduces the dynamics between the inputs and the outputs (exogenous part). The unknowns of the model are the set of coefficients , , the nonlinear part – here represented by highorder polynomials – and the error . A classical way to model the error is to consider it as a movingaverage of unknown whitenoise. The most general strategy of identification based on Eq. 2 is the nonlinear (N) autoregressive (AR) exogenous (X) with moving average (MA) algorithm, usually referred as NARMAX.
A quick overview of the algorithm is provided in A, while for a deeper discussion we refer to (51); (52); (53). For our purposes, it is interesting to focus here on the choice of the inputs; in particular,

if we are interested in revealing the core mechanisms of a dynamical system, we can neglect the input such that all the dynamics is described by the autoregressive part of the model and it is not contained in the inputs; in this case, the model is a discrete representation of an autonomous system.

if a nonautonomous representation is sought, an optimal choice is represented by inputs modelled as stochastic white noise or random binary sequences. In this case, all the dynamics will be described by the identified inputoutput model.

a third choice consists in introducing as input one of the available outputs , thus containing already some of the dynamics. In this case, due to the correlation between input and output, the resulting model will be an optimal transfer function among the elements; this approach is controloriented.
In sec. 5, we use full nonlinear ARMA models; the main goal is to identify the essential dynamics using autonomous models, in order to educe the intrinsic mechanisms at work. In sec. 6, we aim at identifying models that minimise the prediction error in terms of best meansquared errors over a timehorizon, given a measurement as input ; in this case, full nonlinear ARMAX models are computed.
The models discussed in the next sections were computed using the tools developed by Luis A. Aguirre and coworkers (see (50)). The linear models were crossvalidated with inhouse scripts (53). We use polynomial expansions, but alternative approaches are available. For instance, nonlinear combinations of elementary functions can be used. In this case, the system identification problem changes: we do not seek only the coefficients of the model but also the optimal combination of functions fitting the original data; due to the necessity of identifying also the basis of functions, these techniques require full nonlinear optimisation (54) or statisticallearning techniques (see (55), for instance).
5 Nonlinear models: dynamical system analysis from the observables
In this section we consider nonlinear ARMA modelling. With respect to the full equation error model in Eq. 2, we neglect the exogenous terms. The basic model reads
(3) 
The identification process consists in identifying the number of necessary terms and coefficients for an optimal representation of the timeseries over a window of assimilation.
This approach leads to a discrete, autonomous model. The unmodelled part is accounted for in the error term . The final analysis is performed on a surrogate timeseries obtained by running the model over long timehorizons, not driven from external known forcings. Thus, we do not expect a predictive model, but rather a phenomenological model. The aim is to mimic a part of the dynamics of the original system.
We analyse the results obtained by applying the identification approach to the timeseries. In particular, we consider the TSPOD at , nonlinearly filtered (see Sec. 3.1). As previously shown, the TSPOD series are less affected from measurement noise, in the range where the simultaneous data from TSSeries are available; moreover, the nonlinear filter preserves the dynamics across all the frequencies. The chosen location corresponds roughly to the end of the potential core.
In Sec. 5.1, we provide a brief description of the parametric analysis performed during the computation of the nonlinear ARMA models. Results are shown and discussed in Sec. 5.2 for three different models.
5.1 Computation of the chosen models and validation
The analysis involves several parameters, namely: i) the polynomial order; ii) the window of assimilation; iii) the maximum number of coefficients ; iv) the optimisation parameters. Due to the nonlinear terms and the complexity of the system, we expect multiple minima for the optimisation process. In that sense, ranging a large span of parameters is necessary. In the following we briefly motivate and explain our choices.
Polynomial order – nonlinearities We consider as maximum nonlinear order ; the choice is a compromise between numerical and computational issues rising from higherorder models and the physical description of the system. In particular, highorder nonlinearities might lead to inconsistencies due to singular entries in the Hankel matrices (see A), used during the identification process.
Assimilation window The model can depend from the chosen assimilation window. Indeed, considering the spectrogram in Fig. 6, we already observed that the dynamics is characterised by events at given frequencies corresponding to organised behaviour in timedomain. In principle, one should consider more windows of assimilation, stacked over the entire timespan; in practice, for computational and convergence reasons, we consider shorter timewindows of assimilation roughly corresponding to the average length of the highamplitude events observed in the spectrogram. We tested windows of length ; in the present manuscript, we consider , that provided best results.
Coefficient This coefficient dictates the second dimension of the regression matrix (see A), the first being the assimilation window. More importantly, it relates past entries with the actual value during the autoregression process. We consider the range ; the maximum value in convective time corresponds to and it is based on the analysis of the autocorrelation. Larger values would lead to erroneous correlations between past and present dynamics.
Optimisation process parameters For the optimisation process, we consider a maximum of iterations for each model and a maximum elements for the resulting models, whose terms are related to the error part.
As a result of the parametric analysis, we identified a total of models. As said, the models can be biased by the window of assimilation. At a first glance, one might think that the choice of relatively short timewindows of assimilation reduces the validity of the models; however, the main aim of the modelling approach is to explore the dynamics behind the organised behaviour observed in Fig. 6, trying to separate the inherent dynamics observable in small windows of observation from the effects of the noisesensitivity. In that sense, it might not be possible to get endogenous models of the system, unless we were to introduce an exogenous forcing decorrelated from the output . For this reason, we believe that a short window of assimilation results in an identification process consistent with a part of the inherent dynamics of the dynamical system.
Validation process
The simulation of the models is not driven by known forcings; thus, these nonlinear models will quickly deviate from the original timeseries after a short transient roughly corresponding to the initial conditions. The validation process is performed in two steps: i) the minimisation of the error during the system identification; ii) the spectral analysis of the model, compared with the original timeseries. In particular, we can not apply the Welch method due to the short timewindow of assimilation; as an alternative, we apply a power spectral density estimation based on the YuleWalker method (56). This method estimates the spectral content by fitting an autoregressive linear prediction filter to the signal. By applying these two criteria and crossing the results, we were able to identify three different “families” of models. Of all the models identified, around show periodicity, quasiperiodicity, toroidal features or chaotic behaviour. The remaining models are badly conditioned, unstable or decaying to zero after a transient. The selected models are robust with respect to the initial conditions, i.e. they show the same asymptotic behaviour in the presence of uncertainties in the initial conditions. With respect to the position , the optimisation procedure leads to models with different coefficients , but similar dynamics from the phenomenological point of view. Thus, the surrogate models discussed here are still able to qualitatively describe the physics of the underlying system.
5.2 Autonomous models
In this section, we discuss three different models, labelled as follows

Model T – in Fig. 11 the model is characterised by . It is a toroidal (T) example, that settles on the frequencies characterising the windows of learning.

Model S in Fig. 15 the model is characterised by . It is a chaotic model, as it can be seen from the firstreturn map. We believe that it reproduces some short (S) time dynamics related to the jittering.
The coefficient gives to the models distinctive features. The first two models are characterised by a long window of correlation. The third model is obtained considering the shortest, possible correlation length. In the following paragraphs, we comment on the significance of these phenomenological models.
Models and
In Fig. 11, the properties of the model are analysed. As already mentioned, none of these models is predictive: after a short time, the dynamics quickly deviate from the original one (Fig. 11). However, the original timeseries are contaminated by the external disturbances, while the model behaviour is dictated by the intrinsic dynamic educed during the assimilation process. The quality of the model can be assessed by analysing the spectral properties in the window of assimilation: in Fig. 11, the power spectral density is shown based on autoregressive modelling; the current model reproduces the peak frequency ; a second peak appears at . However, we can also observe that model is not properly reproducing the low frequencies . The longtime behaviour is characterized by a toroidal behaviour in phase space, showed with a red solid line in Fig. 11. The Fourier analysis in Fig. 11 and the firstreturn map built with the relative maxima in Fig. 11 confirm the phase space representation. In particular, in Fig. 11, the black dots correspond to the beginning of the trajectory, while the superimposed blue dots the final part.
Model is characterised by a limited match in the autoregressive PSD; for this reason, a second model is presented in Fig. 12, labelled as model . In this case, we observe an excellent agreement at low Strouhal numbers, (see Fig. 12). We note a rather different behaviour when comparing the transient with the asymptotic state. During the transient, the spectrum is not characterized by distinct frequencies (see Fig. 12), although we observe the dominance of frequencies around and (due to the assimilation windows). The overall dynamics is more complex as shown in the phase space portrait in Fig. 12; the corresponding firstreturn map in Fig. 12 can be understood by considering the longtime behaviour of the system reported in Fig. 13. Indeed, after a transient of , the model settles on a limit cycle (Fig. 13); the dominating frequencies in Fig. 13 are the same as those observed in the Welch PSD analysed in Fig. 12. Having the longtime behaviour in mind, it is easier to understand the firstreturn map in Fig. 12. In particular, the blue dots are representative of the limitcycle trajectory: we can observe that the model during the transient progressively settles on the limit cycle, that “attracts” the trajectory.
The transient has a qualitative behaviour similar to the original model also if observed from the spectrogram. In Fig. 14, it is possible to compare the original timeseries (inset ) with the model (), in the same timespan. It is evident that on average the length of the organised events is similar, while a significant difference is given by the frequency range; in particular, as already observed, this feature is an attribute of the extrinsic dynamics at work on the dynamical system. Thus, these two models suggest that if only the endogenous dynamics were considered, the dynamics due to nonlinear interactions among wavepackets would be dominated by the organised behaviour that characterises the “shortliving” events observed in the spectrogram. In particular, the dynamics would evolve after some transient on the surface of a torus (T2) or a limit cycle. In models however, one limit cycle (at least) is stable and attracting. Therefore, the hypothetic underlying chaotic toroidal set is only transient.
Model S
In this section, we analyse a third model, depicted in Fig. 15 and based on a shorter window of correlation. Also in this case, the peak in frequency attains at , as shown in the autoregressive PSD analysis in Fig. 15; this same dominance is observed in the Welch analysis in Fig. 15. The resulting behaviour in phase space and in the firstreturn map is rather distinctive and different compared to the previous models. In particular, this model appears to slowly evolve in time; in this case, we can conjecture that the model captures some inherent mechanisms related to combination of wavepacketspairs being driven by the stochastic turbulence in the real case.
Discussion
Within the limits of the methodology, an analysis of the models provides some hints regarding the possible dynamics underlying the system. The main idea is that the jittering is a mixture of stochastic and deterministic components. Model T and Model LC seems to indicate that frequencies can lock together, leading to nonlinear wavewave interactions. Eventually, without external forcing, they can set on limit cycles or in trajectories embedded on a toroidal surface in phase space.
In this sense, the jitter signature is now more deterministic, but, in the real case, it is ’randomised’ due to continuous forcing by stochastic background turbulence. In fact, because the flow is a convective amplifier, the wavewave interactions can never fully take over and the asymptotic behaviour can never happen because information is continuously convected out of the system. The dynamical systems reproduced by Model T and Model LC evolve according to initial conditions and a dynamic law that cannot evacuate the information convectively.
One might imagine the phase space as a dense ensemble of trajectories evolving on limit cycles and toroidal surfaces, associated with wavewave interactions. Due to the sensitivity of the system to external noise, the system will continuously and erratically deviate from these attracting solutions without setting on any one of them as observed in the spatiotemporal plots and the spectrogram in Sec. 3. The high sensitivity to external disturbances is inherently connected to the convective nature of the flow. This picture can be linked to recent ideas regarding wavepacket theory: the organisation observed is understood as being due to the resolvent operator forced by background turbulence. In this sense, the resolvent operator is an organiser of the flow, while the turbulence that forces wavepackets via this operator is the disorganiser.
Also, note that the resolvent analysis has proved to be particularly efficient for , where the linear wavepackets are in remarkable agreement with the coherent structures observable in the jet flow. On the other hand, at low – and particularly at low Mach number – nonlinear wavewave interactions may be more relevant as the linear analogy fails: this makes our conjecture particularly appealing for explaining this discrepancy.
6 Inputoutput modelling for control: is linear modelling enough?
When the forcing term is represented by an output of the system, used as an input of the model, system identification identifies a transfer function between inputs and outputs. In particular, we use the output at and inputs placed in four positions, . The models include only one input (singleinputsingleoutput modelling, SISO); the performance are shown in Fig. 16, where we compare the model prediction with the original timeseries. The window of assimilation is . The validation is performed by integrating the model in the window .
The application of the full nonlinear ARMAX algorithm results in a polynomial expression: the final models are fully linear and based on a limited number of coefficients and , . The good results can be explained by observing the spectral analysis in Fig. 5. As already noted, the convective nature of the flow implies a progressive drift of the main frequency of amplification and a slight variation of the maximum amplitude in the spectrum as we move further downstream; in this sense, we need only a limited number of degrees of freedom for modulating amplitudes and frequencies, and correctly representing the timedelays of the system. This observation justifies the small dimensions of the models, that here act as filters.
The standard deviation between the resulting model and the original timeseries is generally low. In Fig. 17, the results are summarised for a number of outputs along the streamwise direction. In particular, six different positions along are introduced as input, in the range , while as outputs are considered the locations downstream of the inputs between and . Except the transfer functions with inputs at and output at , characterised by a standard deviation , in all the other cases the standard deviation is below . As already mentioned above, all the shown models are linear: this observation has interesting implications when considering the design of linear controllers.
These results confirm and extend the findings of (38), where linear ARMAX modelling, parabolised stability equation (PSE) and empirical transfer function were compared for the estimation of downstream propagating wavepackets, showing equivalent results. With respect to the work of (38), in our analysis the linearity of the model is a result of the optimisation process and not a working hypothesis: despite the fact that candidate models introduced in the identification procedure were characterised by all the terms of Eq. 2 up to order , nonlinearities appear to be unnecessary for an optimal representation of the transfer functions.
6.1 Robustness and control
The limited amount of dynamical information retained in these models lead to robustness issues. The set of coefficients is obtained as leastsquares solution from the assimilation datasets at given amplitude. The transfer function only retains information regarding the timedelays of the system, the modulation of the amplitude and the frequency bandwidth. Any sources of uncertainty can lead to robustness issues. As an example, in Fig. 17, we illustrate the lack of robustness with respect to timedelays, due to the convective time of propagation of the perturbations. The model assimilated using the input at is applied with different inputs, taken in the range . The validity of the models quickly decays. Values of the standard deviation indicate phaseopposition of the resulting model with respect to the original timeseries.
The robustness issues can be attacked using different strategies. Control design can be conceived by accounting for the robustness issues arising from uncertainties of the system within the robust framework or the adaptive control theory (see (57); (58)). In particular, with respect to the convective nature of the flow, feedforward schemes combined with selftuning, adaptive controllers can be applied. These controllers are characterised by two timescales: a fast timescale, for the realtime estimation of the input , and a slow timescale, based on feedback measurements of the flow necessary for correcting offdesign conditions. Examples of these strategies can be found for weakly nonlinear transitional flows (see (59)) or nonlinear, turbulent cases (see (60)).
7 Conclusions
Coherent largescale structures in turbulent jets – usually referred as wavepackets – have been subject of numerous investigations for the role that such structures might play in sound radiation. Remarkably, it is well established that certain aspects of the wavepacket dynamics can be modelled by means of linear ((1)). In this work, we considered an unforced, isothermal, subsonic jet issuing from a round nozzle with a fully turbulent boundary layers at and . The temporal behaviour of the axisymmetric wavepackets is analysed; starting from the studies by (2), we investigate the experimental, spatiotemporal measurements from a dynamical point of view. To the best of our knowledge, this is the first time such an analysis has been attempted. The key point is the remarkable spatiotemporal coherence of the measurements: despite the high Reynolds number, the lowdimensional projection of the dynamics reveals a wellorganised behaviour. This aspect raises questions regarding the temporal behaviour of the wavepackets, the nature of the nonlinearities involved and the possibility of designing control strategies starting from this knowledge.
In order to answer these questions, a variety of strategies ranging from statistical tools to system identification have been adopted. The estimation of the correlation dimension and the embedding dimension of the dynamical system confirmed that the organised temporal behaviour is lowdimensional and deterministic. The minimal embedding ranges between and , while the correlation dimension is rather low at locations corresponding to the end of the potential core, . These dimensions are nonetheless still too large for useful characterisation of attractor geometry.
For this reason, alternative strategies are necessary for exploring the dynamics of the system. A system identification approach is pursued. The timeseries taken at a streamwise location of the jet is assumed as output; we replace the original dynamics with surrogate, autonomous models using dataassimilation based on this output. These models are limited, as they reproduce only a part of the dynamics and are influenced by the choice of the parameters.
By comparing the results from three different classes of models, we suggest that the jittering dynamics might be described as a combination of extrinsic and intrinsic mechanisms at work. We observe that our surrogate models reproduce the ideal case where the background turbulence is not active; in this limit, the dynamics of the system would be dominated by wavewave interactions leading to stable limit cycles and/or solutions lying on toroidal attractors. Evidence on this behaviour is provided by the spectral analysis of the transient dynamics of these solutions that qualitatively reproduces the local (in time) dynamics of the original system. However, in the real case, the extrinsic dynamics drives the system: wavewave interactions can never fully take over and the asymptotic behaviour of the ideal models can never happen. This interpretation justifies the selforganised features appearing in the timefrequency analysis of the timeseries and it is consistent with the nature of the jets as noiseamplifiers: due to the sensitivity of the system to external noise, the trajectory in phase space will continuously escape from these attracting solutions.
Given these interpretations, control strategies can be conceived. Due to the convective nature of the system and its implications regarding the nonlinear interactions, system identification is used for computing nonautonomous models, where the inputs are represented by local measurements taken upstream of the outputs. These models are filters that i) account for the timedelays and ii) modulate the amplitudes and the frequencies already existing in the inputs. In particular, we found that nonlinearities are unnecessary for an optimal representation of these models. In other words, if the inputs contain a large part of the dynamics, a linear filter is enough to reproduce the temporal behaviour at downstream locations due to the convective nature of the system. The main limit is the robustness of these linear filters; however, in our opinion, this might be more a problem of controllability of the flow, rather than a lack observability and predictability. Thus, we believe that the robustness limitations may be tackled in the control design, by applying strategies such as adaptive control or loopshaping robust control, and verifying the effects of the actuation on the real flow.
Future work should be devoted to further improving the nonlinear models for the dynamical analysis. While accepting the impossibility of capturing all the features of the system in lowdimensional models, a possible path to improve our understanding of the jittering would be the identification of surrogate models where a stochastic excitation is also included; in principle, this strategy would help in separating the instrinsic mechanisms from the extrinsic ones. In this sense, we believe that machine learning and statistical learning techniques may be robust alternatives to classical nonlinear system identification. Indeed, these techniques are quickly spreading in the community as a powerful tool for system identification (see (55)).
The authors wish to acknowledge Luis A. Aguirre for sharing his system identification package and Christophe Letellier for fruitful discussions. This work is supported by the Agence Nationale de la Recherche (ANR) under the ”Cool Jazz” project, grant number ANR12BS090024.
Appendix A System identification using polynomials expansion
We introduced in Sec. 4 the basic ideas behind the application of system identification for the qualitative analysis of dynamical system and the definition of controloriented models. The aim of the section is to provide a quick overview of the algorithms.
We rely on an inputoutput description of timeinvariant systems. The aim of our models is reproducing the behaviour of the output by defining a nonlinear model of the dynamics of the measurement or a transfer function between the input and the output . For sake of conciseness, we consider a singleinputsingleoutput (SISO) case; nonetheless, the extension to the multiinputsmultioutputs cases is rather straightforward from the theoretical point of view.
By assuming a SISO system, we can generalise the equation error models in terms of Volterra series
(4)  
where are the order Volterra kernel. The error can be modelled as a movingaverage of (unknown) whitenoise
(5) 
where are the coefficients for the error modelling. The kernel’s coefficients and the error are the unknowns of our problem. The Volterraseries are an extension of the linear convolution integral and represent nonlinear systems as a series of multisummations (or integrals) of the Volterra kernels and the inputs (61). We can note that the complete model consists of a linear combination of three groups of terms

the autoregressive terms related to the output ;

the exogenous terms due to the input ;

the movingaverage description of the error .
Additionally, the nonlinear combinations are introduced in the full model.
a.1 Solution of the identification problem: a basic ARMAX algorithm
The first step of the identification procedure consists in a the preliminary computation of the coefficients . To this end, we assume that the error is null, such that the leastsquares regression can be performed. Note that the error is modelled as a moving average of white noise , leading to a coloured noise with expected value ; the estimate of the coefficients will be affected by a bias, due to the expected value of the error . In particular, the error is computed as
(6) 
where is a rectangular Hankel matrix including inputs and outputs . The Hankel matrix has the dimensions of the number of the coefficients and the window of assimilation. We can introduce a matrix of the same dimensions as the Hankel matrix such that
(7) 
Thus, with a proper choice , the “debiased” estimate of can be computed as
(8) 
if is invertible. This generalisation of the leastsquares problem is known as Instrumental Variable (IV) method and allows a first estimation of the coefficients in the expansion, while neglecting the error.
The IV method is iterative. For the first iteration, the matrix of instruments is written as ; a rather effective choice consists in choosing , i.e. the delayed input. These values are relaxed till convergence.
The IV method allows to compute the coefficients and , reducing the effects of the bias due to the nonwhiteness of the error . The unknowns of the movingaverage part, namely the coefficients and the whitenoise series can be extracted using appropriate algorithms (see (52) or (53)).
Thus, we can divide the basic algorithm into two steps:

Evaluate the coefficients and using the iterative IV method.

Apply the movingaverage algorithm for defining the coefficients and the whitenoise .
This prototypical algorithm allows a first estimate of the kernel coefficients and the error of the moving average description. We can improve these estimates using the extended leastsquares (ELS) approach; in particular, we can consider the first estimation of the error as an input of the system: in this way, the estimates will be debiased by explicitly accounting the error. Moreover, an iteration can be cast: the difference between the prediction and the original output, allows to refine again the kernels estimate and further improve the prediction.
a.2 Other techniques
The main difference between the linear modelling and the full nonlinear Volterra series is mainly related to the dimensions of the problem and the risk of illconditioness. Already cubic nonlinearities may introduce nearly null entries, leading to rankdeficiency and singularities. To this end, the regression step can be replaced by wellconditioned algorithms, such as the orthogonal leastsquares (OLS); this iterative scheme introduces in the leastsquares process a GramSchmidt orthogonalisation, which ensures that each new column added in the matrices is orthogonal to all previous columns. Moreover, this process allows us to estimate the relative relevance of the kernel’s coefficient in residual terms; thus, elements that are not relevant are discarded during the process. The OLS approach is described by (61) or (50), and allows the identification of minimal nonlinear models by keeping only the relevant terms of the model. This is the strategy implemented in the routines adopted in our manuscript.
References
Footnotes
 journal: Physical Review Fluids (PRF)
References
 P. Jordan, T. Colonius, Wave packets and turbulent jet noise, Annual Review of Fluid Mechanics 45 (2013) 173–195.
 D. E. S. Breakey, P. Jordan, A. V. G. Cavalieri, O. Léon, M. Zhang, G. Lehnasch, T. Colonius, D. Rodriguez, Nearfield wavepackets and the farfield sound of a subsonic jet, in: 19th AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. 2083, 2013, pp. –.
 E. MollöChristensen, R. Narasimha, Sound emission from jets at high subsonic velocities, Journal of Fluid Mechanics 8 (01) (1960) 49–60.
 E. MollöChristensen, M. A. Kolpin, J. R. Martucelli, Experiments on jet flows and jet noise farfield spectra and directivity patterns, Journal of fluid mechanics 18 (02) (1964) 285–301.
 E. MolloChristensen, Jet noise and shear flow instability seen from an experimenter’s viewpoint, Journal of Applied Mechanics 34 (1967) 1–7.
 S. C. Crow, F. H. Champagne, Orderly structure in jet turbulence, Journal of Fluid Mechanics 48 (03) (1971) (547–591).
 J. C. Lau, M. J. Fisher, H. V. Fuchs, The intrinsic structure of turbulent jets, Journal of Sound and Vibration 22 (04) (1972) 379–384.
 A. Michalke, H. V. Fuchs, On turbulence and noise of an axisymmetric shear flow, Journal of Fluid Mechanics 70 (1975) 179–205.
 R. R. Armstrong, A. Michalke, H. Fuchs, Coherent structures in jet turbulence and noise, AIAA Journal 15 (7) (1977) 1011–1017.
 C. J. Moore, The role of shearlayer instability waves in jet exhaust moise, Journal of fluid mechanics 80 (1977) (321–367).
 A. Michalke, Instability of a compressible circular free jet with consideration of the influence of the jet boundary layer thickness, Z. fur Flugwissenschaften,(WestGermany) 19 (8) (1971) (319–328).
 D. Bechert, E. Pfizenmaier, On the amplification of broad band jet noise by a pure tone excitation, Journal of Sound and Vibration 43 (3) (1975) 581–587.
 J. Cohen, I. Wygnanski, The evolution of instabilities in the axisymmetric jet. part 1. the linear growth of disturbances near the nozzle, Journal of Fluid Mechanics 176 (1987) 191–219.
 D. G. Crighton, M. Gaster, Stability of slowly diverging jet flow, Journal of Fluid Mechanics 77 (02) (1976) (397–413).
 P. Plaschko, Helical instabilities of slowly divergent jets, Journal of Fluid Mechanics 92 (02) (1979) 209–215.
 A. Michalke, A wave model for sound generation in circular jets, Zentralstelle für Luftfahrtdokumentation undinformation, 1970.
 A. Michalke, An expansion scheme for the noise from circular jets, Tech. rep., DTIC Document (1971).
 K. A. Bishop, J. E. FfowcsWilliams, W. Smith, On the noise sources of the unsuppressed highspeed jet, Journal of Fluid Mechanics 50 (01) (1971) 21–31.
 S. C. Crow, Acoustic gain of a turbulent jet, Phys. Soc. Meeting, Univ. Colorado, Boulder, paper IE 6.
 D. G. Crighton, Basic principles of aerodynamic noise generation, Prog. Aerospace Sci. 16.
 D. G. Crighton, P. Huerre, Shearlayer pressure fluctuations and superdirective acoustic sources, Journal of Fluid Mechanics 220 (1) (1990) 355–368.
 T. Suzuki, T. Colonius, Instability waves in a subsonic round jet detected using a nearfield phased microphone array, Journal of Fluid Mechanics 565 (2006) 197–226.
 A. V. G. Cavalieri, D. Rodríguez, P. Jordan, T. Colonius, Y. Gervais, Wavepackets in the velocity field of turbulent jets, Journal of Fluid Mechanics 730 (2013) 559–592.
 P. Jordan, T. Colonius, G. A. Bres, M. Zhang, A. Towne, S. K. Lele, Modeling intermittent wavepackets and their radiated sound in a turbulent jet, in: Proceedings of the Summer Program, 2014, p. 241.
 J. W. Nichols, S. K. Lele, Global modes and transient response of a cold supersonic jet, Journal of Fluid Mechanics 669 (2011) 225–241.
 X. Garnaud, L. Lesshafft, P. Schmid, P. Huerre, The preferred mode of incompressible jets: linear frequency response analysis, Journal of Fluid Mechanics 716 (2013) 189–202.
 O. Schmidt, A. Towne, T. Colonius, P. Jordan, V. Jaunet, A. V. Cavalieri, G. A. Brès, Superand multidirective acoustic radiation by linear global modes of a turbulent jet, in: 22nd AIAA/CEAS Aeroacoustics Conference, 2016, p. 2808.
 G. A. Brès, J. Jaunet, M. Le Rallic, P. Jordan, T. Colonius, S. K. Lele, Large eddy simulation for jet noise: the importance of getting the boundary layer right, in: 21st AIAA/CEAS Aeroacoustics Conference, 2015, p. 2535.
 M. Zhang, P. Jordan, G. Lehnasch, A. V. G. Cavalieri, A. Agarwal, Just enough jitter for jet noise, in: 20th AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. 3061, 2014, pp. –.
 Y. B. Baqui, A. Agarwal, A. V. Cavalieri, S. Sinayoko, A coherencematched linear source mechanism for subsonic jet noise, Journal of Fluid Mechanics 776 (2015) 235–267.
 A. V. G. Cavalieri, P. Jordan, A. Agarwal, Y. Gervais, Jittering wavepacket models for subsonic jet noise, Journal of Sound and Vibration 330 (18) (2011) 4474–4492.
 A. V. G. Cavalieri, G. Daviller, P. Comte, P. Jordan, G. Tadmor, Y. Gervais, Using large eddy simulation to explore soundsource mechanisms in jets, Journal of Sound and Vibration 330 (17) (2011) 4098–4113.
 F. Kerhervé, P. Jordan, A. V. G. Cavalieri, J. Delville, C. Bogey, D. Juvé, Educing the source mechanism associated with downstream radiation in subsonic jets, Journal of Fluid Mechanics 710 (2012) 606–640.
 A. Towne, T. Colonius, P. Jordan, A. V. G. Cavalieri, G. A. Brès, Stochastic and nonlinear forcing of wavepackets in a mach 0.9 jet, in: 21st AIAA/CEAS Aeroacoustics Conference, 2015, p. 2217.
 G. Tissot, M. Zhang, F. Lajus Jr., A. V. G. Cavalieri, P. Jordan, T. Colonius, Sensitivity of wavepackets in jets to nonlinear effects: the role of the critical layer, Journal of Fluid Mechanics In review.
 N. D. Sandham, C. L. Morfey, Z. W. Hu, Nonlinear mechanisms of sound generation in a perturbed parallel jet flow, Journal of Fluid Mechanics 565 (2006) 1–23.
 V. Suponitsky, N. D. Sandham, C. L. Morfey, Linear and nonlinear mechanisms of sound radiation by instability waves in subsonic jets, Journal of Fluid Mechanics 658 (2010) 509–538.
 K. Sasaki, S. Piantanida, A. Cavalieri, P. Jordan, Realtime modelling of wavepackets in turbulent jets, in: 21nd AIAA/CEAS Aeroacoustics Conference, AIAA Paper, 2015.
 K. Sasaki, G. Tissot, A. Cavalieri, S. Silvestre, P. Jordan, D. Biau, Closedloop control of wavepackets in a free shear flow, in: 22nd AIAA/CEAS Aeroacoustics Conference, AIAA Paper, Vol. In review, 2016.
 N. Gautier, J.L. Aider, T. Duriez, B. Noack, M. Segond, M. Abel, Closedloop separation control using machine learning, Journal of Fluid Mechanics 770 (2015) 442–457.
 A. V. G. Cavalieri, P. Jordan, T. Colonius, Y. Gervais, Axisymmetric superdirectivity in subsonic jets, Journal of fluid Mechanics 704 (2012) 388–420.
 S. Narayanan, G. H. Gunaratne, F. Hussain, A dynamical systems approach to the control of chaotic dynamics in a spatiotemporal jet flow, Chaos: An Interdisciplinary Journal of Nonlinear Science 23 (3) (2013) 033133.
 F. Takens, Detecting strange attractors in turbulence, in: Dynamical systems and turbulence, Warwick 1980, Springer, 1981, pp. 366–381.
 P. Grassberger, I. Procaccia, Characterization of strange attractors, Physical review letters 50 (5) (1983) 346.
 H. Kantz, T. Schreiber, Nonlinear time series analysis, Vol. 7, Cambridge university press, 2004.
 L. Cao, Practical method for determining the minimum embedding dimension of a scalar time series, Physica D: Nonlinear Phenomena 110 (1) (1997) 43–50.
 R. Hegger, H. Kantz, T. Schreiber, Practical implementation of nonlinear time series methods: The tisean package, Chaos: An Interdisciplinary Journal of Nonlinear Science 9 (2) (1999) 413–435.
 J.P. Eckmann, D. Ruelle, Fundamental limitations for estimating dimensions and lyapunov exponents in dynamical systems, Physica D: Nonlinear Phenomena 56 (23) (1992) 185–187.
 P. Huerre, P. A. Monkewitz, Local and global instabilities in spatially developing flows, Annual review of fluid mechanics 22 (1) (1990) 473–537.
 L. A. Aguirre, C. Letellier, Modeling nonlinear dynamics and chaos: a review, Mathematical Problems in Engineering 2009.
 L. Ljung, System identification, Wiley Online Library, 1999.
 R. Guidorzi, Multivariable system identification: from observations to models, Bononia University Press Bologna, 2003.
 O. Semeraro, L. Mathelin, An opensource toolbox for linear systemidentification., Tech. rep., LIMSICNRS, Orsay, France (2016).
 O. Nelles, Nonlinear system identification: from classical approaches to neural networks and fuzzy models, Springer Science & Business Media, 2001.
 S. L. Brunton, J. L. P. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937.
 P. Stoica, R. L. Moses, Spectral analysis of signals, Vol. 452, Pearson Prentice Hall Upper Saddle River, NJ, 2005.
 D. Sipp, P. J. Schmid, Linear closedloop control of fluid instabilities and noiseinduced perturbations: A review of approaches and tools, Applied Mechanics Reviews 68 (2) (2016) 020801.
 N. Fabbiane, O. Semeraro, S. Bagheri, D. S. Henningson, Adaptive and modelbased control theory applied to convectively unstable flows, Applied Mechanics Reviews 66 (6) (2014) 060801.
 O. Semeraro, S. Bagheri, L. Brandt, D. S. Henningson, Transition delay in a boundary layer flow using active control, J. Fluid Mech. 731 (9) (2013) 288–311.
 R. Rathnasingham, K. S. Breuer, Active control of turbulent boundary layers, Journal of Fluid Mechanics 495 (2003) 209–233.
 S. A. Billings, Nonlinear system identification: NARMAX methods in the time, frequency, and spatiotemporal domains, John Wiley & Sons, 2013.
 L. Ljung, System identification toolbox for use with MATLAB.