domain: nlin spi

text:
P É R I N E T

D A M I R J U R I C 2 A N D L A U R E T T

E S T U C K E R M A N

Numerical simulation of Faraday waves

We simulate numerically the full dynamics of Faraday waves in three dimensions for two incompressible and immiscible viscous fluids. The Navier-Stokes equations are solved using a finite-difference projection method coupled with a front-tracking method for the interface between the two fluids. The critical accelerations and wavenumbers, as well as the temporal behaviour at onset are compared with the results of the linear Floquet analysis of Kumar & Tuckerman (J. Fluid Mech., vol. 279, 1994, p. 49). The finiteamplitude results are compared with the experiments of Kityk et al. (Phys. Rev. E, vol. 72, 2005, p. 036209). In particular, we reproduce the detailed spatio-temporal spectrum of both square and hexagonal patterns within experimental uncertainty. We present the first calculations of a three-dimensional velocity field arising from the Faraday instability for a hexagonal pattern as it varies over its oscillation period.

Historical introduction

The Faraday experiment consists of shaking vertically a container holding two immiscible fluids (the lighter of which can be air) thereby inducing oscillations of the fluids and the interface between them. Beyond a certain threshold, the interface can form many kinds of standing wave patterns, including crystalline patterns and others which are more complex. This phenomenon was first studied by [START_REF] Faraday | On a peculiar class of acoustical figures; and on certain forms assumed by groups of particles upon vibrating elastic surfaces[END_REF] who noticed that the vibration frequency of the interface was half that of the forcing. The results of Faraday were confirmed by Rayleigh (1883a,b). [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF] carried out the first theoretical linear analysis of the Faraday waves, restricted to inviscid fluids. They decomposed the fluid motion into normal modes of the container and showed that the evolution equation of each mode reduced to a Mathieu equation whose stability diagram is well known.

In the 1990s, new behaviours of the interface were discovered, such as quasi-crystalline eight-fold patterns seen by [START_REF] Christiansen | Ordered capillary-wave states: quasicrystals, hexagons, and radial waves[END_REF]. By introducing a forcing which is the sum of two periodic functions with commensurable frequencies, [START_REF] Edwards | Patterns and quasi-patterns in the Faraday experiment[END_REF] were able to produce twelve-fold quasi-patterns. Triangular patterns were observed by [START_REF] Müller | Periodic triangular patterns in the Faraday experiment[END_REF] and superlattice patterns by [START_REF] Kudrolli | Superlattice patterns in surface waves[END_REF], also using two-frequency forcing. Spatio-temporal chaos was studied by [START_REF] Kudrolli | Patterns and spatiotemporal chaos in parametrically forced surface waves: a systematic survey at large aspect ratio[END_REF], who also surveyed the occurrence of lattice patternsstripes, squares or hexagons -as a function of viscosity and frequency. Binks, Westra & van de Water (1997) demonstrated the dependence of the pattern on the depth of the layer. In addition to patterns or quasi-patterns, very localized circular waves called oscillons may occur, as seen by [START_REF] Lioubashevski | Dissipative solitary states in driven surface waves[END_REF]. The Faraday instability is the first macroscopic system in which such structures have been observed. These discoveries endow the Faraday instability with a very great fundamental interest for understanding the natural formation of patterns.

A number of theoretical or semi-numerical analyses were inspired by these experiments. [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] extended the linear stability analysis of [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF] to viscous fluids. This analysis was experimentally confirmed by [START_REF] Bechhoefer | An experimental study of the onset of parametrically pumped surface waves in viscous fluids[END_REF] and used by [START_REF] Kumar | Linear theory of Faraday instability in viscous fluids[END_REF] to predict cases in which the response would be harmonic rather than subharmonic. The method was extended by [START_REF] Besson | Two-frequency parametric excitation of surface waves[END_REF] to calculate the stability tongues in the case of two-frequency forcing. Integral equation formulations of the viscous linear stability problem were derived by [START_REF] Beyer | Faraday instability: linear analysis for viscous fluids[END_REF] and [START_REF] Müller | Analytic stability theory for Faraday waves and the observation of the harmonic surface response[END_REF], who also studied the harmonic response case. [START_REF] Cerda | Faraday's instability in viscous fluid[END_REF] used the lubrication approximation and the WKB method to study shallow viscous layers, obtaining a Mathieu equation that was later used by [START_REF] Huepe | Forcing function control of Faraday wave instabilities in viscous shallow fluids[END_REF] to derive analytic results about the response to multifrequency forcing.

Linear analysis provides no information about the shape of the patterns which appear; other means are necessary to understand the occurence of a given pattern or the amplitude of stabilization. Weakly nonlinear approximations have been derived from the Navier-Stokes equations by Viñals and co-workers, e.g. [START_REF] Zhang | Pattern formation in weakly damped parametric surface waves[END_REF], [START_REF] Chen | Amplitude equation and pattern selection in Faraday waves[END_REF] and by [START_REF] Skeldon | Pattern selection for Faraday waves in an incompressible viscous fluid[END_REF], focusing on the competition between different patterns. [START_REF] Vega | Nearly inviscid Faraday waves in annular containers of moderately large aspect ratio[END_REF] derived equations governing the interaction between Faraday waves and the mean flow. There has been a great deal of analysis of lattices, superlattices and quasi-patterns using equivariant dynamical systems theory, as well as model equations designed to produce specific patterns, e.g. [START_REF] Porter | Pattern control via multi-frequency parametric forcing[END_REF]. The approximation of quasipatterns in spatially periodic domains has also been addressed in [START_REF] Rucklidge | Design of parametrically forced patterns and quasipatterns[END_REF].

Investigation of the full nonlinear viscous problem, however, requires numerical simulations, of which there have been very few up to now, specifically those of [START_REF] Chen | Subcritical bifurcations and nonlinear balloons in Faraday waves[END_REF], [START_REF] Chen | Nonlinear wave dynamics in Faraday instabilities[END_REF], [START_REF] Murakami | Two-dimensional direct numerical simulation of parametrically excited surface waves in viscous fluid[END_REF], [START_REF] Valha | A numerical study of the behaviour of a gas-liquid interface subjected to periodic vertical motion[END_REF], [START_REF] Ubal | A numerical analysis of the influence of the liquid depth on two-dimensional Faraday waves[END_REF] and [START_REF] O'connor | The complex spatiotemporal dynamics of a shallow fluid layer[END_REF]. With the exception of [START_REF] O'connor | The complex spatiotemporal dynamics of a shallow fluid layer[END_REF], all previous simulations have been two-dimensional. The most extensive simulation thus far has been that of [START_REF] Chen | Subcritical bifurcations and nonlinear balloons in Faraday waves[END_REF] and [START_REF] Chen | Nonlinear wave dynamics in Faraday instabilities[END_REF], who used a finite-difference method applied to a boundary-fitted time-dependent coordinate system. At each timestep, the surface is advected and a new 2D grid, adapted to the surface, is recomputed. The amplitude of their numerically computed Faraday waves confirmed the weakly nonlinear analysis of [START_REF] Chen | Amplitude equation and pattern selection in Faraday waves[END_REF], including their prediction of a range of subcriticality. Their calculations also predicted qualitatively new phenomena, such as disconnected solution branches and slow modulated dynamics. [START_REF] Murakami | Two-dimensional direct numerical simulation of parametrically excited surface waves in viscous fluid[END_REF] used a method similar to that of [START_REF] Chen | Subcritical bifurcations and nonlinear balloons in Faraday waves[END_REF] and [START_REF] Chen | Nonlinear wave dynamics in Faraday instabilities[END_REF]. Although they reproduced some features of the experiments by [START_REF] Lioubashevski | Dissipative solitary states in driven surface waves[END_REF], their calculations were limited to accelerations only 0.5 % above critical. The investigation by [START_REF] Ubal | A numerical analysis of the influence of the liquid depth on two-dimensional Faraday waves[END_REF] focused on the influence of liquid depth in two-dimensional simulations using a Galerkin finite-element method in transformed coordinates. In addition to comparing their linear stability predictions with those of [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF] and [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF], they calculated instantaneous surface profiles and velocity fields, as well as the temporal evolution and spectrum. [START_REF] Valha | A numerical study of the behaviour of a gas-liquid interface subjected to periodic vertical motion[END_REF] examined the response of a liquid layer in a vertical cylindrical vessel using the MAC (Marker-and-Cell) method of [START_REF] Harlow | Numerical calculation of time dependent viscous incompressible flow of fluid with free surface[END_REF]. Surface tension was treated by the continuum surface force model of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. [START_REF] O'connor | The complex spatiotemporal dynamics of a shallow fluid layer[END_REF] conducted numerical simulations using an ALE (Arbitrary Lagrangian-Eulerian) spectral-element code in both two and three dimensions; a visualization of a square pattern was presented.

The hexagonal patterns, quasi-patterns and oscillons which motivate our investigation are intrinsically three-dimensional, and have never been calculated numerically from the fluid-dynamical equations. Here we report on the results of fully nonlinear, threedimensional simulations of Faraday waves using a finite-difference front-tracking method. In the classic Faraday problem the lighter fluid is usually taken to be air whose effects can be neglected. However, in contrast to the previously cited investigations, the numerical method described here solves the Navier-Stokes equations for the general case of two distinct superposed fluids. The capability of the method to simulate the motion of both fluids is important in that it permits comparison of numerical results with those of certain experimental configurations, namely those of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] where the lighter fluid cannot be ignored. These experiments were the first to provide quantitative measurements of the complete spatio-temporal Fourier spectrum of Faraday waves and thus form an excellent basis for quantitative comparison with our numerical results in the nonlinear finite amplitude regime, i.e. interfaces with steep slopes.

In the next two sections of this article, we present the hydrodynamic equations that govern the Faraday instability and then describe the computational method. The two sections following these are dedicated to the comparison of our results with the linear theory of [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] and with the experiments of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF][START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF]. After comparing numerical and experimental spatio-temporal spectra for squares and hexagons, we present the three-dimensional velocity field for the hexagonal pattern.

Equations of motion

The mathematical model of the Faraday experiment consists of two incompressible and immiscible viscous fluids in a three-dimensional domain x = (x, y, z) ∈ 2 × [0, h], bounded at z = 0 and z = h by flat walls. The two fluids, each uniform and of densities ρ 1 , ρ 2 and viscosities µ 1 and µ 2 , initially form two superposed horizontal layers with an interface between them. This two-dimensional interface is defined by x = (x, y, ζ(x, y, t)). Within the parameter range we wish to simulate, the height ζ remains a single-valued function of (x, y, t).

The container is shaken vertically in z. In the reference frame of the container the boundary conditions for the fluid velocities u = (u, v, w) are u(x, y, 0, t) = 0, (2.1a) u(x, y, h, t) = 0.

(2.1b)

The gravitational acceleration g is augmented by a temporally periodic inertial acceleration G = (a cos(ωt) -g)e z , (2.2)

where a is the amplitude of the forcing and ω is its frequency. The Navier-Stokes equations for incompressible, Newtonian fluids are

ρ Du Dt = -∇p + ρG+∇ • µ ∇u + ∇u T + s, (2.3a) ∇ • u = 0. (2.3b)
Here p is the pressure and s is the capillary force (per unit volume) and is defined below. Equations (2.3a) and (2.3b) are valid for the entire domain, including the interface, in spite of the fact that the density and viscosity change discontinuously and the surface tension acts only at the interface. In this single-fluid formulation, the density and viscosity fields are defined in terms of the densities and viscosities of the two fluids

ρ = ρ 1 + (ρ 2 -ρ 1 )H, (2.4a) µ = µ 1 + (µ 2 -µ 1 )H, (2.4b)
with the aid of a Heaviside function,

H (x -x ) = 0 if z < ζ(x, y, t) 1 if z ζ(x, y, t) , (2.5)
where we recall that x = (x, y, z) is a point anywhere in the three-dimensional volume and x = (x, y, ζ(x, y, t)) is the vertical projection of x onto the interface. The capillary force is

s = S (t) σκ n δ (x -x ) dS, (2.6)
where σ is the surface tension coefficient, assumed to be constant, n is the unit normal to the interface (directed into the upper fluid) and κ its curvature. δ (x -x ) is a threedimensional Dirac distribution that is nonzero only where x = x . S (t) is the surface defined by the instantaneous position of the interface. To complete the system of equations we need an expression for the motion of the interface. One such expression can be easily derived by noting that mass conservation in an incompressible flow requires Dρ/Dt = 0, which in view of the discontinuous density field (2.4a), is equivalent to DH/Dt = 0.

(2.7)

Thus the interface is represented implicitly by H and advected by material motion of the fluid.

Computational methods

The computational domain is a rectangular parallelepiped, horizontally periodic in x and y and bounded in z by flat walls for which we impose no-slip boundary conditions. The entire domain is discretized by a uniform fixed three-dimensional finite-difference mesh. This mesh has a standard staggered MAC cell arrangement [START_REF] Harlow | Numerical calculation of time dependent viscous incompressible flow of fluid with free surface[END_REF] where the u, v and w velocity nodes are located on the corresponding cell faces and scalar variables are located at the cell centres. Each cell is of dimension ∆x×∆y×∆z.

Within the domain, the two distinct immiscible fluids are separated by a two-dimensional interface which is discretized by a second mesh as sketched in figure 1. This moving and deformable mesh is composed of triangular elements whose motion is treated by a fronttracking/immersed-boundary method [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF][START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]). Because we have assumed that ζ(x, y, t) is single valued, the nodes of the mesh can be fixed in x and y and only their vertical displacements need to be calculated, which is a considerable simplification to the general front-tracking method.

After setting appropriate initial and boundary conditions, the computational solution algorithm for each timestep is composed of three main phases. First, the interface is advected and the density and viscosity fields updated according to the new interface position. The capillary force s is then calculated. Finally, the velocity and pressure are found by means of a standard projection method. Each of these steps is described below.

Advection of the interface

Purely Eulerian interface methods such as Volume of Fluid [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] or Level Set [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]) use a form of (2.7) to advect a scalar field such as H, or a level-set function that implicitly represents the interface. However, in the fronttracking approach that we use here, the interface markers themselves (the nodes of the triangular mesh) are advected. H is then constructed from the position and geometry of the interface. Taking (2.7) as a starting point, we develop an equivalent expression for the vertical displacement of the triangular interface mesh.

The material derivative of (2.5) gives:

DH Dt = ∇H • Dx Dt + ∇ H • Dx Dt , (3.1)
where ∇ = ∂ x and where I is the identity tensor. The specific choice of x = (x, y, ζ(x, y, t)) made here gives for ∇x

∇H = -∇ H = S (t) nδ(x -x )dS. ( 3
    1 0 0 0 1 0 ∂ζ ∂x ∂ζ ∂y 0     . (3.6)
With this, (3.5) leads to the specific displacement relations:

∂x/∂t = 0, (3.7a) ∂y/∂t = 0, (3.7b) ∂ζ ∂t = w -u ∂ζ ∂x -v ∂ζ ∂y . (3.7c)
The application of this advection to the triangular mesh we use for tracking the interface is straightforward. At each vertex the horizontal displacement is zero and for the vertical displacement we compute a first-order approximation to (3.7c):

ζ n+1 -ζ n ∆t = w n (x e ) - ∂ζ n ∂x u n (x e ) - ∂ζ n ∂y v n (x e ).
(3.8)

The superscripts n and n + 1 denote, respectively, the old and new time levels. The derivatives on the right-hand side are evaluated using a simple upwind scheme which requires the usual CFL (Courant-Friedrichs-Lewy) time step restriction. The vertical displacement of the interface mesh requires knowledge of the velocities at the element nodes x e . These in general do not coincide with the Eulerian grid nodes x ijk , whose indices correspond to discretized coordinates along the respective directions x, y and z.

The problem of communicating Eulerian grid velocities to the element nodes is overcome by interpolation between the two grids as is typically done in front-tracking and immersed-boundary methods. Here we use the particular interpolation

u(x e ) = ijk u(x ijk )δ h (x ijk -x e )∆x∆y∆z.
(3.9)

The kernel δ h is a smoothed version of the three-dimensional Dirac delta function with compact support of four grid nodes in each direction (for details of the front-tracking method, see [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF], and for the immersed-boundary method, see [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF]. In (3.9) the weighted information collected from nearby Eulerian grid nodes is interpolated to a given element node.

We now seek to update the density and viscosity fields needed in (2.4), which require H. The equation for H, based on the updated values of x and n, is formulated by taking the divergence of (3.2): The discretized version of the Poisson problem (3.10a) is

∇ 2 H = ∇ • S (t) n δ (x -x ) dS, (3.10a) H(x, y, 0, t) = 0, (3.10b) H(x, y, h, t) = 1. (3.10c)
∇ 2 H ijk = ∇ • e n e δ h (x ijk -x e )∆S e , (3.11)
where standard central differencing is used for the gradient and divergence operators. This numerically calculated Heaviside function is a smoothed transition from 0 to 1 across a distance of 4 grid cells in the direction normal to the interface. In contrast to (3.9), the summation above serves to distribute weighted information from an element node to nearby Eulerian grid nodes. Since an element is triangular, its vertices lie in the same plane, its normal vector is unique and the three tangent vectors are simple to calculate:

n e = t 2 × t 1 ||t 2 × t 1 || (3.12)
where t 1 and t 2 are the tangents on two distinct edges of the triangle (see the sketch in figure 2). We solve (3.11) by fast Fourier transform. Finally, ρ n+1 and µ n+1 are updated using (2.4a) and (2.4b).

Capillary force

From (2.6), the capillary force involves the curvature of the interface and its normal vector. However, from a computational point of view, curvature is a difficult quantity to compute accurately. It is more accurate and physically appealing to calculate the force pulling on the edge of each individual triangular surface element and then sum the contributions for all the elements over the surface. For a given surface element e of surface area δA and perimeter δl, we can write:

s e = σ δA κn dA, = σ δA (n × ∇) × n dA, (3.13) = σ δl t × n dl,
where the last integral represents the sum of the capillary forces exerted around the element perimeter. As sketched in figure 2, the directions of these forces are oriented along the surface and normal to the element's edges. The net force at the shared edge between any two triangles (the sum of the solid and dashed vectors) is directed into the fluid on the concave side of the interface. Following Peskin's immersed boundary method [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF], the discrete version of (2.6) becomes

s ijk = e s e δ h (x ijk -x e ), (3.14)
where we use the same smoothed δ h as in (3.9) and (3.11). Thereby, several interfacial elements contribute to the calculation of the force applied to a single Eulerian node, and a single element influences more than one Eulerian node.

Solution of the Navier-Stokes equations

The Navier-Stokes equations are solved by a projection method [START_REF] Chorin | Numerical simulation of the Navier-Stokes equations[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF]) with incremental pressure correction [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows[END_REF]) applied to a finite-difference scheme which is first order in time and second order in space. In addition a semi-implicit scheme is chosen for the velocities to relax the stability restriction on the time step due to viscous diffusion. All spatial derivative operators are evaluated using standard centred differences, except in the nonlinear term where we use a second-order ENO (Essentially-Non-Oscillatory) scheme [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes[END_REF][START_REF] Sussman | An improved level set method for incompressible two-phase flows[END_REF]. (For an overview of projection methods for the incompressible Navier-Stokes equations, see [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] The time stepping algorithm is thus

u n+1 -u n ∆t = -u n • ∇u n + 1 ρ n+1 ∇ • µ n+1 ∇u + ∇u T n+1 - 1 ρ n+1 ∇p n+1 + s n+1 ρ n+1 + G n+1 , (3.15)
with the boundary conditions on the top and bottom walls

u n+1 Γ = 0. (3.16)
In (3.15), ρ, µ and s depend on x via (2.4-2.6) and have already been updated by (3.8).

We decompose the solution of (3.15) in three steps. The first step is a semi-implicit calculation of an intermediate unprojected velocity ũ, involving only velocities and their gradients:

ũ -u n ∆t = -u n • ∇u n + 1 ρ n+1 ∇ • µ n+1 ∇ũ + ∇ũ T (3.17a) ũ| Γ = 0 (3.17b)
In the second step, we include the capillary, acceleration and old pressure gradient terms to calculate the unprojected velocity u * :

u * - ũ ∆t = G n+1 + s n+1 ρ n+1 - 1 ρ n+1 ∇p n .
(3.18)

Finally we perform a projection step to find the divergence free velocity u n+1 :

u n+1 -u * ∆t = - 1 ρ n+1 ∇ p n+1 -p n , (3.19a) ∇ • u n+1 = 0, (3.19b) u n+1 • n Γ = 0. (3.19c)
Equations (3.19) imply the following elliptic problem for the pressure increment

∇ • u * ∆t = ∇ • 1 ρ n+1 ∇ p n+1 -p n , (3.20a) ∂ p n+1 -p n ∂n Γ = 0, (3.20b)
which we solve with an iterative biconjugate gradient stabilized algorithm [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]).

In the horizontal directions, periodic boundary conditions are imposed on the velocity and pressure, thus excluding a net horizontal pressure gradient. For the simulations we will present here, this choice is consistent with the requirement of no mean horizontal flux in a large but bounded container. We note that in the implicit solution of (3.17), we apply the same biconjugate gradient stabilized solver used for the pressure to each component of ũ = (ũ, ṽ, w) separately. Thus only the diagonal terms of the diffusion operator are treated fully implicitly. The off-diagonal terms are treated quasi-implicitly in that the newest available values of (ũ, ṽ, w) are used in the evaluation of the cross derivatives. To ensure symmetry, we permute the order of solution for each component.

Results: linear analysis

Floquet analysis

In the absence of lateral boundaries, the equations are homogeneous in the horizontal coordinates and the solutions can be represented by a spatial Fourier transform:

ζ(x, y, t) = k ζ(k, t)e ik•x (4.1)
The linear instability of the interface between two fluids is described by (2.2)-(2.7) linearized about a zero velocity field and flat interface ζ = ζ . The linearized equations depend only on the wavenumber k ≡ ||k|| of each wave and not on its orientation and hence the coefficient of e where T = 2π/ω, γ is real and α ∈ [0, 1[. The Floquet modes,

f (k, t mod T ) = ∞ n=-∞ f n (k)e inωt , (4.3)
are not trigonometric, but remain periodic with fundamental frequency ω. Thus, the linearized behavior for a single mode is

ζ(x, y, t) = e ik•x e (γ+iαω)t ∞ n=-∞ f n (k)e inωt . (4.4)
Analogous expressions hold for the velocity u. Equation (4.2) shows that if γ is non-zero or α is irrational, the evolution of the interface motion is not periodic. A non-zero γ indicates that the motion grows or decays according to the sign of γ. An irrational α yields a quasi-periodic evolution function. For the Faraday instability, it can be shown that α can only take two values: 0 and 1/2. As the imposed acceleration a is increased, one encounters regions in the (k, a) plane in which γ > 0 for one dominant temporal frequency, jω/2, where j = 1, 2, 3, ... (see figure 3). Within each instability tongue, the amplitude of the mode grows exponentially. These tongues are called harmonic if α = 0 and subharmonic if α = 1/2. As k is increased, one encounters an alternating sequence of subharmonic and harmonic tongues, which are bounded by neutral curves (k, a c (k)) on which γ = 0. On the neutral curves, the solutions are periodic:

ζ(k, t) = ∞ n=-∞ f n (k)e i(n+ 1 2 )ωt , subharmonic case, (4.5a) ζ(k, t) = ∞ n=-∞ f n (k)e inωt , harmonic case. (4.5b)

Computation of the neutral curves

We first compare our numerically calculated instability thresholds with those found by [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] for the same parameter values. The physical parameters are ρ 1 = 519.933 kg m -3 and µ 1 = 3.908 × 10 -5 Pa s for the lower fluid and ρ 2 = 415.667 kg m -3 and µ 2 = 3.124 × 10 -5 Pa s for the upper fluid. The other parameters are σ = 2.181 × 10 -6 N m -1 and g = 9.8066 m s -2 . The frequency of the forcing is ω/2π = 100 Hz and thus its period is T = 0.01s. The capillary length is defined as

l c = σ/(|ρ 1 -ρ 2 |g).
The container height is taken to be 5l c = 0.231 mm, and the interface, when unperturbed, is equidistant from the top and bottom boundaries. We consider several wavenumbers k and set the x dimension of the box in each case to one expected wavelength λ = 2π/k, i.e. to between 0.074 and 0.224 mm, as listed in Table 1. We can estimate the importance of various physical effects for these parameters by defining dimensionless quantities with length k -1 and forcing period T . The Bond number Bo = (kl c ) -2 = |ρ 1 -ρ 2 | g/(σk 2 ) measures the relative importance of gravitational to capillary effects and ranges between 0.0649 and 0.598. The Reynolds number Re = ρ/(µk 2 T ) is a nondimensional measure of viscous damping and ranges between 0.184 and 1.70 for both fluids.

We have computed the critical acceleration from our fluid-dynamical simulation for the wavenumbers listed in Table 1. Initially, the interface is sinusoidal with wavevector k parallel to the x-axis and the velocity is zero. Moreover, to ensure that the solution corresponds initially to the linear solution, we require the amplitude of the interface displacement to be small compared to λ. In order to maintain a roughly cubic mesh and a minimum x-resolution of about 50 grid cells per wavelength, we vary the resolution in the z direction between 126 and 144 points. Since k points along the x direction, ζ does not depend on y (neither do the velocity nor the pressure) and so the size of the domain and resolution in y are arbitrary. The acceleration a is taken near a c (k), the expected critical acceleration corresponding to each wavenumber. At the threshold, the flow undergoes a pitchfork bifurcation. Since the growth rate is proportional to a -a c close to the neutral curve, it is sufficient to find the growth rates for two values of a and to interpolate linearly between these points.

In figure 3, we plot the values of a c obtained from our fluid-dynamical simulation for several values of k, along with the curves (k, a c (k)) obtained from the method of [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]. Figure 3 shows that these thresholds are in good agreement, despite whatever inaccuracies in a c are introduced by spatial discretization and linear interpolation. The relative error in the critical acceleration at the conditions previously stated is of the order of a few per cent as shown in Table 1. The results suggest that there is a k below which calculation of the growth rates is not possible. Some zones of the diagram are not accessible because a domain of width 2π/k necessarily accommodates all wavenumbers which are integer multiples of k up to the resolution limit π/k∆x. The coefficients of the Fourier expansion of the initial condition differ slightly from zero due to finite-difference spatial approximations and if the growth rate of one of these is greater than that of k itself, then it will quickly come to dominate k. This difficulty is exacerbated by the fact that several forcing periods are required for γ to stabilize. Then the amplitude, whose evolution was expected to be almost periodic, starts to rise before the precise determination of γ is possible, for example in the range of k between 0 and roughly 15 mm -1 . As we see in figure 3, the critical forcing is substantially lower for one of its multiples closer to 32.5 mm -1 . The amplitude corresponding to this wavelength, although initially negligible, increases and rapidly dominates the mode we wish to study, making the calculation of γ unfeasible. In contrast, for k = 48 mm -1 , the growth rate did not vary significantly after having reached a value near zero (relative fluctuations of about 0.1% of the growth rate's limit value were recorded after the stabilization).

Temporal profile of a mode

We recall from section 4.1 that the time dependence of a Floquet mode is not sinusoidal. As a further validation, we can compare the results of our fluid-dynamical simulations to the entire temporal behavior over a period. This is a stronger validation than merely predicting the threshold since it provides a comparison at every time instead of once per period.

In figures 4-6, we plot the deviation ζ -ζ from the flat interface as a function of time at a fixed spatial location from our fluid-dynamical simulation, for values k = 48, 60.9 and 85 mm -1 belonging to the first three tongues. On the same figures, we plot the behavior of (4.5), where the temporal coefficients f n (k) of the Floquet modes have been calculated by the method in [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]. The value of a is set to the interpolated critical acceleration a c (k), so the oscillations approximately retain their initial amplitude as long as they remain small. The comparisons in figures 4-6 show a nearly perfect agreement. The differences observed initially, due to the phase difference between the initial conditions and acceleration, vanish remarkably quickly, in well under one period of oscillation of the container.

Figures 4-6 correspond to tongues jω/2, with j =1, 2, 3, respectively, which show j zero crossings per forcing period T . Odd (even) values of j correspond to subharmonic (harmonic) oscillations, with period 2T (T ). The temporal spectrum f n (k) becomes richer as k increases, leading to increasingly more complex modes, as can be observed by comparing figures 4-6. This strong anharmonicity of the curves is due to the increasing contribution of higher frequency trigonometric functions to the Floquet modes as a increases. The Floquet mode corresponding to k c = 32.5 mm -1 , with the smallest value of a = a c , should be closer to trigonometric, with a fundamental frequency of ω/2.

Results: nonlinear analysis

In the full nonlinear evolution of the interface for a > a c , the amplitude of the interface height grows in time until nonlinear terms in (2.2)-(2.7) become important. After that, the mode whose linear growth rate is maximal gives rise, via nonlinear resonances, to a series of other discrete modes, selected according to the magnitudes and orientations of their wavevector k. This selection is responsible for the formation of patterns that will be the object of our further validations. We seek to compare our calculations with the experimental results of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF][START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] where quantitative data concerning the Fourier spectrum ζ(k, t) are available for squares and hexagons.

We run our numerical simulations with the same experimental parameters as [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF]: ω/2π = 12 Hz (T = 0.0833 s), ρ 1 = 1346 kg m -3 , µ 1 = 7.2 mPa s for the lower fluid and ρ 2 = 949 kg m -3 , µ 2 = 20 mPa s for the upper fluid. The surface tension at the interface is σ = 35 mN m -1 , the total height of the vessel is 1.0 cm and the mean height of the interface, the initial fill height of the heavy fluid, is ζ = 1.6 mm (with some uncertainty; see below). The Floquet analysis for these parameters yields a critical wavelength of λ c = 2π/k c = 13.2 mm and a critical acceleration of a c = 25.8 m s -2 . Here, the Bond number defined in section 4.2 is Bo = |ρ 1 -ρ 2 | g/(σk 2 c ) = 0.49. The Reynolds number Re = ρ/(µk 2 c T ) is Re 1 = 9.9 and Re 2 = 2.52 for the lower and upper fluid, respectively.

Rather than starting from a sinusoidal interface, we chose to add two-dimensional white noise of small amplitude to ζ to define the initial interface height ζ(x, y, t = 0) in order to excite every mode allowed by the box's horizontal dimensions and number of cells. It is thus possible to check that the correct critical mode (that whose growth rate is maximal) emerges from the linear dynamics. In order to reproduce the experimental results in a computational domain of a minimal size, the dimensions in x and y of the box must correspond to the periodicity and symmetries of the expected pattern. The minimal required resolution along these directions has been found to be between 40 and 50 cells per wavelength. The number of triangles used to represent the interface is 16 times the total number of horizontal gridpoints. The number of cells in the z direction is taken so that min S ζ(x, y, t) is greater than about the height of 3-5 cells. The required vertical resolution thus varies with the forcing amplitude. The initial velocity is taken to be zero.

Square patterns

To compare with the experiment of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] for their square patterns, we choose the same forcing acceleration, a = 30.0 m s -2 . Our box has horizontal dimensions which we take both equal to 2π/k c . The timestep is ∆t = 2.78 × 10 -4 s. Figures 7 and8 represent examples of the patterns obtained at saturation under these conditions and are taken from the same simulation at the two instants shown by the two arrows in figure 10. The symmetries characterizing the squares (reflections and π/2 rotation invariance) are clear, showing a first qualitative agreement with [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] where both structures were observed. The pattern oscillates subharmonically, at 2T , where T is the forcing period. Figure 7 is taken when the interface attains its maximum height, while figure 8 is taken at a time 0.24 × 2T later. At this later time, we observe the dominance of a higher wavenumber, which will be discussed below.

Further quantitative investigations of the patterns involve the spatial Fourier transform of the interface height. In the case of square patterns, the distribution of the spatial modes is shown in figure 9. The modes with non-negligible amplitude are ±k c e x and ±k c e y , with |k| = k c and amplitude A(k c ); ±2k c e x and ±2k c e y , with |k| = 2k c and amplitude A(2k c ); and k c (±e x ± e y ), with |k| = √ 2k c and amplitude A( √ 2k c ). (For a square pattern, the amplitude of each mode is identical to that of each of its images through rotation by any integer multiple of π/2.) The interface height is written as:

ζ(x, t) = ζ + A(k c , t) 4 j=1 e ikcej •x + A(2k c , t) 4 j=1 e i2kcej •x + A( √ 2k c , t)
4 j=1 e i √ 2kce j •x + higher order terms, (5.1) where e j ≡ e x cos(πj/2) + e y sin(πj/2) and e j ≡ e x cos(π/4 + πj/2) + e y sin(π/4 + πj/2) for j = 1, . . . 4. We have chosen this notation, rather than ζ(k, t) as used in equation (4.1), to facilitate comparison with [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF][START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF].

We have compared the evolution of the three principal spatial modes (figure 10) and their temporal Fourier transform (figure 11) with the experimental results [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF]). Here we turn the reader's attention to the recent erratum by [START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] for correct quantitative comparisons of the spectra.

Figure 10 compares the experimental evolution of each spatial wavenumber to numer- The main features found in [START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] are recognized in figure 10. In particular, both the fundamental periodicity of each mode (harmonic or subharmonic) and the form of each numerical curve in figure 10 are very similar to the experimental data. The amplitudes and the phases are also quite close. Most of the experimental amplitudes are bracketed by the numerical ones. Thus, they lie in the interval of amplitudes allowed by the range of uncertainties which is surely underestimated since only the uncertainty in ζ has been taken into account. A(k c ) crosses zero at times different from the two higher wavenumbers, A(2k c ) and A( √ 2k c ). At these instants, the higher wavenumbers dominate the pattern. In particular, the pattern of figure 8, taken near the second arrow in figure 10,when A(k c ) is low, contains more peaks than that of figure 7, taken when A(k c ) is high. The large ratio between the amplitude of k c and the others makes this phenomenon very short-lived.

Figure 11 shows the temporal Fourier decomposition of the curves in figure 10. These spectra for the experiment [START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] and for the computation are quite similar too. All of the square patterns that we have observed, once saturation is attained, remain so for the entire duration of the calculation. We present a brief numerical grid convergence study in figure 12. All qualitative features, such as the square symmetry, were observed with each of the three resolutions chosen, despite the coarseness of the 20 × 20 × 40 and 40 × 40 × 80 grids. With increasing resolution, the principal spatial modes converge to the experimental curves shown in figure 10, with only a small difference between the curves with the two highest resolutions. The order of numerical convergence of the maximum and minimum of the amplitudes of each of the three modes in figure 12 shows that the convergence is between first and second order, which is expected to be the case with the immersed-boundary method. In particular, we would expect that a further doubling of the resolution would change the results by at most 4 % for the principal k c mode.

Hexagonal patterns

When the amplitude of the forcing acceleration a is further increased, the modes can reorganize. The symmetries change and, in the experiments of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF], the initial square pattern becomes hexagonal. Though k c remains constant, the horizontal dimensions of the minimal computational box necessary to support the periodic pattern must change too. These dimensions become 4π/k c in y and 4π/( √ 3k c ) in x, as shown in figure 13. The wavevector lattice for hexagonal patterns is shown in figure 14. The principal modes are again of three amplitudes: k c , 2k c and √ 3k c . When a pattern is hexagonal, a mode will have the same amplitude and temporal behaviour as each of its images through rotations by any integer multiple of π/3. The interface height is thus where e j ≡ e x cos(πj/3) + e y sin(πj/3) and e j ≡ e x cos(π/6 + πj/3) + e y sin(π/6 + πj/3). for j = 1, . . . , 6. Our simulations are carried out at acceleration a = 38.0 m s -2 and mean height ζ = 1.6 mm. We have used two different initial conditions: a rectangular pattern, and also white noise, as in our previous simulations of the square patterns. In both cases, hexagons emerge and saturate. The results shown below are those that emerge from the white noise. The time step varies during the calculation, depending on the viscous diffusion limit and the CFL. The spatial resolution is 58x100x180 in the x, y, z, directions, respectively.

ζ(x, t) = ζ + A(k c , t) 6 j=1 e ikcej •x + A(2k c , t) 6 j=1 e i2kcej •x
In figures 15-18, we show visualizations of the patterns at four instances in time. A movie of the temporal evolution of the hexagon pattern over one subharmonic oscillation is available in the online version of this article. The π/3 rotational symmetry confirms that the rectangular numerical grid does not forbid the formation of hexagonal patterns, which are not aligned with this grid. The patterns reproduce several prominent features from the visual observations of hexagons in the experiments. For example, one can observe in figures 15 and 17 the up and down hexagons shown in the experimental snapshots (figure 10 of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF]. The pattern in figure 18, when the surface elevation is minimal, is dominated by wavenumbers higher than k c , as is also the case in figure 10 of [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF]. This is reflected by the disappearance of A(k c ) and the resulting The spectra from experiments and simulations are represented in figures 19 and 20. Given that experimental uncertainties concerning the hexagons are greater than for the squares (A. Kityk & C. Wagner, private communication), the agreement is remarkable. The principal mode is well reproduced while the other two modes show rough agreement. It is striking that, in contrast to square patterns, every wavevector is a superposition of harmonic and subharmonic temporal modes, so that each has temporal period 2T . This phenomenon was explained by [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] as a spatio-temporal resonance as follows. In the case of the square lattice, two critical subharmonic modes (e.g. k c e x and k c e y ) interact to yield a higher wavenumber harmonic mode (e.g. k c (e x + e y )). In the hexagonal case, two critical subharmonic modes (e.g. -k c e y and k c ( √ 3e

x -e y )/2) interact to yield a higher wavenumber harmonic mode (k c ( √ 3e x -3e y )/2), as in triangle I of figure 14. Further interaction of this mode with a critical subharmonic mode (k c ( √ 3e

x + e y )/2) yields subharmonic contributions to the higher spatial wavenumber mode (k c ( √ 3e

x -e y )), as shown in triangle II. Other quadratic interactions between critical subharmonic modes can contribute to a third harmonic mode of wavenumber k c (triangle III).

In addition to the interface height, our simulations also produce the entire velocity field, which is the focus of figures 21-23. These figures show the velocity fields on horizontal planes at three instants spanning the oscillation period of a hexagonal pattern, as well as the vertical velocity on the interface. Figures 21,22 been decreased to 36.0 m s -2 , and the number of triangles used to represent the interface has been increased to 64 times the total number of horizontal gridpoints. Figure 21 is taken at t = 0.07×(2T ), just after the interface reaches its maximum height (at t = 0), when the peaks are beginning to descend. Consequently, the fluid converges horizontally towards the interface peaks, then descends dramatically below them. The that the rims of the wide flat craters seen in figure 22 have in turn collapsed inwards, forming circular waves which invade the craters, whose remnants are visible as dimples. The velocity field of figure 23 shows fluid converging horizontally below these dimples. These are erupting at velocities which are the largest in the cycle, and will eventually reconstitute the high peaks seen in figure 21.

Figures 21-23, as well as figures 15-18, show that these cases pose great computational difficulties. The interface periodically forms a very thin film (approximately 0.1 mm; see wide crater in figures 17 and 22) over large portions of the lower boundary, within which the velocity may be significant. These features make it difficult to adequately resolve the flow in this layer. For the time being we use a uniform grid spacing; however, in this case an adaptive grid would be more efficient and is under development. We have also simulated hexagons with resolutions of 70 × 120 × 100 and 70 × 120 × 50. Although we do not show these, the two highest resolutions lead to very similar spatial spectra with a maximum difference in amplitudes of the principal modes of about 5% between the 70 × 120 × 100 and 58 × 100 × 180 resolutions. The case resolved by only 50 cells in the z direction shows differences mainly in the 2k c mode where the difference between the 70 × 120 × 50 and 58 × 100 × 180 resolution is about 25%; for the two other modes the difference is about 10%. Hexagonal motifs were observed for all of the resolutions.

The calculation for the hexagon case, for the resolution of 58 × 100 × 180 takes about 7 h per subharmonic oscillation on a 2.16 GHz Intel processor. This corresponds to 42 h of calculation time for 1 s of physical time.

In contrast to the square patterns, all of the hexagonal patterns that we have observed are transient. In our calculations, they last for several seconds, i.e. about 15-20 subharmonic oscillation periods, over which time the amplitudes and periods of the principal modes remain constant. This is also the case for the experimental observations (A. Kityk & C. Wagner, private communication), although the experimental lifetimes are longer. In our simulations, hexagonal patterns alternated with patterns with other symmetries, whose lifetimes were long (on the order of several seconds) but irregular. This behaviour suggests that the hexagonal state may belong to a heteroclinic orbit. A more extensive examination of the hexagonal regime will be the subject of a future investigation.

Conclusion

We have carried out full nonlinear three-dimensional simulations of Faraday waves. The incompressible Navier-Stokes equations for two fluid layers of different densities and viscosities are solved using a finite-difference method. The interface motion and surface tension are treated using a front-tracking/immersed-boundary technique. The simulations are validated in several ways. First, for small oscillation amplitudes, our computations match the solution of [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] to the Floquet problem which results from the linearized evolution equations. The boundaries of the instability tongues, i.e. the critical amplitude as a function of horizontal wavenumber are calculated for several wavenumbers on several tongues and are in good agreement with the theoretical values. The temporal dependence of the Floquet modes is also well reproduced by our numerical results, an even more quantitatively significant validation.

For finite oscillation amplitude, our computations reproduce the square and hexagonal patterns observed by [START_REF] Kityk | Spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF][START_REF] Kityk | Erratum: spatiotemporal characterization of interfacial Faraday waves by means of a light absorption technique[END_REF] at moderate and high-oscillation amplitudes, respectively. Although the domains shown in figure 13 were chosen to accommodate square and hexagonal patterns respectively, we consider the emergence of these patterns at the appropriate parameter values a non-trivial test of our program, since these domains can also accommodate rectangles and stripes. Quantitative comparisons were made between experiment and simulation of the spatio-temporal spectra. Our numerical results lie well within the experimental uncertainty. The hexagonal patterns are long-lived transients and show intriguing dynamical behavior. Our direct numerical simulations provide velocity fields and pressure throughout the entire domain of calculation. Thus, we have been able to ascertain precisely the fluid motion for the Faraday waves, both above and below the interface between the two fluids.

Our future studies of Faraday waves will include a more detailed investigation of the dynamics of the hexagonal patterns, and the simulation and interpretation of oscillons.

Figure 1 .

 1 Figure 1. Spatial discretization of the domain.

Figure 2 .

 2 Figure 2. A triangular element of the interface mesh illustrating the action of the capillary forces according to (3.13). For the shaded triangle the forces act perpendicular to the triangle's edges (unlabelled solid arrows), the dashed arrows are the corresponding forces on the edges of neighboring triangles. The net capillary force at the shared edge between any two triangles (the sum of the solid and dashed vectors) is directed into the fluid on the concave side of the interface.

Figure 3 .

 3 Figure 3. Critical acceleration ac/g as a function of the wavenumber k. The solid curves represent the neutral curves obtained by[START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]. The ac found with the simulation are indicated by the circles.

Figure 4 .

 4 Figure 4. Linear evolution of the surface height deviation ζ(t) -ζ for k = 48 mm -1 , in the first instability tongue. Our simulation results are plotted with symbols and those derived from a Floquet analysis with the solid line. The height and time are nondimensionalized by the wavelength λ = 2π/k and forcing period T , respectively.

Figure 5 .

 5 Figure 5. Linear evolution of the surface height deviation ζ(t) -ζ for k = 60.9 mm -1 , in the second instability tongue. Same conventions used as figure 4.

Figure 6 .

 6 Figure 6. Linear evolution of the surface height deviation ζ(t) -ζ for k = 85 mm -1 , in the third instability tongue. Same conventions used as figure 4.

Figure 7 .

 7 Figure 7. Example of square pattern. Height of interface as a function of the horizontal coordinates, at the instant corresponding to first arrow of figure 10, when height is maximal. Resolution in x, y, z directions: 80 × 80 × 160. Note that the vertical scale is stretched with respect to the horizontal scale. Each horizontal direction in the figure is twice that of the calculation domain.

Figure 8 .

 8 Figure 8. Example of square pattern, at the instant corresponding to second arrow of figure 10, time 0.24 × 2T after figure 7.

Figure 9 .

 9 Figure 9. Lattice formed by the spatial modes comprising a square pattern. The principal modes, with wavenumbers kc, 2kc and √ 2 kc, whose evolution will be studied in figures 10 and 11 are indicated by hollow black circles.

Figure 10 .

 10 Figure 10. Temporal evolution of the amplitudes of the spatial modes with wavenumbers kc, 2kc and √ 2kc. Solid curves represent the experimental results of Kityk et al. (2009), dashed curves and crosses represent numerical results for ζ = 1.6 mm and ζ = 1.7 mm, respectively. Resolution in x, y, z directions: 80 × 80 × 160. Arrows, from left to right, show the time at which figures 7 and 8 have been plotted.

Figure 11 .

 11 Figure 11. Temporal Fourier transform of the amplitudes in figure 10. Circles indicate experimental results of Kityk et al. (2009), while crosses and plus signs indicate numerical data with ζ = 1.6 mm and ζ = 1.7 mm respectively.

Figure 12 .

 12 Figure 12. Temporal evolution of the amplitudes of the spatial modes with wavenumbers kc, 2kc and √ 2kc for square patterns. Study of the convergence with three different spatial resolutions. Circles indicate a resolution (in x, y, z directions) of 20 × 20 × 40, dashed curves 40 × 40 × 80 and continuous curves 80 × 80 × 160. The timestep is the same, ∆t = 2.78 × 10 -4 s, for all curves shown.

Figure 13 .

 13 Figure 13. Boxes supporting the periodic patterns in the square and hexagonal cases. In black, the borders of the box. Light lines, pattern contained by each box; λ = 2π/kc.

Figure 14 .

 14 Figure 14. Lattice formed by the spatial modes comprising a hexagonal pattern. The principal modes, with wavenumbers kc, 2kc and √ 3 kc, involved in later quantitative investigations are indicated by hollow black circles. The labelled triangles illustrate resonance mechanisms leading to harmonic contributions to higher wavenumbers.

  and 23 correspond approximately to the visualizations of figures 15, 16 and 18, where the structures are more visible since the interface has been repeated periodically in the horizontal directions for clarity. The parameters are the same as those given previously, except that the acceleration a has

Figure 15 .

 15 Figure15. Snapshot of hexagonal pattern, taken when height of the interface peaks is maximal. Each horizontal direction is twice that of the calculation domain.

Figure 16 .

 16 Figure 16. Snapshot of hexagonal pattern taken t = 0.3 × 2T after instant of maximal interface height.

Figure 17 .

 17 Figure 17. Snapshot of hexagonal pattern taken t = 0.48 × 2T after instant of maximal interface height.

Figure 18 .

 18 Figure 18. Snapshot of hexagonal pattern taken t = 0.68 × 2T after instant of maximal interface height.

Figure 19 .

 19 Figure 19. Temporal evolution of the amplitudes of the spatial modes with wavenumbers kc, 2kc and √ 3kc. Solid curves represent experimental results (A. Kityk & C. Wagner, private communication) at a ≈ 38.5 m s -2 . Dashed curves represent the simulation for a = 38.0 m s -2 at resolution (in x, y, z directions) of 58 × 100 × 180. Arrows indicate times corresponding to figures 15-18.

Figure 20 .

 20 Figure 20. Temporal Fourier transform of the amplitudes in figure 19. Circles represent experimental results (A. Kityk & C. Wagner, private communication) for a ≈ 38.5 m s -2 . Crosses represent numerical results for a = 38.0 m s -2 at resolution 58 × 100 × 180.

Figure 21 .

 21 Figure 21. Velocity field at time t = 0.07 × (2T ) after the instant of maximum height. Interface is colored according to the vertical velocity w. Arrows show velocity field at z = 0.53 mm and z = 6.08 mm. (Total height is 10 mm, average interface height is 1.6 mm.) For clarity, velocity vectors are plotted only at every fourth gridpoint in each direction. Note that the vertical and horizontal scales are different. One computational domain is shown.

Figure 22 .

 22 Figure 22. Velocity field at time t = 0.41 × (2T ) after the instant of maximum height. Vectors shown at z = 0.083 mm and z = 6.25 mm. Vector and color scales differ from those of figure 21.

Figure 23 .

 23 Figure 23. Velocity field at time t = 0.73 × (2T ) after the instant of maximum height. Arrows show velocity field at z = 1.58 mm. Vector and color scales differ from those of figures 21 and 22.

Table 1 .

 1 Comparison of the computed ac with Floquet theory for various wavenumbers k.

	)

k( mm -1 ) λ( mm) No. of gridpoints in x/z ac/g (Theor.) ac/g (Comp.) Error(%

Acknowledgments

The authors acknowledge J. Fineberg, E. Knobloch and A. Rucklidge for insights on theoretical aspects of the Faraday phenomenon, J. Chergui, M. Firdaouss and K. Borońska for advice regarding implementation of certain numerical algorithms and, especially, A. Kityk and C. Wagner for extensive discussions and for sharing their experimental data.


L Kahouadji

N Périnet

L S Tuckerman

S Shin

J Chergui

D Juric

Numerical simulation of super-square patterns in Faraday waves

We report the first simulations of the Faraday instability using the full three-dimensional Navier-Stokes equations in domains much larger than the characteristic wavelength of the pattern. We use a massively parallel code based on a hybrid Front-Tracking/Level-set algorithm for Lagrangian tracking of arbitrarily deformable phase interfaces. Simulations performed in square and cylindrical domains yield complex patterns. In particular, a superlattice-like pattern similar to those of [Douady & Fauve, Europhys. Lett. 6, 221-226 (1988); Douady, J. Fluid Mech. 221, 383-409 (1990)] is observed. The pattern consists of the superposition of two square superlattices. We conjecture that such patterns are widespread if the square container is large compared to the critical wavelength. In the cylinder, pentagonal cells near the outer wall allow a square-wave pattern to be accommodated in the center.

Introduction

In the [START_REF] Faraday | On a peculiar class of acoustical figures; and on certain forms assumed by groups of particles upon vibrating elastic surfaces[END_REF] problem, the interface between two superposed fluid layers, subjected to periodic vertical vibration of sufficient amplitude, forms sustained standing wave patterns, called Faraday waves. Many universal dynamical-systems phenomena were first discovered through the Faraday experiment. In particular, a number of experimental investigations have focused on the variety of patterns formed in square (e.g. [START_REF] Douady | Pattern selection in Faraday instability[END_REF][START_REF] Douady | Experimental study of the Faraday instability[END_REF][START_REF] Simonelli | Surface wave mode interactions: effects of symmetry and degeneracy[END_REF] and cylindrical (e.g. [START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF][START_REF] Das | Parametrically forced gravity waves in a circular cylinder and finite-time singularity[END_REF][START_REF] Batson | The Faraday threshold in small cylinders and the sidewall non-ideality[END_REF][START_REF] Rajchenbach | Observation of star-shaped surface gravity waves[END_REF] containers, which have led to important theoretical developments (e.g. [START_REF] Meron | Parametric excitation of multimode dissipative systems[END_REF][START_REF] Silber | Parametrically excited surface waves in square geometry[END_REF][START_REF] Crawford | Period-doubling mode interactions with circular symmetry[END_REF][START_REF] Crawford | Normal forms for driven surface waves: Boundary conditions, symmetry, and genericity[END_REF][START_REF] Gomes | Steady PDEs on generalized rectangles: a change of genericity in mode interactions[END_REF] in the fields of nonlinear physics and pattern formation.

Recently, [START_REF] Périnet | Numerical simulation of Faraday waves[END_REF][START_REF] Périnet | Alternating hexagonal and striped patterns in Faraday waves[END_REF] performed the first three-dimensional direct numerical simulation of Faraday waves in horizontally periodic domains containing a few basic cells. Here, we describe numerical simulations of Faraday waves in much larger domains which have yielded more complicated patterns. For this purpose, we have used the free-surface code developed by [START_REF] Shin | A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows[END_REF] for simulations of two-phase flows on massively parallel computer architectures. In a square domain (i.e. with square horizontal cross-section), we have obtained patterns with two different length scales, very similar to those reported in experiments by [START_REF] Douady | Pattern selection in Faraday instability[END_REF] and [START_REF] Douady | Experimental study of the Faraday instability[END_REF]. Patterns of this type were later termed superlattices by [START_REF] Kudrolli | Superlattice patterns in surface waves[END_REF] and [START_REF] Arbell | Pattern formation in two-frequency forced parametric waves[END_REF], who observed a large variety of such patterns in their experiments, but by using a temporal forcing with two frequencies. The spectral analysis of our patterns, produced by single-frequency forcing, shows that they result primarily from the superposition of two superlattices constructed from very similar wavelengths. We have also simulated the same conditions in a cylindrical container. The resulting interface contains squares and pentagons.

The paper is organized as follows. First, we briefly present the problem and explain the methods involved in the code. We then present the numerical linear stability analysis that is validated by comparison with the theoretical treatment of [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]. We demonstrate the robustness of the obtained superlattice-like pattern by testing several boundary conditions. This pattern is described and analysed spectrally. Then, we show the results of simulations with the same conditions but in a cylindrical domain.

Problem formulation, governing equations and numerical scheme

The governing equations for an incompressible two-phase flow can be expressed by a single field formulation:

ρ ∂u ∂t + u • ∇u = -∇P + ρG + ∇ • µ ∇u + ∇u T + F, ∇ • u = 0 (2.1)
where u is the velocity, P is the pressure, ρ is the density, µ is the dynamic viscosity, G represents the homogeneous volume forces, and F is the local surface tension force at the interface. Here, G contains the gravitational term supplemented by a time-dependent force accounting for the vibrations of the domain:

G = -(g + γ cos(ωt)) e z (2.2)
where g is the gravitational acceleration, e z is the upward vertical unit vector, γ is the amplitude of the inertial forcing and ω = 2π/T its frequency.

The code essentially consists of two modules: one, which solves the three-dimensional incompressible Navier-Stokes equations, and the other, which treats the free surface using a parallel Lagrangian Tracking method and is only active if the flow is a two-phase flow. Material properties such as density or viscosity are defined in the entire domain:

ρ(x, t) = ρ 1 + (ρ 2 -ρ 1 ) I(x, t) µ(x, t) = µ 1 + (µ 2 -µ 1 ) I(x, t).
(2.

3)

The indicator function, I, is a numerical Heaviside function, ideally zero in one phase and one in the other phase. I is resolved with a sharp but smooth transition across 3 to 4 grid cells and is generated using a vector distance function computed directly from the tracked interface (Shin & Juric (2009a)). The fluid variables u and P are calculated by a projection method [START_REF] Chorin | Numerical simulation of the Navier-Stokes equations[END_REF]. The temporal scheme is first order, with implicit time integration used for the viscous terms. For spatial discretization we use the staggered-mesh marker-in-cell (MAC) method [START_REF] Harlow | Numerical calculation of time dependent viscous incompressible flow of fluid with free surface[END_REF]) on a uniform finite-difference grid with second-order essentially non-oscillatory (ENO) advection [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes[END_REF]. The pressure and distance function are located at cell centers while the x, y and z components of velocity are located at the faces. All spatial derivatives are approximated by standard second-order centered differences. The treatment of the free surface uses a hybrid Front-Tracking/Level-Set technique which defines the interface both by a discontinuous density field on the Eulerian grid as well as by triangles on the Lagrangian interface mesh.

The surface tension F is implemented by the hybrid formulation

F = σκ H ∇I, κ H = F L • N N • N (2.4)
where σ is the surface tension coefficient assumed to be constant and κ H is twice the mean interface curvature field calculated on the Eulerian grid, with

F L = Γ(t) κ f n f δ f (x -x f ) ds, N = Γ(t) n f δ f (x -x f ) ds (2.5)
Here, x f is a parameterization of the time-dependent interface, Γ(t), and δ f (x -x f ) is a Dirac distribution that is non-zero only where x = x f ; n f stands for the unit normal vector to the interface and ds is the length of an interface element; κ f is twice the mean interface curvature obtained on the Lagrangian interface. The geometric information, unit normal, n f , and interface element length, ds in N are computed directly from the Lagrangian interface and then distributed onto the Eulerian grid using the discrete delta function and the immersed boundary method of [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF]. A detailed description of the procedure for calculating F, N and I can be found in [START_REF] Shin | Computation of the curvature field in numerical simulation of multiphase flow[END_REF]. The Lagrangian interface is advected by integrating dx f /dt = V with a second-order Runge-Kutta method where the interface velocity, V, is interpolated from the Eulerian velocity.

The time step ∆t is chosen at each iteration in order to satisfy a criterion based on

{∆t CFL , ∆t int , ∆t vis , ∆t cap } (2.6)
which ensures stability of the calculations. These bounds are defined by:

∆t CFL ≡ min j min domain ∆x j u j ∆t int ≡ min j min Γ(t) ∆x j V ∆t vis ≡ min ρ 2 µ 2 , ρ 1 µ 1 ∆x 2 min 6 ∆t cap ≡ 1 2 (ρ 1 + ρ 2 )∆x 3 min πσ 1/2 (2.7)
where ∆x min = min j (∆x j ). With the periodic volume force term G of (2.2), a supplementary upper bound ∆t ω = 2π/(50ω) is required. Fortran 2003 allows the definition of a set of dynamically allocated derived types and generic procedures associated with the grids, scalar and vector fields, operators as well as the various solvers used in the Navier-Stokes and Lagrangian tracking modules. The parallelization of the code is based on algebraic domain decomposition, where the velocity field is solved by a parallel generalized minimum residual (GMRES) method for the implicit viscous terms and the pressure by a parallel multigrid method motivated by the algorithm of [START_REF] Kwak | Multigrid algorithm for the cell-centered finitedifference method II: Discontinuous coefficient case[END_REF]. Communication across process threads is handled by message passing interface (MPI) procedures.

Finally, the code also contains a module for the definition of immersed solid objects and their interaction with the flow, which we have used to simulate Faraday waves in a cylindrical container. A Navier-slip dynamic contact line model is implemented where hysteresis is accounted for by fixing the contact angle limits to prescribed advancing or receding angles as in Shin & Juric (2009b). Further details are available in [START_REF] Shin | A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows[END_REF].

Parameters and linear stability

We consider a layer of silicone oil of thickness h = 14.5 mm with density ρ 1 = 965 kg m -3 and dynamic viscosity µ 1 = 0.02 kg m -1 s -1 . This layer is overlaid with 29.5 mm of air of density ρ 2 = 1.205 kg m -3 and dynamic viscosity µ 2 = 1.82 × 10 -5 kg m -1 s -1 . The surface tension at the interface between the oil and the air is σ = 0.02 kg s -2 . This choice of parameters was originally motivated by experiments on Faraday waves in a rotating cylindrical container (Clément, Pucci & Couder, unpublished), which will be the subject of a future investigation.

Before simulating the patterns, we computed the critical acceleration γ c from direct numerical simulations (DNS) and compared it with that of the Floquet linear stability analysis of [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF]. These simulations are conducted in a small periodic domain of length λ c and width λ c /2, where the critical wavelength is λ c = 2π/k c = 11.3 mm with k c the critical wavenumber. We compute the growth rate of the total kinetic energy of the system (figure 1(a)) after the growth becomes exponential. Near the threshold γ c , the growth rate varies linearly with γ, so that γ c can be deduced by linear interpolation. Table 1 compares the Floquet predictions of γ c to the results from DNS for three vibration frequencies. The discrepancy is under 3% in all three cases. The ratio of gravitational to capillary forces is measured by the Bond number, whose value here is g∆ρ/(σk 2 c ) = 1.5, indicating that both forces are at work.

Results and discussion

Visualisations

The top and bottom of the domain are taken to be rigid; we impose no-slip boundary conditions on the velocity field. Three types of boundary conditions have been tested on the lateral walls, (a) no-slip, (b) free-slip, and (c) periodic. Although our parameters were originally chosen for a different purpose, we find that the patterns that we observe are strikingly similar to those observed in earlier experiments by [START_REF] Douady | Pattern selection in Faraday instability[END_REF], figure 1(b) and [START_REF] Douady | Experimental study of the Faraday instability[END_REF], figures 10(b), 11(b) and 15(b), with different parameters: water, with ρ 1 = 1000 kg m -3 , µ 1 = 9 × 10 -4 kg m -1 s -1 , and σ = 0.03 kg s -2 , in a container of dimensions 80 mm × 80 mm × 5 mm oscillated with f = 37 Hz or f = 70 Hz. For all three boundary conditions, we observe superlattice-like patterns shown in figure 3. The interface in the x-y plane consists of four large square sub-blocks, each composed of smaller squares whose sides have approximate length λ c . In each case, the overall pattern has symmetry D 2 , i.e., it is invariant under reflections through the two diagonals. Each of the four sub-blocks is in phase opposition with its two adjacent neighbors. The patterns in (a) and (c) have periodicity length L while (b) has a period of 2L. Although one would expect the periodic case (c) to be the simplest to analyze, we observe that the blocks in the free-slip case (b) are more homogeneous than those in cases (a) and (c); we will therefore focus on this case.

Figure 4 displays the temporal evolution of the interface at saturation. (A movie of this evolution is available as supplementary material to this article.) After one period of forcing oscillation T , troughs are replaced by crests (see figure 4(a) and (c)) as a consequence of the subharmonic behavior of the interface. This subharmonic behavior is also displayed by the large sub-blocks, since the appearance of each block is transformed into that of its two adjacent neighbors after time T and then returns to its initial shape after time 2T has passed. The pattern is thus invariant under the combined operations of rotation by π/2 and time-translation by T , a spatio-temporal symmetry.

Figure 5 shows the interface and the velocity on a vertical slice at the same three instants. The contrasting behavior of the velocity in the left and right halves of the slice is displayed. The interface of figure 5(a) has a maximum (minimum) at the left (right) boundary and vice versa for (c), illustrating again that the field is not periodic in a domain of length L.

Fourier spectra

For free-slip boundary conditions, the pattern is periodic when it is doubled horizontally by reflection at the boundaries. The doubled pattern is invariant under the symmetry group D 4 of the square, consisting of reflections through the x and y coordinate axes and rotations of multiples of π/2. We decompose the doubled interface height profile ζ(x, y, t) into spatial Fourier modes, and then into spatio-temporal modes as follows:

ζ(x, y, t) = ∞ l,m=-∞ ζl,m (t)e ik l,m •x = ∞ l,m=-∞ ∞ n=-∞ ζl,m,n e inΩt/2 e ik l,m •x (4.1)
where k l,m ≡ πl L , πm L , x = (x, y) and k l,m = |k l,m |. The square symmetry of the doubled interface implies that the eight modes in an octet ζ±l,±m and ζ±m,±l all have the same amplitude. Figure 6 shows two such octets, with representative elements ζ23,1 and ζ23,3 , located on circles which we will call C 1 and C 3 , respectively, i.e. with l = 23 and m = 1 or m = 3.

The circle C 1 has radius k c , the critical wavenumber. The dominant wavevector octet of our pattern, represented by k 23,1 , is in accordance with the experimental observations of [START_REF] Douady | Pattern selection in Faraday instability[END_REF], who scanned over the forcing frequency in order to vary the selected wavenumbers and patterns. They noted that the range of existence was greatest for pairs (l, m) with a small ratio m/l and was maximum for m = 1, as in our case, i.e. for wavevectors which are almost parallel to the walls.

The modes on C 3 can be seen as the result of cubic interactions between modes of C 1 :

k 23,3 = k 23,1 + k 23,1 + k -23,1 (4.2)
Modes on C 3 are fed continually and damped slowly, if at all, since |k 23,3 -k c |/k c = 0.75%. Generalizing (4.2) to any m l shows that such cubic interactions lead to k l,3m , with |k l,3m -k c |/k c ∼ 4(m/l) 2 . Due to its proximity to the critical circle, this mode is weakly damped and should often be present when Faraday waves are excited in a square domain which is much larger than the critical wavelength. A second harmonic of C 1 , arising from the quadratic interaction

k 46,0 = k 23,1 + k 23,-1 (4.3)
is seen in figure 6(b) along with its image under π/2 rotation. These modes occur in a quartet rather than an octet because of their reflection symmetry in the coordinate axes and are of smaller amplitude than the modes on C 1 or C 3 . The experimental techniques of [START_REF] Douady | Pattern selection in Faraday instability[END_REF] did not give them access to amplitudes of Fourier modes other than the dominant one.

According to the terminology of [START_REF] Simonelli | Surface wave mode interactions: effects of symmetry and degeneracy[END_REF], the four components ζ±l,±m comprise one pure mode, while ζ±m,±l comprise another pure mode. [START_REF] Douady | Pattern selection in Faraday instability[END_REF] observed mixed modes (which they called symmetric) for small m/l, and pure modes (which they called dissymmetric) for larger m/l. Pure modes may be combined to form two types of mixed modes, the in-phase mode in which all of the components have the same sign (as in our figure 6), and the out-of-phase mode in which the components of the two pure modes are of opposite sign. [START_REF] Silber | Parametrically excited surface waves in square geometry[END_REF] showed that these two types of mixed modes are equivalent if l and m have opposite parities. [START_REF] Crawford | Normal forms for driven surface waves: Boundary conditions, symmetry, and genericity[END_REF] showed that this was also true for l and m with the same parity (as in our case) by invoking the concept of hidden symmetries. In a square container with Neumann boundary conditions, the doubled pattern formed by reflecting at each boundary, as we have done above, satisfies periodic boundary conditions, whose inherited (hidden) translation invariances are incorporated into the normal form.

For each of C 1 or C 3 , the resulting pattern has a superlattice structure, in which squares of different sizes coexist, as shown in figure 7 The temporal evolution of ζ23,1 , ζ23,3 and ζ46,0 is displayed in figure 8(a) and the spatio-temporal spectrum is shown in figure 8(b). The mode amplitudes ζ23,1 and ζ23,3 are related by a multiplicative constant and contain only subharmonic temporal components (i.e. n odd in (4.1)). The temporal component n = 1 corresponding to Ω = ω/2 dominates the temporal evolution. Both ζ23,1 and ζ23,3 evolve in phase opposition, crossing zero at the same time and yielding patterns such as that in figure 4(b). In contrast, ζ46,0 is harmonic and dominated by two temporal components n = 0 and n = 2, as ex- pected from the quadratic resonance of ζ23,1 with itself or with ζ23,3 . Each spatial mode has a finite n = 0 mean component; the offset of ζ46,0 is especially noticeable.

Simulation in a cylindrical container

We performed simulations under the same conditions, but in a cylindrical container of diameter D = L and height h = L/3 (see figure 2(b)). The lateral boundary conditions are free-slip and the advancing and receding contact angles are fixed at 100 • and 60 • respectively. Figure 9 shows snapshots of the interface every T /2. (A movie of this evolution is available as supplementary material to this article.) The pattern has the spatio-temporal symmetry of invariance under combined time-evolution by T /2 and rotation by π/2, while each instantaneous pattern is reflection-symmetric about both diagonal axes. At the center, we observe squares whose characteristic wavelength remains close to λ c of Table 1 for f = 30 Hz. Closer to the boundary, five-sided cells can be seen as the pattern organizes itself to fit in the circular domain. Aside from this, the pattern is relatively insensitive to the circular shape of the domain because L/λ c = 11.7 1, in contrast with small-radius experiments and theory. In most of these (e.g. [START_REF] Ciliberto | Chaotic mode competition in parametrically forced surface waves[END_REF][START_REF] Crawford | Period-doubling mode interactions with circular symmetry[END_REF][START_REF] Das | Parametrically forced gravity waves in a circular cylinder and finite-time singularity[END_REF][START_REF] Batson | The Faraday threshold in small cylinders and the sidewall non-ideality[END_REF], the interface takes the form of Bessel functions, including axisymmetric modes, while [START_REF] Rajchenbach | Observation of star-shaped surface gravity waves[END_REF] report intriguing patterns of the form of stars and pentagons.

Conclusion

We have presented the first three-dimensional numerical results of the Faraday instability in domains much larger than the critical wavelength in both square and cylindrical containers. In the square domain, the interface displays patterns that are very similar to those found in the work of [START_REF] Douady | Pattern selection in Faraday instability[END_REF] and [START_REF] Douady | Experimental study of the Faraday instability[END_REF]. By means of a spatial spectral decomposition, we have described these patterns quantitatively and have found that they have a complex superlattice-like structure due to resonances between modes of very similar wavelengths. We conjecture that this pattern arises from two effects: (i) In a square container which is large compared to the critical wavelength, Faraday wave patterns tend to be mixed modes whose wavevectors are almost parallel to the boundaries of the box, as stated by [START_REF] Douady | Pattern selection in Faraday instability[END_REF]. (ii) In this case, cubic nonlinearities generate another wavelength which is very close to the dominant one. The combination of the two effects leads to the superlattice pattern that we observe in our simulations and we infer that such patterns should be observed in a wide variety of systems. The present work demonstrates that it is now possible to numerically compute all the fields of Faraday waves in domains sufficiently large to produce complex patterns such as superlattices, quasicrystalline patterns or oscillons. The hydrodynamic mechanisms responsible for the formation of these patterns can then be explored numerically.

the Agence Nationale de la Recherche (ANR) for the TRANSFLOW project. This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and future planning (NRF-2014R1A2A1A11051346).

Figure 1 .

 1 Figure 1. a) Evolution of the total kinetic energy for ω/2π = 30 Hz and various accelerations in a small domain. The fact that the peaks for a given acceleration lie along a line shows that the energy growth is very close to exponential. b) Evolution of the various bounds on the time step during the course of the simulation. The smallest, ∆tcap, is used in the simulations.

Figure 2 .

 2 Figure 2. Calculation domain of size L × L × L/3 (L = 132 mm), divided into 8 × 8 × 4 subdomains. a) Interface profile. b) an immersed solid cylinder of diameter D = L and height h = L/3 surrounding the flow. The flow outside of the cylinder is not computed. For both (a) and (b) the Cartesian resolution per subdomain is 64 3 , which gives a global resolution of 512 × 512 × 256.

Figure 3 .

 3 Figure 3. Snapshots of the top view of the interface colored by the vertical velocity for (a) no-slip, (b) free-slip, and (c) periodic boundary conditions.

Figure 4 .

 4 Figure4. Interface profiles at times separated by intervals of T /2, where T is the forcing period and 2T is the subharmonic response period. In (a) and (c), the small squares consist of wells surrounded by ridges with peaks at each corner, while in (b), low central peaks are surrounded by four higher peaks. In (a), a weak diagonal "bridge" connects the large squares on the bottom left and top right; in (c), the bridge links the top left and bottom right.

Figure 5 .

 5 Figure 5. Velocity profiles on the plane y = 3L/4, corresponding to each of the three instants shown in figure 4. The vectors are colored according to their vertical component.

Figure 6 .

 6 Figure 6. Spatial spectrum of a snapshot of the interface in the doubled domain. Circles C1 (dark blue) and C3 (dark brown) have slightly different radii. C1 contains the fundamental mode ζ23,1 and its images under reflection and rotation (dark blue dots surrounded by circles). C3 contains the highest third-harmonic mode ζ23,3 and its images (dark brown dots surrounded by squares). Only these modes and, to a lesser extent, the weaker second harmonic ζ46,0 (orange) have significant amplitude. a) Restriction to [-25, 25] 2 highlights the square symmetry of the pattern. b) Restriction to the first quadrant [0, 50] 2 includes ζ46,0.

  (a) and (b). Although shown in domains of length L, the periodicity length is 2L. The smallest squares have length scale 2π/k 23,m for m = 1, 3. The patterns formed by modes on C 3 have an additional intermediate length scale 2L/3. Close examination of figure 7(b) shows that adjacent medium squares are almost, but not exactly identical. Figure 7(c) shows a superposition of the two patterns of figure 7(a) and (b), weighted by the amplitudes of their respective Fourier components. During most of the time, the interface resembles this superposition; compare with figure 4. When the amplitudes of modes on C 1 and C 3 approach zero, which occurs twice during one subharmonic period (see figure 8(a) and figure 4(b)), it is necessary to take modes like ζ46,0 into account.

Figure 7 .

 7 Figure 7. Patterns constructed from Fourier modes (a) in C1, (b) in C3, and (c) a superposition of the two. Lighter (darker) zones correspond to the peaks (troughs) of the interface. Superlattices, consisting of smaller and larger squares, are observed in all of the plots, Smallest boxes are approximately of size λc. b) Medium squares of size 2L/3 are present (long-dashed box). c) Reconstructed pattern, made by combining patterns in (a) and (b) weighted by the corresponding Fourier coefficients ζ23,1 and ζ23,3. Compare with figure 4(c).

Figure 8 .

 8 Figure 8. a) Temporal evolution of the three main spatial modes. Modes ζ23,1 (blue curve and circles) and ζ23,3 (red curve and + symbols) are subharmonic: their period of oscillation is twice that of the forcing. Mode ζ46,0 (black curve and × symbols) is harmonic. b) Temporal spectrum of the main modes. ζ23,1, ζ23,3 contain only odd temporal harmonics while ζ46,0 has mainly even harmonics. Each spatial mode has a finite n = 0 component and hence a non-zero temporal mean.

Figure 9 .

 9 Figure 9. Interface profile in the cylinder at time intervals of T /2. The central part of the domain is still occupied by squares, but five-sided cells are present near the boundaries, allowing the pattern to fit inside a circle.

Table 1 .

 1 Comparison of Floquet and DNS instability thresholds for various frequencies.

	30 11.3 0.733 0.713	2.73 %
	60	5.9	2.65	2.61	1.50 %
	90	4.3	5.32	5.3	0.38 %

ω/2π λc Floquet DNS relative error Hz mm γc/g γc/g Floquet/DNS

Acknowledgements

We thank S. Douady, S. Fauve, and G. Pucci for helpful discussions. This work was performed using high performance computing resources provided by the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Équipement National de Calcul Intensif). N. Périnet was supported by a grant from CONICYT -FONDECYT/postdoctoral research project 3140522. L.S.T. acknowledges support from