Transition, intermittency and phase interference effects in airfoil secondary tones and acoustic feedback loop

A large eddy simulation is performed to study secondary tones generated by a NACA0012 airfoil at angle of attack of $\alpha = 3^{\circ}$ with freestream Mach number of $M_{\infty} = 0.3$ and Reynolds number of $Re = 5 \times 10^4$. Laminar separation bubbles are observed over the suction side and near the trailing edge, on the pressure side. Flow visualization and spectral analysis are employed to investigate vortex shedding aft of the suction side separation bubble. Vortex interaction results in merging or bursting such that coherent structures or turbulent packets are advected towards the trailing edge leading to different levels of noise emission. Despite the intermittent occurrence of laminar-turbulent transition, the noise spectrum depicts a main tone with multiple equidistant secondary tones. To understand the role of flow instabilities on the tones, the linearized Navier-Stokes equations are examined in its operator form through bi-global stability and resolvent analyses, and by time evolution of disturbances using a matrix-free method. These linear global analyses reveal amplification of disturbances over the suction side separation bubble. Non-normality of the linear operator leads to further transient amplification due to modal interaction among eigenvectors. Two-point, one time autocovariance calculations of pressure along the spanwise direction elucidate aspects of the acoustic feedback loop mechanism in the non-linear solutions. This feedback process is self-sustained by acoustic waves radiated from the trailing edge, which reach the most sensitive flow location between the leading edge and the separation bubble, as identified by the resolvent analysis. Leading edge disturbances arising from secondary diffraction and phase interference among the most unstable frequencies computed in the eigenspectrum are also shown to have an important role in the feedback loop.


Introduction
Trailing-edge noise is an overriding concern for design of quiet air vehicles. At low and moderate Reynolds numbers, tonal noise becomes an important component of the acoustic spectrum. Several studies on trailing-edge aeroacoustics were conducted starting in the 1970s to examine the tonal noise generation by airfoils (Paterson et al. 1973;Tam 1974;Fink 1975;Longhouse 1977;Arbey & Bataille 1983). Noise measurements were performed by Paterson et al. (1973) for symmetric NACA airfoils at various angles of attack over a Reynolds number range between 10 5 and 10 6 . Their results showed the existence of multiple tones in a ladder-like structure in terms of frequency and freestream velocity. These authors also found strong two-point correlation of surface pressure along the spanwise direction, which indicated that the flow phenomenon associated with tonal noise generation could be modeled as two-dimensional. Tam (1974) suggested that the ladder-like structure is due to a self-excited feedback loop between disturbances in the boundary layer and the airfoil wake. Fink (1975) assumed that the discrete tonal frequencies are related to the laminar boundary layer on the pressure side. In order to elucidate aspects of airfoil noise, Arbey & Bataille (1983) performed experiments in an open wind tunnel for different NACA airfoils at α = 0 • for 1 × 10 5 Re 7 × 10 5 . The aforementioned studies showed that the noise spectrum has a broadband component with a main tonal peak plus a set of equidistant secondary tones due to a feedback mechanism closing at the point of maximum flow velocity along the airfoil. The broadband component was assumed to appear due to scattering of Tollmien-Schlichting (TS) instabilities. For airfoils at incidence, Lowson et al. (1994) found that the presence of secondary tones was related to a separation bubble developed on the airfoil pressure side. In this case, TS instabilities developing along the laminar boundary layer would lead to acoustic scattering on the trailing edge and acoustic waves would then propagate upstream closing the feedback loop, with the separation bubble acting as an amplifier of acoustic disturbances. Nash et al. (1999) performed experimental studies of airfoil noise for a NACA0012 profile up to a Reynolds number of 1.45 × 10 6 and several angles of attack. A closedsection wind tunnel, with and without acoustic-absorbing lining on its walls, was used in the experiments and results from the hard-wall tunnel revealed multiple frequency peaks. However, the authors argued that these tonal peaks were correlated to resonant frequencies of the wind tunnel. Thus, they carried out measurements with lined walls simulating anechoic conditions and a single dominant tone was observed instead of several peaks. Furthermore, no ladder-like structure of tonal frequency was observed, in disagreement with previous studies of Paterson et al. (1973), Fink (1975) and Arbey & Bataille (1983). It is important to mention that secondary tones were often observed in experiments conducted in open-jet facilities.
More recently, Plogmann et al. (2013) also found multiple tones in their experiments (NACA0012, 3.1×10 5 Re 1.5×10 6 , 0 • α 9 • ) and demonstrated that tripping the pressure side boundary layer leads to turbulent flow, eliminating the separation bubble and the secondary tones. These authors emphasize that the feedback loop is extremely sensitive to small variations in the flow conditions what, in turn, lead to changes in the tonal components. This occurs particularly at higher Reynolds numbers, where the flow is more prone to transition. The dependency of angle of attack and Reynolds number on tonal noise emission is highlighted by . These authors performed experiments for a NACA0012 airfoil for 0.3 × 10 5 Re 2.3 × 10 5 and effective angles of attack 0 • α 6.3 • . Tripping devices were applied separately on each side of the airfoil to identify their respective role in the noise generation and it was found that suction-side (pressure-side) events dominate at lower (higher) Reynolds numbers. Moreover, it was observed that, at low angles of attack, interactions between the two sides of the airfoil become increasingly important.
As discussed above, early experiments report the presence of a separation bubble on the airfoil pressure side. Most of these investigations were conducted at higher Reynolds numbers, where the flow was likely turbulent on the suction side. These observations are in agreement with the large eddy simulations of a NACA0012 at α = 5 • and Re = 4×10 5 performed by Wolf et al. (2012a) and Ricciardi et al. (2019). Flow visualization and proper orthogonal decomposition were used in these references to identify coherent structures shed from the pressure side near the trailing edge. Such flow structures were found responsible for the intense tonal noise generation despite the fact that a turbulent boundary layer developed on the airfoil suction side. In this case, the pressure side boundary layer was laminar due to the favorable pressure gradient. In agreement with the experiments of Plogmann et al. (2013), Wolf et al. (2012a) showed that the tonal component vanishes when both boundary layers are tripped.
At lower Reynolds numbers, a laminar separation bubble (LSB) exists on the airfoil suction side and is responsible for the overall flow dynamics and noise generation. In the context of airfoil flows, this flow feature has been studied for different purposes. For instance, direct numerical simulations were performed by Jones et al. (2008) for a NACA0012 airfoil at Reynolds number Re = 5×10 4 , M ∞ = 0.4 and α = 5 • . It was shown that, despite being absolutely stable by means of linear stability analysis, turbulence is self-sustained even with the absence of forcing. Later, experimental investigations by Pröbsting & Yarusevych (2015) with a NACA0012 airfoil at moderate Reynolds numbers 0.65×10 5 Re 4.5×10 5 and α = 2 • present intermittent laminar-turbulent transition that affects the advection of coherent structures from the LSB towards the trailing edge. Further analyses were also presented by Kurelek et al. (2016Kurelek et al. ( , 2018Kurelek et al. ( , 2019, Michelis et al. (2018) and Pröbsting & Yarusevych (2021) to study the impact of acoustic excitation, three-dimensional effects and shedding/merging of vortices from the suction side LSB.
The first numerical simulations on airfoil secondary tones were conducted by Desquesnes et al. (2007) considering two-dimensional flows over a NACA0012 airfoil for a Reynolds number Re = 1 × 10 5 at an angle of attack α = 5 • , and for Re = 2 × 10 5 at α = 2 • . In agreement with most experimental observations, multiple tonal peaks were observed. The previous authors also performed a local stability analysis assuming parallel flow and showed that the main tone frequency radiated to the far-field was close to that most amplified along the pressure side boundary layer. An assessment of the linear dynamics of wavepackets driving the feedback loop mechanism and the flow receptivity to acoustic forcing was investigated by time-marching the linearized Navier-Stokes equations by Jones et al. (2010) for a NACA0012 airfoil at Re = 5 × 10 4 , M ∞ = 0.4 and α = 5 • . Fosas de Pando et al. (2014b) applied bi-global stability analysis to investigate the dynamics coupling the boundary layers on both airfoil sides and the wake. The authors studied the 2D flow over a NACA0012 airfoil for Re = 2 × 10 5 , M ∞ = 0.4 and α = 2 • and found multiple frequencies in the eigenspectrum related not only to the main tonal peak but also to the secondary tones. Prior to this work, the stability analyses were limited a parallel flow assumption, and only a discussion of the dominant tonal frequency was presented. By means of adjoint and resolvent analyses, Fosas de Pando et al. (2014a also identified sensitive regions of the flow which are prone to close the feedback loop mechanism. Hence, linear stability theory has proven to be an important methodology to investigate the generation of airfoil secondary tones and the feedback loop mechanism. Another important observation from Desquesnes et al. (2007) regards the amplitude modulation of velocity fluctuations computed near the trailing edge. These authors discuss that such modulation is caused by interference of vortical structures from both sides of the airfoil which, combined with the feedback loop mechanism, would lead to the presence of multiple tones. Following the discussion on modulation of flow structures, Pröbsting et al. (2014) employed particle image velocimetry to study airfoils at Reynolds numbers 1 × 10 5 Re 2.7 × 10 5 and 2 • α 4 • . They investigated the mechanisms associated with tonal noise generation and the interference effects between suction and pressure sides of the airfoil. For Re ≈ 1.5 × 10 5 , at the lowest angle of attack, the authors showed that the amplitude modulation discussed by Desquesnes et al. (2007) was related to destructive interference of boundary layer instabilities occurring on both sides of the trailing edge. More recently, Ricciardi et al. (2020) performed 2D simulations of a NACA0012 at Re = 1.0 × 10 5 at α = 3 • and showed that the multiple tones are related to modulation of the vortical structures developing on the suction side. In this case, the instantaneous main frequency alternates in time due to phase modulation of the flow structures shed by the laminar separation bubble. As a consequence, multiple equidistant frequencies must appear in the Fourier transform to reconstruct the modulated signals. However, as pointed by Pröbsting et al. (2014), flow transition on the suction side occurs in experiments and their conclusions may be different from those observed in 2D numerical simulations. In this regard, the simulation of Sanjose et al. (2019) for a thincambered airfoil at Re = 1.5 × 10 5 exhibits flow transition where intermittency played a key role in the flow dynamics and noise emission for the configuration investigated. Thus, experimental (Pröbsting et al. 2014) and numerical (Sanjose et al. 2019;Nguyen et al. 2021) findings contradict previous assumptions that 2D simulations are sufficient to explain all the mechanisms of secondary tones.
In this study, a large eddy simulation (LES) is performed to investigate a NACA0012 airfoil at α = 3 • immersed in a freestream flow with Mach number M ∞ = 0.3 and Reynolds number Re = 5 × 10 4 . This condition is chosen based on experimental results from  and Pröbsting & Yarusevych (2015). These authors show that, for this flow configuration, suction side events are responsible for the noise generation without a contribution from the pressure side. Moreover, the previous references show that the flow behavior at this particular angle of attack is somehow independent of the Reynolds number up to a point where the relevant vortex dynamics switches to the pressure side. For the current flow setup, some key aspects are investigated including the dynamics of the suction side separation bubble and its vortex shedding, flow intermittency and its impact on amplitude modulation. Linear stability theory is applied to investigate the presence of sensitive regions and amplification mechanisms within the flow.
The present work is organized as follows: section 3.1 presents the laminar separation bubbles which are observed over the suction side and near the trailing edge, on the pressure side. Then, flow visualization is employed in section 3.2 to investigate vortex shedding aft of the suction side separation bubble. It is shown that vortex interaction results in merging or bursting such that coherent structures or turbulent packets are advected towards the trailing edge, leading to different levels of noise emission. Despite the intermittent occurrence of laminar-turbulent transition, the noise spectrum presents multiple equidistant secondary tones as discussed in section 3.3. To understand the role of flow instabilities on tonal noise generation, in section 3.4, the linearized Navier-Stokes equations are examined in its modal form using bi-global stability and resolvent analyses, and by the time evolution of disturbances which results in periodic wavepackets. Aspects of the acoustic feedback loop mechanism in both the linear and non-linear solutions are examined in section 3.5 and we highlight the role of leading edge disturbances from secondary diffraction in the closure of the acoustic feedback loop mechanism. The intermittent transition which leads to amplitude modulation of noise generation is analyzed in section 3.6 in terms of phase interference of the dominant frequencies from hydrodynamic instabilities. Then, the main findings and conclusions are presented in section 4.

Large eddy simulation
Large eddy simulation is performed to solve the compressible Navier-Stokes equations in general curvilinear coordinates. The spatial discretization of the governing equations employs a sixth-order accurate compact scheme for derivatives and interpolations on a staggered grid (Nagarajan et al. 2003). The time integration is carried out by a hybrid implicit-explicit method. The implicit second-order scheme of Beam & Warming (1978) is applied in the near-wall region to overcome the stiffness problem due to a fine boundary layer grid, while a third-order Runge-Kutta scheme is used for time advancement of the equations in flow regions away from the solid boundary. For the communication across the different methods, information is exchanged at overlapping points. An explicit sub-grid scale model is not applied. However, outside the boundary layer, a sixth-order compact filter (Lele 1992) is applied to control high-wavenumber numerical instabilities arising from grid stretching and interpolation between staggered grids. The transfer function associated with such filters has been shown to provide an approximation to sub-grid scale models (Mathew et al. 2003).
No-slip adiabatic wall boundary conditions are enforced along the airfoil surface and characteristic plus sponge boundary conditions are applied in the far-field locations to minimize wave reflections (Wolf 2011). Periodic boundary conditions are used in the spanwise direction. Length scales, velocity components, density, pressure and temperature are non-dimensionalized as , respectively. Here, L * is the airfoil chord, a * ∞ is the freestream speed of sound, ρ * ∞ is the freestream density, T * ∞ is the freestream temperature and γ is the ratio of specific heats. Variables with superscript * are given in dimensional units. Herein, time and frequency (Strouhal number) are presented non-dimensionalized by freestream velocity as t = t * U * ∞ /L * and St = f * L * /U * ∞ , respectively. The present numerical methodology has been extensively validated for various 2D and 3D simulations of compressible airfoil flows at different configurations (Wolf et al. 2012a(Wolf et al. ,b, 2013Ramos et al. 2019).
The O-grid employed for the current LES is shown in gray lines for every 3 points in figure 1(a). For a smooth O-grid generation, the original airfoil trailing edge is truncated at 98% of the chord and it has a curvature radius r/L * = 0.4%. The reference length scale L * is the unit chord from the original NACA0012 airfoil. The leading edge is placed at (x, y) = (0, 0) and the airfoil is pivoted about this point. The total wingspan has 0.4L * , being larger compared to previously reported studies of airfoil flows for this range of Reynolds number (Jones et al. 2008;Ducoin et al. 2016). The domain extends 37 chord-lengths outwards and a circular sponge, given by A[(r − r 0 )/∆r] n , occupies the last ∆r = 10 chords. In this function, r is the radial distance, r 0 is the sponge starting position, and its coefficients are A = 20 and n = 4.
The mesh has a distribution of points in the streamwise, wall-normal and spanwise directions given by n x = 660, n y = 600 and n z = 192, respectively, which results in approximately 76×10 6 grid points. The ratio of grid points along the suction side relative to the pressure side is approximately 5:3. The wall-normal distance of the first grid point is ∆n = 0.0001 and the stretching ratio is 1.5%. As discussed by Desquesnes et al. (2007), Sandberg (2011) andFosas de Pando et al. (2014b), the important mechanisms for tonal noise generation arise from two-dimensional flow instabilities in the laminar region of the flow. Motivated by this observation, a grid refinement study in terms of mean and fluctuation properties was conducted for 2D simulations and shown by the present authors in Ricciardi et al. (2021). However, in the present study, the flow  is three-dimensional and transitions to turbulence on the airfoil suction side, near the trailing edge. Considering only the turbulent flow region, the estimated mesh resolution in terms of wall units is given by ∆x + < 10, ∆y + ≈ 0.3 and ∆z + < 5. These values follow the best practices for wall-resolved LES. The time step is ∆t = 1.5 × 10 −4 and 75 convective time units are employed for post-processing and analysis of results. The 3D simulation starts from a 2D flow superposed with random noise and more than 30 convective time units are discarded before collecting flow statistics.

Linear stability and resolvent analyses
The mechanisms of generation and amplification of flow instabilities can be explained by linear stability theory. Hence, bi-global linear stability and resolvent analyses are employed in this work to investigate the most unstable frequencies and their sensitivity to disturbances, besides their role in the tonal noise generation. Through a Reynolds decomposition, it is possible to split the unsteady flow q(x, t) in a mean base flowq(x) plus a time-dependent fluctuation component q (x, t). If the fluctuations are sufficiently small, the Navier-Stokes equations can be linearized about the (mean) base flow. Although it is possible to consider the turbulent mean flowq as the base flow, the linear stability analysis would not hold because such state is not an equilibrium point of the NS equations. Nonetheless, the use of a time-averaged base flow may provide some insights as a model . With addition of a time-dependent external forcing f (x, t), the linearized Navier-Stokes equations (LNS) can be written as is an operator which serves as a spatial window where forcing is applied. The matrix L = L(q) is the bi-global LNS operator with the spanwise time-averaged flow propertiesq(x) as base flow. The evolution of linear disturbances by equation 2.1 can be performed by a direct analysis of the operator L or time-integrating the disturbances in a linearized version of the CFD code. In the former case, with a transformation where the wavenumber β ∈ R and the frequency ω ∈ C, it is possible to write the forced LNS equations in discrete form as In this case, the linear operator becomes a function of the spanwise wavenumber L = L(q, β) and B is the discrete version of B. If forcingf is absent, the LNS equations can be analyzed separately for each wavenumber β as an eigenvalue problem (linear stability analysis) as Here, Q holds the eigenvectors of L and the eigenvalues appear in the diagonal matrix Λ = −iωI, where the frequency and growth rate are the real and imaginary parts of ω, respectively. Due to the velocity gradients appearing off-diagonal, the linear operator becomes a non-normal matrix such that L H L = LL H . As a consequence, the eigenvectors are also non-normal. In this equation, the superscript H denotes the complex conjugate transpose (Hermitian). The matrix L H is the adjoint operator, being related to regions of sensitivity within the flow. Its eigenvalue decomposition yields The eigenvalues of the adjoint operator are the complex conjugate of those from direct analysis. The eigenvectors Q † represent the region of flow sensitivity and they are equal to Q −1 for normal systems. On the other hand, this is not true for non-normal systems. A thorough review on the significance of adjoint operators is presented by Luchini & Bottaro (2014) and references therein. In case forcing is applied at a frequencyω, equation 2.1 yieldŝ where the matrix H = H(q, β,ω) is the resolvent operator (Reddy & Henningson 1993;Schmid 2007;McKeon & Sharma 2010;Taira et al. 2017Taira et al. , 2020. For non-normal operators, eigenvalue sensitivity and energy amplification are related to the induced L 2 norm of operator H, i.e., the leading singular value σ obtained from singular value decomposition In this equation, V and U are unitary matrices holding right and left singular vectors and Σ is a diagonal matrix containing the singular values σ, such that The first column of V contains the forcing term that produces the largest response in the flow (first column in matrix U) with the amplification ratio given by Σ.
A transformation using matrix W is applied to convert variablesq andf with an appropriate energy norm prior to performing the SVD. In the case of compressible flows, the Chu norm which relates density, velocity and temperature (Chu 1965) is typically used. A spatial window C = C(x) is applied to limit the domain of analysis and the system response in Fourier domainŷ is given bŷ (2.9) Combining equations 2.6 and 2.9 leads to the modified resolvent operator (2.10) The amplification mechanisms of flow disturbances can be identified by using the eigenvalue decomposition of the linearized NS operator, equation 2.4, in the resolvent, yielding where Λ and Q are the eigenvalues and eigenvectors from the solution of equation 2.4. The bounds of the induced L 2 norm are where the lower bound is the case of an operator with orthonormal eigenvectors. In this case, the norm depends only of the resonances, following a 1 /R decay based on the distance R in complex plane with respect to the eigenvalues. For non-normal systems, the norm also depends on the pseudoresonances, measured by the product of eigenvectors Q and their inverse Q −1 . The weighting W and the two windowing matrices B and C are included in the norm calculation (see Schmid & Henningson (2001); Symon et al. (2018) for more details).
To perform the global stability (spectral) and resolvent (pseudospectral) analyses, the discrete linear operator is computed using the second-order accurate methodology from Sun et al. (2017) and . The base flow is the time-spanwised averaged solution from the LES,q = [ρ,ū,v,w,T ]. For the far-field and airfoil surface, Dirichlet boundary conditions are set for [ρ , u , v , w ] = [0, 0, 0, 0] and a Neumann boundary condition is set for T such that the wall-normal derivative ∂T /∂n = 0. At the outlet boundary, the same Neumann boundary condition is set for all flow variables. With these boundary conditions and the mean base flowq, the linear operator is computed in its discrete form L(q, β) for a prescribed spanwise wavenumber β.
The mean flowq obtained is interpolated from the O-grid in figure 1(a) to the H-mesh shown in figure 1(b) for the modal analysis. This is necessary to improve spatial accuracy of the linear operator downstream and upstream of the airfoil, where direct and adjoint eigenvectors are supported, respectively . The present two-dimensional H-mesh has an extent of x ∈ [−5, +7], y ∈ [−5, +5] and it is composed of approximately 520 × 10 3 grid points. The mesh refinement targets a Strouhal number cut-off St = 8.0. The sponge region comprises a circle with parabolic growth f = a(r −r 0 ) 2 for r > r 0 . The sponge is centered at (x, y) = (0.5, 0) with a = 1.0 and r 0 = 1.5. The mesh refinement and sponge placement are assessed by convergence of physically meaningful eigenvalues while suppressing spurious eigenvalues. If the mesh is coarse or the sponge is placed too far from the airfoil, it was observed that the meaningful eigenvectors also exhibit spatial support in the same region of the spurious eigenvectors.
For the pseudospectrum calculation, the dashed red rectangle in figure 1(b) depicts the region considered for the present resolvent analysis in terms of operators B and C in equation 2.10. This windowing is important not only from a physical point of view, where the noise sources close to the airfoil surface are most important, but it also eliminates spurious eigenvectors, which are mostly observed far away from the body. Finally, the SVD in equation 2.10 is performed with the randomized algorithm presented by Ribeiro et al. (2020).

Mean flow perturbation
The properties of the resolvent operator can be studied as an initial value problem where large transient energy amplification is observed, even in stable problems, due to the non-normality (Trefethen et al. 1993;McKeon & Sharma 2010). The time evolution of perturbations with respect to the mean flow can be performed with a linearized version of the CFD code according to equation 2.1. However, here we employ an approach that uses the present LES code with the addition of a force term to the right-hand side of the equations (Touber & Sandham 2009;Bhaumik et al. 2015;Ranjan et al. 2018). In this approach, known as mean flow perturbation, the Reynolds decomposition q =q + q is applied to the non-linear NS equations which are expressed as The turbulent mean flowq is not an equilibrium state of the NS equations and its time derivative ∂ tq is non-null. Hence, in order to use it as a base flow, an additional term R must be added in the numerical procedure to keep the base flow stationary. A Taylor series expansion about the base stateq yields (2.14) Grouping together the base flow terms results in where the base flow is stationary when Considering that the imposed perturbation amplitude q is sufficiently small leads to a linearized solution of the NS equations as where the Jacobian ∂N (q) ∂q is the linearized operator L. In contrast to the previous approach in section 2.2, this methodology is matrix-free and the linear operator is never explicitly computed. Despite this, the method still accounts for the modal interaction that leads to transient energy amplification within a non-normal analysis of the linear operator. The system response is obtained by the time evolution of the disturbances and the initial condition is set as an impulsive excitation that triggers the dominant response. In the current analysis, the grid and numerical methodology employed are the same used in the LES.

Results
The physical mechanisms responsible for the generation of secondary tones in airfoil flows are investigated by post-processing the non-linear results from LES and also employing linear stability theory. The linear analysis provides insights on the onset and growth of disturbances in the flow besides its receptivity. On the other hand, LES results allow the investigation of transition, intermittency and phase interference effects in the context of the multiple secondary tones and the acoustic feedback loop.

Mean flow
The mean flow is presented with contours of u-velocity normalized by the freestream speed of sound in figure 2. Important flow features can be observed including a long recirculation bubble that extends over the suction side, from the detachment at x = 0.34 until its reattachment at x = 0.82. In the magnified view, a small bubble can be observed on the pressure side, followed by a recirculation region at the trailing edge. The blue contours enclosed by magenta lines indicate regions of reversed flow.
Distributions of root-mean-square (RMS) for kinetic energy k and pressure p are presented in figures 3(a) and (b), respectively. As it can be seen from the plots, fluctuations start amplifying along the bubble and the highest values appear at the reattachment region, on the suction side. The green and blue dashed lines in figures 3(a) and (b), respectively, depict the locations of maximum k and p along the shear layer forming on the suction side. In the following analyses, these lines will be used as reference locations for data extraction to track flow disturbances. Again, the solid magenta line shows the recirculation bubble.
The negative mean pressure coefficient −C p is shown in figure 4(a) with blue and red solid lines for suction and pressure sides, respectively. To highlight the regions with intense fluctuations, the pressure coefficient is presented as −C p ∓ C pRMS , where the blue dashed and red dotted lines correspond to the suction and pressure sides, respectively. The peak value of pressure coefficient on the suction side is observed near the leading edge at x = 0.02. From this point onward, a pressure increase leads to a drop in −C p until reaching a plateau along the laminar separation bubble. The plateau extends up to x = 0.7, where the pressure further increases towards the trailing edge, causing another drop in −C p . Pressure fluctuations are shown in figure 4(d) for both suction and pressure sides. The RMS values of C p are low along the entire pressure side, and upstream of the bubble (x 0.5) on the suction side. However, they increase towards the peak at x ≈ 0.78 on the suction side. The mean skin friction C f distribution over the airfoil is presented in figure 4(b). The RMS values of skin friction C fRMS , shown in figure 4(e), are also added to the mean values. A magnified view highlights the presence of the laminar separation bubble which starts at x = 0.34 and reattaches at x = 0.67. The reattachment is followed by a second detachment which reaches a minimum at x = 0.75 and reattaches at x = 0.82. The nature of the double detachment profile is discussed in details by Duck et al. (2000). The RMS peaks at x ≈ 0.76 due to the high value of k shown in figure 3(a).
Mean tangential velocity profiles extracted along the wall normal direction ∆n are presented in figure 4(c) at different locations along the chord. At the detachment point, x = 0.34, the red curve indicates the zero wall-normal derivative in the velocity profile. At x = 0.73, the velocity profile shown as a grey line exhibits the maximum reversed flow of −0.13U ∞ . The orange line shows a velocity profile at x = 0.95, where threedimensional effects are important due to turbulent transition. In figure 4(f), RMS values of the tangential velocity are presented and a single peak is observed for the profile computed at x = 0.34. Velocity fluctuations increase inside the bubble at x = 0.64 and become dominant at x = 0.73. In both these locations, the velocity fluctuations present triple peak profiles as also observed by Nash et al. (1999), Desquesnes et al. (2007) and Duck et al. (2000). Near the trailing edge at x = 0.95, where a turbulent regime may occur, mixing results in a smoother velocity fluctuation profile.

Vortex dynamics
Instantaneous three-dimensional flow fields over the airfoil are presented in figure 5 with iso-surfaces of Q-criterion colored by u-velocity. A magenta shade along the airfoil surface depicts the region of reversed flow along the separation bubble. This figure shows two-dimensional rolls that are visible along the separation bubble and amplify on the suction side, leading to vortex shedding. Moreover, different flow regimes are observed at the trailing edge, where coherent structures alternates with periods of turbulent packets. For example, figure 5(a) shows a single laminar-like roll reaching the trailing edge at t = 20.82 while figure 5(b) shows that, at t = 21.78, vortex breakdown leads to a turbulent regime near the trailing edge. Readers are referred to the movie of the flow field provided as supplemental material (movie 1) for detailed inspection of the present flow dynamics.
The analysis of flow snapshots is important to understand the vortex dynamics over the airfoil suction side. Based on the acoustic analogy of Ffowcs-Williams & Hall (1970), efficient sound generation is achieved when pressure fluctuations are aligned with an edge, such that it will be maximum for 2D-like perturbations. In this sense, figure 6 shows the spanwise averaged ω z vorticity, where blue and red contours represent negative and positive values, respectively. A dark blue color indicates strong spanwise coherence of a 2D vortex, while light blue contours are representative of uncorrelated turbulence. A thin black line is used to identify individual vorticity packets, based on the iso-contour of an entropy measure given by p /ρ γ − ( p∞ /ρ γ ∞ ) = 0.1%. This separates the region where vorticity is important from where the flow is irrotational, outside the boundary layer. A magenta dashed line represents the border of averaged negativeū velocity, as discussed in figure 2. The simulation time is shown in convective time units on the lower left corner of all plots. A movie of figure 6 is provided as supplemental material (movie 2) to aid the dynamic visualization of the flow features.
As shown in figure 6, the flow dynamics is dominated by events on the suction side, where vortices are shed from the laminar separation bubble. During the shedding process, there is a possibility that the 2D laminar vortices undergo a process of vortex pairing. In case this process is not successful, vortex bursting leads to turbulent transition and low coherence along the span, represented by dotted lines. These turbulent packets may interact with other vortices and further reduce their coherence, as indicated by the dashed-dotted line. On the other hand, if the vortices merge or only a solitary vortex is shed from the bubble, the trend is for either structure to keep its high coherence up to the trailing edge. Both of these processes are represented by the dashed (merging) and solid (single vortex) lines in figure 6. A detailed discussion on the consequence and cause of this behavior is presented in sections 3.3 and 3.6, respectively. Similar dynamics of vortex shedding from LSBs are observed experimentally by Kurelek et al. (2016).
As shown in the experimental investigations of Pröbsting & Yarusevych (2015) and , for this Reynolds number and angle of attack, it is expected that the pressure side does not play an important role in terms of the flow dynamics. In order to confirm this observation, a temporal signal of pressure fluctuation is extracted at x = 0.98, marked by the blue dashed line in figure 3(b). A second signal is obtained for the velocity derivative in the wall-normal direction on the pressure side at x = 0.95. Both signals are presented in figure 7(a) as well as the time averaged dū/dn, indicating that the flow is, on average, recirculating (dū/dn < 0). In the velocity derivatives, the wall-normal direction is positive pointing outwards from the airfoil surface. Snapshots of the flow dynamics are presented in 7(b), where the red-blue contours are the ω z vorticity, similarly to the previous figures, the green hatched colors show instantaneous regions of negative u-velocity while the dashed magenta line indicates the boundaries of the timeaveraged pressure side bubble. In this figure, it is possible to see that when a vortex arrive at the trailing edge from the suction side at t = 20.59 (negative p peak), the flow reattaches on the bottom side (positive du/dn at the wall). As the vortex leaves the airfoil surface on the suction side at t = 20.73, the pressure increases on the suction side (positive p ) and, as a consequence, the flow detaches on the pressure side (negative du/dn). Thus, the magenta dashed line in figure 7(a) shows that the recirculation bubble exists on an time-averaged sense. Moreover, these figures show that the pressure side dynamics seems to be a reaction to the suction side events in the present flow.

Intermittency and sound generation
The acoustic analogy of Ffowcs-Williams & Hall (1970) states that efficient acoustic scattering occurs due to flow fluctuations near a trailing edge. In this sense, the hydrodynamic nearfield with pressure fluctuations is presented in figures 8(a,c) with black-white contours, while the acoustic field is presented in figures 8(b,d). In the latter, the levels are presented 10 times lower compared to the former. The figures also display ω z vorticity in red-blue colors. To highlight the spanwise coherence, regions of |ω z | > 24 are plotted as filled contours. The vortex cores are related to lower pressures (white background contours) while the gaps between the cores display higher pressures (black contours). A coherent vortex at the trailing edge is presented in figure 8(a) at a time instant t = 20.82. The highly coherent hydrodynamic structure leads to efficient acoustic scattering at the trailing edge which, in turn, induces the emission of an intense acoustic wave. The wave is scattered with a π phase opposition relative to the incident fluctuations and it reaches    Temporal signals are acquired at different positions along the blue dashed line in figure 3(b) in terms of spanwise averaged non-dimensional pressure. The signal at the trailing edge is presented in figure 9(a) and displays intermittent deep valleys related to the passage of coherent structures as shown in figure 8(a). The signal also shows low amplitude oscillations related to the turbulent packets that carry a negative pressure envelope, as described in figure 8(b). Such behavior can be characterized as an amplitude modulation of the signal. The same trends are also observed in the acoustic pressure signal, which is extracted one chord above the trailing edge, at (x, y) = (1c, 1c), and it is presented in figure 9(c). It is possible to see that the low pressure valleys from figure 9(a) result in high pressure peaks that reach the observer position at a retarded time.
The magnitude of the Fourier coefficients is presented in figures 9(b,d) for the hydrodynamic and acoustic fields, respectively, displaying equidistant multiple tones. The Fourier transforms employ a frequency domain averaging with an overlap of 67% amongst 4 bins covering the entire period. The frequency resolution is ∆f = 0.03 and a Hanning window function with energy correction is applied. All tonal peaks are integer multiples of the lowest frequency tone, St ≈ 0.5. As it can be seen in figure 9(b), the dominant peak at the trailing edge is at St ≈ 3.5 for probes located close to the leading edge (x = 0.05), at the detachment point (x = 0.34) and the trailing edge (x = 0.98). This frequency is related to the passage of low pressure packets, either turbulent or coherent. At x = 0.64, in the shear-layer roll-up region, the dominant peak moves to a higher frequency, at St ≈ 4.5, which may be related to the smaller scale instabilities yet to be merged into larger vortices. In the acoustic field, figure 9(d) presents results extracted at 1 and 5 chords above the trailing edge. The same frequencies are observed where St = 3.5 is still the most prominent peak.
The identification of time instants when coherent vortical structures pass by the trailing edge is performed using the two-point, one-time autocovariance of pressure fluctuations along the spanwise direction ∆z at each time t according to R(x, y, ∆z, t) = p (x, y, z, t) p (x, y, z + ∆z, t) . (3.1) Calculations are performed near the trailing edge at x = 0.98 and results are presented in figure 10. The dark lines display higher spanwise coherence which are almost constant and independent of ∆z. These corresponds to the two-dimensional vortex rolls illustrated in figure 5(a). On the other hand, in moments where the flow transitions and only small scale eddies are observed at the trailing edge, such as shown in figure 5(b), the spanwise coherence presents very low levels. Furthermore, in these time instants, the pressure fluctuations are considerably lower compared to the low pressure-core of the vortices, which are responsible for large pressure fluctuations. Hence, the products of small p in equation 3.1 results in even lower values of autocovariance. As it can be seen from figure 10, coherent structures often reach the trailing edge with a slow time-scale of ∆t ≈ 2.0 convective time units, which corresponds to a frequency St = 0.5. However, due to the chaotic motion, this process is not perfectly periodic and multiple vortices may reach the trailing edge within this period, as indicated by the time differences ∆t = 0.22, 0.25, 0.29, 0.33 and 0.40 in the figure. In this faster time-scale, the different time intervals are related to frequencies St = 4.5, 4.0, 3.5, 3.0 and 2.5, respectively, which are integer multiples of the low frequency tone at St = 0.5.
To better highlight this behavior, a time-frequency analysis is also performed and results are presented in figure 11. The top row shows the continuous wavelet transform, CWT, using the complex Morlet wavelet (Farge 1992). This function is the product of a monochromatic complex exponential with a Gaussian envelope. Its spectral response is an exponentially decaying band-pass filter centered around the frequency of the complex exponential. The filter width depends on the Gaussian standard deviation. The scalogram presents how the correlation of the wavelet with a reference signal changes in time based on the wavelet frequency. Hence, higher values indicate strong fluctuations of the signal at a particular frequency and time. Results are obtained for the same probe locations considered in figure 9(b). The scalograms from figures 11(a) and (c) exhibit longer periods of intense pressure fluctuations at St = 3.5, while in figure 11(b), a broader range of frequencies is excited more intensely. As it can be seen from the probe at x = 0.64, frequencies higher than St = 3.5 are excited, particularly St = 4.5, 5.0 and 5.5. One should note that the color levels are different for each location, with lower and higher amplitudes for probes at x = 0.34 and 0.98, respectively. These figures show that the secondary tones present intermittent behavior, with each frequency experiencing periods of silence (blue contours) and loud noise emission (green, yellow and red contours). Another important aspect captured by the scalograms is the amplitude modulation of The bottom row shows the two-point, one-time pressure autocovariance along the span, computed for each time step over the entire simulation period. It is possible to see that the dark lines in the covariance plots are associated to the excitation of the secondary tones in the wavelet scalogram, indicating that specific frequencies are related to advection of quasi-2D flow structures. The coherence level and the time interval between successive vortices impacts not only the instantaneous magnitude of the wavelet transform but also may change its dominant frequency.

Linear analysis and amplification mechanisms
In order to investigate the growth of velocity and pressure fluctuations at the tonal frequencies, figure 12 presents Fourier transforms of the temporal signals extracted along the airfoil suction side, following the green (velocities) and blue (pressure) dashed lines from figures 3(a) and (b), respectively. In this analysis, the temporal signals are computed for the spanwise averaged normal and tangential velocity components u n and u t , respectively, besides pressure p. Results are plotted together with the RMS values (solid black line) extracted from figure 3. Bothû t andp show strong amplification downstream of x = 0.5, which coincides with the location of the recirculation bubble. Furthermore, a plateau is observed upstream of x = 0.5 for pressure fluctuations, which indicates that acoustic waves dominate over hydrodynamic disturbances along this region. Pressure fluctuations at St ≈ 3.5 exhibit higher amplitudes along the entire chord and similar observations can be made for the velocity components. As discussed by Kurelek et al. (2019), it is possible that these fluctuations act as a sub-harmonic forcing that stimulates the vortex merging process.
The mechanisms of generation and amplification of flow instabilities is explained by linear stability theory. In this regard, bi-global linear stability and resolvent analyses are employed to understand the most sensitive frequencies in terms of optimal forcing and response characteristics for 2D perturbations, i.e., only for spanwise wavenumber β = 0. This is justified since 2D structures are expected to radiate noise more efficiently (Sano et al. 2019). In the linearized NS operator, the off-diagonal terms, related to velocity gradients and shear, lead to a non-normal operator. Hence, orthogonality of the stability  modes is not expected and interaction among modes within the system may lead to transient response that affect its short-term dynamics (Schmid 2007). This behavior is properly captured by the resolvent analysis, which provides the information on eigenvalue sensitivity and transient energy amplification rate. Following the theoretical and numerical methodology described in section 2.2, the eigenvalue spectrum from a bi-global stability analysis is presented in figure 13. The black circles represent the physically meaningful eigenvalues while the empty grey symbols are spurious eigenvalues. The latter could be identified due to their large displacement in the complex plane during the grid/sponge convergence study (not shown) as well as by visualization of the associated eigenvectors (modes). Sensitivity of the operator to external disturbances is computed by the resolvent analysis along both real (frequency) and imaginary (growth rate) axes and it is measured by the pseudospectrum, presented as the contour plot in figure 13. The region of analysis in the resolvent operator is limited by means of matrices B and C in equation 2.10. Only the near-field grid points are considered, i.e., those inside the dashed red rectangle in figure 1. This is physically justified since the most intense quadrupole noise sources which are responsible for the incident field in the acoustic scattering process are close to the trailing edge (Wolf et al. 2012a). As an additional positive side effect, this improved results by reducing spatial support of spurious eigenvectors.
Results from the bi-global stability analysis in terms of the eigenspectrum are presented in figure 13. It is possible to see multiple eigenvalues at frequencies close to those from the LES, previously discussed in figures 9(b,d). For instance, the maximum relative deviation is less than 5%, at St ≈ 2.0, and around 2% for the most unstable frequency St ≈ 5.0. Thus, the vortex dynamics and the multiple tones seem to be triggered by linear mechanisms. The pseudospectrum, obtained from the resolvent analysis, is shown as a contour plot in logarithmic scale of the leading singular value from the SVD according to equation 2.7. The input-output amplification ratio depends on both resonance and pseudoresonance (see equation 2.12) such that peaks appear when approaching an eigenvalue. The pseudoresonance levels depend on the mutual excitation of eigenvectors with more sensitive eigenvalues presenting higher peaks. In this sense, it is possible to see in figure 13 that higher frequencies are more sensitive compared to the lower ones. Moreover, the stable eigenvalues may impact the flow dynamics in terms of transient energy amplification given their sensitivities to disturbances.
The most unstable frequencies of St ≈ 4.5, 5.0 and 5.5 are observed in the LES on the shear-layer roll-up region at x = 0.64, as illustrated by figure 11(b). The interaction among multiple frequencies of the flow instabilities and the vortex merging reduce the shedding frequency in the LES, where the dominant tonal peak becomes St ≈ 3.5. This non-linear phenomenon is captured by the LES but not by the linear analysis. A more detailed discussion regarding the necessary conditions for vortex merging is presented in section 3.6.
The modal shape of each frequency is crucial to understand the system dynamics and to identify amplification regions within the flow. In this sense, figures 14(a), (b) and (c) present the real part of the global stability modes for the most unstable frequency St = 5.08 in terms of u , v and p , respectively. The modes are normalized with respect to the maximum value for each variable and the magnitude doubles every level to better represent lower values. The modal shapes computed for other relevant frequencies in the spectrum (resonances) show similar patterns compared to the one described here and differences are only related to their characteristic wavelengths. Along the separation bubble, the modal structures of u and v are tilted against the mean shear, suggesting an Orr mechanism for transient growth of energy and amplification of disturbances (Schmid & Henningson 2001). After flow reattachment, the structures become aligned with the wall-normal direction. The pressure disturbances p display upstream traveling acoustic waves radiated from the trailing edge due to scattering mechanism. The modes displayed in the figure originate as shear layer instabilities on the airfoil suction side and extends along the wake. Hence, wake-boundary layer coupling mechanisms cannot be discarded.
The dominant resolvent modes are computed along the neutral stability axis (ω i = 0.0) for St = 5.08 in terms of kinetic energy k and presented in figure 15 normalized with respect to their maximum value. The response mode, shown in white-purple contours, depicts the location where disturbances are more energetically relevant. Here, its modal shape magnitude is presented such that it doubles every 2 levels. It is possible to see negligible levels from the leading-edge until the recirculation bubble. Then, for x > 0.4, its magnitude increases quickly, indicating that the suction side bubble is an amplifier of disturbances. The forcing mode, upper-left plot with white-blue contours, highlights the most receptive location, where minimal disturbances result in larger transient growth by the system. This mode is related to the adjoint operator of the linearized NS equations and it is introduced in the context of airfoil secondary tones by Fosas de Pando et al. (2017). In the present results, the leading edge region on the suction side shows the higher values of forcing modes. Furthermore, these modes are absent in the pressure side which corroborates to the absence of pressure side-driven events for this flow configuration.

Acoustic feedback loop mechanism
As discussed by several authors (Tam 1974;Lowson et al. 1994;Desquesnes et al. 2007;Jones & Sandberg 2011;Fosas de Pando et al. 2014b), airfoil flows at transitional Reynolds numbers are subject to an acoustic feedback loop mechanism that is related to the multiple tones. In this feedback process, hydrodynamic perturbations are excited in the vicinity of the leading edge and convected downstream. The disturbances amplify along the recirculation bubble on the airfoil surface and, upon reaching the trailing edge, scatter acoustic waves that travel upstream and trigger new disturbances. In this context, figures 16(a-h) present results from the resolvent analysis for different frequencies. The real part of the pressure forcing modes is shown in red-blue contours while black-white contours display the real component of the leading pressure response modes. The modes have been normalized with respect to their maximum values. To highlight the most receptive locations for each individual frequency, the yellow spots in the forcing modes are shown only for regions where the magnitude reaches higher values than 98% of the maximum. As presented in figure 16, the locations of maximum receptivity change with frequency and vary along the interval 0.10 < x < 0.18.
Aspects such as the amplitude and phase between the acoustic waves and flow receptivity at different frequencies impact the onset of hydrodynamic disturbances. To understand such complex interrelation, the acoustic feedback is analyzed using the mean flow perturbation technique to solve the linearized NS equations around the 3D mean flow in time-domain (Jones & Sandberg 2011;Fosas de Pando et al. 2014b). The method assumes that the base flow is steady and that it is not modified by the linear perturbations. The dynamics of the pressure fluctuations is presented in figure 17. For the initial condition, a pressure disturbance is applied at the leading edge, from x = 0.0 to 0.005 over a short duration of ∆t = 1 × 10 −4 to excite all frequencies (impulse response). This triggers a wavepacket that is advected by the mean flow along the suction side. Figures 17(a,b) present the growth of disturbances for x > 0.4, where spatial amplification described by the response modes becomes important, as shown in figure 15. When the wavepacket reaches the trailing edge, acoustic waves are generated and travel upstream as shown in figure 17(c). After the wavepacket reaches the wake, hydrodynamic disturbances are no longer observed on the suction side, as displayed by figure 17(d). Afterwards, disturbances in the wake slowly decay while a new wavepacket amplifies on the suction side (figures 17(e,f)). Similarly to the low frequency observed in the LES calculation, the wavepackets have a period of ∆t = 2.0 and the next cycle is presented in figures 17(g-j). The combined response from the suction side wavepacket, wake instability and upstream traveling acoustic waves keeps the fluctuation field active, amplifying the globally unstable frequencies in the eigenspectrum of figure 13, which ultimately become dominant. Readers are referred to the movie provided as supplemental material (movie 4) presenting the evolution of the linearized flow dynamics. The feedback loop mechanism is investigated here by analyzing the generation of flow instabilities by upstream propagating acoustic waves, as presented in figure 18. In this p : figure, the entropy measure p /ρ γ is given by red-blue contours while pressure fluctuations are given by blue-green-red lines, where solid and dashed lines are positive and negative values, respectively. The entropy measure is used to partially filter out fluctuations associated to acoustic waves while pressure is used to highlight acoustic waves. In figure  18(a), it is possible to see an incoming acoustic wave at t = 5.888 that, upon reaching the leading edge at t = 5.920, triggers an entropy fluctuation due to the secondary diffraction. This latter phenomenon is typical of airfoil noise as discussed by Miotto et al. (2017). The entropy fluctuation at the leading edge, presented in red color and highlighted by a dotted line box, is propagated downstream in figures 18(b,c,d) at t = 5.920, 5.952 and 5.984, respectively. This disturbance overlaps with an incoming entropy wave moving upstream, carried by an acoustic wave depicted with blue dashed lines. Then, in figures 18(e,f), respectively at t = 6.016 and 6.048, the downstream propagating entropy fluctuation encounters a second positive acoustic wave in the region where the forcing modes are maximum, i.e., for 0.10 < x < 0.18. This upstream propagating acoustic wave carries a positive entropy fluctuation which impinges on that generated at the leading edge. Thus, although the onset of instabilities occurs very close to the leading edge, the interaction between upstream traveling acoustic waves and downstream propagating boundary layer instabilities cannot be discarded since it occurs along the entire sensitive region of the flow. Albeit the analysis is conducted for pressure and entropy, local velocity fluctuations also occur and contribute to the incoming perturbations seen by the bubble. This overall process of the feedback loop mechanism is better visualized in movie 5 submitted as supplemental material.
Flow disturbances are tracked along the black dashed line in figure 18 (also shown as blue dashed line in figure 3(b)) and presented as a function of time and space in figures 19(a) and (b). The former shows pressure fluctuations p while the latter displays disturbances of an entropy measure ( p /ρ γ ) . In these plots, the values are normalized with respect to the maximum and presented in a logarithmic scale. The positive and negative values are presented as red and blue colors, respectively. The first wavepacket shown in figures 17(a,b) is observed traveling downstream with a local mean velocityū as bottom-up contours starting at t = 0.0. Upon reaching the trailing edge, at t > 1.2, the wavepacket generates acoustic waves propagating upstream withū − a speed. These waves are observed as top-down contours with lower magnitudes. The acoustic pressure has a phase opposition of π compared to the hydrodynamic fluctuations on the trailing edge and, hence, an incident negative hydrodynamic pressure fluctuation (blue) creates a scattered positive acoustic wave (red).
Upstream of x ≈ 0.4, the acoustic waves dominate and it is not possible to see the onset of instabilities. Hence, a measure of entropy ( p /ρ γ ) is shown in figure 19(b) to partially filter the upstream acoustic disturbances and highlight the hydrodynamic content. In this figure, the contours of hydrodynamic fluctuations from bottom-up are better visualized. In both plots, a green region covering the entire period is placed at 0.1 x 0.18, illustrating the maximum forcing region from figure 16. Along this region, disturbances generated in the vicinity of the leading-edge (x ≈ 0.0) travel downstream and appear to interact with the incoming acoustic waves, as depicted in figures 18(e,f). Such interaction results in a chess board-like interference pattern where the entropy levels carried by upstream propagating acoustic waves (inside the boundary layer) are of the same order of the entropy disturbances generated at the leading edge. We conjecture that the simultaneous combination of acoustic and hydrodynamic interactions at this location results in the maximum forcing. Finally, black lines are drawn on the plots using the convective velocity information for both pressure and entropy disturbances. These lines are used to track the feedback dynamics and show that, if the feedback indeed closes along the most receptive region, the period of ∆t ≈ 2.0 is satisfied.
The analysis of the feedback mechanism for the LES results is presented based on the RMS values of either the pressure p or the entropy measure p /ρ γ averaged along the zdirection (spanwise direction). Similarly to the linear analysis, time signals are extracted along the blue line in figure 3(b) and the values of R 0 are presented in a logarithmic scale in figure 20(a) for pressure and (b) for the entropy measure. In both cases, it is possible to see high fluctuation values for x > 0.4 traveling downstream towards the trailing edge. For x < 0.4, smaller values of p propagate upstream towards the leading edge. Similarly to the linear analysis, the former and latter disturbances represent hydrodynamic and acoustic waves, respectively. The entropy measure partially filters the acoustic waves moving upstream. Based on linear analysis results, the onset of disturbances occurs very close to the leading edge in a position that coincides with the small hump observed in figure 4 for pressure and friction coefficients in terms of RMS values, respectively C pRMS and C fRMS . Thus, the wave secondary diffraction at the leading edge excites boundary layer disturbances which are also observed in the LES. The entropy fluctuations generated at the leading edge appear to interact with the upstream-propagating acoustic waves in the region of maximum receptivity, similarly to the phenomenon observed in the linear analysis. In figure 20, black lines are used to track disturbances assuming that the feedback mechanism takes place at 0.1 < x < 0.18, indicated by the green shaded rectangle which represents the region of maximum receptivity from figure 16. As a final remark, the repeating patterns observed in the simulations demonstrate that acoustic feedback occurs in the present flow. Results indicate that the loop depends on the leading edge secondary diffraction and the region of maximum sensitivity, located downstream of the leading edge and upstream of the separation bubble. However, the intermittent interaction of vortex merging and bursting on the suction side may eventually disrupt the cycle by altering the phase and magnitude of the tonal frequencies with a temporary impact on the feedback loop.

On vortex merging and phase interference
From the previous sections, it is observed that the acoustic feedback loop is selfsustained by strong acoustic waves scattered from coherent flow structures advected past the trailing edge. These upstream-propagating acoustic waves create new flow instabilities which are amplified by the separation bubble and undergo a process of vortex pairing. Thus, it is important to understand which circumstances lead to generation of coherent structures or small-scale eddies, as discussed in section 3.3. For this task, the dynamics of wavepackets in the linear simulation is first investigated, followed by a similar analysis in terms of the LES data.
The amplification of the wavepacket along the recirculation bubble results in a maximum value of the perturbations evaluated by the MFP at x = 0.72. The temporal signal of p at this location is presented in figure 21(a) depicting three envelopes. These are evaluated using the Hilbert transform to obtain the imaginary part and thus the magnitude of a real signal. Then, using the continuous wavelet transform, it is possible to extract the phase information of the previous signal. To highlight the relevant periods of the dynamics, the phase is weighted by the magnitude of the CWT coefficients as presented in figure 21(b). Comparing the temporal signal with the weighted phase, it is possible to see that the wavepacket envelopes peak when there is phase alignment across the multiple frequencies of the eigenspectrum shown in figure 13. Hence, constructive interference of multiple modes is crucial for maximum energy amplification. When the waves are out of phase, the wavepacket magnitude drops.
The phase information from the LES data is presented together with the spanwise pressure autocovariance in figures 22(a) and (b) at x = 0.64. This location is represen- tative of the vortex pairing region. These figures show that the high coherence depends on the phase alignment of frequencies 4 St 6. This frequency range corresponds to that where the most amplified modes from the linear analysis are observed. If more (fewer) frequencies are aligned, a stronger (weaker) signal coherence is achieved. If the disturbances are in phase, the result is vortex merging and strong coherence. On the other hand, if they are out of phase, the vortex pairing fails and coherence is reduced. The repeating patterns of high coherence show that the feedback mechanism has a period of ∆t ≈ 2.0, related to the lower tone in the spectrum at St = 0.5.
Destructive and constructive interference of vortical structures during the pairing process can be interpreted as an amplitude modulation of the signal (Desquesnes et al. 2007;Ricciardi et al. 2020). The modulation mechanism from the 2D simulations performed in these references differs from the current 3D case, where turbulent transition occurs. Here, the destructive interference leads to vortex breakdown and smaller pressure fluctuations. On the other hand, constructive interference of multiple frequencies creates spikes in the pressure signal as shown for certain time instants in figures 9(a) and (c).

Conclusions
A study of trailing edge tonal noise is conducted through a combination of large eddy simulation and linear stability theory. The focus of this investigation is on the secondary tones and acoustic feedback loop which arise in airfoil flows at moderate Reynolds numbers. A compressible LES is performed for the flow over a NACA 0012 airfoil at α = 3 • immersed in a freestream flow with Mach number M ∞ = 0.3 and Reynolds number Re = 5 × 10 4 . Despite intermittent flow transition to turbulent regime, the noise spectrum depicts a main tone with multiple equidistant secondary tones, and all tonal peaks observed are integer multiples of the lowest frequency peak.
The mean flow displays two recirculation regions, one on each side of the airfoil, where the separation bubble on the suction side is longer than that on the pressure side. The flow dynamics on the airfoil is dominated by events on the suction side, aft of the separation bubble, as indicated by the RMS values of kinetic energy and pressure. The mean friction coefficient shows a double detachment pattern, typical of large fluctuations downstream the recirculation bubble. Mean velocity profiles show that the maximum reversed velocity inside the separation bubble is 13% of the freestream value.
Linear stability analysis of the mean flow using a bi-global operator shows marginally stable and unstable eigenvalues at frequencies close to those of the tonal peaks observed in the non-linear simulation, indicating that the overall flow dynamics originates from linear mechanisms. Results from resolvent analysis of the linear operator are displayed as the pseudospectrum, showing that the higher frequencies are more susceptible to be excited by pseudoresonances. The response modes show that the suction side recirculation bubble acts as an amplifier of disturbances. On the other hand, the forcing modes show a high sensitivity region within the flow, starting downstream of the leading edge and extending until upstream of the bubble.
In the large eddy simulation, the amplification of fluctuations along the laminar separation bubble results ultimately in vortex shedding on the suction side. The vortices undergo a pairing process that may lead to either merging or bursting. This process moves the dominant tonal peak in the LES to a lower frequency than that observed in the unstable eigenvalues of the linear spectrum. Another result of the vortex pairing is that either coherent vortices or turbulent spots are advected towards the trailing edge. The success of vortex merging aft of the bubble depends on the phase alignment between the unstable frequencies observed in the linear eigenspectrum. When such frequencies are in phase, it is likely that the vortices will successfully merge. On the other hand, destructive interference results in vortex bursting. The different regimes where coherent structures alternate with uncorrelated eddies at the trailing edge act as an amplitude modulation of the signal. The magnitude of the radiated acoustic waves depends on the coherence level of the vortical structures at the trailing edge, measured by two-point, onetime autocovariance of spanwise pressure fluctuations. Together with the autocovariance, a time-frequency analysis using continuous wavelet transforms show that the flow has two different time scales. One is related to the lowest tonal frequency, which drives the acoustic feedback loop and depends not only on the coherent structures but also on the uncorrelated turbulent packets. The other is a faster time scale related to the passage of multiple coherent structures at the trailing edge. It has been found that the intermittent dynamics of this faster time scale impacts the magnitude of tonal peaks with the possibility of changing the instantaneous dominant tonal frequency.
The feedback loop mechanism is investigated by the time integration of the linearized Navier-Stokes equations using the mean flow perturbation technique. The system response to an impulsive actuation at the leading edge displays a wavepacket advected towards the trailing edge with subsequent radiation of acoustic waves due to a scattering mechanism. The wavepacket is composed of multiple frequencies in the eigenspectrum and its magnitude is maximum when the phases of the unstable frequencies from the eigenspectrum are aligned. Visualization of the linearized solution shows that entropy disturbances are triggered by the secondary diffraction and propagate downstream along the boundary layer. These disturbances play an important role since they interact with upstream-propagating waves at the region of maximum forcing. Thus, it is likely that the feedback loop closes at the most sensitive region, described by the forcing modes, in order to satisfy the low frequency period of the wavepacket cycle.
The spanwise pressure covariance is analyzed along the airfoil chord and shows an influence of the feedback mechanism, where specific patterns repeat in time. Intermittency due to transition to turbulence may shift the phase and magnitude of all frequencies, disrupting the feedback loop. Despite this, the flow reaches a new equilibrium and reestablish the periodic merging/bursting process, indicating that the periodicity of multiple stability modes may be the self-sustaining mechanism for multiple tones and flow transition.