Three-dimensional fluid motion in Faraday waves: creation of vorticity and generation of two-dimensional turbulence

We study the generation of 2D turbulence in Faraday waves by investigating the creation of spatially periodic vortices in this system. Measurements which couple a diffusing light imaging technique and particle tracking algorithms allow the simultaneous observation of the three-dimensional fluid motion and of the temporal changes in the wave field topography. Quasi-standing waves are found to coexist with a spatially extended fluid transport. More specifically, the destruction of regular patterns of oscillons coincides with the emergence of a complex fluid motion whose statistics are similar to that of two-dimensional turbulence. We reveal that a lattice of oscillons generates vorticity at the oscillon scale in the horizontal flow. The interaction of these vortices explain how 2D turbulence is fueled by almost standing waves. Remarkably, the curvature of Lagrangian trajectories reveals a"footprint"of the forcing scale vortices in fully developed turbulence. 2D Navier-Stokes turbulence should be considered a source of disorder in Faraday waves. These findings also provide a new paradigm for vorticity creation in 2D flows.

Parametrically excited waves are a ubiquitous phenomenon observed in a variety of physical contexts. They span from Faraday waves on the water surface to spin waves in magnetics, electrostatic waves in plasma and second sound waves in liquid helium. Parametrically excited Faraday waves on the surface of vertically vibrated liquids quickly become nonlinear. In dissipative liquids, or in granular media, these nonlinear waves form regular lattices of oscillating solitons (oscillons), resembling in some aspects 2D crystals. If the vertical acceleration is increased, the oscillons do not solely grow in amplitude, their horizontal mobility is also greatly enhanced, and ultimately the lattice melts and becomes disordered. Until recently, the physics of these self-organized waves and their transition to disorder have been studied almost exclusively based on the analysis of the wave motion rather than the motion of their constitutive components, whether they are solid grains or fluid particles.
It has recently been discovered that the fluid motion on a liquid surface perturbed by Faraday waves reproduces in detail the statistics of two-dimensional turbulence. This unexpected discovery shifts the current paradigm of order to disorder transition in this system: instead of considering complex wave fields, or wave turbulence, it is conceivable that the 2D Navier-Stokes turbulence, generated by Faraday waves, feedbacks on the wave crystal and disorders it in a statistically predictable fashion. To date, the very mechanism behind the turbulence generation in such waves remains unknown. A better understanding of this phenomenon is important for a wide spectrum of physics applications involving parametric waves.
In this paper, we visualize 3D trajectories of floating tracers and reveal that the fluid particles motion injects 2D vortices into the horizontal flow. This is an unexpected and new paradigm for vorticity creation in a 2D flow. The horizontal energy is then spread over the broad range of scales by the turbulent inverse energy cascade. Two-dimensional turbulence destroys the geometrical order of the underlying lattice. The crystal order, however, can be restored by increasing * Nicolas.Francois@anu.edu.au viscous dissipation in the fluid which hinders vorticity creation and thus the development of turbulence.

Abstract:
We study the generation of 2D turbulence in Faraday waves by investigating the creation of spatially periodic vortices in this system. Measurements which couple a diffusing light imaging technique and particle tracking algorithms allow the simultaneous observation of the threedimensional fluid motion and of the temporal changes in the wave field topography.
Quasi-standing waves are found to coexist with a spatially extended fluid transport. More specifically, the destruction of regular patterns of oscillons coincides with the emergence of a complex fluid motion whose statistics are similar to that of two-dimensional turbulence. We reveal that a lattice of oscillons generates vorticity at the oscillon scale in the horizontal flow. The interaction of these vortices explain how 2D turbulence is fueled by almost standing waves. Remarkably, the curvature of Lagrangian trajectories reveals a "footprint" of the forcing scale vortices in fully developed turbulence. 2D Navier-Stokes turbulence should be considered a source of disorder in Faraday waves. These findings also provide a new paradigm for vorticity creation in 2D flows.
Parametric surface waves can self-organize into various motifs and they have been the focus of the pattern formation physics for many years [11][12][13]. Stable patterns, crystals, or even quasi-crystals are produced on the surface of dissipative liquids [14][15][16] or granular media [17,18]. At higher wave amplitudes such patterns break down into disordered lattices. Different pathways leading to disorder have been described such as: defectmediated turbulence, oscillatory transition phase, or the lattice melting scenario [18][19][20][21][22].
Until recently the physics of these parametrically excited waves, and their transition to disorder, have been studied almost exclusively based on the analysis of the wave motion rather than the motion of their constitutive solid grains or fluid particles. Nevertheless recent studies of the motion of floaters on the surface of Faraday waves (FW) revealed an unexpected physics picture: the particle motion, in these three-dimensional waves, reproduces in detail the statistics of two-dimensional turbulence. Both the inverse energy cascade and the spectral condensation have been found, suggesting that Navier-Stokes two-dimensional turbulence could be a source of disorder at work in low dissipation systems [23,24].
Although FW are quasi-standing waves, early studies already pointed out that they can generate a random fluid particle motion [25][26][27]. It was suggested that this Lagrangian motion can be modelled by an extension of the Stokes drift phenomenology for random wave field [26,28,29]. However, such a phenomenology shares almost no connection with the recent observation of FW driven two-dimensional turbulence [23,24,30]. Indeed this modelling approach describes poorly the presence of a spectrally localized horizontal forcing mechanism for 2D turbulence. Moreover, it contains no ingredient to capture the existence of an inertial range. Thus, what remains unknown is how and why waves produce 2D turbulence. In particular, how the energy is injected into horizontal flow from vertical oscillations and why it is injected in a narrow range of scales leaving it to turbulence to spread this energy over a broad inertial interval.
In this work, we address some of these shortcomings and describe previously unexplored Lagrangian features of particle trajectories in disordered and ordered Faraday waves. We use a combination of a high speed camera, a diffusing light imaging technique and tracking algorithms to reconstruct the three-dimensional (3D) trajectories of floating tracers as well as the topography of the wave field. This approach gives new insights on how an essentially vertical energy injection is converted into spatially extended horizontal fluid motion. We confirm a clear distinction between the wave and the fluid motion in disordered FW [24]. The waves motion is equivalent to that of an assembly of oscillons confined to the sites of a weakly irregular lattice. In this sense, even disordered FW remain close to standing waves. Conversely, particles motion is highly erratic and spatially extended. We show that disordered quasi-standing waves generate vorticity in the horizontal fluid motion. Those vortices are determined by a unique characteristic length scale, and are the fuel of the 2D turbulent flow. As demonstrated in this paper, a substantial increase in viscosity results in the recovery of a perfect standing wave crystal while the particle motion becomes restricted to closed orbits. The emerging physical model of Faraday ripples consists of a spatially periodic 2D array of oscillating solitons [17,32] which contain vertical oscillatory energy but also inject horizontal vortices into the flow at the scale of the oscillon size. This horizontal energy is then spread over the broad range of scales by the inverse energy cascade and at the same time 2D turbulence disorders the underlying lattice.

II. RESULTS
In these experiments, the three-dimensional motion of floating tracers in Faraday waves and the topography of the wave field have been measured simultaneously.
Faraday waves are formed in a circular container (178 mm diameter, 30 mm deep) filled with water. The liquid depth is larger than the wavelength of the perturbations at the surface (deep water approximation). The container is vertically vibrated by an electrodynamic shaker. The forcing is monochromatic and is set to either f 0 = 30 Hz or f 0 = 60 Hz. Beyond a certain vertical acceleration threshold (a > a th ), the surface of the liquid becomes unstable. This is the Faraday instability. The waves that appear on the fluid surface, are parametrically forced; the dominant frequency of the excited surface ripples is at the first subharmonic of the excitation frequency, f = f 0 /2. In the range of frequencies f 0 = (30 − 60) Hz, the Faraday wavelength is λ ≈ (8 − 14) mm. We define the supercriticality as = (a−a th )/a th , where a th is the threshold for the parametric generation of Faraday waves.
The liquid surface is seeded with floating tracers to visualize the horizontal fluid motion. The use of surfactant and plasma treatment ensures that tracers do not aggregate (see the Appendix).
We use a diffusing light imaging technique to measure the topography of the wave field. A few percent of milk added to water provides sufficient contrast to compute high resolution reconstruction of the parametrically ex-   Figure 1(a) shows the 3D reconstruction of a trajectory recorded for 1.6 second in a disordered Faraday wave field forced at f 0 = 60 Hz and a vertical acceleration of a = 1.6g ( ≈ 1.7). This trajectory is superimposed on a reconstruction of the wave field which is built upon two images corresponding to two successive maxima of the parametric excitation (i.e: the blue and pink wave fields are separated by a time interval dt = 1/f 0 = 0.016 s). Figure 1 illustrates qualitatively two important features of the particle motion in disordered FW: i) the particle motion can be extremely convoluted, often forming twisted patterns or even cusps which coincide with a strong orientational randomization of the trajectory ( Fig. 1(b)). These intricate loops are usually confined within a circle whose radius is ≈ λ/4. The lower the vertical acceleration, the longer the time a particle spends doing this complex motion. A pronounced cusp can be seen in Fig. 1(b). This is a singular point in the trajectory geometry which corresponds to a reversal of the particle velocity. ii) a particle can escape this trapping and it experiences intense velocity peaks as shown in Fig. 1(c). From the top view ( Fig. 1(b)), the particle trajectory during such event looks like a long ballistic flight. However, this ballistic behaviour is a characteristic of the horizontal motion only. Indeed, Fig. 1(a) clearly shows that the particle follows the complex z -topography of the wave field during such a flight. In the following, we refer to 'traps' and 'flights' to characterize these geometrical features observed in all particle trajectories in disordered FW [27]. We also note that the particle trajectory shows poor resemblance to the classical Stokes drift picture [28], where a trajectory is composed of slowly drifting vertical orbits, called trochoids.
Following the approach described in [27], we now give a quantitative definition of traps and flights along a trajectory. First, we compute the Lagrangian correlation time T L = ∞ 0 ρ(t)dt where ρ(t) is the temporal autocorrelation velocity function [33], as well as the spatial Lagrangian correlation scale L L = ∞ 0 ρ(r)dr, where ρ(r) is the spatial autocorrelation velocity function defined in [30]. L L is the characteristic length scale of the particle dispersion. It was shown that L L ≈ λ/2 for developed 2D turbulence [30,31]. Then we consider, along the particle trajectory, the 'jumping function' defined as [27]: Finally we refer to as traps, any sequences of consecutive stagnant points which correspond to a time span larger than 5T L . On the other hand, a flight is a succession of jumps which lasts at least 3T L . Points located at the trajectory ends (δt < T L /2) are not considered for the trap/flight analysis. Figure 2(a) shows the nature of traps and flights detected by these definitions. As expected, traps are characterized by extremely convoluted twists in all trajectories. On the other hand, flights correspond usually to the 'straightest' parts of the trajectories.
To assess correlations existing between the wave field topography and the horizontal features of the trajectory, the autocorrelation of the particle elevation was calculated in both traps and flights. It is defined as: where denotes the statistical average, z p (t) is the particle elevation at time t and z p (t 0 ) is its initial elevation which is obtained either at the beginning of a trap or a flight.
We note that both autocorrelation functions in Fig. 2(b) show an exponentially decaying envelope with a decay time T dec longer than the Faraday period T F . For traps we found T dec = 2.9 T F , while T dec = 1.9 T F for the flights. Figure 2(b) also shows that traps are clearly correlated with the subharmonic of the frequency excitation. This is because trapped particles rarely cross over the nodal lines around one wave peak. Conversely, the elevation autocorrelation function computed on flights reveals a substantial phase shift which can be interpreted as a Doppler shift. Indeed, in Fig. 1(a) the direction of the ballistic flight seems to be insensitive to the wave elevation and the flight itself corresponds to a 'straight' path probing different wave phases. In Fig. 1(a), we also note that a trap is clearly correlated with an area delimited by the nodal lines around a wave peak. It has recently been shown that the Faraday wave properties are consistent with those of an ensemble of oscillating solitons, or oscillons [32]. The oscillon peaks ( Fig. 3(a)) can be tracked and these positions are used to construct oscillon trajectories. Figure 3(b) shows a visualization of these oscillon tracks along with six particle trajectories; both were followed over 4 s at a supercriticality of = 1.67. The oscillons and the particles motions occur on different length scales. An oscillon wanders erratically within a cage whose typical size is roughly λ/2 [34]. Consequently, their confined tracks mark the sites of an irregular lattice, as confirmed by a peaked wave number spectrum of the surface elevation (see Fig. 3(c)). But particles move substantially in the wave field, similarly to what was observed in Fig. 1(a). Although the presence of these flights is connected with wave disorder [24], Fig. 3(b) emphasizes a clear distinction between the disorder-induced oscillon motion and the fluid motion. Figure 1(c) suggests a possible connection between the particles horizontal motions and high velocity events, or flights. The key role played by inertia in particle behavior is even better illustrated in the Eulerian frame. spectrum was computed from Eulerian horizontal velocity fluctuations measured using PIV technique (see the Appendix) at f 0 = 60 Hz and a vertical acceleration a = 1.6g ( = 1.67). Two scaling laws can be identified: k −5/3 at wave numbers k < 1500 m −1 and k −3 at k > 1500 m −1 . As reported in [23,24], these features are consistent with Kolmogorov-Kraichnan predictions for an inverse energy cascade and a direct enstrophy cascade respectively, which were derived for an incompressible fluid [35]. In these experiments, the vertical motion of the particles might result in non-zero values for the instantaneous divergence of the horizontal velocity field, however the divergence becomes close to zero when it is averaged over space or in time over several Faraday periods [23]. We also note that particles stay homogeneously distributed over the water surface. These observations confirm that the 2D flow generated by FW can be considered as incompressible. In contrast, floaters driven by an underlying 3D turbulence form massive ribbon-like clusters and their motion at the fluid surface exhibits properties of a compressible fluid [36,37].
The kink observed in the spectrum of the horizontal kinetic energy at k ≈ 1500 m −1 gives a forcing wavenumber corresponding roughly to half the Faraday wavelength λ/2 ≈ 4 mm (see Fig. 3(c-d)). This observation raises the following question: how do Faraday waves generate a spectrally localized forcing mechanism of 2D turbulence?
Wave induced vorticity. Figure 4 shows details of a three dimensional trajectory followed for 2 seconds in FW at relatively low drive ( = 0.67). By decreasing the vertical acceleration, the average 'trapping time' increases [27]. Thus, some trajectories can consist almost exclusively of consecutive trapping sequences as seen in Fig. 4(c). Figure 4 also reveals an interesting feature of particle trajectories that is somewhat hidden in the complex Lagrangian picture of Figure 1(a): traps can look like almost perfectly circular loops around an oscillon. It suggests that waves may generate local vorticity in the horizontal flow. In Figures 4(a-c), these horizontal loops have a diameter comparable to the size of an oscillon (≈ λ/2).
Moreover, Figure 4 emphasizes that this circular particle motion at the oscillon scale is somehow independent of the Faraday period but seems related to a much longer time scale. In Fig. 4(a) for instance, a particle describes an almost circular loop which takes 24 Faraday periods corresponding to the peaks in the green trajectory.
We now consider the horizontal projection of the particle trajectory and test if there is any resilient feature of those circular loops when the vertical acceleration is increased. The geometry of the horizontal projection of particle trajectory can be characterized by its local curvature κ. This quantity reveals the amount by which a trajectory deviates from a straight line and thus indicate the presence of loops or cusps: low values of κ are related to locally straight flights while higher values of κ occur in traps. The instantaneous curvature κ(t) along a trajectory is given by: where → v xy and → a xy stand for the horizontal velocity and acceleration along a trajectory, → z is the unit vector in the vertical direction. As the computation of → a xy is quite sensitive to the noise level, the curvature κ(t) was smoothed over one Faraday period T f , such that: x y On a technical note, the local curvature diverges when the velocity → v xy is zero (i.e. the displacement is smaller than the spatial resolution which is 100 µm for the data plotted in Fig. 5). These events correspond to the presence of sharp cusps along the trajectory (see Fig. 1(b)). In such rare cases, locally diverging points in κ(t) are removed before the quantity is smoothed over T f .
In Figure 5, we show the probability density function (PDF) of the curvature for different excitation frequencies f 0 and vertical accelerations. First, these PDFs are wide with distinct tails. The PDFs are found independent of the vertical acceleration for a fixed f 0 , Fig. 5(a). However, the PDF becomes more peaked and narrows around the origin with the decrease in the forcing frequency f 0 . It keeps a power law tail but high κ T f events get less probable, Fig. 5(b). Although the entire PDF follows a q-gaussian distribution, it can be fitted by an exponential decay law P (κ T f ) = αexp(−l c κ T f ) for the curvature values within the range κ T f ≤ 1000 m −1 .
Since the curvature is the inverse of the curvature radius, l c has the dimension of a length scale. This characteristic length scale is roughly a quarter of the Faraday wavelength as seen from the fitting parameters for two different cases at f 0 = 30 Hz and f 0 = 60 Hz in Figure 5(b).
In Figure 5(c), the absolute value of the horizontal ve- At the excitation frequency of f 0 =10 Hz, a quasi-regular square lattice of 7 × 7 oscillons is excited. Shortly after the parametric waves are generated, a matrix of vortices is formed which roughly corresponds to the underlying quasi-regular lattice of oscillons Fig. 6(a). The characteristic diameter of the vortices is close to the size of an oscillon. This configuration however is not stable as the vortices interact with each other. The flow evolves as illustrated in Figure 6(c), the vortex interaction results in the emergence of a disordered flow. The evolution of the vertical vorticity field → ω z = (δu y /δx − δu x /δy) → z computed using PIV velocity field (see details in the Appendix) and averaged over 4 seconds is shown in Fig. 6(b) and (d). In Fig. 6(b), the vorticity field consists of an almost regular square pattern made of counter-rotating vortices. These vortices ultimately aggregate into larger eddies (Fig. 6(d)). The total vorticity remains zero.

III. DISCUSSION AND CONCLUSIONS
How do Faraday waves generate a spectrally localized forcing mechanism of 2D turbulence? Figure 6 demonstrates that Faraday waves create horizontal vortices on the liquid surface at the early stage of the flow development. Initially a periodic array of vortices can be detected. These vortices have a diameter given by the size of an oscillon. As the flow develops, these vortices merge into larger eddies in the process of the inverse energy cascade. The initial vortices are the fuel of 2D turbulence.
The creation of vorticity by Faraday waves bears striking similarity with that in electromagnetically driven turbulence. Though the forcing scale vortices in electromagnetic turbulence are generated by periodic arrays of permanent magnets, the electromagnetically driven vortices injected into the flow are very similar to those illustrated in Fig.6(a). This explains why the 2D turbulence produced by these two different methods is very similar [30]. In both methods, the forcing scale vortices are only visible in the Eulerian frame in the initial stage of the turbulence development [43].
The experimental data show that the shape of the curvature PDF is a resilient feature of particles trajectories in disordered FW. The most probable values for the curvature can be fitted by an exponential decay law. The characteristic decay diameter is determined by the period of the oscillon lattice [24,34] which is also the forcing scale of the kinetic energy spectrum, Fig. 3(d). We note that a connection between a statistical property y (m) x (m) of the Lagrangian curvature (l c in Fig 5(b)) and the forcing scale detected in the kinetic energy spectrum (Fig. 3(d)) is not trivial. Indeed the Lagrangian and Eulerian descriptions of a flow are difficult to connect. Hence it is remarkable that the curvature of the Lagrangian trajectories still reveals a "footprint" of the forcing scale vortices in fully developed turbulence.
The Lagrangian motion in a 3D Faraday wave crystal. Figure 7 is to be compared with Figure 1. It shows the particle trajectories observed in a crystal-like wave field. An almost perfect wave crystal is realized in a viscous solution of glycerol (73%), milk (10%), and water (viscosity = 40 cP) [16]. In such a crystal-like array of waves, Fig. 7(a), particles no longer experience spatially extended random walks. Instead, they move along closed orbits of different geometries, Fig. 7(b). The orbital mo-tion is almost vertical near a wave antinode while it shows small horizontal oscillations close to a nodal line. Along a nodal line, particle trajectories are U-shaped. The crossover of two nodal lines marks the location of a saddle point in the wave topography where there is no particle motion.
As the vertical acceleration is increased (or if dissipation is decreased), this ordered wave structure breaks down into a disordered lattice of steep waves. Although Faraday waves remain almost standing waves, the presence of disorder in the wave field coincides with an increase in the particle mobility and the emergence of 2D Navier-Stokes turbulence. These observations suggest that 2D turbulence should be considered as a potential source of disorder in FW.

Conclusions.
These experiments prove that the deformation of a fluid interface can be considered as a new paradigm for vorticity creation in a 2D flow [38]. Besides, this finding sheds new light on the wave induced fluid motion which is commonly considered in textbooks as an irrotational flow. As a consequence, the influence of vorticity in surface wave phenomena was only studied relatively recently compared to the long standing history of surface wave physics. To our knowledge, the influence of vorticity has been mainly addressed in viscous boundary layers [39,40], or in coastal shear flows in oceanography [41]. In both cases, the presence of vorticity can have significant effects on the flow. In the latter case, the presence of strong nonlinearities in the governing equations has already stressed the importance of interacting vortices.
In Faraday ripples, a new physical model emerges: a 2D lattice of oscillating solitons contains vertical oscillatory energy, but it also produces horizontal vortices at the oscillon size: the 'turbulent fuel'. This horizontal energy is ultimately spread over a wide range of scales by the inverse energy cascade. This unexpected generation of 2D turbulence is an interesting twist in our understanding of the order-disorder transition in Faraday waves [19][20][21].

APPENDIX: MATERIALS AND METHODS
Faraday waves are formed in a circular container (178 mm diameter, 30 mm deep) filled with a liquid whose depth is larger than the wavelength of the perturbations at the surface (deep water approximation). The container is vertically vibrated by an electrodynamic shaker (Bruel&Kjaer). The forcing is monochromatic and set to either f 0 = 30 Hz or f 0 = 60 Hz. The amplitude a of the vertical acceleration imposed by the shaker is measured by an accelerometer. An accurate control of the acceleration amplitude is performed by a proportionalintegral-derivative controller.
We use a diffusing light imaging technique to measure the topography of the wave field, see Fig. 8(a) and 8(d).
The fluid surface is illuminated by a LED panel placed underneath the transparent bottom of the container. A few percent of milk (from 2% to 10%) added to water provides sufficient contrast to obtain a high resolution reconstruction of the parametrically excited wave field. The absorption coefficient is calibrated before each experiment, by gradually increasing the liquid depth and by measuring the change in the light intensity transmitted when the liquid surface is flat (i.e. when the container is not vibrated). This procedure allows to calibrate the wave field elevation in millimeters. In these experiments, the dynamic range of our images (16 bits) allows to resolve a 20 µm change in the fluid elevation.
Videos are recorded at a 16 bit resolution and at a high frame rate using the Andor Neo sCMOS camera which was mounted above the tank. The typical field of view is either a 8 × 8 cm 2 domain, imaged at 120 Hz with a resolution of 100 µm or a 3 × 3 cm 2 area, imaged at 587 Hz with a resolution of 200 µm. Black floating particles spread on the fluid surface are easy to observe and allow to visualize the horizontal motion of the fluid. We use particles with a diameter within a range of 150-300 µm (300 µm is the mesh size of the sieve used). Particles are made of carbon glass and have been plasma treated to reduce their intrinsic hydrophobicity. The use of surfactant and plasma treatment ensures that particles do not aggregate on the surface. This is illustrated in Fig. 8(a) with particles homogeneously distributed at the surface of the wave field.
To test if the particle size affects the observation of turbulence generation in Faraday waves, additional measurements were performed using 50 µm diameter polyamid particles. No effect of the particle size was detected, in agreement with observations reported in [23,30].
3D-PTV techniques: Three-dimensional Lagrangian trajectories are retraced using a combination of two-dimensional PTV technique and a subsequent estimation of the local elevation along the trajectory (see Fig. 8 3) and f0 = 60 Hz. Dark and light blobs give the location of wave peaks and troughs respectively. The tiny red dots are tracer particles (diameter=150-300 µm), the locations of which are extracted from the wave field using a 'rolling ball' filter. A green circle of radius λ/4 = 2 mm is indicated. (b) Slight zoom-in view on the horizontal projection of an individual trajectory followed for 4 s at 120 frames per second. This trajectory is reconstructed using PTV technique. The blue line (barely visible) corresponds to raw data points, the red line is the trajectory smoothed over one Faraday period (T = 1/f ≈ 33 ms). (c) 3D reconstruction of the trajectory shown in b). (d) 3D Visualization of the wave field topography shown in a). Particles have been removed from the image by using a 'rolling ball' filter, the original image (80 × 80 mm 2 ) was smoothed over a 0.8 × 0.8 mm 2 window.
horizontal projection (x-y coordinates) of each point on a trajectory is tracked using a nearest neighbor algorithm [42]. In such algorithm, the maximum displacement allowed for a particle in consecutive frames is set to be smaller than the minimal distance separating particle pairs in the field of view. Then the particle elevation (z coordinate) is estimated as the mean of the wave elevation over a local window (400 µm radius) which is centered on the x-y particle coordinates at a given time. In all our experiments, the localness of the z coordinate estimation is ensured by choosing a sufficiently large Faraday wavelength (at f 0 = 60 Hz, λ >> 400 µm). For each vertical acceleration a at a given forcing frequency f 0 , thousands of particle trajectories were followed for 4 s at 120 or 587 frames per second. The 3D trajectories of the particle and the wave field are visualized using the Houdini (T M ) 3D animation tools (by Side Effects Software). For the study of the horizontal features of particle trajectories, x-y particle coordinates are first smoothed over one Faraday period T f . This filtering limits the noise induced by the projection on the horizontal plane of the vertical oscillatory motion. Similarly, the determination of the curvature κ(t) requires the computation of the horizontal acceleration → a xy , a quantity quite sensitive to the noise level. To reduce the noise, the measured curvature κ(t) was also smoothed over one Faraday period T f .
2D-PIV techniques: The particle image velocimetry (PIV) technique is used to obtain the velocity field of the horizontal motion of the flow in Figure 3(d) and Figures 6(b) and (d). Images are split into a large number of interrogation windows which form a regular grid. A velocity vector is calculated for each interrogation area by temporal cross-correlation of the particle intensity distribution within each interrogation tile. The PIV technique allows to resolve sub-pixel displacement and thus gives a high degree of accuracy in the estimation of the velocity field.
In Figure 3(d), the flow is recorded at a frame rate twice the shaker frequency (i.e. equal to four times the parametric wave frequency) for 4 seconds. The field of view is 80 × 80 mm 2 with a spatial resolution of 100 µm (pixel size). The PIV velocity fields are computed on a 90×90 spatial grid (grid mesh size is 0.89 mm), with a 2.6 × 2.6 mm 2 interrogation window size (the interrogation windows are overlapping). The measurements of the instantaneous displacement are accurate down to 10 µm. The energy spectrum in Figure 3(d) is averaged over 400 snapshots of the velocity field.
In Figures 6(b) and (d), the flow is recorded at a frame rate of 5 fps for more than 10 minutes. The field of view is 400 × 400 mm 2 with a spatial resolution of 200 µm (pixel size). The PIV velocity fields are computed on a 70×70 spatial grid (grid mesh size is 5.7 mm), with a 10×10 mm 2 interrogation window size (the interrogation windows are overlapping). The measurements of the instantaneous displacement are accurate down to 20 µm. The vertical vorticity fields in Figures 6(b) and (d) are averaged over a 4 second time interval (= 20 T f ). The average vorticity measurement is accurate down to 10 −3 s −1 .

ACKNOWLEDGMENTS.
This work was supported by the Australian Research Council's Discovery Projects funding scheme (DP110101525). HX would like to acknowledge the support by the Australian Research Council's Discovery Early Career Research Award (DE120100364).