Wave-based liquid-interface metamaterials

The control of matter motion at liquid–gas interfaces opens an opportunity to create two-dimensional materials with remotely tunable properties. In analogy with optical lattices used in ultra-cold atom physics, such materials can be created by a wave field capable of dynamically guiding matter into periodic spatial structures. Here we show experimentally that such structures can be realized at the macroscopic scale on a liquid surface by using rotating waves. The wave angular momentum is transferred to floating micro-particles, guiding them along closed trajectories. These orbits form stable spatially periodic patterns, the unit cells of a two-dimensional wave-based material. Such dynamic patterns, a mirror image of the concept of metamaterials, are scalable and biocompatible. They can be used in assembly applications, conversion of wave energy into mean two-dimensional flows and for organising motion of active swimmers.

C ontrolled manipulation of particles on a surface is important in a variety of material science and bioengineering applications. Accordingly self-assembly, or the autonomous organization of components into patterns or structures, has become a fast growing interdisciplinary research area [1][2][3] . The controlled arrangement of particles on a liquid surface would offer a flexible and tunable method of engineering surface properties, such as, for example, electrical or thermal conductivity.
Recently new approaches to the manipulation of particles at a fluid surface have been proposed. They rely on the generation of surface waves to control the motion of particles at the liquid-gas interface. For instance, parametrically excited waves are generated when a liquid surface is vertically vibrated beyond a certain acceleration threshold. Such waves, often referred to as Faraday waves 4 , are modulationally unstable 5 and are readily broken into ensembles of localised oscillating solitons, or oscillons 6,7 . In viscous liquids, oscillons can create spatially periodic patterns and it has been proposed that such patterns can be viewed as metamaterials 8 . The idea of employing Faraday wave patterns as templates for micro-scale assembly applications has recently been discussed 9 . Current understanding of what is achievable with those waves often relies on properties of particles such as their wettability and their density 10,11 .
At lower fluid viscosity, steep oscillons form disordered wave fields which are coupled to a random motion of fluid particles at the surface 12,13 . This motion represents a macroscopic Brownian walk, where the diffusion of particles at the surface can be modified by changing the wave height and the wave length 14 . It has been shown that the properties of this interface are consistent with diffusion at thermal equilibrium, and as such, it can be viewed as a tunable thermal metafluid 15 .
Another approach relies on propagating waves originating from a localized source to control the motion of fluid particles on the surface. In this case particles can be guided away from the wave source, as well as towards the source in the so-called 'tractor beam regime' 16 . In this case, as for Faraday waves, the waves are essentially nonlinear.
Here we show that there is another way of organising particles on a liquid-gas interface which is analogous to a method used in low-temperature physics to trap atoms in spatially periodic optical lattices 17 . We propose a method of remotely shaping the particle trajectories. It uses linear three-dimensional (3D) surface waves and relies solely on hydrodynamic forces. The main idea can be described as a combination of the concept of rotating waves, created by a superposition of two small-amplitude standing waves, and the Lagrangian drift of particles along closed paths in such waves. Though rotating waves (electromagnetic or acoustic) are usually generated in cylindrical resonators 18,19 , we show below that periodic patterns of rotating waves can be created in a square geometry, similarly to optical lattices. On a fluid surface, such waves rotate within subwavelength cells and possess local angular momentum which is transferred to the matter. This mechanism produces particle trajectories in the form of a spatially periodic lattice of nested orbitals. This method offers a high degree of control over the particle motion and the ability to confine particles to spatially periodic cells. We present experimental results on the creation and control of such dynamical liquid metamaterials, develop a theoretical model of particle trajectories and characterize the transport properties of a multi-unit cell lattice.

Results
Experimental set-up. Stationary surface waves are studied in a square container (400 Â 400 mm) filled with water to a depth of d ¼ 81 mm. The waves are generated by two orthogonal horizontally oscillating paddles whose motion is driven by two electrodynamic shakers. The motion of the paddles (amplitude, acceleration) is accurately controlled via an accelerometer-based feedback loop (see 'Methods' section). In these experiments, we study linear waves of small amplitude (HE1 mm o ol, where H is the wave crest/trough amplitude and l is the wavelength). These waves obey the gravity-capillary wave dispersion relation: where o is the angular frequency of the waves, K ¼ 2p/l is the wave number, g is the gravitational acceleration, a is the surface tension and r is the density of the liquid. Figure 1 shows schematically the experimental set-up (a,b) and a photo of the laboratory set-up (c,d).
In these experiments the frequencies of the paddle oscillations are chosen to fit an integer number of wavelengths into the square paddle-wall cavity (312 Â 312 mm 2 ). The relative temporal phase of the paddle oscillations can be tuned in the range of ±180°w ith an accuracy of ± 0.1°. This set-up allows the superposition of two planar standing surface waves to create a periodic wave field for which the relative temporal phase is controlled.
Wave-driven fluid motion. We study the motion of floating micro-particles on the water surface perturbed by surface waves. Two orthogonal plane standing waves create a 3D wave field as the one shown in Fig. 2a. First, we investigate trajectories of surface fluid particles tracked for one wave period T ¼ 2p/o at the nodal points. Nodal points are places on the surface where the local amplitude of the standing wave is zero at every instant in time. If the wave frequencies are equal, o 1 ¼ o 2 , the shape of the projection of the particle trajectory on the horizontal plane varies depending on the phase shift f between the waves, as shown in Fig. 2b-d. A straight line corresponds to f ¼ 0, a circle to f ¼ p/2 and an ellipse to f ¼ p/4. When o 1 ¼ 2o 2 , the trajectory is represented by a figure of eight (Fig. 2e). The projections of these trajectories on the horizontal plane are reminiscent of the Lissajous figures in a two-dimensional (2D) harmonic potential. In the case of an ideal irrotational fluid, the velocity field u P ¼ ( u P x , u P y ; u P z ) associated with such trajectories can be represented by the gradient of a velocity potential F(x, y, z, t): u P ¼ rF. The velocity potential in this wave field is given by: where K is the wave number, d is the fluid depth and A is the potential amplitude related to the wave amplitude H. The z-direction is upwards with z ¼ 0 being the level of the undisturbed liquid surface.
The changes in trajectories observed in Fig. 2b-d highlight that both the wave dynamics and the trajectories of the fluid particles at the surface depend on the phase f. This can be seen in the temporal evolution of the normal n f to the surface at a nodal point. Figure 2f,g shows the evolution of the horizontal projection of n f for f ¼ 0 and f ¼ p/2 over one wave period. At zero phase shift, the temporal trace is a straight line, while it is a circle for the p/2 phase shift. This means that two orthogonal waves, phase shifted by f ¼ p/2, possess local angular momentum. This angular momentum originates from the rotating character of standing waves, which is null at f ¼ 0 and increases with f. This effect is somewhat surprising since equation (1) for standing waves produced in a square geometry shows decoupled temporal and spatial evolution. Indeed, this decoupling is manifest when f ¼ 0. However, when f ¼ p/2, a progressive rotating phase exists locally in the system. This can be described through a Taylor series expansion of F near a nodal point (x n , y n , 0), which reads: where dx ¼ dr cos y, dy ¼ dr sin y are expressed in polar coordinates (r, y) centred at the nodal point and we note that sin(Kx n ) ¼ ± sin(Ky n ) ¼ ± 1. In this case a rotating phase (y±ot) appears. This rotation is clearly seen in the experimental visualizations of the wave (Supplementary Fig. 1, Supplementary Movie 1). Now we focus on the wave motion and particle trajectories when f ¼ p/2. Figure 3a shows a snapshot of the wave topography measured at half the wavelength and marks the location of some remarkable points (peak, trough, nodal point, saddle points). First, the motion of the wave peak traces a square path, as illustrated by experimental data in Fig. 3b. The size of such a unit cell is L c ¼ l/2. The motion of the wave extrema is discontinuous along the edges of this cell. In contrast, the z ¼ 0 wave isoline rotates continuously at the frequency o around the nodal point within the unit cell, as shown in Fig. 3b and in Supplementary Figs 1 and 3. The conservation of mass underpins the small-scale orbital motion of fluid particles at the time scale of the wave period T (Fig. 2c). Here we show that a rotating wave can also transfer momentum to fluid particles. A slow drift of the orbits is observed in the direction of the wave rotation. This drift occurs along closed loops with a larger characteristic size (BL c ) and a large time scale (about 50T) as seen in Fig. 3c. The direction of the orbital drift is opposite in adjacent unit cells and it follows the rotation of the wave (Fig. 3d).
3D visualisation of travelling waves and the particle drift. From a geometrical viewpoint, the existence of a small-scale gyroscopic motion coexisting with a slow drift motion is reminiscent of the Stokes drift observed in a planar progressive wave 20 (see the Theoretical Model section). This drift is revealed when a water wave is described from the Lagrangian perspective, or from the point of view of the fluid particle paths.
Experimentally, the Lagrangian nature of the particle drift can be illustrated by visualizing 3D trajectories of particles. Here we use a recently-developed method of 3D particle tracking on the surface perturbed by waves 16,21 with high spatial resolution (E10 À 3 l, where l is the wavelength) and high temporal resolution (E 0.05T). An example of a reconstructed trajectory is shown in Fig. 4a (red) along with its projection on the horizontal plane (green). Figure 4b and the Supplementary Movie 2 show instantaneous surface elevations at times t i superimposed on the trajectory with the position of the particle indicated by a small sphere. During half the wave period (high wave amplitude), the particle progresses in the direction of the wave rotation, while it moves backward during the second half. The particle's speed when the wave crest hits it, is higher than during the wave trough moment. This results in a small displacement of the particle in the direction of the wave propagation when a wave cycle is completed. Note that the physics of the Lagrangian circular drift revealed here is intrinsically different from a recent Eulerian theory of vorticity generation on a surface perturbed by waves, which considers bulk viscosity as the essential ingredient of the mechanism 22 .
It should also be mentioned that we are able to generate such drift trajectories at the surface of various liquids with different viscosities (like glycerol-water solutions with viscosities in the range of (1-10) mPa?s).
Self-organised periodic lattice of rotating waves. We now analyse the emergent properties of the fluid motion in a multi-unit cell lattice of such wave-driven metamaterial. Figure 5a-c shows particle streaks in flows produced at different phase shifts f between two orthogonal standing waves. At f ¼ 0 the flow is stationary but disordered. The flow becomes more ordered as f increases. ARTICLE First, we compute the compressibility of the flow in the horizontal (x À y) plane as: Here u x,y denote fluid velocity components in the horizontal plane and Á Á Á h i Tav denotes averaging over time T av . Since the parameter C is computed only on the horizontal components of the velocity, it characterizes the dimensionality of the flow. Quantitatively, C can take on a value between 0 and 2. The lower the value, the more 2D is the flow. The value CE0.5 marks the onset of 3D effects 23,24 . Figure 5d shows that C is a decreasing function of T av . For time scales comparable to the drift characteristic time (B50T), the compressibility C is actually much smaller than 0.5 independently of f (Fig. 5d, inset). While the drift mechanism is a 3D phenomenon, the slow flow produced as a result of this drift is essentially 2D.
To characterize the fluid transport properties of the lattice we use the structure function r r where E is the total horizontal kinetic energy of the surface flow. The angular brackets denote averaging over different positions r 0 in the flow of the products of horizontal speeds u separated by a distance r. Figure 5e shows the Fourier transform of the structure function r(k) (where k ¼ 2p/r) for different relative phases f. At f ¼ 0, the spectrum r(k) spreads over low wave numbers (ko100 m À 1 ) corresponding to the large-scale disordered streams seen in Fig. 5a. As the phase f is increased, the broad distribution at low wave number is replaced with a strong peak at k w ¼ 2p=L C $ 160 m À 1 . This peak characterizes the order emerging in the wave-guided transport at the spatial scale corresponding to the unit cell. While the phase f is akin to a control parameter of the rotating wave momentum, the magnitude of the peak r(k w ) can be viewed as a structure factor of the liquid-interface metamaterial. Figure 5f shows that r(k w ) grows exponentially in the range f ¼ [0°, 40°] by a factor of about 20 and then saturates.
The emergence of order on the liquid surface with the increase in the phase f is not the only strong effect. The total kinetic energy E of the horizontal fluid motion increases with f, Fig. 5g, and for fZ60°it is more than twice as large as that at f ¼ 0°. Importantly, in the experimental data shown in Fig. 5, the wave amplitude produced by the oscillating paddles (that is, the energy injected in the system) is kept constant for all f. Thus the emergence of rotating waves at the unit cell scale increases the conversion rate of the wave vertical oscillatory energy into spatially ordered horizontal kinetic energy.
Theoretical model of the 3D rotating drift. As mentioned above, it is tempting to draw an analogy between the 3D rotating drift and the 2D Stokes drift for planar progressive waves 20 . In the classical picture uncovered by Stokes, the Eulerian velocity at a point averages to zero over one wave period, Fig. 6a, but the Lagrangian velocity of a fluid particle gyrating in the neighbourhood of that point does not, such that the fluid particle drifts in the direction of the wave propagation (Fig. 6b).
To test the relevance of such an analogy, we use the incompressible Euler equations to model the flow u ¼ (u x , u y , u z ) and the liquid surface Z ¼ Z(x, y, t) of an ideal fluid perturbed by orthogonal standing waves. Such equations read: @u @t þ u Á ru ¼ À 1 r rp ð3Þ with the following boundary conditions: where z is measured upwards from the unperturbed liquid surface and r is the fluid density. The modified pressure p includes the gravitational pressure rgz (constant atmospheric pressure is assumed at z ¼ Z). Equation (5) imposes zero flow through the rigid bottom of the container. Equation (6) expresses horizontal momentum balance at the surface ('dynamical free surface boundary condition'). Equation (7) relates changes in the surface displacement to the vertical component of the flow ('kinematic free surface boundary condition').  Under these conditions, the velocity field u P ¼ rF (equation (1)) is a solution. Assuming the phase difference f ¼ p/2, we obtain: where L Ao cosh(Kd)/g has the dimension of a length. Lagrangian particle trajectories can be computed by numerical integration of equation (8) with z ¼ Z. The resulting trochoid-like trajectories circulate in nested orbits about the unit cell in the direction of the rotating wave (Fig. 6c,d) in agreement with the experimental observations. A given trajectory exhibits smallradius gyrations at frequency o superimposed with a circular drift about the unit cell at a much lower frequency. Analytically, the drift velocity U d can be computed as:   where Á Á Á ð Þ ¼ 1 Þdt is a gyro-average over one gyration period T ¼ 2p/o (note that U d is a time-averaged Lagrangian property and u P is an Eulerian velocity). The drift speed U d approximately matches the experimentally measured one. Counter-intuitively, a quasi-linear irrotational model can produce vortex-like structure at a fluid interface 25,26 . For more quantitative agreement, we include an additional steady rotary motion at the surface to take into account the finite steepness of the waves (see 'Methods' section). The existence of such steady rotation has also been demonstrated in numerical simulations of non-linear waves at the surface of an ideal fluid 27,28 .

Discussion
Our results show that a horizontal flow on a liquid surface can be ordered into a perfect lattice of counter-rotating vortex-like structures by adjusting the temporal phase between two orthogonal surface waves to p/2. Such vortices confine floating particles within the unit cells. Manipulation of matter using waves is well known in a range of physical contexts [29][30][31][32][33] . It is particularly interesting to compare trapping of fluid particles within unit cells of the surface-wave metamaterial with the case of optical lattices in which atoms are trapped in the potential landscape created by two standing optical waves 29,30 . In optical lattices, the formation of radiation pressure vortices has been reported 17,31 . Radiation pressure is related to the momentum of electromagnetic waves, and the vortices are generated when two orthogonal standing waves have their temporal phases shifted by p/2, similar to the case of the surface waves. The wave energy flows along closed paths that resemble fluid vortices and the radiation force acting on the atoms is non-conservative 31,32 .
Although such an analogy between optical lattices and the surface-wave-based metamaterial is remarkable, the two systems are quite different.
Surface waves are strongly dispersive and they obey mass conservation in the fluid. In addition, waves in vacuum and waves in a medium have other fundamental differences. There is usually no clear connection between a mechanical wave and the momentum it may generate in a fluid [34][35][36] . Examples include an acoustic linear wave which can propagate in a fluid without any momentum present, or a small gravity wave packet which can generate pressure disturbances (that is, momentum) in a fluid far away from the packet spatial location 34 .
In this context, the Stokes drift has a special role: it can be viewed as a wave momentum 34 . Indeed, the phenomenon intimately links the transport of surface fluid particles to a propagating planar surface wave. Here we uncover a wave configuration where the Stokes drift exists along a closed path: a rotating Stokes drift. This path is created by a locally rotating wave. In this sense, the results presented here point to the existence of a radiation-pressure-like force guiding particles at the liquid surface.
It would also be interesting to find an analogue of a particle polarizability by light waves in the context of water waves. In optics, the particle polarizability plays a major role in both the generation of optical vortices and in trapping of atoms. Similarly to the case of optical traps, it was shown that in the standing surface waves, particles' inertia and wettability conspire to trap particles either at the wave extrema or at the nodes 10 . This suggests that the analogue of the polarizability in the surface waves could be related to the wettability and inertia of particles. In our experiments, the rotating Stokes drift clearly dominates the trapping effects. Beyond these analogies, the scope of the recently introduced concept of a wave-driven metafluid 15 can be broadened to that of a wave-based liquid-interface metamaterial. As we show, by dynamically shaping a fluid interface using rotating waves, one can produce an effective 2D material endowed with prescribed transport properties. An important feature of such a system is the establishment of unit cells confining nested spatio-temporal structures. These structures would allow remote control of particles on multiple length scales.

Methods
Experimental set-up. Figure 1 shows schematically the experimental set-up (a,b) and a photo of the laboratory set-up (c). Computer controlled electrodynamic shakers (TIRA TV51140) are used to drive synchronized motion of the two wave paddles. The paddle accelerations are measured using two accelerometers (B&K 4507, 1,000 mV g À 1 ) that provide feedback to the system's motion controller (Vibration Research, VR9500). The phase delay f between the paddles is adjustable in the range of ±180°via an arbitrary waveform generator (HP 33120 A).
Particular attention was paid to boundary conditions at both the oscillating paddles and the fixed wall facing them. It was noticed that the presence of a meniscus strongly influences our ability to produce a well-controlled standing wave field. Indeed we observe that a meniscus affects the spatial control on the phase f in the cavity, and it was recently shown to affect the reflection coefficient of gravitycapillary waves on the wall 37 . Therefore we machine specific grooves in the container walls and paddles such that the contact line is pinned to the wall edge with no meniscus. The wave fields produced in this cavity match well the numerically modelled waves (see Supplementary Figs 1 and 3).
Flow measurements. The horizontal fluid flow is visualised using buoyant tracer particles (Polyamid, 50 mm) illuminated by light-emitting diode strip lights surrounding the transparent acrylic fluid tank. A high-resolution video camera (Andor Zyla X5.5; 2,560 Â 2,160 pixel; 100 fps) is used to film the motion of imaging particles. With 16-bit resolution, the camera provides sufficient pixel intensity and spatial resolution (B50-200 mm per pixel; Nikon f1.4/50 mm lens) for quantitative data analysis using well-developed particle imaging velocimetry and particle tracking velocimetry (PTV) algorithms.
The particle imaging velocimetry technique is used to obtain velocity fields of the horizontal fluid motion. During an experimental run, the particle motion (corresponding to different f) is recorded twice: once at a frame rate equal to the wave frequency f ¼ o/2p and then at (10 Â f) fps. This allows us to characterize accurately both the fast gyration motion and the slow rotating drift. For the data presented in Fig. 5, the field of view used for the analysis is 234 Â 234 mm 2 corresponding to a 6 Â 6 unit cell lattice (for waves generated at f ¼ o/2p ¼ 4.58 Hz, l/2 ¼ p/K ¼ 39 mm). The spatial resolution is 125 mm per pixel. The velocity fields are computed on a 40 Â 40 spatial grid (grid mesh size is E5.35 mm), with a 10.7 Â 10.7 mm 2 interrogation window size (the interrogation windows are overlapping). The measurement resolution of the instantaneous displacement is subpixel. The wave number spectra in Fig. 5e are averaged over the field of view and 200 snapshots of the velocity field.
We use diffusive light imaging to measure the wave motion. The fluid surface is illuminated using a light-emitting diode panel placed underneath the transparent bottom of the container. A few per cent of milk added to water provides sufficient contrast to obtain a high-resolution reconstruction of the wave field. The absorption coefficient is measured before each experiment, which allows calibrating the surface elevation with a vertical resolution of 80 mm. For 3D PTV measurements, floating black carbon glass particles are used to visualise simultaneously the fluid and the wave motions. A few drops of non-ionic surfactant are added to ensure that these particles are homogeneously distributed at the fluid surface before starting an experimental run. For a given set of parameters (paddles acceleration, frequency), no difference could be measured between particle horizontal motion in these experiments and the trajectories of Polyamid particles at the surface of water with no milk and no surfactant.
3D Lagrangian trajectories (see Fig. 4, Supplementary Movie 2 and Supplementary Fig. 2) are retraced using a combination of a 2D PTV technique and a subsequent estimation of the local elevation along the trajectory 21 . First, the horizontal motion (x À y coordinates) of a particle is tracked using threshold filters and a nearest-neighbor algorithm 38 . The wave field evolution is then obtained by removing particles from the movies with local filter techniques. Then the particle elevation (z coordinate) is measured as the wave elevation over a local window (300 mm radius), which is centred on the x À y particle coordinates at a given time.
Finally the 3D trajectories of the particle and the wave field are visualized using the Houdini 3D animation tools (Side Effects Software).
Theoretical effect of the finite wave steepness. Analytically, the drift U d (equation 9) associated with the Eulerian velocity (equation 8) is where Á Á Á ð Þ ¼ 1 Þdt is a gyro-average over one gyration period T ¼ 2p/o. This expression agrees quantitatively with the gyro-averaged velocity observed by numerical integration of the Eulerian equation (8). The parameters L, K, o and d are taken from the experiment.
The small parameter KZ is a measure of the steepness of the wave and is much less than 1 for surface vertical displacements of several millimetres across a 4 cmwide unit cell as considered here. In our modelling, we have also considered the wave/flow coupling to first order in O KZ ð Þ. The modelled surface displacement Z comes from equation (6) with u ¼ u P . To find Z, we substitute equation (8) into equation (6), invert the gradient and retain terms of order O KZ ð Þ. The result is: where Z 1 ¼ L sin ot cos Kx; Z 2 ¼ À L cos ot cos Ky Z 3 ¼ L 2 gK 2 4o 2 cosh 2 Kd ð Þ cos 2 ot cos 2Kx; Z 4 ¼ L 2 gK 2 4o 2 cosh 2 Kd ð Þ sin 2 ot cos 2Ky Z 5 ¼ À L 2 gK 2 2o 2 sin 2ot cos Kx cos Ky: Defining H as the maximum surface displacement, we find that HE ffiffi ffi 2 p L when KLo o1. For a time sequence of Z for parameters relevant to the experiment (L ¼ 1.3 mm, K ¼ 58.1 rad m À 1 , o ¼ 23.9 rad s À 1 , d ¼ 8 cm) and comparison with experimental data, see the Supplementary Figs 1 and 3 and Supplementary Movies 1 and 3. There is a strong match between the experimentally measured wave field dynamics and the predictions of the model.
To take into account higher-order effects of the finite steepness of the waves on the flow, we add to the model a rotational flow u R in addition to the intrinsic Lagrangian drift. Specifically: where b is a rotational velocity amplitude. The horizontal components of u R are chosen to match the form of the drift velocity U d . The introduction of such a term is supported on physical grounds. Indeed we have noted that the model (11) at order O KZ ð Þ predicts a time-averaged surface gradient outward from the centre of each unit cell, that is, a stationary well, providing a centripetal acceleration within the unit cell. A time-averaged remanent level profile has been discussed in the context of the effect of radiation pressure on 2D standing waves 35 . On a related vein, we note that a horizontal drift of a similar form was mentioned in a nonlinear analysis of non-viscous standing surface waves 27 . When the rotation amplitude b is chosen to match the experimental data, u R is of the order of magnitude of the drift U d (see Supplementary Fig. 3 and Supplementary Movie 4).
The particle trajectories shown in Fig. 6c,d were produced by integrating the velocity u ¼ u P þ u R using a 3D, fourth-order Runge-Kutta algorithm. The model parameters used were L ¼ 1.3 mm and b ¼ 0.8456 cm s À 1 for c and L ¼ 0.433 mm and b ¼ 0.4882 cm s À 1 for d. The trajectories were integrated for 2,500 points in time over 37 gyration cycles (c) and 25 gyration cycles (d).
Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.