Moulding hydrodynamic 2D-crystals upon parametric Faraday waves in shear-functionalized water surfaces

Faraday waves, or surface waves oscillating at half of the natural frequency when a liquid is vertically vibrated, are archetypes of ordering transitions on liquid surfaces. Although unbounded Faraday waves patterns sustained upon bulk frictional stresses have been reported in highly viscous fluids, the role of surface rigidity has not been investigated so far. Here, we demonstrate that dynamically frozen Faraday waves—that we call 2D-hydrodynamic crystals—do appear as ordered patterns of nonlinear gravity-capillary modes in water surfaces functionalized with soluble (bio)surfactants endowing in-plane shear stiffness. The phase coherence in conjunction with the increased surface rigidity bears the Faraday waves ordering transition, upon which the hydrodynamic crystals were reversibly molded under parametric control of their degree of order, unit cell size and symmetry. The hydrodynamic crystals here discovered could be exploited in touchless strategies of soft matter and biological scaffolding ameliorated under external control of Faraday waves coherence.


Supplementary
. Linear response at low driving amplitudes. Surface wave spectrum in the linear response regime (black curve), obtained for very low vertical excitation amplitude Γ = 0.03 ≪ Γ , at 0 = 47 Hz (capillary regime).The observed resonance occurs exactly at the pumping frequency induced by the lauder (green curve). The observed spectral broadening is compatible with the natural bandwidth of the excitation channel (green straight line).

Supplementary Figure 3. Nonlinear surface waves regimes in water surface.
Images and spectral response of nonlinear surface waves driven on pure water under parametric excitation ( 0 = 12 and 0 = 47 , upper and down panels, respectively). For increasing vertical acceleration determining the driving force (from a to d) different wave regimes were observed, as described in the main paper: (a) linear regime exhibiting a single peak at 0 , (b) Kolmogorov-Zakharov spectrum of weak wave turbulence exhibiting container-shape surface waves and a −17/6intensity decay typical of cascades of nonlinear capillary waves, (c) Faraday regime exhibiting disordered surface roughening and a spectrum with subharmonic response at 0 /2, and an intensity decay cascade that scales as −5 ; (d) Chaotic regime featuring strongly disordered surface and a continuous spectrum (without any resonance) and a intensity decay proportional to −2 .

a) b) c) d)
Supplementary Figure 4. Hybrid spectrum of turbulence for pure water in the capillary regime ( 0 = 47 ) where both Faraday (F: ~ −5 ) and Kolmogorov-Zakharov (KZ: ~ −17/6 ) spectral decay coexist. The former was found in the low frequency range, while the latter in the high frequency spectrum. For increasing Γ (0.07 − 2.5) the hybrid spectrum converts into a pure Faraday wave cascade. As discussed in the main manuscript, the KZ spectrum is restricted to the higher frequencies range due to its lower intensity in comparison with Faraday wave cascade. ⁄ . This physics compatible with resonance between the liquid inertia and the external driving force is explained in Supplementary Note 1. Figure 6. Hybrid spectrum of turbulence obtained at the transition from Faraday to continuous spectrum for bare water surface in the capillary regime ( 0 = 47 ) by increasing Γ in the range (0.8 − 2). At higher values of Γ continuous spectrum emerges with characteristic Brownian-like energy cascade (~− 2 ). This cascade firstly develops in the lower frequency range and, for increasing Γ, propagates to higher harmonics. Also, it is noteworthy the fact that for increasing excitation amplitudes the spectral peak broad until they disappear at higher Γ.  aescin covered (right panels) surface waves driven in the gravity regime ( 0 = 10 ) for increasing vertical acceleration (from top to bottom panels). The aescin concentration was 1 mM. Only the introduction of aescin was able to induce an ordered surface pattern on the Faraday waves (magenta surrounded box) b) Comparison of images of bare water (left panels) and aescin covered (right panels) surface waves driven in the in the capillary regime ( 0 = 47 ) by increasing the vessel vertical acceleration (from top to bottom panels). By increasing the vertical acceleration, surfaces waves emerged due to the parametric resonance and organized in a regular packing leading to the formation of a hydrodynamic crystals which achieved their maximum packing order at intermediate accelerations (magenta surrounded boxes). By further increasing the driving acceleration a crystal breakup occurred.

Supplementary
Supplementary Figure 9. Representative log-log spectra of nonlinear surface waves turbulence in aescin covered water surface. Below the capillary frequency gravity waves (left panels) are found, while capillary waves (right panels) are observed above. The sequence of spectra spans from top with the linear regime (only the fundamental response is present at 0 ) to bottom; then, the incoherent KZ-turbulence characterized by a spectral decay exponent that depends on the dispersion regime ( = −4 for gravity regime and = − 17 6 ⁄ for capillary regime), through of the Faraday instability (where the subharmonic response appears at 1 2 ⁄ = 0 2 ⁄ and the Faraday waves turbulent cascades decay as = −5), up to the chaotic regime in the lower panels ( = −2). Crystal formation is only observed in the Faraday regime ( = −5).

Supplementary Tables
Supplementary Table 1 Critical acceleration thresholds for the non-linear surface-wave states We considered bare water surface with a high surface tension, low-tension aqueous surface with a CTAB-adsorption layer (light gray), and high-viscosity water/glycerol mixture (dark gray). We report the experimental values of critical acceleration (aX, where X stands for L: linear; KZ: Kolmogorov-Zakharov; F: Faraday; C: Chaotic) on the onset of different parametric wave states as quantitatively identified by the slope of the turbulence cascade detected by means of laser Doppler vibrometry (LDV). The onset for Faraday instability (aF; marked in blue) measured as a function of the driving frequency (ω 0 ) follows the theoretical prediction as given by Eq. (1) of the main text for the quoted values of surface tension (σ) and the kinematic viscosity (η) (see Supplementary Figure 1 and Fig. 3c and the main text for details). The Faraday wave (FW)-dispersion status is determined by the value of ω 0 with respect to the capillary frequency ω (for water, ω = 14 ). Thus, gravity waves occur below or close or below ω , whilst capillary waves emerge for ω 0 > ω . To study the FWs in a defined dispersion regime, the experiments reported in the main manuscript corresponded to the reference frequency ω = 47 , i.e. in the capillary regime. By analysing the results reported in this Supplementary Supplementary Table 2 Critical acceleration thresholds in water covered with aescin To study the appearance of NLSWs in the different regimes of wave propagation, we focused on the amplitude transitions at variable acceleration in three well-delimited dispersion regimes; these are: hybrid gravity-capillary waves close to the capillary frequency (GCW at ≈ = 14 ), pure gravity waves (GW at ≪ ) and pure capillary waves (CW at ≫ ). We detected the acceleration thresholds as reported in Supplementary

1-1. Nonlinear surface waves: resonant interactions
The nonlinear hydrodynamics of interacting water waves has been the subject of intensive research as a first principle of fluid motion in two dimensions [1][2][3]. The notion of resonant interaction has been central on all those advances not only as constitutes the hydrodynamic skeleton for harmonic ensembles of nonlinear surface waves (NLSW) [3][4][5], but specially for building Faraday waves (FWs) at subharmonic resonance [6,7]. Fluid surfaces respond to periodic forces at their natural frequencies as far as they behave as harmonic oscillators when coupled to their bulk counterpart [8]. By forcing a fluid to oscillate, different classes of NLSWs are generated with a hydrodynamic structure that depends on two key factors [5][6][7][8][9]: a) The specific nature of the driving force; either random or monochromatic, eventually leading harmonic resonance under extrinsic forcing.
b) The underlying hydrodynamic skeleton; inducing wavefield interactions and resonances between coupled surface modes.

1-2. Nonlinear stiffening as a NLSW-FW organizer
As refers to the new class of surface hydrodynamic skeletonization involved in this work (as a stiffness-induced coherent ordering of the Faraday wavefield at a fluid surface), the ordered NLSWs appeared as FW patterns (hydrodynamic crystals) are envisaged to emerge as ordered condensates of harmonic and subharmonic resonances with the external monochromatic source of driving [5][6][7].
Supplementary Figure N1. Mechanical equivalent of nonlinear surface wave motion in a vertically vibrated liquid. We argue around a minimal dynamical depiction as a Duffing Non-Linear Oscillator (DNLO) that qualitatively explains the main features observed in experiments. A) Schematics on NLSW-DNLO dynamics as a transverse surface wavefield externally excited upon liquid vibration that is globally exerted in the container by an oscillating AC driver (of force amplitude 0 and monochromatic frequency 0 ). This waving motion is restored by liquid inertia and surface elasticity (linear and nonlinear), and frictional dissipation operated under viscous drag. Mechanical ingredients: effective inertial mass defining linear momentum ( ); surface elastic response as a linear spring constant ( ; which recapitulates gravity and surface tension); Duffing-like wavefield self-interaction as produced upon nonlinear surface elasticity ( is a nonlinear rigidity). B) Elastic free-energy landscape. A Hookean harmonic response characterized by the linear spring constant is assumed at a first instance. Nonlinear rigidity causes the quadratic symmetry to be broken in the Faraday wave selffocusing characteristic of the hydrodynamic crystals experimentally observed upon surface stiffness functionalization with adsorbing aescin (FW-freezing). The nonlinear coefficient represents either surface stiffening leading wave self-focussing ( > 0), or surface softening causing wave defocussing ( < 0). C) Mechanical equivalent corresponding to the force ingredients depicted in A) as coupled in series; from left to right: inertial mass ( ); viscous damper with Newtonian friction coefficient ( ); elastic spring ( is the Hooke-spring constant; the nonlinearity rigidity coefficient describing the strength for wavefield self-interaction); time dependent external force ( ) = 0 ( 0 ), which excites the DNLO at its natural frequency 0 = ( ⁄ ) 1  As a minimal implementation we propose a mechanical model with subharmonic resonances (FW alike) at the core of the field self-focussing interaction (FW-ordering).
Supplementary Figure

1-3. Towards a minimal NLSW theory
Harmonic resonance is assumed the organizational principle that guides NLSW ordering in self-focussed wavefields. We invoke the essential ingredients of harmonic resonance to build a mechanical equivalent system at extrinsic resonance with the driving force. Our theoretical construction captures the emergence of FWs as parametric resonances with liquid inertia, the so-called mass-force resonances. We then suggest hydrodynamic crystals as FW-moulded wavefields with a standing wave pattern coherently imposed by a structure of harmonic (and subharmonic) resonances with the external source. Surface stiffness is assumed to be the key factor leading to enough phase coherence upon wavefield self-focusing. Chaos is ultimately viewed as a frictional dominance of decoherence at high wavefield velocities. The weak turbulence regime is finally envisaged as the natural response at low amplitude when mass-force resonance does not occur yet. In this sense, the classical Kolmogorov-Zakharov (KZ) picture of dispersion-governed weak wave turbulence will be recovered [3,10]. A one-dimensional wavefield ( , ) is assumed as the mechanical response in a single surface point upon external driving with a monochromatic AC structure ⁄ for the FWs [13,14]).
The minimal prerequisite for harmonic resonance is some spring system that when deflected from some rest state experiences a restoring force that pushes it back toward the equilibrium state. We assume the surface Hookean response embedded into a spring constant ( , ), as composed of the gravity and capillary restoring forces; this fixes the natural frequency at 0 = ( ⁄ ) 1 2 ⁄ .
Also required is bulk inertia, or momentum ingredient, that makes the wavefield to overshoot the equilibrium point and pass on through. As an indispensable ingredient for FWs, we consider liquid inertia embedded in a fictitious inertial mass that enables the elastic surface to oscillate between a mass-force system that allows parametric resonance as a time-dependent restoring force. The system is driven back and forth the equilibrium at frequency 1 2 ⁄ [5,8]. Bulk inertia is usually neglected in the models of weak NLSW [3,[9][10][11][12], which focus on statistical kinetic equations for the random wavefield velocities as governed by interactions between the dispersive surface modes (see Supplementary Figure 1). The FW regime has been shown as dominated by bulk inertia [13][14][15][16], exhibiting a very characteristic spectral density decay, at −5 rate, independently of the wave dispersion domain (see Supplementary Figure 3).
Furthermore, bulk friction is considered as a sink for dissipation. The viscous drag ingredient is described by a friction coefficient . We expect a leading role for friction in decohering resonant correlations at increasing velocity, thus being necessary to describe frictional death at low forcing (weak turbulence), and the dynamic pathway to chaos observed at high forcing (Landau's catastrophe).
Finally, as the distinctive component representing stiffening nonlinearity, we consider a self-interacting wavefield term imposed by surface shear rigidity (promoting wave freezing) [15][16][17]. In a minimal symmetry-preserved, nonlinear form, the wavefield response would be perturbed as (1 + 2 + ⋯ ), with being a stiffening constant leading to field self-focussing (see Supplementary Figure N1  factor for the phase locking required for stable FW patterning (see Figure 4).

1-4. Mechanical equivalent: Duffing nonlinear oscillator
We consider the vertical displacement at a single surface point (see Supplementary Figure N1). The wavefield is folded in ( , ) as captures a discretized NLSW ensemble = { }, which is extrinsically determined under a highly narrow distribution of normal stresses driven by the monochromatic source. Because the injected energy is assumed to continuously flow from the driving source through the vibrating liquid up to the surface (where it spreads over the whole NLSW ensemble), we expect the above dynamical components as coupled in series. Once the fluctuating parts have been separated from the mean flow, our deterministic description of the wavefield adopts the form as a dynamical system with the mathematical formulation of a resonantly forced, nonlinear oscillator (NLO) [8,15,16]; this is written as the Duffing oscillator in a generalized form:  and its spectral signature, in the corresponding scaling-law decaying cascades of power (with the observed subharmonic response at 1 2 ⁄ , and the = −5 decay rate of the power spectral densities; see Figure 2).

1-5. Faraday waves as a parametric resonator
Parametric resonance is a phenomenon caused by periodic changes in the spring parameters [15,17]. If the frequency of this modulation is an integer multiple of the natural frequency 0 , a subharmonic resonance can occur at doubled period ( 1 2 ⁄ = 0 2 ⁄ ), causing the newly created mode to become amplified [17]. Faraday waves (FWs) constitute a paradigm of subharmonic NLSW resonance at extrinsic coupling with an external driving force causing vertical liquid vibration [6,2,17]. Because energy is cyclically injected, the FW response has been described as a forced oscillator at parametric resonance [18]. Resonantly forced oscillators like the DNLO, undergo a bifurcated response leading inertia-restored subharmonic response in addition to nonlinear resonance [19]. In quantitative terms, DNLO solutions split in two superposed cascades of harmonics and their counterpart subharmonics as = 0 + (1), one gets: Assuming the harmonic cascade dominate over the subharmonics ( < 1), the stiffening Duffing term can be expanded as (up to second leading order): so that, the time-dependent spring term writes as (up to second perturbative order): 1) The main, generating wavefield 0 composed by the fundamental mode (and their harmonics) as excited by the external source: Under oscillatory forcing, it exhibits a natural response at extrinsic resonance with the external source; so that, the fundamental frequency is 0 = √ ⁄ (= Ω). Because of the self-focusing structure of the Duffing-like nonlinearity chosen (cubic spring of elastic restoring force = − 0 (1 + 0 2 )), only odd harmonics are generated, i.e. 0 , 3 0 , 5 0 , … A complete sequence including even harmonics can be predicted from more complex NLO structures (from possible symmetry breakings, e.g. = − 0 (1 + 0 + 0 2 )). For the sake of simplicity, we restrict to the minimal DNLO construction for the FWs: 2) The subharmonic response at parametric resonance 1 2 ⁄ , with time-dependent modulation of the spring parameter imposed with the same selffocusing structure as the nonlinear Duffing term [8,17]; this is: Assuming a parametrically resonant response imposed by the driving force, with leading term in the form of 0 = 0 (Ω ), Eq. (6-7) correspond the damped Mathieu equation usually invoked to describe the FW subharmonic instability [13,14]: with a natural response at resonant frequencies Ω 's, and Mathieu parameters: If = 0 ( = 0), then = 1 and Ω = 0 , so that the harmonic oscillator is settled.
Very pertinently, the subharmonic wavefield 1 2 ⁄ elicits the DNLO system to resonate at half of the natural frequency, i.e. at 1 2 ⁄ = 0 2 ⁄ (as observed in our NLSW experiments as an idiosyncratic FW distinctive; see Figure 2). In the DNLOdescription, the wavefield generator depicts the subharmonic resonance as emerging from the intrinsic coupling between liquid inertia and the natural surface response (in the absence of external driving). The subharmonic field 1 2 ⁄ is coupled with the main wavefield 0 via parametric modulation, see Eq. (6-7). The wavefield decomposition enlightens the link between the extrinsic harmonic resonance (the periodically forced DNLO field 0 ), and the intrinsic sub-harmonic resonance due to field self-focusing (the Faraday wavefield 1 2 ⁄ ). Physically, the subharmonic resonance emerges from the periodic modulation of the Hookean spring response, which yields parametric resonance under frequency ( 1 2 ⁄ = Ω 2 ⁄ ) and amplitude conditions ( 0 2 > 1) [14][15][16][17][18].
The DNLO structure imposes a commons inertia-governed decay within the hybrid harmonic and sub-harmonic cascades governed by the effective mass , independently of the specific natural response as embedded in the spring constant (and the friction constant ) [15,16,19]. The appearance of these spectral features at specific FW-signature is analysed below.

1-6. Power spectral density
To describe the DNLO spectral responses in frequency domain, we follow the prevalent formalism considered in the surface wave turbulence (SWT) theory [3,9,10].
Likely, we compare theory and DNLO-predictions with measurements from NLSW experiments, which are harvested in terms of surface heights considered as wave amplitudes (more often used for analysis than the velocities directly measured by LDV). Here, we use the same definitions as introduced in the main text. This convenient representation deals with wavefield amplitude spectra as determined by the density of power per unit of frequency [3,9]: This Fourier transform is called the power spectral density (PSD), which phenomenologically represents the hydrodynamic structure underlying the energetic cascade ~ − (as described by a Kolmogorov-like scaling law characterized by the decaying slope , which describes the envelope of the spectral peaks). The PSD function measures the amount of covariance of the wavefield by each frequency , being determined by the height-to-height wavefield autocorrelation function calculated at periods spaced a lag time between events oscillating at frequencies = 2 ⁄ ; for periodic wavefunctions: which deals with the time average 〈 〉( ) calculated for each lag time .
In the experiments, we obtain the PSDs as the Fourier transform of the time series tracked for the surface roughness ( , ), as observed at a given surface displacement by the local LDV vibrometer; this is ( ) = ℱ{〈 ( , ) ( , + )〉}( ) [20].

1-7. Model prediction: DNLO spectral slopes
We get DNLO-wavefield solution obtained upon numerical integration of Eq. Eq. (10-11) for trajectory tracks ( ) obtained from the above scheme as implemented in Wolfram Mathematica [21]. We numerically perform the integration by invoking the autocorrelation theorem based on the Parseval's identity [22]: where ℱ{ ( )}( ) is calculated as the Fourier transform of the time series; this is: From a completely analytic approach below, we get further insight on the inertial character of the spectral FW-signature ( = −5).

N1-8. DNLO harmonic balance: inertial and frictional signatures
Analytic expressions can be predicted for the PSD's envelopes in the theoretical DNLO setting. We exploit the harmonic balance to perform homotopy analysis between the driving force ( ) = 0 Ω and the expected solution ( ). Assuming oscillatory solutions with the oscillatory form ( ) = 0 , from Eq. (1) we get: and multiplying this equation by its complex conjugate, one obtains: For = 0, it reduces to the energy spectrum of the forced harmonic oscillator [8]: By expanding Eq. (15), for relatively small nonlinearity ( ≈ 0 + ( 0 2 ) at > 0), we get the approximate formula: i.e. although strengthened by the stiffening term ( 0 2 ), the DNLO energy is essentially determined by the frequency dependence of the function 0 2 (forced oscillator), which measures which frequencies contain the signal's energy (in units 2 ). This is somewhat different to the power spectral density (PSD) as defined by Eq. (13), which is the measure of signal energy content per frequency; the PSD can be looked upon as a frequency-domain plot of power per unit Hz vs. frequency (in units 2 ⋅ ). While the energy spectrum in Eq. (15) estimate the area under the signal plot, the PSD assigns units of power to each unit of frequency and thus, enhances periodicities.
Consequently, the DNLO power spectral densities can be calculated from the energy spectrum in Eq. (17) as: which represents the spectral envelope under stiffening perturbation modulated by a form factor centred at the natural frequency 0 : Supplementary Figure N3. DNLO spectral envelopes at harmonic balance. PSD profile for the Duffing oscillator at inertial resonance ( = 1) with a driving source at 0 , Eq.(18) in black; linear resonance ( = 0; thin line); resonance focussing under Duffing stiffness, = 1, under form factor Ξ as given by Eq. (19); thick line. A tail envelope with a characteristic -5 slope (FWs) is predicted for the discrete cascade of nonlinear harmonics: ordinary; ≡ { 0 } (yellow); extraordinary subharmonics; ≡ { 1 2 ⁄ } (blue). In the non-inertial, pure frictional case ( = 0; ≠ 0)) a -3 slope is predicted, in agreement with experiments at highly frictional conditions. As shown by Supplementary Figure N3, two limiting behaviours are clearly discernible in the spectral tails of the PSD cascades: a) Pure inertial ( ⁄ ≫ 0 ), with a PSD of amplitude proportional to the imposed acceleration (Γ = 0 ⁄ , in reduced units), decaying with frequency as: corresponding to the dominance of the extrinsic resonance between the driving source and the liquid inertia at the natural frequency 0 (with characteristic inertial signature

1-10. Loosing Faraday wavefield self-focussing: frictional death
The −5 spectral slope has been theoretically identified as a DNLO-wavefield attribute controlled by inertia (Phillips-like; at resonant balance with the driving source). This spectral behaviour is equivalent to the FW resonance observed in experiments (see Fig. 2 and Fig. 4 in the main text and Supplementary Figures 3-6). In addition, a transition to a −3 spectral slope has been realized as a progressive dominance of friction at loosing resonance. This behaviour is known in the DNLO's literature as a reverse bifurcation upon frictional defocusing [28]. That KZ-regime will be discussed in Supplementary Note 2 as a special case of noninertial (incoherent) DNLO dynamics governed upon the nonlinear Schrödinger (NLS) equation [3].

1-11. Transition to unsteady chaotic regime: Landau's spectrum
Our experiments showed a transition from a steady FW-regime of discretized modes (Γ ≤ Γ ≤ Γ ℎ ), to an unsteady chaotic regime of unstable NLSW motion above a force threshold at Γ ℎ ≈ 1. Whereas the hybrid cascades of harmonics { } and subharmonics { 2 ⁄ } coexist superposed at isolated frequencies in the FW-regime, at Γ ≥ Γ ℎ we observe a continuous spectrum that contains many modes filling the spectral density of states (see Fig. 2 and Fig. 4; lower panels). As envisaged by L.D.
Landau on his Paper #52 On the Problem of Turbulence [29,30], this class of continuum spectrum does appear as a signature for the deterministic initiation of a progressive unsteadiness of the principal motion composed of discrete modes at isolated frequencies [29].
Building upon the steady motion that forms the discrete spectrum for the principal stationary motion, Landau supposed that the frequencies of the continuum spectrum correspond to unstable modes created as new degrees of freedom from the former discrete modes, i.e. the turbulent flow field can be expressed as a Fourier series = ∑ , with time-varying components of the form ~ ( is a rate of unsteadiness). As turbulence increases, the number of degrees of freedom becomes likewise progressively large, thus leading chaotic turbulence once superposed [29].
As a matter of fact, the modulus 2 of the nonstationary chaotic motion should not increase infinitely, but rather remains bounded to a certain upper limit for stability, i.e.

1-12. DNLO chaotic regime: Landau's mode superposition
We have checked the above Landau's turbulence in the DNLO setting. Figure SN6 shows a series of simulations at subharmonic resonance. A transitional scenario is evidenced with increasing external driving force even at low nonlinearity ( ≈ 0.01); we detect two characteristic features: a) New and more offspring modes observed to emerge as satellites around the discrete modes that constitute the principal response (Supplementary Figure N6 a)

2-1. Discrete KZ-cascades in the massless DNLO wavefield: One-dimensional MMT turbulence.
Weak NLSW turbulence is known to take place in systems of nonlinear dispersive waves [3,9,10]. The energy transfer between waves compatible with their dispersion structure occurs mostly among resonant sets of waves in a random mixing of phases.
Zakharov showed the weak turbulent dynamics to satisfy linear evolution equations for the mode amplitudes; building upon the Hamilton equations through of a perturbed interaction ℋ = ℋ 2 + ℋ at their kernel, the dynamic equation primarily describes kinetic driving upon wave action ℋ 2 = ∑ 2 , plus mode-mode interactions ℋ ; in Fourier space [3,9,10]: For weak and random wave interaction, the conservative interaction term is calculated  [31], which is calculated for a continuum of coupled modes across the inertial interval [3,10]. The dynamic status of such a system is called turbulent as far the large wavelength greatly differs from the smaller dissipation scale (as occurred in turbulent flows [32]). The cascade's concept fixes the macroscopic manifestation of turbulence: the mean rate of the energy pumping does not depend on viscosity along the inertial domain of turbulence [32].
Assuming weak interaction at lowest order perturbative expansion (ℋ = ℋ 3 + ℋ 4 + ⋯ ), Zakharov described a transport of the integral of motion propagating wave action and wave collision interactions from the driving source to small scales. For the cascade picture to be valid, the collision integral assumes locality in k-space [3]. Under random forcing (ℱ chosen as a Gaussian distribution of noise), at the inertial interval  [3,9]. In this case the equation of motion takes the same form as a forced, damped oscillator devoid of bulk inertia. In terms of field variables, the noninteracting version of Eq. (24) becomes: with a friction coefficient that recapitulates the damping term .
Since we are here interested in a unifying framework, we refocus the simplified description based on NLO-equations at one dimension. A family of dispersive equations of the NLO-class was developed by Majda, McLaughlin and Tabak, namely the MMT model [33]. They assessed the validity of the statistical theory of weak turbulence for one-dimensional waves in the KZ-regime at random forcing; the generalized MMT equation is structured as [33,34,35]: where the exponent describes the spatial structure of the dispersive NLO. The nondispersive case is = 0, corresponding to constant propagation frequency (no spatial structure, e.g. the harmonic oscillator). The parameter describes stiffening nonlinearity with a spatial structure embedded in the exponent ; the relative = 0 and = 0 corresponds to the DNLO without inertia, = 0 in Eq. (1). Interestingly, = 2 describes Fickian-like dynamics.
An essential behavioral difference arose between MMT and KZ turbulence; the PSDscaling predictions differed very significantly with large scale forcing and dissipation [33][34][35]. A source of failure of KZ-conditions in MMT-NLOs was identified as a lack for the locality assumption for wave interactions inside resonant sets [35]; indeed, basic MMT does not consider intrinsic couplings explicitly.

2-2. Massless MMT-DNLO at strong turbulence: self-focused Schrödinger regime
When nonlinearity becomes dominant, or exceeds dispersion, different organized structures are predicted to occur because of strong turbulence (even at no bulk inertia; = 0) [4,9]. They belong, however, to the MMT class ( = 2 and = 0) of nonlinear wave packets described by the forced NLS equation [10]: in which the nonlinear term accounts for wave attraction (wavefield self-focussing under stiffening at > 0). Wave repulsion is described under softening (at < 0).

2-3. Mass and stiffness are at the core of the ordering wave interaction: a reliable FW evolver for hydrodynamic crystal condensation
The forced NLS equation describing KZ-developed turbulence recapitulates into a DNLO devoid of inertia. From the evidence raised, one can hypothesize that typical turbulent states (phase incoherent) correspond to a balance between dispersion and localized nonlinearity in -space [3][4][5][6][7][8][9]. Induced wave ordering (phase coherence) has been evidenced as a consequence of the detailed resonant balance between external forcing and inertia. The DNLO momentum term has been revealed as the crucial ingredient for predicting the −5 -feature typical of FW-spectra at parametric resonance under external forcing; lack of inertia causes, however, the extrinsic resonance quickly to dissipate into an overdamped −3 -regime (see Supplementary Figure N3) At insufficient self-focusing below a threshold for wavefield condensation ( ≤ ), FWs should be unable to spatially organize as observed in experiments with weaker surface stiffening agents that aescin (see Fig. 5c). With the strong surface stiffener aescin, however, hydrodynamic crystals do appear at sufficient rigidization ( ≥ ). This condensation is the natural consequence of self-focused wavefield conditions causing coherence to spatially organize the wave pattern with a permanent form, which is compatible with the dispersion structure of the piloting wave.
Supplementary Figure N9. Unification schematics for DNLO-NLSWs in the generalized MMT framework. The fundamental dynamic DNLO skeletonization is determined as a trade-off between the external driving force ( ), systematic momentum (as determined by the inertial mass ), and frictional stress (fixed by viscous drag ). A transverse timeline of NLSW state evolution is envisaged as a DNLO dynamics at modulated inertial response and non-linearity; from top to bottom: 1) Linear response ( = 0; > 0; = 0). 2) NLS regime from weak to strong turbulence ( = 0; ≠ 0). 3) Forced inertial DNLO for parametric resonance leading FW formation ( = 1): whereas surface stiffening causes wavefield self-focusing (the necessary condition for hydrodynamic crystal formation; > 0), surface softening elicits phase decoherence ( < 0). 4) Further forcing causes the system to enter the Landau's regime of chaotic turbulence, which resumes into a non-inertial inverted DNLO ( = 0 and < 0) at unsteady chaotic evolution in a repulsive barrier regulated by friction; the unsteadiness rate is determined by the inverted oscillator constant referred to the friction coefficient ( = − ⁄ ); surface Duffing stiffening refocuses the wavefield thus relaxing chaotic unsteadiness ( = − ). The spatial structure of the NLSW states is envisaged as a transversal infrastructure along the dynamic DNLO skeleton. Wave dispersion appears at the core of this superposed spatial skeletonization as encoded in the MMT exponents for the spatial derivatives ( for the linear term and for the nonlinear wavefield self-interaction; see Eq. S26). The genuine non-dispersive DNLO is identified as the simplest MMT class ( = = 0). As a well identified spatially organized status, the NLS states of strong turbulence are stablished as a diffusive structure ( = 2; = 0). Other classes of spatial organization are available among the family of MMT relatives.

2-4. NLSW unification into the DNLO model: dynamic NLSW states
In a step forward the generalized wisdom of the physics underlying the discrete NLSW cascades here involved, the discretized wavefield in Eq. (1) should fulfil the Navier-Stokes equations subjected to appropriate boundary conditions. These are primarily given by the extrinsic spatiotemporal structure of the driving force ( , ). In our NLSW setting this force is oscillatory, spatially homogenous and perfectly monochromatic (as The different NLSW states are classified in terms of the evolution of their dynamic organization and the spatial structure as endowed by wave dispersion. B) Phase diagram at the disorder-order transition induced by surface stiffening ( > 0). The appearance of hydrodynamic crystals is a natural consequence of stiffening-induced Faraday wavefield self-focusing causing the strong wave phase locking observed in experiments; FW-molded patterns (see Figure 4 in the main text).

3-2. Hydrodynamic crystalline state: FW-freezing under surface stiffening
As a first cause for the FW-distinctive spectral slope ( = −5), bulk inertia has been

3-4. Wavefield discretization: elemental quantum of action
The equivalent string response appears classically quantized in the given sense that which takes the same form as the Klein-Gordon (KG) equation for the quantized version of the relativistic energy-momentum relation [40]. The action appears here "quantized", or discretized in units of 0 , corresponding to the "quantum of action" for the mass-spring equivalent particle that generates the field; in terms of the elemental pieces of mass, length and time of the system, the quantum writes:

3-5. Linear Lagrangian: effective mass interaction
Since each term in Eq. (30) is linear in the field , the forced wave equation represents a wavefield under discretized action from a quadratic Lagrangian density: which is composed by kinetic energy (the wave action) and two potential terms appearing in 2 . The first one, of "quantized" strength 2 = 2 0 2 , is the Hookean field on the surface strain ( 2 = 2 ⁄ ). The term proportional to 2 is known as a mass term, or inertial potential due to its subsequent interpretation in terms of particle mass.
This quadratic mass term is now at the core of the interaction as the first cause for inertial resonance (as in the DNLO theory; see Supplementary Note 1).
The linear Lagrangian with the KG-structure can be interpreted as an effective field equation for the discrete wave-particle oscillations created upon monochromatic excitation. Because the quantum of action is fixed by the driving frequency 0 ≡ 2 0 , the harmonic frequencies correspond to multiple integers of this quantum, i.e. = 0 = 0 2 ⁄ . This naturally discretized response has been cartooned in Supplementary Figure N11, which brings to light why only discrete states are available from the elemental quanta of action 0 . So that responding linearly at the natural frequency 0 requires the system with an action 0 . A two-folded action 2 = 2 0 is required for a nonlinear response at the first harmonics 2 = 2 0 , 3 = 3 0 for 3 = 3 0 , … and so on, up to the highest level = 0 . The discrete solutions are obtained as a quantized classical field = { }, whose discretized quanta at the harmonic frequencies , are system-representing quasi-particles with an effective potential mass considered as the wavefield generator for given driving strength and frequency, this is ≡ 0 0 2 ⁄ .

3-6. Accounting for DNLO-stiffening: quartic self-interaction
In order to account for nonlinear Duffing-like stiffening interaction, the quadratic Lagrangian can be complemented with a scalar potential for the quartic self-interaction 4 = 4 4 4 ⁄ . in terms of the nonlinear Lagrangian density, ℒ 4 = ℒ 2 − 4 , the global action for the discretized DNLO field can be written as: from which the equation of motion is obtained as the Euler-Lagrange minimizer: with nonlinear strengthening proportional to the Duffing parameter ( 4~ 0 2 ), which is also subjected to discretization imposed by the driving frequency.
More synthetically, we have constructed a scalar field theory for a nonlinear (DNLOlike) self-interaction in addition to a mass term. They both recapitulate the generating role of inertia and wavefield self-focusing as the fundamental components at the core of a scalar potential, which defines the Lagrangian connection between potential energy and momentum in a system-equivalent quasiparticle at wavefield interaction (see Lifshitz's field below). As an outlook on the future development of this novel physics, we have preliminarily formulated the essentials of the classical field theory the undergrounds the new material object here discovered, the hydrodynamic 2Dcrystal. Aligned with previous arguments and evidences on the analogy of Faraday waves and De Broglie matter waves [41,42], we are further explaining how the physics here described could represent a unique realization of a classical matter wave able to freeze a crystalline pattern in two dimensions.

3-7. Faraday matter waves
The inertia-generated FWs exhibit a material nature as a momentum equivalent in the moving liquid surface at parametric resonance with the bulk fluid: the DNLO-model has shown the parametric FW's supported in the subharmonic wavefield as generated in an elastic body (the fluid surface), at resonance with a material content that endows strong inertia (the bulky fluid). Both entities exhibit an intrinsic wave-particle duality.
Such reciprocity can be stated as matter waves in an analogous formulation as De Broglie duality, but in the same classical sense stated above. Because we have encoded the discretized structure of the normal DNLO-modes of a structured surface as an infinity of coupled oscillators that represent its elastic response, the KG equation which describes the systemic spectral blue-shift of the massive system with respect to the natural frequency; for the massless oscillator ( = 0), the non-dispersive harmonic oscillator is recovered as = .
By expressing Eq. (35) in terms of the linear momentum = , the quantized energymomentum relation for the DNLO is: which represents an amplified wave with an exponential growth. Note that these unstable high-modes would be killed by friction.

3-8. Self-focused wavefield spatial structure: Lifshitz's field for Faraday waves
Under the inertial resonant conditions captured by the nonlinear KG-equation at high enough self-focusing (induced e.g. under the surface "stiffening" endowed by the adsorbed layer of aescin), the Faraday wavefield should be prone to condense with a spatial organization mastered by the mother wavelength 1 2 ⁄ , which is inversely proportional to a piloting subharmonic frequency one expects harmonic solutions as: with dispersion law giving the natural frequency of the mode as: By adopting the field theory perspective, these behavior can be described by a Lagrangian density of the Lifshitz class for a FW-equivalent particle of inertial mass taking part of the self-interacting field characterized by the spatial structure operator Through of the exponent > 0, the Lagrangian in Eq. (43) contains self-interaction terms at higher order derivatives than the pure harmonic oscillator without spatial structure ( = 0). Consequently, the equation of motion can be expressed as: with an effective natural frequency: where the MMT-exponent could recapitulate the level of wavefield self-interaction in terms of stiffness as a breakdown of the discretized symmetry of the harmonic Hamiltonian. The case = 1 describes the well-known fields at unitary self-interaction, mass and diffusive-like wave dispersion = 2 .
The Lifshitz's action generator in Eq. (43) is envisaged as a minimal kernel for the development of a forthcoming theory of hydrodynamic crystals as frozen Faraday waves under a high amount of surface stiffness. This theory is sufficiently comprehensive of the essential material contents, i.e. bulky inertial mass and surface elastic response both discretized upon the monochromatic forcing.

Supplementary Note 4
Synthesis of aqueous suspensions of microparticles PS-MAA particle suspension.
Aqueous suspensions of spherical particles of poly(styrene-co-methacrylic acid) (PS-MAA) were prepared from surfactant-free emulsion by a conventional polymerization process of both monomers [44]. The polymerization reaction was carried out in a 500 mL round-bottomed, five-necked flask. In the flask outlets, a water-cooled reflux condenser, a T-shaped stirrer, a gas inlet and a contact thermometer were fitted. One outlet was used for introducing the chemicals of the reaction. The flask was introduced in a thermostat water bath which controls the reaction temperature to ±2℃. The synthetic process starts by introducing 253 mL of water, 18 g of styrene and 1.58 g of methacrylic acid. The styrene was washed with a NaOH 0.1M water solution to remove the monomer inhibitors. The stirring was fitted to 350 rpm and nitrogen gas was bubble during all reaction time to remove any oxygen traces. When the temperature arrived to 85 ºC, 0.22 g of potassium persulfate (KPS) was added as initiator of the polymerization reaction. After eight hours of reaction, the flask is removed from the bath and cooled in a water-ice mixture. The resulting polymer suspension is dialyzed against water over three weeks.
Scanning electron microscopy (SEM) was performed with a JEOL-JSM-6330F electron microscope operating at 10kV. By analysing SEM-images using ImageJ, the average diameter was estimated into 190 nm. Measurements of dynamic light scattering (DLS) with a 90-Plus/BI-MAS apparatus was performed to characterize particle size distribution. An average particle diameter of 200 nm was obtained, with a relative standard deviation of 0.030 indicating an extremely low size polydispersity.
Electrophoretic mobilities were estimated by using the ZETA-PALS mode of the BI-MAS equipment. The measures were performed at room temperature and an ionic strength of 1mM of KNO3. The analysis of the data gave a zeta-potential of 47.3 mV on negative charged particles, with charges come from the persulfate fragments.

PS particle suspension
We also prepared Polystyrene (PS) microparticles by following the same synthetic route as with PS-MMA, without adding MAA monomer and introducing a small monomer charge of the reactor (9 g). The slow polymerization rate of styrene required to perform the polymerization during 24h.
A dilute PS suspension was measured by DLS and the average particle diameter was 313 nm, with a relative standard deviation of 0,037.

Supplementary Note 5 Dispersion relation of GC-waves in the presence of a surface viscoelastic film
We are interested in the correction to GC-wave dispersion due to surface elasticity arising from a viscoelastic film adsorbed to the surface. The classical treatment of Hansen and Mann (HM) considers linear hydrodynamics for the propagation of capillary waves on viscoelastic surfaces [45]. To account for the case of surface waves propagating on aescin covered surfaces as the data represented in the Fig. 3b) of the main text, we considered the HM-theory for a highly rigid surface film, this is at the material limit = ⁄ ≫ 1. Essentially, the viscoelastic film is known to affect capillary motions inducing weak changes on the propagation characteristics ( and ) companioned by strong wave damping ( ). For the propagation frequency of the capillary waves, the HM-theory predicts a non-monotonic variation with the capillary wavenumber, from the ideal Kelvin value corresponding to the bare surface, this is 0 = ( 3 ⁄ ) 1 where we have substituted the factor = ⁄ = ∞ ⁄ in the HM-theory as written in terms of the overdamped frequency ∞ = ⁄ . This equation represents the highrigidity perturbation limit to the GCW propagation in a free surface covered with a viscoelastic film. It was numerically resolved to obtain the limiting solution →∞ ( ) plotted in Fig. 3b).

Supplementary Note 6 Experimental setup for the interfacial rheology measurement
We used a commercial hybrid rheometer (HDR100, TA Instruments) working in the surface shear mode for measuring the viscoelastic modulus of surface films at the air/water interface. The rationale is explained in Supplementary Figure N12.
Supplementary Figure N12. a) Du Nöuy ring configuration with the thermostatic vessel. b) The ring is in contact with the aqueous solution surface with surfactant, which develops a surface monolayer (here simplified by drawing a schematic head-tail amphiphilic molecule configuration). c) The rheological measurement is carried out by in-plane shear force stressed on the surface monolayer by the ring tool. d) Experimental measurements of the stress-strain relationship in the same systems considered in this work. Data in the Figure 5 of the main text were measured at 0.1% strain, corresponding to the linear regime.