Aditi Chakrabarti

Serge Mora

Franck Richard

Ty Phou

Jean-Marc Fromental

Yves Pomeau

Basile Audoly

Selection of hexagonal buckling patterns by the elastic Rayleigh-Taylor instability

Keywords: A. Buckling, B. Elastic material, B. Finite strain, B. Plates, C. Stability and bifurcation

Selection of hexagonal buckling patterns by the elastic Rayleigh-Taylor instability 1. Introduction Elastic buckling phenomena have been classically investigated in thin or slender structures such as shells, plates and rods, that are effectively compliant at large scale owing to their thin or slender geometry. More recently, the attention has been extended to buckling instabilities arising in elastic solids made of soft materials, such as gels or soft polymers. Soft elastic solids are attractive for applications as they can be actuated by forces that are otherwise too weak to induce significant elastic strains such as their own weight [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF], electric forces [START_REF] Arun | Electric-field-controlled surface instabilities in soft elastic films[END_REF][START_REF] Wang | Creasing to cratering instability in polymers under ultrahigh electric fields[END_REF][START_REF] Huang | Giant, voltage-actuated deformation of a dielectric elastomer under dead load[END_REF][START_REF] Bense | Buckling of elastomer sheets under non-uniform electro-actuation[END_REF], magnetic forces [START_REF] Danas | Instability of a magnetoelastic layer resting on a non-magnetic substrate[END_REF], adhesive forces [START_REF] Ghatak | Meniscus instability in a thin elastic film[END_REF][START_REF] Mönch | Elastic instability of rubber films between solid bodies[END_REF], or even the capillary forces present at a curved solid-fluid interface [START_REF] Mora | Capillarity driven instability of a soft solid[END_REF][START_REF] Mora | Solid drops: Large capillary deformations of immersed elastic rods[END_REF]. As they undergo large strains, soft elastic solids display a non-linear response and are prone to a variety of buckling instabilities [START_REF] Biot | Surface instability of rubber in compression[END_REF][START_REF] Tanaka | Mechanical instability of gels at the phase transition[END_REF][START_REF] Mora | Surface instability of soft solids under strain[END_REF][START_REF] Ciarletta | Shear instability in skin tissue[END_REF][START_REF] Lagrange | From wrinkling to global buckling of a ring on a curved substrate[END_REF]; some of these instabilities are discontinuous and are therefore difficult to approach analytically: this is the case of the creasing instability for example [START_REF] Hong | Formation of creases on the surfaces of elastomers and gels[END_REF][START_REF] Cao | From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling[END_REF][START_REF] Hohlfeld | Scale and nature of sulcification patterns[END_REF][START_REF] Ciarletta | A semi-analytical approach to Biot instability in a growing layer: Strain gradient correction, weakly non-linear analysis and imperfection sensitivity[END_REF]. This is also the case of the elastic Rayleigh-Taylor instability, which we investigate in this paper. This instability is obtained when a thick slab of elastic material is hung below a rigid plane, see figure 1. A competition takes place between the elasticity of the slab (which tends to keep the slab undeformed) and its weight (which acts as a destabilizing force). When elasticity wins, the slab remains undeformed and its lower free interface appears as planar when observed from below. When gravity overcomes elasticity, however, the free surface of the slab buckles into a hexagonal pattern [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF], as shown in figure 1c. The main goal of the present paper is to explain the selection of the hexagonal buckling pattern and the discontinuous character of the transition.

The elastic Rayleigh-Taylor instability is the counterpart of the classical Rayleigh-Taylor instability in fluids, observed when a layer of a denser viscous fluid is placed above a layer of lighter viscous fluid and tries to 'push into it'. In typical experiments, the difference in density results from the presence of a vertical temperature gradient and from thermal expansion. When the density mismatch is large enough, fluctuations at the interface between the two layers grow in time, and the heavier fluid eventually forms droplets that detach as they fall into the lighter fluid [START_REF] Rayleigh | Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density[END_REF][START_REF] Chandrasekhar | The character of the equilibrium of an incompressible heavy viscous fluid of variable density[END_REF][START_REF] Fermigier | Two-dimensional patterns in Rayleigh-Taylor instability of a thin layer[END_REF]. By contrast the elastic Rayleigh-Taylor instability is reversible, and the formation of dropletlike patterns halts before detachment occurs [START_REF] Riccobelli | Rayleigh-Taylor instability in soft elastic layers[END_REF]. Like its elastic counterpart, the fluid Rayleigh-Taylor instability can produce a hexagonal networks of drops: this happens for instance when a rectangular plate is coated with a viscous fluid from below [START_REF] Fermigier | Two-dimensional patterns in Rayleigh-Taylor instability of a thin layer[END_REF]. The underlying selection mechanism is well known, and the Rayleigh-Taylor instability is now a classical example of pattern formation in fluid mechanics, together with the Rayleigh-Bénard instability which can produce hexagons by a similar mechanism.

By contrast, the selection of the patterns produced by the elastic Rayleigh-Taylor instability has not been documented. The present work is an attempt to bridge this gap, building up on ideas originally developed in the context of pattern formation in fluids that have later been combined with the classical methods of elastic bifurcation and stability.

In previous work, some of us have carried out a linear bifurcation analysis of the elastic Rayleigh-Taylor instability [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF]. This linear analysis selects critical values of the load and of the wavenumber that have been verified to approximately agree with the experiments. Due to the symmetries of the system, however, a continuum of critical wavevectors appears at the bifurcation: they have all the same norm and point in all possible directions in a horizontal plane. The corresponding modes can be combined arbitrarily, and each combination yields a different pattern: the linear bifurcation analysis is not sufficient to account for the pattern selection. A non-linear analysis is required to address the cooperation between the linear modes. This is the goal of the present work. Specifically, we derive equilibrium solutions for the different buckling patterns, and characterize their stability close to the bifurcation threshold by means of a post-bifurcation expansion [START_REF] Lindstedt | Beitrag zur integration der differentialgleichungen der strungtheorie[END_REF][START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste: Méthodes de MM[END_REF].

The unbuckled configuration of the slab is transversely isotropic, i.e. it is invariant by any translation in the horizontal direction and with respect to any rotation about the vertical axis. This symmetry plays a key role in the analysis of the elastic Rayleigh-Taylor instability. In the linear bifurcation analysis, it makes a continuum of critical wavevectors available at bifurcation, as noted already. In the non-linear bifurcation analysis, the transverse symmetry allows for resonances involving three linear modes. These resonances produce the hexagonal buckling pattern. The selection mechanism is geometric in essence, and our analysis of the Rayleigh-Taylor instability uses ideas that have been developed for the analysis of bifurcations in the presence of symmetries [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF][START_REF] Golubitsky | Singularities and groups in bifurcation theory[END_REF], which can similarly apply to other transversely isotropic structures. This includes swelling gel layers [START_REF] Tanaka | Mechanical instability of gels at the phase transition[END_REF], thin films under equi-biaxial compression on an elastic foundation [START_REF] Breid | Surface wrinkling behavior of finite circular plates[END_REF][START_REF] Cai | Periodic patterns and energy states of buckled films on compliant substrates[END_REF], thin films deformed by high electric fields [START_REF] Wang | Creasing to cratering instability in polymers under ultrahigh electric fields[END_REF] and spherical shells under external pressure [START_REF] Carlson | Experimental studies of the buckling of complete spherical shells[END_REF]-at least on length-scales smaller than the sphere radius where defects in the hexagonal pattern can be ignored. All these systems have been reported to display a preference towards hexagonal patterns close to buckling threshold; this is at odds with the conclusions of most of the buckling analyses from prior work, which suggested that other patterns such as one-dimensional stripes or networks of square should also be observed. Here, we show that a generic selection mechanism applicable to transversely isotropic structures is applicable and favors hexagons over other buckling patterns; this restores consistency with the experimental observations.

The main ideas explaining the formation of hexagonal patterns have emerged in the context of Rayleigh-Bénard convection, and we start by reviewing them briefly. A century ago, [START_REF] Bénard | Les tourbillons cellulaires dans une nappe liquide propageant de la chaleur par convection: en régime permanent[END_REF] reported the occurrence of hexagonal convection cells in a horizontal layer of fluid heated from below. When the top interface of the layer is free, gradients of surface tension induced by temperature are involved in the formation of the patterns; this is known as the Bénard-Marangoni instability [START_REF] Scriven | On cellular convection driven by surface-tension gradients: effects of mean surface tension and surface viscosity[END_REF]. When the upper interface is in contact with a cool plate instead, the flow patterns are due to gradients in buoyancy, again effected by temperature; this is known as the Rayleigh-Bénard instability [START_REF] Rayleigh | On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side[END_REF]. [START_REF] Palm | On the tendency towards hexagonal cells in steady convection[END_REF] was the first to explain the formation of hexagons in Rayleigh-Bénard convection by a qualitative argument based on the resonance of three linear modes, whose wavevectors form the sides of an equilateral triangle. In the fluid experiments as well, the hexagons appear by a discontinuous transition. This feature has been explained by [START_REF] Segel | On the question of the preferred mode in cellular thermal convection[END_REF], [START_REF] Busse | Non-linear properties of thermal convection[END_REF] and [START_REF] Joseph | Bifurcating time periodic solutions and their stability[END_REF] who showed that the hexagons result from a trans-critical bifurcation. For a more complete survey on Rayleigh-Bénard instabilities, the reader is referred to the works of [START_REF] Palm | Nonlinear thermal convection[END_REF]; [START_REF] Normand | Convective instability: a physicist's approach[END_REF]; [START_REF] Koschmieder | Bénard cells and Taylor vortices[END_REF]; Manneville (2006); [START_REF] Pomeau | Front motion, metastability and subcritical bifurcations in hydrodynamics[END_REF].

In the context of elastic buckling, the role of three-modes resonances in the formation of hexagonal patterns has been appreciated by some authors. [START_REF] Tanaka | Mechanical instability of gels at the phase transition[END_REF] invokes these resonances in connection with his observations of hexagons at the surface of swollen gels. [START_REF] Hutchinson | Imperfection sensitivity of externally pressurized spherical shells[END_REF] carries a post-bifurcation expansion of spherical shells under external pressure based on a combination of three linear modes producing a hexagonal pattern, which indeed captures the resonance phenomenon. Jia and Ben Amar (2013) identify a branch of hexagonal solutions that appear by a transcritical bifurcation in the post-bifurcation analysis of a swollen elastic layer resting on a rigid support. Here, we build upon these previous works and propose a systematic analysis of the patterns, addressing not only hexagons but also other patterns such as stripes and squares. As outlined throughout this paper, our findings are fully consistent the predictions of the general theory of bifurcations in the presence of symmetries due to [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF]; [START_REF] Buzano | Bifurcation on the hexagonal lattice and the planar bénard problem[END_REF]; [START_REF] Golubitsky | Singularities and groups in bifurcation theory[END_REF]; our presentation is self-consistent and no prior knowledge of this theory is assumed.

This paper is organized as follows. In Section 2 we report our experimental observations on the patterns produced by the elastic Rayleigh-Taylor instability. In Section 3 we provide a complete list of equations governing the equilibrium of an elastic slab undergoing finite strain. In Section 4 we recall the linear stability analysis from our previous work, and provide explicit expressions for the linear modes. Next, we analyze the small-amplitude patterns, starting with the non-symmetric hexagons ( §5), and proceeding next to the symmetric patterns such as stripes and squares ( §6). This yields a complete bifurcation diagram in the neighborhood of the bifurcation point, as summarized in Section 7. In Section 8, we characterize the finite-amplitude buckling patterns that appear further away from the bifurcation point, by running of finite-element simulations and by analyzing a simplified energy functional.

Experiments

A parallelepipedic container with dimensions × × h = (40 cm) 2 × 2.5cm is filled with the reagents generating a cross-linked polyacrylamide gel. The cured gel behaves as an incompressible, isotropic and purely elastic solid for strains up to several hundreds percent. Its shear modulus µ can be finely tuned by varying the concentrations in monomers and cross-linkers or the temperature, and can be as low as a few tens of Pascal [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF]; its density ρ is very close to that of water. After the gel has cured, the container is turned upside down, so that the free surface of gel is facing downwards, see figure 1a. The gel is then observed from below and the deformation of its lower interface is monitored as a mean to track buckling instabilities, see figure 1b-c.

Several experiments have been performed with different concentrations and temperatures, corresponding to different values of the shear modulus µ. The shear modulus µ is measured by indentation of a nonadhesive rigid sphere (Teflon) at the surface of a control sample prepared simultaneously and with the same concentrations: the control sample is fully covered with pure water to effectively remove the capillary forces at the fluid-air interface, and the indentation force is measured along with the indentation depth, leading to the determination of the shear modulus [START_REF] Czerner | Determination of elastic modulus of gelatin gels by indentation experiments[END_REF][START_REF] Tong | Comparison of spherical and flat tips for indentation of hydrogels[END_REF]. Returning to the experiments in figure 1, the downward-facing free surface is found to remain perfectly flat in all the experiments where the gel's shear modulus is larger than µ = 40.5 Pa, as in figure 1b. By contrast, a pattern of dimples appear for all the tested gels having a shear modulus lower than µ = 40.5 Pa, as shown in figure 1c; once formed, the pattern remains in place permanently. These dimples are separated by vertical The vertical line is the buckling threshold αc predicted by the linear bifurcation analysis in §4 and [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF]. Hexagonal patterns are observed well below the critical load αc predicted by the linear bifurcation analysis (dashed vertical line). This discrepancy will be explained by our non-linear analysis.

walls of highly deformed gel, and are organized in close-to-hexagonal patterns. The hexagonal pattern is perturbed by the container's walls and is more regular toward the center of the container.

To quantify the buckling amplitude, we have measured the difference w b of the maximum thickness h + w b after buckling, minus the initial thickness h, see figure 1a. The amplitude w b is also the displacement of the points that end up on the bottom of the walls separating the hexagonal cells, measured positively in the downward direction. To avoid any bias caused by edge effects, we have ignored any wall adjacent to the sides of the container in the determination of w b . The buckling amplitude w b has been measured as a function of the shear modulus of the specimen. In figure 2, w b is rescaled by the initial thickness h, and plotted as a function of the dimensionless loading parameter α = (ρgh)/µ, where g is the acceleration of gravity. At a critical value of the control parameter, α = (ρgh)/µ = 5.67 ± 0.1, the observed pattern goes from flat (unbuckled) to hexagonal. This transition is clearly discontinuous, as shown by the jump in amplitude in the figure. This discontinuity can be observed in containers whose size is significantly larger than the gel thickness, such as that used in our experiments, for which /h = 16 1. The discontinuity was not present, however, in our previous set of experiments [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF] where we used a smaller aspect-ratio /h = (18 cm) / (2.75 cm) = 6.5.

Model

With the aim to explain the selection of the hexagonal buckling patterns and the discontinuous nature of the bifurcation, we analyze the buckling of an elastic slab hanging from below a rigid horizontal plane. We start in this section by deriving the equations governing the non-linear equilibrium of the slab for finite strain, which will serve as the foundation for our buckling analysis later on.

Geometry and kinematics

We choose the reference configuration of the gel to be its natural configuration when the container has been flipped upside down, i.e. with its unconstrained face facing downwards. We introduce Cartesian coordinates (X, Y, Z), with the third axis Z oriented in the upward vertical direction, i.e. from the free interface to the clamped interface, with the force due to gravity acting in the -Z direction. The coordinates X = (X, Y, Z) in reference configuration are used as Lagrangian coordinates that label material points in the gel. The free interface at the bottom, and the upper interface bound to the rigid support correspond to Z = 0 and Z = h, respectively.

In the deformed configuration, the position of a generic point in the gel is denoted as x(X), where x = (x, y, z) are the Cartesian coordinates. Let F = ∇ x(X) denote the deformation gradient, C = F T • F Green's deformation tensor and I 1 and I 2 its invariants,

I 1 = tr C -3 I 2 = tr 2 C-C : C 2 -3,
where the double contraction has been denoted as

C : C = tr(C • C).
The attachment to horizontal support is written as

x(X, Y, h) = (X, Y, h).

(1)

Generic constitutive model for the gel

The gel is modeled as an incompressible and isotropic hyper-elastic material. The incompressibility is expressed by the constraint det F = 1.

(2)

We consider a generic constitutive law for an isotropic material: the dependence of the strain energy density W (I 1 , I 2 ) on the deformation invariants occurs through a function W that will remain unspecified for the moment. This allows us to carry the buckling analysis in full generality: as we shall see, the only input to the bifurcation analysis will be the derivatives of W evaluated in the natural configuration (I 1 , I 2 ) = (0, 0). The particular cases of a neo-Hookean material (NH) or of Gent material, for instance, correspond to the following strain energy densities,

W NH (I 1 , I 2 ) = µ 2 I 1 W Gent (I 1 , I 2 ) = -µJm 2 ln (1 -I 1 /J m ) , (3)
where µ is the shear modulus and J m is an additional material parameter specific to the Gent model.

Returning to a general constitutive model W (I 1 , I 2 ), we obtain the (initial) shear modulus as

µ = 2 ∂W ∂I 1 (0, 0) + ∂W ∂I 2 (0, 0) , (4)
where the gradients are evaluated in the natural configuration (I 1 , I 2 ) = (0, 0). Of particular importance is the dependence of the constitutive law on the second invariant I 2 . It is measured by a material parameter τ defined as

τ = 1 µ ∂W ∂I 2 (0, 0). ( 5
)
This parameter cancels for any constitutive law such that the strain energy is independent of the second invariant I 2 : this is the case for both the neo-Hookean and Gent models, see (3),

τ NH = 0 τ Gent = 0. ( 6
)
reference configuration actual configuration

Figure 3: Reference and actual configurations of the elastic slab. We use dimensionless units, with a scaled thickness h = 1, a scaled shear modulus µ = 1, and a scaled weight α defined in (8). These drawings are shown in 2 dimensions for the sake of legibility; the analysis, however, is done in 3 dimensions: the direction Y transverse to the place of the figure is accounted for.

Energy

The total potential energy of the system is the sum of the elastic strain energy, and of the potential energy due to gravity. Including an additional constraint term to account for the incompressibility, we obtain the augmented energy as

E = h 0 dZ [0, ] 2 dXdY W (I 1 , I 2 ) + ρgz(X) -p(X) (det F -1) , (7)
where p(X, Y, Z) is a Lagrange multiplier proportional to the hydrostatic pressure in the gel. Non-linear equilibria of the system can be obtained from the condition that E is stationary with respect to virtual motions that are compatible with the boundary condition (1).

Rescaled energy

The dimensionless parameters governing the equilibrium of the gel are the aspect-ratio h/ and the dimensionless gravity number α,

α = ρgh µ . (8)
In the forthcoming analysis we ignore any effect related to the finite size of the container, i.e. consider the limit h/ 1. The only dimensionless parameter left in the problem is then α. To discuss the selection of patterns, it will be convenient to work in units such that both the initial shear modulus and the initial height are equal to one,

µ = 1, h = 1. (9)
Choosing µ = 1 > 0 is sufficient to warrant material stability in the neighborhood of the undeformed configuration. The reference and actual geometries are sketched in figure 3 in terms of the rescaled quantities. By a small abuse of notation, we continue to employ the same symbols when working in rescaled variables: for instance, the rescaled augmented energy (7) reads

E = 1 0 dZ dXdY W (I 1 , I 2 ) + αz(X) -p(X) (det F -1) , (10)
where the double integral now involves the entire plane (X, Y ) ∈ R 2 . The boundary condition on the upper interface reads x(X, Y, 1) = (X, Y, 1). ( 11)

In rescaled variables, (4) and (5) yield the gradient of the strain energy function in terms of the material parameter τ as

∂W ∂I 1 , ∂W ∂I 2 (I1,I2)=(0,0) = 1 2 -τ, τ (12)

Non-linear equilibrium

All unknowns are collected in a state vector t = (x, y, z, p): a configuration of the gel is then parameterized as t(X). For a given value of the load parameter α, a configuration t(X) is an equilibrium if it satisfies the incompressibility condition (2), the boundary condition (11) and makes the augmented energy E stationary,

(∀ t(X) / x(X, Y, 1) = 0) DE(α, t) • [ t] = 0. (13a)
Here, t(X) denotes a virtual motion, t(X) = (x(X), ŷ(X), ẑ(X), p(X)), i.e. an infinitesimal increment of all the unknowns, and DE(α, t) • [ t] stands for the first variation of the energy functional E in the direction of the virtual motion t(X) evaluated in the configuration t(X), namely

DE(α, t) • [ t] = lim u→0 E(α, t + u t) -E(α, t) u .
Equation ( 13a) is the non-linear principle of virtual work for the incompressible hyper-elastic slab. This variational equality is restricted to virtual motions t(X) that satisfy the incremental kinematic boundary condition x = (x, ŷ, ẑ) = 0 at the upper interface Z = 1. By contrast, the stress-free boundary condition on the lower interface will be a consequence of the equilibrium (13a); there is no need to impose this condition on the virtual motions. There is no need to enforce a condition of volume preservation on the virtual motions of (13a) either, as the incompressibility is enforced by the Lagrange multiplier.

The first variation of the energy in (13a) can be calculated as

(∀ t(X) / x(X, Y, 1) = 0) 1 0 dZ dXdY ∂W 1 ∂I 1 Î1 + ∂W 1 ∂I 2 Î2 + αẑ -p det F F -T : F -p (det F -1) = 0 (13b)
where the increments are given by F = ∇ x(X), Î1 = F T : F and Î2 = Î1 tr C -C : (F T • F ). By means of integrations by parts, this non-linear principle of virtual work can be transformed into a set of non-linear equations of equilibrium [START_REF] Biot | Mechanics of incremental deformations[END_REF]. We skip the derivation of this equivalent form of the equilibrium, as both our weakly non-linear expansions and our numerical simulations use the weak form (13) of the equilibrium.

Linear bifurcation analysis

In this section, we recall the linear bifurcation analysis of an hyper-elastic slab as done by [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF]. Doing so, we extend it to an arbitrary constitutive law and introduce a set of notation that will be useful for the subsequent non-linear analysis.

Base solution

As the gel is incompressible, the unbuckled configuration t 0 (X) remains undeformed. The only change with respect to the reference configuration is the hydrostatic pressure p 0 accounting for the weight of the material below the current point,

t 0 (α, X) = (x 0 = X, y 0 = Y, z 0 = Z, p 0 = -α Z) .
This solution satisfies the principle of virtual work (13), as well as the incompressibility and boundary conditions ( 2) and (11).

Linear bifurcation analysis

The linear bifurcation analysis considers an infinitesimal perturbation in the form

t(X) = t 0 (α, X) + t 1 (X). ( 14
)
In view of the invariance of the unbuckled solution with respect to the horizontal coordinates X and Y , it is sufficient to consider a perturbation that depends harmonically on the in-plane coordinates,

t 1 (X, Y, Z) = T 1 (Z) e ık•(X,Y ) , (15)
where ı = √ -1 is the pure imaginary number, stands for the real part, k = (k x , k y ) is a 2D wavevector, and T 1 (Z) = (x 1 (Z), y 1 (Z), z 1 (Z), p 1 (Z)) are the four complex amplitudes of the perturbations that depend on the vertical coordinate Z only.

Inserting the perturbation ( 14) into the non-linear equilibrium, we obtain the linearized equations of equilibrium in weak form as

(∀ t(X, Y, Z) / x(X, Y, 1) = 0), D 2 E(α, t 0 (α)) • [t 1 , t] = 0, ( 16
)
subject to the kinematic boundary condition

x 1 (X, Y, 1) = (x 1 , y 1 , z 1 ) (X,Y,1) = 0. ( 17
)
In equation ( 16), D 2 E stands for the second variation of the energy: this is a symmetric bilinear form with respect to its arguments in the square brackets, namely the real increment t 1 and the virtual increment t.

Inserting the special form (15) of the perturbation and integrating by parts, the problem (16-17) is transformed into a set of 12 homogeneous linear differential equations with 12 boundary conditions. There exist non-trivial solutions to these equations when the critical condition α = α * (|k|) is met, and only then; the expression of the critical wavenumber has been derived in [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF] as

α * (k) = 2k 1 + 2k 2 + cosh(2k) sinh(2k) -2k . ( 18
)
The only dependence of the critical load α * (k) on the strain energy function W is through the initial shear modulus, which is hidden in our dimensionless units, as µ = 1. The present linear bifurcation analysis therefore applies to any incompressible isotropic constitutive model.

The wavenumber k = |k| that minimizes α * is denoted by k c . It is the root of

dα * dk (k c ) = 0. ( 19a
)
In terms of k c , the critical load is obtained as

α c = α * (k c ). ( 19b
)
The numerical values are found by numerical root-finding as [START_REF] Mora | Gravity driven instability in elastic solid layers[END_REF])

α c 6.223, k c 2.120. ( 19c
)
As a result of the transverse invariance of the system, which entails rotational invariance in the (X, Y )plane, the critical condition (18) depends on the norm k = |k| of the 2-d wavevector but not on its direction. Let φ denote the direction of a critical wavevector in the (Oxy) plane, which is a free parameter for the moment,

k = (k x , k y ) = k c (cos φ, sin φ). (20)
When written in strong form, the linear bifurcation problem (16-17) can be solved analytically. This yields the expression of the marginal mode associated with the wavevector k = (k x , k y ) as

x 1 (X) = -ıξe ık•(X,Y ) f k (Z) cos φ y 1 (X) = -ıξe ık•(X,Y ) f k (Z) sin φ z 1 (X) = (ξe ık•(X,Y ) )f z (Z) p 1 (X) = (ξe ık•(X,Y ) )f p (Z) (21)
where stands for the real part, ξ is a complex amplitude, e ık•(X,Y ) = e ıkc(X cos φ+Y sin φ) is a complex phase, and f k , f z and f p are known functions capturing the dependence of the linear mode on the Z coordinate,

f k (Z) = Z sinh(kcZ)(αc sinh kc cosh kc-kc(αc+2k 2 c ))+kc(2kc(Z-1)-Z sinh(2kc)) cosh(kcZ) 2kc-sinh(2kc) f z (Z) = (sinh(2kc)-2kc)(αc+2k 2 c Z) sinh(kcZ)+kc cosh(kcZ)(2kc(2k 2 c Z+αcZ+2)-sinh(2kc)(αcZ+2) 4kc(kc-sinh kc cosh kc) f p (Z) = - kc cosh(kcZ)(4k 3 c (αcZ-2)+2α 2 c kcZ-α 2 c Z sinh(2kc))+(sinh(2kc)-2kc)(α 2 c +2k 2 c (αcZ-2)) sinh(kcZ) 4kc(kc-sinh kc cosh kc)
.

There is an arbitrary scaling factor in the functions f k , f z and f p , as the norm of the complex amplitude ξ remains unspecified in the linear bifurcation analysis. We have used the normalization convention that f z (0) = 1, which warrants that ξ can be interpreted as a complex amplitude measuring the deflection of the free surface,

z 1 (X, Y, 0) = (ξe ıkc(X cos φ+Y sin φ) ). ( 22
)
The linear mode ( 21) can be written in compact form as

t 1 (φ, ξ, X) = ξ T 1 (φ, Z)e ıkc(X cos φ+Y sin φ) (23)
where

T 1 (φ, Z) = (-ıf k (Z) cos φ, -ıf k (Z) sin φ, f z (Z), f p (Z)
) are the complex amplitudes associated with a critical wavenumber pointing in direction φ.

As the unbuckled solution is transversely isotropic, the linear modes with different directions φ in the horizontal plane all appear concurrently. Therefore, a general solution of the linear bifurcation problem is obtained by superposing linearly an arbitrary number n of pure-Fourier modes t 1 (φ j , ξ j , X), each having its own complex amplitude ξ j and wavevector direction φ j ,

t 1 (ξ 1 , • • • , ξ n , X) = n j t 1 (φ j , ξ j , X) = n j ξ j T 1 (φ j , Z)e ıkc(X cos φj +Y sin φj ) . (24)
In the following, we construct post-buckled solutions. This entails identifying the combinations of the parameters n, ξ j and φ j that satisfy the equilibrium equations at higher orders. The stability of these solutions is analyzed.

Weakly non-linear analysis of hexagonal patterns

With the aim to predict the buckling patterns in the vicinity of the bifurcation point, we push the perturbation method initiated in the previous section to second order. This expansion method is known under different names, such as 'Koiter's method'-although asymptotic expansion method of this kind have been developed well before Koiter by Lyapunov, Schmidt and others. It has been successfully used to analyze post-bifurcation in a variety of elastic structures in the past [START_REF] Koiter | On the stability of an elastic equilibrium[END_REF][START_REF] Hutchinson | Imperfection sensitivity of externally pressurized spherical shells[END_REF][START_REF] Hutchinson | Postbuckling theory[END_REF][START_REF] Budiansky | Theory of buckling and post-buckling behavior of elastic structures[END_REF][START_REF] Peek | Worst shapes of imperfections for space trusses with many simultaneously buckling members[END_REF][START_REF] Peek | Postbuckling behavior and imperfection sensitivity of elastic structures by the Lyapunov-Schmidt-Koiter approach[END_REF]. There exist several variants of Koiter's method: we use the one without imperfection and with multiple linear modes. Koiter's expansion is rederived in a self-contained way and no prior knowledge of this method is assumed. We emphasize that this method is asymptotically exact close to the bifurcation threshold, and is not simply an ad hoc upper-bound construction.

The only pattern that we will find at this quadratic order are the hexagons, hence the title of this section. Different patterns will be derived later in section 6, when we push the expansion of the energy further, to cubic order. The fact that the hexagonal patterns appear first the Koiter expansion sets them apart from the other patterns, and ultimately explains that they are the only stable patterns near the buckling threshold.

The analysis of the hexagonal pattern proceeds as follows. We first outline the principle of the expansion in §5.1. Next, we carry out the calculation explicitly, by inserting the specific form of the energy E relevant to the elastic Rayleigh-Taylor instability ( §5.2). This yields an equation for the buckling amplitude as a function of the distance to threshold, whose solutions are shown in §5.3 to describe the hexagonal pattern. In §5.4 we proceed to analyze the stability of this branch.

Outline of the Koiter expansion for a non-symmetric pattern

Following Koiter's method (van der Heijden, 2008), we derive a solution t(X) to the non-linear equilibrium equation ( 13) in the vicinity of the bifurcation point t=t 0 and α = α c , by continuing the expansion ( 14) to higher orders. For improved clarity, we introduce a small parameter η explicitly in the expansion: we rename the small perturbation from t 1 to ηt 1 , where the new t 1 is no longer viewed as a small quantity. Both the load parameter α and the solution are perturbed in terms of a small parameter η as

α = α c + ηα 1 + η 2 α 2 + • • • t(X) = t 0 (α, X) + ηt 1 (X) + η 2 t 2 (X) + • • • (25)
The quantity η is an arc-length parameter along the bifurcated branch. In the non-linear expansion for the solution t(X), the linear term is the linear bifurcation mode t 1 obtained in (24). In (25), we prescribed the zero-order term to be the unbuckled solution t 0 (α, X) evaluated at the current value of the load α, and not at the critical value of the load α c : this trick simplifies the calculation of the non-linear term t 2 considerably (van der [START_REF] Van Der Heijden | Koiter's Elastic stability of solids and structures[END_REF][START_REF] Triantafyllidis | Stability of solids: from structures to materials[END_REF]. The Koiter expansion involves inserting the expansions (25) into the non-linear equations of equilibrium ( 13) and reading it off order by order. Doing so, we obtain a hierarchy of equations, which we simplify systematically using the fact that the trivial branch, corresponding to t 0 (α, X) and t 1 (X) = t 2 (X) = • • • = 0, is a solution for any value of α.

At linear order η, we recover the linear bifurcation problem ( 16) at the critical load α c ,

(∀ t(X) / x(X, Y, 1) = 0), D 2 E(α c , t 0 (α c )) • [t 1 , t] = 0, (26)
whose general solution t 1 has been expressed in (24) as a combination of pure Fourier modes with complex amplitudes ξ 1 , . . . , ξ n . At order η 2 , we find

(∀ t(X) / x(X, Y, 1) = 0) D 2 E(α c , t 0 (α c )) • [t 2 , t] • • • • • • + 1 2 D 3 E(α c , t 0 (α c )) • [t 1 , t 1 , t] + α 1 d(D 2 E(α, t 0 (α))) dα αc • [t 1 , t] = 0. ( 27
)
This equation for the second-order correction t 2 involves the same operator D 2 E(α c , t 0 (α c )) as in the linear bifurcation problem; this is confirmed by comparing with (26). By construction, this operator is singular at the critical load α = α c : a solvability condition must be enforced before attempting to solve (27) for t 2 . It is found by considering virtual displacements t(X) matching the form of the bifurcation modes in ( 27) with arbitrary virtual amplitudes ξj :

t 1 ( ξ1 , • • • , ξn , X) = n j ξj T 1 (φ j , Z)e ıkc(X cos φj +Y sin φj ) . (28)
Inserting this virtual displacement into (27), we find that the first term vanishes by ( 26). This eliminates the unknown t 2 and we are then left with the solvability condition

(∀ ξ 1 , • • • , ξn ) 1 2 D 3 E(α c , t 0 (α c )) • [t 1 , t 1 , t 1 ] + α 1 d(D 2 E(α, t 0 (α))) dα αc • [t 1 , t 1 ] = 0. ( 29
)
In this equation, ξ 1 , • • • , ξn are arbitrary complex amplitudes, the first term involves the third variation of the energy D 3 E evaluated at the bifurcation point with 2 increments taken as the marginal mode (24) and the last increment taken as the virtual increment t1 in (28), and the last term involves the second variation D 2 E of the energy which, as implied by the notation, must be-in the following order -evaluated in the fundamental branch t 0 (α), differentiated with respect to α, evaluated in the marginally stable state α = α c , and finally evaluated with a combination of the real and virtual increments.

The left-hand side of ( 29) is linear with respect to the complex numbers ξj 's: eliminating them we find n equations that are quadratic with respect to the real amplitudes ξ j 's (recall that t 1 is linear with respect to the ξ j 's). By solving these equations for α 1 and for the ξ j 's, one determines the weakly non-linear solutions (small-amplitude patterns); using (25), the coefficient α 1 yields the slope of the bifurcated branch in the bifurcation diagram. This is the standard recipe of Koiter, and in the forthcoming sections we apply it to our particular system.

In our system, the upward and downward directions are not equivalent as the slab is not up-down symmetric: the displacement is restrained on the surface in contact with the rigid plane, and unconstrained on the opposite surface. Consistent with this observation, hexagonal patterns are analyzed based on equation (29), which is the classical equation applicable to non-symmetric systems by the Koiter method (van der [START_REF] Van Der Heijden | Koiter's Elastic stability of solids and structures[END_REF][START_REF] Triantafyllidis | Stability of solids: from structures to materials[END_REF].

Bifurcation equation

We proceed to carry out Koiter's method, as outlined above. To do so, we need to insert the energy functional E relevant to the hyper-elastic slab into (29), and to calculate the second and third variations of the energy explicitly.

Expanding each of the three modes t 1 , t 1 and t 1 , in Fourier components using ( 24) and ( 28), one can rewrite the first term in (29) as

1 2 D 3 E(α c , t 0 (α c )) • ξ j1 T 1 (φ j1 , Z)e ık j 1 (X,Y ) , ξ j2 T 1 (φ j2 , Z)e ık j 2 (X,Y ) , ξj3 T 1 (φ j3 , Z)e ık j 3 (X,Y ) ,
where k j = k c (cos φ j , sin φ j ) are the wavenumbers. A key remark is that the term in the equation above vanishes upon integration over the in-plane variables (X, Y ) ∈ R 2 unless the following resonance condition is satisfied,

k j1 + k j2 + k j3 = 0.
As the three wavenumbers have the same norm |k ji | = k c , this requires that the wavevectors form the sides of an equilateral triangle,

|φ j2 -φ j1 | = |φ j3 -φ j2 | = |φ j1 -φ j3 | = 2π 3 . ( 30
)
Therefore, the only type of mode coupling arising from equation ( 29) involves triples of wavevectors forming an equilateral triangle [START_REF] Busse | Non-linear properties of thermal convection[END_REF][START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF]. Furthermore wavevectors belonging to different equilateral triangles do not satisfy the condition (30) and therefore remain uncoupled: given the invariance of the system by rotation in the (X, Y ) plane, it is sufficient to analyze any particular triple of wavevectors, say

φ 1 = 0 φ 2 = 2π 3 φ 3 = 4π 3 n = 3. ( 31
)
We proceed to apply equation ( 29) governing the non-linear buckling amplitudes to the particular case of three modes oriented according to (31). Inserting the expansion (24) of the bifurcation modes and (28) of the virtual displacement, and evaluating the second and third variations of the energy functional (10) explicitly, we find

∀( ξ1 , ξ2 , ξ3 ) ∈ C 3 -γ(τ ) ξ 1 ξ 2 ξ3 + ξ 2 ξ 3 ξ1 + ξ 1 ξ 3 ξ2 - α 1 2 ξ1 ξ1 + ξ2 ξ2 + ξ3 ξ3 = 0, ( 32
)
where ξ denotes the complex conjugate of ξ and γ(τ ) is a known real function of the material parameter τ introduced in (5). Equation (32) has a particularly simple form, which is in fact dictated by the transverse isotropy of the unbuckled configuration-this is most easily realized by inspecting the reduced energy F, the stationarity of which implies (32), see the discussion below equation ( 41). With the help of the symbolic calculation language Wolfram Mathematica (Wolfram Research, Inc., 2014), we have calculated the function γ(τ ) relevant to the elastic Rayleigh-Taylor problem as

γ(τ ) = - k 2 c 16(sinh(2k c ) -2k c ) 3 96(4τ -1)k 3 c cosh(k c ) + 32(4τ -1)k 3 c cosh(3k c ) • • • + 12 -8τ + (8τ -1)k 2 c + 2 sinh(4k c ) -48k c -8τ + (8τ -1)k 2 c + 2 cosh(2k c ) • • • + 4k c 39(4τ -1) + 6(7 -24τ )k 4 c + (68 -320τ )k 2 c • • • + 3 -68τ + 4(24τ -7)k 4 c + 8(64τ -15)k 2 c + 17 sinh(2k c ) • • • + (4τ -1) sinh(6k c ) + 60(1 -4τ )k c cosh(4k c ) . (33)
This expression has been obtained by integrating symbolically over the depth Z a number of products of the functions f k (Z), f z (Z) and f p (Z) and of their derivatives.

With k c given by equation (19c), the numerical value of γ is found as

γ(τ ) = 1.1916 -9.1397τ. ( 34
)
This result holds for any isotropic incompressible constitutive law, as specified by the strain density function W (I 1 , I 2 ) and the material parameter τ defined in (12). In fact, it is remarkable that the only dependence on the constitutive law in (34) occurs through a linear dependence on τ = ∂W ∂I2 (0, 0): the second derivatives

∂ 2 W ∂I 2 1 (0, 0), ∂ 2 W ∂I1 ∂I2 (0, 0) and ∂ 2 W ∂I 2 2
(0, 0) do enter into the intermediate expressions of γ(τ ) but we found that they ultimately cancel out.

Having identified three-modes resonances as the primary coupling mechanism between linear modes, we proceed in the following section to solve the bifurcation equation ( 32) for the complex amplitudes ξ j . This will allow us to identify the buckling pattern as hexagons.

Amplitude of the hexagonal patterns

Equation ( 32) depends linearly on the virtual amplitudes ξj . Upon elimination of the latter, we obtain three coupled quadratic equations for the complex amplitudes (ξ 1 , ξ 2 , ξ 3 ) of the pattern,

γ ξ 1 ξ 2 = -α1 2 ξ3 γ ξ 2 ξ 3 = -α1 2 ξ1 γ ξ 3 ξ 1 = -α1 2 ξ2 ( 35
)
where the dependence of γ on τ will be implicit from now on for the sake of legibility. A similar equation has been derived by [START_REF] Hutchinson | Imperfection sensitivity of externally pressurized spherical shells[END_REF] in the analysis of spherical shells under pressure, using real variables and including the additional effect of imperfections. The form of the bifurcation equation ( 35) can be established without calculations based on the symmetries of the problem [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF]; however, the explicit expression of the coefficient γ obtained in (33), which is needed for a detailed comparison to experiments and simulations, is available from a detailed Koiter expansion only. Equation ( 35) is solved by a standard argument, given for instance by [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF]. Combining the first and second equation, we have

γα 1 |ξ 1 | 2 ξ 2 = γα 1 |ξ 3 | 2 ξ 2 , which implies that α 1 = 0, ξ 2 = 0 or |ξ 1 | = |ξ 3 |.
Repeating the argument with the two other pairs of indices, we find that

α 1 = 0, ξ 1 = ξ 2 = ξ 3 = 0 or |ξ 1 | = |ξ 2 | = |ξ 3 |.
In the first special case, α 1 = 0, two amplitudes ξ j 's out of three cancel, as can be seen from ( 35): this corresponds to the stripe pattern, which is studied later in §6. The second special case ξ 1 = ξ 2 = ξ 3 = 0 corresponds to the unbuckled branch (i.e. to the undeformed configuration in our experiments).

In this section, we focus on the last (and most generic) case, |ξ 1 | = |ξ 2 | = |ξ 3 | = 0 and α 1 = 0. It can be checked that the arguments of the complex numbers ξ 1 , ξ 2 and ξ 3 can be made identical by an appropriate change of origin of the coordinate system (X, Y ) in the horizontal plane. Therefore, we can assume without loss of generality that ξ 1 = ξ 2 = ξ 3 . It is convenient to seek these complex phases in rescaled form as where the rescaled complex amplitude ξ ∈ C is our new unknown. Inserting into (35), we obtain ( ξ) 2 = ( ξ).

ξ 1 = ξ 2 = ξ 3 = - α 1 2γ ξ, (
The solutions of this complex equation are the cube roots of 1, namely ξ ∈ {1, e 2ıπ/3 , e 4ıπ/3 }. One can check that multiplying the three amplitudes ξ j by e 2ıπ/3 amounts again to a translation of the pattern in the (X, Y ) plane, so we can assume without loss of generality that ξ = 1. We have just shown that the solution to the bifurcation equation ( 35) can be expressed in a particular orthonormal coordinate system as

ξ 1 = ξ 2 = ξ 3 = - α 1 2γ . ( 36
)
The corresponding deflection z i (X, Y ) = z(X, Y, 0) at the lower interface can be calculated with the help of the normalization condition (22) as

z i (X, Y ) = - ηα 1 2γ e ı kc X + e ı kc (-X 2 + Y √ 3 2 ) + e ı kc (-X 2 -Y √ 3 2 ) . (37)
It is plotted in figure 4. As observed from below, the pattern (37) can describe both a hexagonal network of dimples or a hexagonal network of bumps. Specifically, dimples are obtained if ηα1 γ < 0, see figure 4a, while bumps are observed if ηα1 γ > 0, see figure 4b. The sign of ηα1 γ and, thus, the type of hexagons is determined (i ) by the function γ(τ ) given in (34) in terms of the material parameter τ , and (ii ) by the signed distance to threshold ηα 1 ≈ α -α c , see (25). Both for the neo-Hookean and for Gent's constitutive models, we have τ = 0 and so γ = 1.1916 > 0, see ( 6) and ( 34): for such materials, dimples are obtained below threshold (α < α c ) and bumps are obtained above threshold (α > α c ), as depicted in figure 4c.

Of particular interest is the deflection z * of the lower surface at the center of any of the hexagonal cells, as denoted by the black dot in figure 4a-b. As can be seen from ( 37), the material point initially at the origin of the coordinate system (0, 0, 0) is mapped to the center of a hexagonal cell, and so the signed buckling amplitude reads

z * = z i (0, 0) = - 3 2 ηα 1 γ = - 3 2γ (α -α c ) . ( 38
)
With this definition, z * > 0 corresponds to a network of dimples, and z * to a network of bumps, as observed from below. Equation ( 38) has been used to sketch the bifurcation diagram in figure 4c. Based on the bifurcation equation ( 32), we have just characterized the type of pattern and the buckling amplitude in terms of the distance to threshold α -α c and of the material parameter γ = γ(τ ) defined in (34). To compare with the numerical simulations, it is convenient to use the root mean square displacement of the lower interface z rms i as a measure of the buckling amplitude. It is defined as

z rms i = z 2 i (X, Y ) = z 2 (X, Y, 0) , ( 39
)
where the angular bracket denote average over the horizontal coordinates (X, Y ). The prediction (37) of the buckling analysis can be rewritten in terms of z rms i as

z rms i = 3 1 2 α -α c 2γ 2 1/2 = √ 3 2 √ 2 |γ| |α -α c | (hexagons) ( 40
)
This prediction is verified against numerical simulations in figure 6 (dashed pink lines merging at the point of bifurcation). We have just found that amplitude z rms i of the hexagons is proportional to the load increment |α -α c | near the buckling threshold, see (38); this is a key property of a trans-critical (non-symmetric) bifurcation.

Anticipating that the bifurcated branches appearing in the bifurcation diagram in figure 4c are unstable, we will not attempt to compare them to the experimental results for the moment.

Stability of the hexagonal pattern

We now proceed to analyze the stability of the hexagonal patterns. To this end, we carry out the classical approach in stability theory outlined in Appendix A; see also §7 in the the article of [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF], and the book of [START_REF] Iooss | Elementary stability and bifurcation theory[END_REF].

For the stability analysis, we need the values of the coefficients U ijk and V ij that are defined in terms of the second and third variations of the energy of the system in equation (A.2) in the Appendix, and then the reduced energy F(ξ 1 , ξ 2 , ξ 3 ) defined in (A.5). We can save ourselves the effort of calculating them, by noticing that the quantities U ijk and V ij are the coefficients entering in the bifurcation equation ( 32), which we have already derived in explicit form. The tensors U ijk and V ij can therefore be obtained by identifying the coefficients appearing in (27) to the definition (A.2) of U ijk and V ij .

Inserting the values of U ijk and V ij identified in this way into the definition (A.5) of the reduced energy, one obtains an explicit expression of F(ξ 1 , ξ 2 , ξ 3 ) as

F(ξ 1 , ξ 2 , ξ 3 ) = - α 1 4 (|ξ 1 | 2 + |ξ 2 | 2 + |ξ 3 | 2 ) -γ (ξ 1 ξ 2 ξ 3 ). ( 41
)
The unknown buckling amplitudes (ξ 1 , ξ 2 , ξ 3 ) are identical to the generic amplitude Ξ used in the Appendix: we refer collectively to them as a complex vector Ξ = (ξ 1 , ξ 2 , ξ 3 ). The condition that F is stationary with respect to arbitrary increments Ξ = ( ξ1 , ξ2 , ξ3 ) of its arguments yields the bifurcation equation ( 32) again, as shown by the general method outlined in the Appendix, see equation (A.4).

The simple form of F in (41) reflects the symmetries of the problem, namely its invariance by permutation of the amplitudes (ξ 1 , ξ 2 , ξ 3 ) (corresponding to rotations by 2π/3 in the plane and to mirror-symmetries about a vertical plane), by a change of the phases of the ξ i 's leaving the total phase unchanged (corresponding to translations in the plane). Any elastic structure featuring transverse isotropy is therefore governed by a reduced energy of the same form as (41) near threshold [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF][START_REF] Sattinger | Group theoretic methods in bifurcation theory[END_REF]. The only element of the bifurcation equation ( 41) that is specific to the elastic Rayleigh-Taylor instability is the value of the constant γ, given in terms of the material parameter τ by (34) for our hyper-elastic slab. Different expressions of γ can be derived for different transversely isotropic systems, such as an elastic shell subjected to external pressure or a thin plate on an elastic foundation with isotropic pre-stress. Note that the sign of γ captures the preference for inward versus outward hexagonal buckling, as we will show later.

We return to the stability analysis of the different branches near the bifurcation point. We start by the fundamental branch. As recalled in Appendix A.2, its stability is governed by the second variation of F evaluated in the unbuckled state (ξ 1 , ξ 2 , ξ 3 ) = (0, 0, 0). This second variation reads

D 2 F(0, 0, 0) • [(δξ 1 , δξ 2 , δξ 3 ), (δξ 1 , δξ 2 , δξ 3 )] = d 2 F (uδξ1,uδξ2,uδξ3) du 2 u=0 = -α1 2 (|δξ 1 | 2 + |δξ 2 | 2 + |δξ 3 | 2 ).
Note that the cubic term in F has disappeared: it does not contribute to the stiffness in the unbuckled state. As a result, the quadratic form (α 1 ) -1 D 2 F(0, 0, 0) is negative definite. By the rule given at the end of Appendix A.2, the unbuckled configuration is then stable below the critical load α c , and unstable beyond α c , as sketched in figure 4c.

To address the stability of the bifurcated branch, one must analyze the eigenvalues of the second variation of F, now evaluated with ξ j given by ( 36). The second variation in the direction (δξ 1 , δξ 2 , δξ 3 ) is calculated as

D 2 F -α1 2γ , -α1 2γ , -α1 2γ • [(δξ 1 , δξ 2 , δξ 3 ), (δξ 1 , δξ 2 , δξ 3 )] = d 2 F (-α 1 2γ +uδξ1,- α 1 2γ +uδξ2,- α 1 2γ +uδξ3) du 2 u=0 = -α1 2 (|δξ 1 | 2 + |δξ 2 | 2 + |δξ 3 | 2 ) + α 1 (δξ 1 δξ 2 + δξ 2 δξ 3 + δξ 3 δξ 1 ).
When tested in different directions (δξ 1 , δξ 2 , δξ 3 ), the quadratic form D 2 F α1 appears to be either positive, zero or negative:

• with (δξ 1 , δξ 2 , δξ 3 ) = (ı, -ı/2, -ı/2) and (δξ 1 , δξ 2 , δξ 3 ) = 0, ı √ 3/2, -ı √ 3/2 , we find 1 α1 D 2 F -α1 2γ , • • • • [(δξ 1 ,
δξ 2 , δξ 3 ), (δξ 1 , δξ 2 , δξ 3 )] = 0 (these increments correspond to an infinitesimal translation of the pattern in the x or y directions, which indeed leaves the energy invariant);

• with (δξ 1 , δξ 2 , δξ 3 ) = (1, 1, 1), we have

1 α1 D 2 F -α1 2γ , • • • • [(1, 1, 1)
, (1, 1, 1)] = 3 2 ;

• with (δξ 1 , δξ 2 , δξ 3 ) = (ı, ı, ı), we have

1 α1 D 2 F -α1 2γ , • • • • [(ı, ı, ı), (ı, ı, ı)] = -9 2 .
When evaluated in the bifurcated branch, D 2 F α1 has thus both positive and negative eigenvalues. By the rule given at the end of Appendix A.2, this implies that the bifurcated branch is unstable on both sides in the neighborhood of the bifurcation point, see figure 4c. More accurately, the bifurcated branch is unstable by an amplitude mode (δξ 1 , δξ 2 , δξ 3 ) = (1, 1, 1) below the critical load α c , and by a phase mode (δξ 1 , δξ 2 , δξ 3 ) = (ı, ı, ı) beyond the critical load α c (recall that ı is our notation for the pure imaginary number).

To sum up, we have characterized the stability of the fundamental branch and of the hexagonal branch in the neighborhood of the bifurcation point: as summarized in figure 4c, the stable branch is stable below the critical load α and unstable above it, while the hexagonal branch is always unstable.

These results are now extended in two ways: in §6 we derive additional patterns, such as stripes and squares, and in §8 we discuss the patterns and their stability far from the bifurcation point.

Weakly non-linear analysis of symmetric patterns

Outline of the expansion for stripes

In this section, we return to the bifurcation equation ( 35) and consider the stripe pattern (known as the roll pattern in Benard convection). It is possible to re-use the previous analysis by viewing the stripe pattern as a special case of a three-modes pattern for which two out of the three complex amplitudes cancel, say ξ 2 = ξ 3 = 0. The linear increment α 1 of the load along the stripe branch can then be calculated simply by inserting ξ 2 = ξ 3 = 0 into the bifurcation equation ( 35) obtained previously. The result is

α str 1 = 0. ( 42
)
where the label 'str' refers to 'stripe'. A consequence of ( 42) is that the bifurcation equation ( 29) classically used for the analysis of non-symmetric systems is identically zero, and cannot be used to determine the amplitude of the stripe pattern; to do so, we will need to push the Koiter expansion to the following order. Another consequence is that the stripe branch is parabolic in the bifurcation diagram, as the expansion of the load in terms of the arc-length parameter η along the stripe branch reads

α = α c + η 2 α 2 + • • • , see (25).
All these features point to the fact that the stripe pattern is effectively a symmetric bifurcation, even though it takes place in a non-symmetric system. Choosing the orientation of the X axis so that it is aligned with the wavevector of the stripe mode, we can describe the stripe mode as

t 1 (ξ 1 , X) = ξ 1 T 1 (φ 1 = 0, Z)e ıkcX .
The lower interface displays parallel sinusoidal stripes (undulations) that are invariant in the Y direction. This pattern is cylindrically symmetric in the direction Y parallel to the stripes. We proceed to solve the Koiter expansion order by order for the stripe pattern. To do so, we set α 1 = 0, ξ 2 = ξ 3 = 0 in the equations derived previously for a generic perturbation with three modes. The solvability condition (32) or ( 35) is automatically satisfied. Equation ( 27) can therefore be solved for the second-order correction to the solution t 2 . With α 1 = 0, this yields

(∀ t(X) / x(X, Y, 1) = 0) D 2 E(α c , t 0 (α c )) • [t 2 , t] + 1 2 D 3 E(α c , t 0 (α c )) • [t 1 , t 1 , t] = 0.
Inserting the detailed expression of the linear mode T 1 given in ( 21), the solution t 2 of this equation is obtained explicitly as the sum of three contributions: (i ) the squared pattern amplitude (ξ 1 ) 2 times known functions of Z, (ii ) (ξ 1 ) 2 times the second harmonic e 2ıkcX times known functions of Z, and (iii ) arbitrary combinations of linear modes T 1 (φ, Z), which play no role in the following as they represent a small correction to the sinusoidal pattern described by ξ 1 .

As mentioned above, the bifurcation equation ( 29) is satisfied automatically. To derive the amplitude equation relating the pattern amplitude ξ 1 and the load increment α 2 , the non-linear equilibrium ( 26) must be expanded one order further than earlier, namely to order η 3 . This yields an equation for t 3 , whose solvability condition is again obtained by considering a special virtual motion of the form t(X) = t1 ( ξ1 , • • • , ξn , X) = ξ 1 T 1 (φ 1 = 0, Z)e ıkcX . This eliminates t 3 from the equation, and we obtain

(∀ ξ 1 ∈ C) 1 6 D 4 E(α c , t 0 (α c ))•[t 1 , t 1 , t 1 , t1 ]+D 3 E(α c , t 0 (α c ))•[t 2 , t 1 , t1 ]+α 2 d(D 2 E(α, t 0 (α))) dα αc •[t 1 , t1 ] = 0, ( 43
) This bifurcation equation is usually applicable to symmetric systems [START_REF] Van Der Heijden | Koiter's Elastic stability of solids and structures[END_REF][START_REF] Triantafyllidis | Stability of solids: from structures to materials[END_REF]. Even though our system is not symmetric, as noted at the end of §5.1, this equation is applicable to the stripe pattern, which is effectively symmetric: changing the sign of the amplitude ξ 1 amounts to a translation of the stripe pattern by a half-period, and this leaves the energy invariant. By contrast, changing the signs of the three amplitudes that make up a hexagonal pattern transforms a network of bumps into a network of dimples and vice versa, and this does not leave the energy invariant. This is what explains the special role played by hexagons: being non-symmetric, they are described by a different bifurcation equation than the stripe patterns, and they stand out from the other (symmetric) patterns.

Amplitude of the stripe pattern

An explicit form of the amplitude equation can now be obtained calculating the second-order solution t 2 in terms of (ξ 1 ) 2 , e 2ıkcX and Z, and then inserting it into (43). This calculation is not difficult, but it would be tedious to do by hand. By a similar approach as that used for hexagons earlier, see (32), and with the help of the symbolic calculation language Wolfram Mathematica (Wolfram Research, Inc., 2014), we could rewrite (43) as

(∀ ξ 1 ∈ C) ω(τ )|ξ 1 | 2 ( ξ1 ξ1 ) -α 2 2 ( ξ1 ξ1 ) = 0, ( 44
)
where ω(τ ) is a real function depending solely on the material parameter τ . An explicit expression of the function ω(τ ) is available as the outcome of the symbolic calculation, but it is too complicated to be included here. It can easily be evaluated numerically for any particular value of the material parameter, however.

For materials such that τ =0, for example, which includes both neo-Hookean and Gent materials,

ω(τ = 0) = 0.22450.
Eliminating the virtual quantity ξ1 , we can rewrite the amplitude equation above as

|ξ 1 | 2 - α 2 2ω(τ ) ξ 1 = 0. ( 45
)
The root ξ 1 = 0 corresponds to the unbuckled branch, and the other root |ξ 1 | 2 = α2 2ω(τ ) corresponds to the stripe pattern.

In view of the normalization condition ( 22), the stripes are associated with a deflection at the lower interface z i (X, Y ) = z(X, Y, 0) given by z i (X, Y ) = 0 + ηξ 1 e ı kc X , hence the amplitude of the stripes

|ηξ 1 | = η 2 |ξ 1 | 2 = (η 2 α2) 2ω(τ0) 1 √ 2ω(τ0) (α -α c ) 1/2
. The square-root dependence on the distance to threshold is typical of a symmetric (pitchfork) bifurcation.

Using the numerical value ω(τ = 0) = 0.2245 relevant to neo-Hookean and Gent materials, the buckling amplitude of the stripe reads

z rms i = |ηξ 1 | √ 2 = 1.055 (α -α c ) 1/2 (τ = 0, stripes). ( 46
)
Given that ω(τ = 0) > 0, the square root is defined for α > α c , which means that the bifurcated branch lies above the critical load α c , as in a super-critical pitchfork bifurcation. Unlike in a typical pitchfork bifurcation, however, this branch is unstable, as explained below. This important fact has not been appreciated in earlier work on the elastic Rayleigh-Taylor instability, to the best of our knowledge. The analytical prediction of the buckling amplitude (46) has been verified using finite-element simulations (see the pink dashed curve in figure B.8b).

Stripes are unstable by hexagonal perturbations

The stability analysis of the stripes follows the general approach recalled in Appendix A.2. Based on the arguments given in §5, the most effective perturbations are those involving three Fourier modes whose wavevectors form the sides of an equilateral triangle. If none of these three wavevectors are aligned with the wavevector of the stripe pattern, no coupling occurs: we focus on the interesting case where one of the three wavevectors that makes up the perturbations coincides with the wavevector of the stripe pattern. Effectively, this amounts to taking a linear combination involving three modes as earlier in (31), with ξ 1 given by the amplitude equation ( 45) plus a small perturbation, while ξ 2 and ξ 3 are pure perturbations. We can therefore re-use the stability analysis done in §5.4, now setting α 1 = 0, as well as ξ 2 = ξ 3 = 0 in the base configuration.

In particular, we can reuse the reduced energy (A.5) relevant to three-modes configurations, now dropping the quadratic term proportional to α 1 = 0. In view of (41), we are left with

F str (ξ 1 , ξ 2 , ξ 3 ) = -γ (ξ 1 ξ 2 ξ 3 ). ( 47
)
In view of the bifurcation analysis of §6.2, the stripe branch is described by ξ str 1 = α2 2ω , ξ str 2 = ξ str 3 = 0. The second variation of the reduced energy F str along this branch as

D 2 F str (ξ str 1 , 0, 0) • [(δξ 1 , δξ 2 , δξ 3 ), (δξ 1 , δξ 2 , δξ 3 )] = -2γ (ξ str 1 δξ 2 δξ 3 ) = -2γ α 2 2ω (δξ 2 δξ 3 )
One can check that the perturbation χ = (δξ 1 , δξ 2 , δξ 3 ) = (0, 1, 1) is an eigenvector of this symmetric operator associated with a negative eigenvalue σ 1 = -γ α2 2ω . As explained in the Appendix, any such eigenvalue σ 1 corresponds to an eigenvalue σ of the stiffness operator D 2 E that varies along the bifurcated branch as

σ = ησ 1 = -γ (η 2 α2) 2ω = -γ √ 2ω √ α -α c .
For both neo-Hookean and Gent materials, this eigenvalue reads, from ( 34) and ( 6

.2), σ = -1.779 √ α -α c (τ = 0)
and it is negative (unstable case). Note that the square root is defined on the side α > α c of the bifurcation point, which is indeed where the strip branch lives, see ( 46). We have just shown that the stripe pattern is unstable in the vicinity on the bifurcation point, as there exists a negative eigenvalue σ of the stiffness operator D 2 E on this branch. The corresponding perturbation tends to produce a hexagonal pattern: when added to the base solution (ξ str 1 > 0, 0, 0), the perturbation (0, 1, 1) tends to make all three amplitudes (ξ 1 , ξ 2 , ξ 3 ) equal, and this corresponds to hexagons, see (36).

Squares and other symmetric patterns are unstable as well

In the above stability analysis, the only property of the stripes which we have actually used is α 1 = 0. This property warrants that the reduced energy F, relevant to the vicinity of the bifurcation point, contains a cubic term only, see (47), and the cubic term, unsurprisingly, gives rise to unstable patterns. Other symmetric patterns are possible, such as an array of squares combining two buckling modes having orthogonal wavenumbers. Squares are symmetric, as a change of sign of the amplitude of the pattern amounts to a combined translation by a half-wavelength in both the directions of the wavenumbers. By a similar argument as for the stripes, see §6.3, squares and other symmetric patterns satisfy α 1 = 0, their reduced energy is cubic, and these patterns are unstable close to bifurcation.

Summary of the weakly non-linear bifurcation analysis

The weakly non-linear analysis of the hexagonal and straight patterns is summarized in the bifurcation diagram in figure 5. The branches corresponding to symmetric patterns such as stripes and squares have α 1 = 0 and therefore have an initially parabolic shape in the vicinity of the bifurcation point, see (25). The parabola fitting the stripes branch has been calculated in (46), and as an outcome of this calculation we found that the stripes exist above the critical load, for α > α c . We have carried out a similar postbuckling analysis for the square pattern, but it is omitted as this branch is unstable as well; the outcome of the analysis is that the squares live on the other side of the critical load, for α < α c , as sketched in the bifurcation diagram in figure 5.

The hexagonal pattern stands out from the other patterns, due to the three-modes resonances analyzed in §5.2: hexagonal perturbations make all patterns unstable, and the hexagonal branch has a finite slope α 1 = 0 near the bifurcation point, see (38), while the other branches have a parabolic shape. Note that when one plots the unsigned buckling amplitude as in figure 5, the hexagonal branch from figure 4c folds up into two half-lines meeting at the bifurcation point.

The main conclusion of the weakly post-buckling analysis is that no small-amplitude bifurcated solution exists near the bifurcation point. This explains why in the experiments the system jumps to a finiteamplitude buckled solution, as reported in figure 2. These finite-amplitude patterns cannot be characterized based on a perturbative method such as a Koiter expansion, and different approaches are used. All the other branches are unstable (dashed/dotted curves). An unsigned buckling amplitude is used on the vertical axis, which makes the hexagonal branch appear as a broken line.

Finite-amplitude patterns

In this section, we analyze patterns of finite amplitudes. To this end, we set up numerical simulations of the hexagonal pattern using the finite-element method first ( §8.1) and we use a simplified mathematical model to gain insights into the competition of the various patterns far from threshold ( §8.2).

Numerical simulations of the finite-amplitude hexagons

With the aim to tackle the fully post-buckled regime, we have set up numerical simulations of the equilibria of hyper-elastic slab (7) using the non-linear finite element method. In this section, we focus on simulations of the hexagonal pattern. Additional numerical results are given in Appendix B for the stripe and square patterns.

Our simulations make use of the neo-Hookean strain energy W NH in (3). The incompressibility is enforced by means of a mixed formulation, and the unknowns are the components of the displacement as well as the pressure p. The implementation is based on the open-source finite-element library FEniCS [START_REF] Logg | Rayleigh-bénard convection: thirty years of experimental, theoretical, and modeling work[END_REF].

Both the shear modulus and the thickness are set to one, µ = 1, h = 1, following the same conventions as in the analysis, see (9). The simulation domain is a hexagonal prism with thickness h = 1, as shown in figure 6a. The width (edge-to-edge distance) of the hexagonal base of the prism is set to = 4π √ 3k , with periodic boundary conditions applied on the opposite edges of the hexagon. The simulation domain is the generating cell of a periodic hexagonal pattern. The cell size depends on a parameter k, which we can freely adjust in our simulations. In most of the forthcoming simulations, k has been set to the critical value k c = 2.120 predicted by the linear bifurcation analysis, see (19c), but we have varied k in some simulations to capture any evolution of the wavelength of the hexagonal network in the post-buckled regime. Numerical convergence to hexagonal patterns is achieved by introducing a bias featuring hexagonal symmetry in the initial guess passed to the non-linear root-finding algorithm.

A set of simulations of the hexagonal pattern obtained by varying the load parameter α, while locking the wavelength to its critical value k = k c are shown in figure 6. The root mean square buckling amplitude on the lower face z rms i , defined in (39), has been extracted from the simulation, and is used to construct the bifurcation diagram shown in figure 6b. Near threshold, the simulations confirm the results of the bifurcation analyses ( §4-5): the critical load α c and the weakly non-linear amplitude of the hexagons (dashed segments meeting at the bifurcation point in the figure) are correctly predicted. The small-amplitude hexagonal pattern is a network of dimples below α c , and a network of bumps above α c , in accord with the weakly non-linear expansion, see figure 5.

The numerical simulations reveal the presence of a fold point F along the branch of hexagonal dimples, at α = α F (k c ) 5.79. As a result, there are two solutions corresponding to a hexagonal network of dimples above α F (k c ) 5.79: a small amplitude one, shown as C in the figure, whose amplitude goes to zero at α c , and a large amplitude one, shown as D in the figure. The finite-amplitude pattern of dimples is the only one that is stable near the bifurcation threshold, as we will learn from the next section: this is the one that 2) with the finiteelement simulations with k = kc = 2.12 (grey curve: solid part is stable, dashed part is unstable according to §8.2), with no adjustable parameter. Note that the experimental data-points accumulate close to the stable part of the numerical curve.

is observed in the experiments. The existence of this finite-amplitude hexagonal solution is fully consistent with the fact that all the small-amplitude patterns are unstable with respect to hexagonal perturbations in the Koiter post-bifurcation expansion, see §5.4 and §6.4.

A comparison of the finite-element simulations with the buckling amplitude w b measured in experiments is shown in figure 6d (recall that w b is defined as the largest downward displacement measured on the lower face of the slab). The comparison involves no adjustable parameter and shows qualitative agreement. The simulations predict that no hexagonal pattern exists below the load value α F (k c ) = 5.79 corresponding the fold point, and this value is indeed close to the load α = 5.67 at which the transition to hexagons takes place in the experiments. Even though the simulations successfully captures the salient features of the buckling instability, the agreement in figure 6d is not perfect; the discrepancy can likely be attributed to the finite size of the experimental slab (while the simulations assume an infinite domain), and the difficulty of reliably measuring the constitutive law of the gel (especially for the large strains of several hundreds of percent predicted by the simulations).

In an additional set of simulations (data not shown), we have relaxed the constraint k = k c , i.e. we have varied the width = 4π/ √ 3k of the unit cell. We have found that the optimal value of k does not depart from k c by more than 10% for the range of loads shown here, suggesting that the variations of k have a limited influence on the bifurcation diagram.

Analysis based on a toy energy functional

In this section, we introduce a simple energy functional that captures finite-amplitude patterns. This alternative approach to the finite-element method allows for a simple, albeit approximate, analysis the stability of the patterns and of their competition. The energy functional is obtained by extrapolating to finite amplitudes the perturbative analyses of §5-6. The extrapolation cannot be justified rigorously, but it is expected to capture the main features of the bifurcation diagram qualitatively.

To start with, we rewrite the reduced energy F governing hexagons and stripes in non-scaled form:

F * = -α-αc 4 (|ζ 1 | 2 + |ζ 2 | 2 + |ζ 3 | 2 ) -γ (ζ 1 ζ 2 ζ 3 ) (hexagons) -α-αc 4 |ζ 1 | 2 + ω 4 |ζ 1 | 4 (stripes).
These expressions have been derived by restoring the factor η 3 that had been scaled out from the energy F for hexagons in (41), i.e. F * = η 3 F, and the factor η 4 that had been scaled out from the energy F for stripes in (44), i.e. F * = η 4 F; we have also identified ηα 1 ≈ α -α c as the load increment for hexagons, η 2 α 2 ≈ α-α c as the load increment for stripes (for which α 1 = 0), and ζ j = ηξ j as the non-scaled amplitudes of the Fourier components of the pattern, see (25).

The terms in F * that are quadratic in the buckling amplitudes ζ j agree exactly for hexagons and stripes, set aside the obvious difference that ζ 2 = ζ 3 = 0 for stripes. Indeed, the different Fourier components are uncoupled at this order and these contributions are just the sum of the quadratic energy -α-αc 4 |ζ j | 2 of each of the relevant Fourier modes.

Limiting attention to patterns made up of three Fourier modes with complex amplitudes (ζ 1 , ζ 2 , ζ 3 ), and whose wavevectors all have the same norm k c and form an equilateral triangle (this includes hexagons but also stripes as a particular case), we postulate the following energy, based on the expression of F just derived,

G(α, ζ 1 , ζ 2 , ζ 3 ) = - α -α c 4 (|ζ 1 | 2 + |ζ 2 | 2 + |ζ 3 | 2 ) -γ (ζ 1 ζ 2 ζ 3 ) . . . + ω 4 (|ζ 1 | 4 + |ζ 2 | 4 + |ζ 3 | 4 ) + σ(|ζ 1 | 2 |ζ 2 | 2 + |ζ 2 | 2 |ζ 3 | 2 + |ζ 3 | 2 |ζ 1 | 2 ).
This postulated energy has the following properties. It is consistent with the reduced energy F * for hexagons derived by a rigorous expansion, up to order |ζ j | 3 included. It yields exactly the reduced energy F * for stripes when we set ζ 2 = ζ 3 = 0. It is invariant by the symmetries of the system, namely by rotations of ±2π/3 and mirror-symmetries (which amount to a permutation of the amplitudes ζ j ) and by translations by a vector (τ x , τ y ) in the horizontal plane (which amount to incrementing the complex phases of the ζ j 's by

k c τ x , k c -τx 2 + τ y √ 3 3 and k c -τx 2 + τ y √ 3 3
, respectively). In fact, it can be shown that the proposed form of G(ζ 1 , ζ 2 , ζ 3 ) is the most general one that is polynomial of order 4 in terms of the real and imaginary parts of the complex amplitudes ζ j , and is consistent with the symmetries of the system. In this expression of G, we use the values of the parameters determined earlier by the rigorous buckling analysis: focusing on a neo-Hookean or Gent material, we set α c = 6.223 from (19c), γ = γ(τ = 0) = 1.1916 from (34) and ω = ω(τ = 0) = 0.2245 from (6.2). Then σ is the only parameter left. By trial and errors, we found that the value σ = 0.4 reproduces well the branch of hexagons found the by finite-element method in §8.1: this is the value we use in the following.

The benefit of introducing the quartic term (proportional to σ) in our toy-energy G is that the energy becomes bounded from below for large buckling amplitudes (provided σ > -ω/2, an inequality which is indeed satisfied with our set numerical values). The energy F * for hexagons, by contrast, was not bounded from below due to the cubic term. The new term proportional to σ which we have introduced by hand in the definition of G can be viewed as an simple and approximate way to capture the next orders in the Koiter expansion-in practice deriving these higher orders would be a daunting task.

With the aim to capture the finite-amplitude buckling patterns, we proceed to derive a bifurcation diagram from the toy-energy G. To do so, we solve the equilibrium ∂G ∂ζj (α, ζ 1 , ζ 2 , ζ 3 ) = 0, with 1 j 3, with respect to the amplitudes ζ j 's for fixed values of the load α. The analysis is very similar to that done earlier in §5.3. The equilibrium writes, in explicit form,

γ ζ 1 ζ 2 = -α-αc 2 + ω|ζ 3 | 2 + 2σ(|ζ 1 | 2 + |ζ 2 | 2 ) ζ3 γ ζ 2 ζ 3 = -α-αc 2 + ω|ζ 1 | 2 + 2σ(|ζ 2 | 2 + |ζ 3 | 2 ) ζ1 (48) γ ζ 3 ζ 1 = -α-αc 2 + ω|ζ 2 | 2 + 2σ(|ζ 3 | 2 + |ζ 1 | 2 ) ζ2
Multiplying these three equations sides by sides, one can prove that any solution of ( 48) is such that -γζ j + (ω + 4σ)ζ 2 j = 0, which is consistent with the amplitude equation derived earlier in (36) in the limit where the amplitudes ζ j are small, i.e. when the quadratic term (ω + 4σ)ζ 2 j is dropped out. For the stripes, (48) yields exactly the same amplitude equation -α-αc 2 + ωζ 2 1 = 0 as earlier in (45). A systematic analysis of (48) shows that yet another pattern is possible, which we call the modulated stripes: this pattern is such that ζ 1 = ζ 2 = ζ 3 and looks like a hybrid of stripes and hexagons, as shown in the inset in the upper-right hand side corner of figure 7.

(ζ 1 ζ 2 ζ 3 ) 3 is real. A similar
By calculating the Hessian of G, which is a 3 × 3 complex matrix or, equivalently, a 6 × 6 real symmetric matrix, one can determine the stability of the different branches of solutions. The result is shown in figure 7, where thick, solid lines denote stable solutions and thin, dashed lines denote unstable solutions. As the energy G is bounded from below, for any value of the load α there is at least one, and sometimes two, patterns that are stable.

In the neighborhood of the critical point (α c , 0), the bifurcation diagram is consistent with the non-linear buckling analysis: all branches are unstable in this neighborhood, except for the part of the unbuckled branch located below the critical load, which is stable. This is in full agreement with figure 5, as could be expected from the fact that the toy-energy G extrapolates by design the perturbative buckling analysis.

Zooming out and considering now on the lower-left hand side corner of the diagram, we see that its predictions are fully consistent both with the experimental findings reported in figure 1 and with the finite element simulations shown in figure 6: thanks to the quartic term included in G, the branch of hexagons corresponding to a network of dimples that emanates from the bifurcation point (α c , 0) passes through a fold point F , where the hexagons branch become stable again. As a result, the only stable pattern present immediately above α c is the hexagonal network of dimples, and the latter has a finite amplitude at α = α c : this explains the discontinuous transition observed in the experiments. This pattern continues to exist for values of the load parameter α significantly below α c , down to the fold point F . As the actual value α F depend on the value of free parameter σ, the present model cannot make a prediction on the value of α F ; this value happens to be α F = 5.83 for the particular value σ = 0.4 chosen here. In the interval α F α α c both the planar and hexagonal patterns are stable. Further analysis (data not shown) reveals that the pattern with lowest energy switches from being the planar pattern (for lower values of α) to the hexagonal pattern (for larger values of α), and this happens at a specific value of the load in this interval. A similar exchange of the configurations having lowest elastic energy takes place in the finite-element simulations as well, see figure 6c.

When the load α is increased, the hexagons remain the only stable pattern until a much larger value α ≈ 8 is reached. At this point, stripes (of finite amplitude) become stable, along with the hexagons. It is unclear whether the predictions of the toy-energy are accurate this far from the bifurcation threshold, however, and these findings would need to be confirmed by additional numerical simulations. In any case, the toy model highlights the possibility that branches which are unstable near the bifurcation point re-stabilize far from threshold.

Note that the square pattern does not appear from this bifurcation diagram, as we have limited attention to patterns made up of three waves when we introduced the toy-energy G. For the same reason, the stability of the branches in figure 7 is with respect to perturbations comprising three waves only.

Overall, the bifurcation diagram in figure 7 bears a striking resemblance with the bifurcation diagram derived figure 4 of [START_REF] Busse | Non-linear properties of thermal convection[END_REF] in the context of convective instabilities in fluids. This is not surprising, as the form of our toy energy is severely constrained by the symmetries of the system, which are identical in the two systems, see for instance [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF].

The toy energy analyzed in this section has been proposed based on a crude extrapolation of the reduced energy governing the Koiter expansion. Yet, its predictions are fully consistent with all the features of the bifurcation that have been gathered so far from the experiments, from the bifurcation analysis and from the numerical simulations. Most notably, this includes the formation of finite-amplitude hexagons at the onset of bifurcation, and the presence of hexagonal patterns for values of the load parameter α below the critical value α c predicted by the linear bifurcation analysis.

Discussion and conclusion

We have addressed the pattern selection of the elastic Rayleigh-Taylor instability based on a different methods, including a weakly non-linear Koiter expansion. Owing to the transverse isotropy of the unbuckled configuration, the dominant coupling mechanism between linear modes takes place through three-modes resonances. This singles out the hexagonal pattern, which ends up being the only stable pattern such that the load remains close to the bifurcation load: stripes and squares, for example, are unstable with respect to hexagonal perturbations.

The hexagonal pattern appears by a trans-critical bifurcation, as illustrated in the bifurcation diagram in figures 5 and 7. It is qualitatively different from the pitchfork bifurcations commonly found in systems that are up-down symmetric. The branch of hexagons is unstable in the vicinity of the bifurcation point on both sides, implying that the transition from unbuckled to hexagons is discontinuous. A stable branch of hexagons of finite amplitude exists at the critical load, see figure 7, and persists down to a value of the load α F significantly smaller than the critical load α c predicted by the linear bifurcation analysis.

The hexagonal patterns have a preferred direction: for both the neo-Hookean and Gent constitutive laws, we have shown that the centers of the hexagons tend to buckle inwards (i.e. upwards), forming a hexagonal network of dimples when observed from below. This preference is consistent with the absence of up-down symmetry in the system: the slab is hanging below the rigid plane. As far as symmetries are concerned, there are two requirements for three-modes resonances to appear (and, hence, for hexagons to be selected): the presence of transverse isotropy and the absence of up-down symmetry.

Overall, these results agree well with the experimental findings reported in §2.

Since these features of the bifurcation diagram are mostly a consequence of the symmetries, other transversely isotropic elastic structures can be expected to produce hexagons by a similar mechanism, and to display a similar bifurcation diagram. We provide a quick critical review of the literature on the selection of buckling patterns in transversely isotropy structures in the light of our results.

In his post-bifurcation analysis of an elastic slab under equi-biaxial compression, Ciarletta (2014) derives a branch of hexagonal solutions but the expansion is based on scaling assumptions that are relevant to symmetric (pitchfork) bifurcations. As a result, the three-modes resonances are suppressed, the transcritical character of the bifurcation is lost, i.e. the buckling amplitude is no longer proportional to the load increment as it is in (40), and the preference for a hexagonal network of dimples over a hexagonal network of bumps is not accounted for. Squares and stripe patterns are derived as well in this work but their unstable character, which involves three-waves resonances, is not mentioned. In a recent work on the same system, Jia and Ben Amar (2013) discuss the hexagonal branch and nicely point out the existence of a transcritical bifurcation; solutions for the stripes and square patterns are obtained as well. Our results show that the latter are actually unstable close to bifurcation.

A thin film under equi-biaxial compression resting on an elastic foundation is another well-studied buckling problem featuring transverse isotropy. In their post-buckling analysis, [START_REF] Audoly | Buckling of a thin film bound to a compliant substrate (part 1). Formulation, linear stability of cylindrical patterns, secondary bifurcations[END_REF] address the hexagonal pattern but they overlook the three-modes resonances and end up with a pitchfork (symmetric) bifurcation; as a result of this, the energy of the hexagonal pattern is erroneous and they incorrectly identify the squares as the pattern of lowest energy near threshold. The same system has been analyzed by [START_REF] Cai | Periodic patterns and energy states of buckled films on compliant substrates[END_REF]; [START_REF] Breid | Effect of stress state on wrinkle morphology[END_REF], who report that some predictions of their buckling analysis are at odds with experimental findings: in their theory the lowest energy pattern appears to be the squares at low overstress, although the experiments reported in the paper show an 'overwhelming' preference towards hexagons; there is no predicted preference for inward over outward buckling in the hexagonal patterns although the experiments show a marked preference for inward buckling. These unexplained features seen in these experiments are qualitatively similar to those analyzed here in the context of the elastic Rayleigh-Taylor instability, and the two phenomena may have a common explanation.

The stability of hexagonal, stripe and square patterns arising in transversely isotropic structures has been discussed based on a general and powerful set of methods for analyzing bifurcations in the presence of symmetries [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF]; [START_REF] Golubitsky | Singularities and groups in bifurcation theory[END_REF]; [START_REF] Buzano | Bifurcation on the hexagonal lattice and the planar bénard problem[END_REF]. As in the more recent work on thin films under in-plane compression, these authors restrict their stability analyses to perturbations featuring the same set of eigenmodes as the base pattern. Roll patterns (stripes), for instance, are viewed as a degenerate form of square patterns (with one of the complex amplitudes set to zero), and tested against perturbations featuring square symmetry only: it is concluded that stripe patterns could be stable near the bifurcation point-see for instance table 6.1 in the paper by [START_REF] Sattinger | Group representation theory, bifurcation theory and pattern formation[END_REF]-, although the argument presented in §6.3 and 6.4 shows that they are unstable with respect to hexagonal perturbations.

The following limitations of our analysis will hopefully be addressed in future work. Based on input from the experiments, we have focused on the hexagonal patterns; a careful analysis of the symmetry group associated with transverse isotropy shows that other bifurcation branches are possible, corresponding to non-hexagonal pattern [START_REF] Golubitsky | Symmetries and pattern selection in rayleigh-bénard convection[END_REF]. A complete classification of these branches together with their stability could be addressed in future work. Note that the hexagonal pattern will still appear at the quadratic order in the Koiter expansion, i.e. before any other pattern, including those we have left out from the present analysis: this suggests that the pattern that we have left out from the present analysis will be dominated by hexagons, much like the stripes. Besides, we have not taken into account the finite extent of the slab in the horizontal direction: in a slab of large but finite size, the perturbations caused by the boundaries may be addressed by deriving an amplitude equation, as has been done for fluid convection [START_REF] Pomeau | Front motion, metastability and subcritical bifurcations in hydrodynamics[END_REF]. Finally, the strongly post-buckled regime α 1 could be investigated in future work: in this regime, it may be possible to obtain an asymptotic description of the pattern consisting of thin and severely stretched walls surrounding deep dimples.

Acknowledgments

This paper was prepared using the T E X macs typesetting system [START_REF] Van Der Hoeven | Gnu texmacs: a scientific editing platform[END_REF]; we are grateful to the developers for making this outstanding tool freely available. We are also grateful to the anonymous reviewers for their valuable comments on the initial version of this manuscript.

Figure 1 :Figure 2 :

 12 Figure 1: (a) Principle of the experiment: a gel is prepared in a container having a square base, which is turned upside down once the gel is cross-linked. The side of the square is = 40 cm and gel thickness is h = 2.5 cm. (b,c) Pictures, taken from below, of the lower surface of the gel. For a higher shear modulus µ = 45 Pa as in (b), the surface remains flat. For a lower shear modulus µ = 37 Pa as in (c), a close-to-hexagonal pattern appears spontaneously as soon as the sample is flipped upside down.

Figure 4 :

 4 Figure 4: (a,b) Hexagonal pattern predicted by the post-bifurcation analysis, see equation (37): observing the lower face from below, one sees a hexagonal pattern made up (a) of dimples if ηα 1 γ < 0 or (b) of bumps if ηα 1 γ > 0. (c) Bifurcation diagram showing the signed buckling amplitude z as a function of the load α for a material such that γ = 1.1916 -9.1397τ > 0 (which includes the case of a neo-Hookean or a Gent material, for which τ = 0), based on equation (38). The stability of the branches in the neighborhood of the bifurcation point is based on the results of §5.4: solid lines correspond to stable equilibria, and dotted lines to unstable equilibria.

Figure 5 :

 5 Figure 5: Sketch of the bifurcation diagram in the neighborhood of the bifurcation point αc, as obtained by the weakly nonlinear analysis, see §5 and §6. Only the part of the unbuckled pattern below the first critical load (α < αc) is stable (solid line).All the other branches are unstable (dashed/dotted curves). An unsigned buckling amplitude is used on the vertical axis, which makes the hexagonal branch appear as a broken line.

Figure 6 :

 6 Figure 6: Non-linear finite-element simulation of the hexagonal pattern. (a) The simulation domain, shown in light grey, is a hexagonal prism with periodic boundary conditions; its width = 4π √ 3k is constrained as k = kc = 2.12 here, in accord with the prediction of the linear bifurcation analysis. (b) Bifurcation diagram showing the unsigned buckling amplitude z rms i rescaled by the initial thickness h as a function of the load parameter α. The dashed lines emanating from the bifurcation point α = αc show the hexagonal branch predicted by the weakly non-linear analysis, see equation (40). Insets: visualization of numerical solutions, flipped upside down to aid visualization (the free surface is visible while the surface with blocked displacement is hidden at the back). These images have been obtained by replicating the simulation cell multiple times to aid visualization. The colors encode the magnitude of the vertical displacement, as indicated in the scale bar. (c) Total potential energy of the pattern per unit area as a function of the load parameter α. The zero of energy is chosen as the unbuckled configuration. (d) Comparison of the experimental buckling amplitude w b of the gel surface (symbols, same as in figure2) with the finiteelement simulations with k = kc = 2.12 (grey curve: solid part is stable, dashed part is unstable according to §8.2), with no adjustable parameter. Note that the experimental data-points accumulate close to the stable part of the numerical curve.

  argument as that given at the beginning of §5.3 then shows that each of the complex amplitudes ζ i 's can be assumed to be real without loss of generality, (ζ 1 , ζ 2 , ζ 3 ) ∈ R 3 . Equation (48) admits the following particular solutions: the unbuckled pattern (ζ 1 = ζ 2 = ζ 3 = 0), the hexagons (ζ 1 = ζ 2 = ζ 3 ), and the stripes (ζ 2 = ζ 3 = 0). For the hexagonal pattern, (48) yields an amplitude equation -α-αc 2

FFigure 7 :

 7 Figure 7: Bifurcation diagram predicted by the toy model (48) for patterns comprising at most three waves. Parameters are αc = 6.223, γ = 1.1916, ω =0.2245 (this set being relevant to both neo-Hookean and Gent materials), and σ = 0.4. Stable and unstable parts of the branches are shown using solid and dashed curves, respectively. This diagram is both consistent with experiments (compare the lower-left corner of this figure with figure 2) and with the non-linear bifurcation analysis (compare the vicinity of the bifurcation point α = αc in this figure with figure 5).

  Figure B.8: Finite-element simulation of the square pattern. (a) The domain, shown in light grey, is parallelepipedic box with periodic boundary conditions with size × × 1. Its side = 2π/kc is taken to match the critical wavenumber kc = 2.12 here. (b) Unsigned buckling amplitude z rms i rescaled by the initial thickness h as a function of the load parameter α. The dashed curve near the stripe branch shows the predictions of the weakly non-linear post-buckling analyses in equation (46). Inset: numerical snapshots, turned upside down to aid visualization (these images have been obtained by replicating the simulation cell multiple times). The colors encode the magnitude of the vertical displacement, as indicated in the scale bar. (c) Energy of the patterns per unit area, with the zero of energy chosen as the unbuckled configuration.

  Figure B.9: Comparison of the energy per unit area of the hexagonal and square patterns, based on the finite-element simulations. For both patterns, the wavenumber of the Fourier components have been constrained to k = kc. The hexagonal pattern of finite amplitude appears to be the one having the lowest energy in the entire range of values of α considered here.

Appendix A. A compendium on the theory of elastic bifurcations and stability Here, we recall the classical method that derives the branches near a bifurcation point from a weakly non-linear Koiter expansion (non-linear bifurcation analysis) and analyze their stability (stability analysis). We present the case relevant to the Rayleigh-Taylor instability, namely that of a non-symmetric bifurcation with multiple bifurcation modes. We refer to the reader to the book of [START_REF] Iooss | Elementary stability and bifurcation theory[END_REF] for a general introduction to bifurcation and stability theory, and to the lectures notes of [START_REF] Triantafyllidis | Stability of solids: from structures to materials[END_REF] for an exposition in the context of continuum mechanics.

Appendix A.1. Bifurcation equation

We consider the general case of a conservative system whose total potential energy E(α, t) is a function of a load parameter α and of the current configuration t. In weak form, non-linear equilibria t are found by solving (∀ t) DE(α, t)

see equation ( 26), where t is referred to as a virtual motion. Let t = t 0 (α) denotes a known branch of solutions, called the trivial branch later: (α, t 0 (α)) is by assumption a solution of the non-linear equilibrium above for any value of the load α.

Consider a critical value of the load α = α c , and denote by m the number of independent linearized equilibrium about the trivial branch ( T 26), where m is called the multiplicity.

Near the critical point (α c , t 0 (α c )), one seeks an expansion of a bifurcated branch in terms of an arclength parameter η as

These expansions are inserted into the non-linear equilibrium (A.1), and solved order by order. At linear order in η, we find that t 1 is a combination of the marginal modes, 23), where (Ξ 1 , . . . , Ξ m ) are real amplitudes which will be determined later.

For the rest of the analysis, it is useful to introduce the tensor representation of the energy gradients in the space of marginal modes,

Note that these operators are symmetric with respect to all their indices, by construction. At quadratic order η 2 , equation (A.1) yields a solvability condition, called the bifurcation equation, which reads (∀i,

This equation appeared in an equivalent form (32) in our analysis of the hexagonal patterns: the amplitudes of the linear modes were then searched in the form of complex amplitudes, which means that the present Ξ i 's can be identified with the real and imaginary parts of the former ξ j 's, i.e. m = 2n and

The bifurcation equation (A.3) can be viewed as a stationarity condition

with respect to the unknown amplitudes Ξ = (Ξ 1 , . . . , Ξ m ), of a reduced energy defined by

Each solution Ξ of (A.4) corresponds to a different bifurcated branch.

Appendix A.2. Stability of the branches Next, we present the classical method for addressing the stability of each of the bifurcated branches. The problem is to calculate the eigenpairs (σ(α), w(α)) of the rigidity operator D 2 E(α, t(α)) evaluated on each of the branches, in perturbation with respect to the parameter η: the existence of negative eigenvalues σ(α) implies that the corresponding part of the branch is unstable. Mathematically, this is done by analyzing the signs of the eigenvalues of the rigidity operator D 2 E(α, t(α)), after inserting the non-linear expansion t(α) = t 0 (α) + ηt 1 + η 2 t 2 + • • • relevant to any particular equilibrium branch.

Each of the eigenpairs satisfies the eigenvalue problem

where •|• is a bilinear dot product. We focus on the eigenvalues that cancel at the bifurcation point, as they are the only ones that can change sign upon crossing the bifurcation point. We use the following expansion for the corresponding eigenpair,

and insert it into the eigenvalue problem just stated. At order η 0 , we find

which shows that w 0 is a linear combination of the marginal modes,

with coefficients χ = (χ 1 , . . . , χ m ) to be determined next.

At order η 1 , the eigenvalue problem yields

This equations yield both the eigenvalue along the bifurcated branch, σ(η) ≈ ησ 1 , and the unknown components χ j of the eigenvector w 0 on the marginal modes T i . In terms of the reduced energy F introduced in (A.5), this equation can be rewritten as

In this equation, Ξ has to be the solution of the bifurcation equation (A.4) corresponding to the branch of interest. The fundamental branch, for instance, is analyzed by setting Ξ = 0. The stability of any branch is then concluded as follows. If the quadratic form (α 1

• positive definite, then all the rescaled eigenvalues are such that σ1 α1 > 0, implying that the original eigenvalues σ ησ 1 have the same sign as ηα 1 α -α c : the corresponding branch is unstable below the bifurcation point (α < α c ) and stable beyond the bifurcation point (α > α c );

• negative definite, then all the rescaled eigenvalues are such that σ1 α1 < 0, implying that the original eigenvalues σ ησ 1 have the opposite sign to ηα 1 α -α c : the corresponding branch is stable below the bifurcation point (α < α c ) and unstable beyond the bifurcation point (α > α c );

• has both strictly positive and negative eigenvalues σ1 α1 , then the corresponding branch is unstable on both sides of the bifurcation point.

Carrying this entire procedure is actually quite simple: once the reduced energy F(α 1 , Ξ) has been identified, the branches are found by canceling the gradient of F(α 1 , Ξ) with respect to the amplitudes Ξ; next, their stability is found by calculating the eigenvalues σ 1 of the second gradient of F(α 1 , Ξ) on these different branches. Therefore, the branches and their stability are entirely governed by the reduced energy F(α 1 , Ξ) in the vicinity of the bifurcation point.

Appendix B. Additional finite-element results

In figure B.8, additional numerical results concerning the square and the stripe patterns are shown. The same material parameters as in §8.1 are used, but the simulation domain is now a prism with a square base, as sketched in part (a) of the figure. Much like the hexagonal branch, the square branch features a fold F , implying that a pattern made of squares of finite amplitude exist at the critical load α c .

In figure B.9, we compare the energies of the square and hexagonal patterns. It appears that the branch corresponding to hexagons of finite amplitude has significantly lower energy, over the range of loads α considered here, than both the finite-amplitude and small-amplitude square patterns. This observation is consistent with the fact that squares have not been observed in the experiments.