High-Acceleration Patterns in Thin Vibrated Granular Layers

Theoretical and experimental study of high-acceleration patterns in vibrated granular layers is presented. The order parameter model based on parametric Ginzburg-Landau equation is used to describe strongly nonlinear excitations including hexagons, interface between flat anti-phase domains and new localized objects, super-oscillons. The experiments confirmed the existence of super-oscillons and bound states of super-oscillons and interfaces. On the basis of order parameter model we predict analytically and confirm experimentally that the additional subharmonic driving results in controlled motion of the interfaces.

The collective dynamics of granular materials is a subject of current interest [1][2][3][4][5].The intrinsic dissipative nature of the interactions between the constituent macroscopic particles give rise to several basic properties specific to granular substances and setting granular matter apart from the conventional gaseous, liquid, or solid states.
Driven granular systems manifest collective fluid-like behavior: convection, surface waves, and pattern formation (see e.g.[1]).One of the most fascinating examples of this collective dynamics is the appearance of long-range coherent patterns and localized excitations in vertically-vibrated thin granular layers [2][3][4][5][6][7].The particular pattern is determined by the interplay between driving frequency f and acceleration of the container Γ = 4π 2 Af 2 /g (A is the amplitude of oscillations, g is the gravity acceleration) [2,3].
Patterns appear at Γ ≈ 2.4 almost independently of the driving frequency f .At small frequencies f < f * [3,4] the transition is subcritical (hysteretic), leading to formation of squares.In the hysteretic region, localized excitations such as individual oscillons and various bound states of oscillons appear as Γ is decreased.For higher frequencies f > f * the onset pattern is stripes, and at the frequency slightly higher than f * the transition becomes supercritical.Both squares and stripes, as well as oscillons, oscillate at half of the driving frequency f /2.At higher acceleration (Γ > 4), stripes and squares become unstable, and hexagons appear instead.Further increase of acceleration at Γ ≈ 4.5 converts hexagons into a domain-like structure of flat layers oscillating with frequency f /2 with opposite phases.Depending on parameters, interfaces separating flat domains, are either smooth or "decorated" by periodic undulations.For Γ > 5.7 various quarter-harmonic patters emerge.
The pattern formation in thin layers of granular material was studied theoretically by several groups.Direct molecular dynamics simulations [8,9] (see also [10]) reproduced a majority of patterns observed in experiments and many features of the bifurcation diagram, although until now have not yielded oscillons and interfaces.Hydrodynamic and phenomenological models [11,12] reproduced certain experimental features, however neither of them offered a systematic description of the whole rich variety of the observed phenomena.In Refs.[13,14] we introduced the order parameter characterizing the complex amplitude of sub-harmonic oscillations.The equations of motion following from the symmetry arguments and mass conservation reproduced essential phenomenology of patterns near the threshold of primary bifurcation: stripes, squares, and, localized objects, oscillons.
In this paper we describe high acceleration patterns on the basis of the same order parameter model and compare it with experimental observations.Our preliminary results were published earlier in Refs.[5,14,15].Here we show that at large amplitude of driving both hexagons and interfaces emerge.We find morphological instability leading to the formation of "decorated" interfaces.We study the motion of the interface under the influence of a small subharmonic component in the driving acceleration.We also find a new localized structure, "super-oscillon", which exists for high acceleration values.We discuss possible mechanisms of saturation of the interface instability.Our experimental results demonstrate the existence of super-oscillons and bound states of super-oscillons and interfaces.They also confirm theoretical predictions of the externally controlled interface motion.
The structure of the article is following.In Sec.I we introduce our phenomenological order-parameter model.In Sec.II we develop the weakly-nonlinear theory for the flat period-doubled state of the vibrated layer.In Sec.III we analyze the interface solution and study its stability with respect to transverse perturbations.In Sec.IV we discuss new localized solutions.In Sec.V we present the combined theoretical phase diagram of the model.In Sec.VI we demonstrate that in a certain limit our model can be reduced to the real Swift-Hohenberg equation.This equation also possesses a similar variety of patterns including stripes, hexagons, stable and unstable interface solutions, as well as localized solutions.In Sec.VII we demonstrate that additional subharmonic driving results in drift of the interface with the velocity determined by the amplitude of the driving and the direction determined by the relative phase.In Sec.VIII we present our experimental results.It includes the study phase diagram and the effect of additional subharmonic driving.Section IX summarizes our conclusions.

I. PARAMETRIC GINZBURG-LANDAU EQUATION
The essence of the model [13,14] is the order parameter equation for the complex amplitude ψ of parametric layer oscillations h = ψ exp(iπf t) + c.c. at the frequency f /2 coupled with the equation for the thickness of the layer ρ averaged over the period of vibrations, Here γ is the normalized amplitude of forcing at the resonant frequency f .Linear terms in Eq. ( 1) can be obtained from the complex growth rate for infinitesimal periodic layer perturbations h ∼ exp[Λ(k)t + ikx].Expanding Λ(k) for small k, and keeping only two leading terms in the expansion Λ(k) = −Λ 0 −Λ 1 k 2 we obtain the linear terms in Eq. ( 1), where b = ImΛ 1 /ReΛ 1 characterizes ratio of dispersion to diffusion and parameter ω = (Ω 0 −πf )/ReΛ 0 , Ω 0 = −ImΛ 0 characterizes the frequency of the driving.In Eq. ( 2) ζ, β are transport coefficients for ρ.Slowly-varying thickness of the layer ρ controls the dissipation rate (the last term in Eq.( 1)).The second equation ( 2) describes re-distribution of the averaged thickness due to the diffusive flux ∝ −∇ρ, and an additional flux ∝ −ρ∇|ψ| 2 caused by the spatially nonuniform vibrations of the granular material.This coupled model was used in Refs.[13,14] to describe the pattern selection near the threshold of the primary bifurcation.It was shown that at small ζρ 0 β −1 (which corresponds to low frequencies and thick layers), the primary bifurcation is subcritical and leads to the emergence of square patterns.For higher frequencies and/or thinner layers, transition is supercritical and leads to roll patterns.At intermediate frequencies (f ≈ f * ), the stable localized solutions of Eqs.( 1),(2) corresponding to isolated oscillons and variety of bound states were found in agreement with experiment.
In this paper we focus on high-acceleration patterns at high frequencies.In Ref. [13] it was indicated that the density transport coefficient β is proportional to the energy of plate vibration (∝ f 2 ), and, therefore, it should increase with the driving frequency f .As a result, for high frequencies the coupling between ρ and ψ becomes less relevant, and one can assume ρ =const and exclude it from Eq.(1) by rescaling.Then the model can be reduced to a single order-parameter equation Eq.( 3) also describes the evolution of the order parameter for the parametric instability in vertically oscillating fluid layers (see [17,18]).
It is convenient to shift the phase of the complex order parameter via ψ = ψ exp(iφ) with sin 2φ = ω/γ.The equations for real and imaginary part ψ = A + iB are: where s 2 = γ 2 − ω 2 .

II. STABILITY OF UNIFORM STATES
At s < 1, Eqs. (4),(5) has only one trivial uniform state A = 0, B = 0, At s > 1, two new uniform states appear, A = ± √ s − 1, B = 0.The onset of these states corresponds to the period doubling of the layer flights sequence, observed in experiments [2] and predicted by the simple inelastic ball model [2,19].Signs ± reflect two relative phases of layer flights with respect to container vibrations [16].
The trivial flat state A = B = 0 loses stability if the following condition fulfilled (compare [13]): Because of the symmetry of Eq.( 3), small perturbations from the trivial state lead to the formation of a periodic sequence of rolls or stripes.
Let us analyze the stability of the non-trivial states A = ± √ s − 1, B = 0 with respect to small perturbations with wavenumber k At the threshold the uniform state loses its stability with respect to periodic modulations with the critical wavenumber k c at s < s c (correspondingly, γ < γ c ), where Small perturbations with every direction of the wavevector grow at the same rate.The resultant selected pattern is determined by the nonlinear competition between the modes.In the presence of the reflection symmetry ψ → −ψ, quadratic nonlinearity is absent, and cubic nonlinearity near the trivial state favors stripes corresponding to a single mode.Near the fixed points A = ± √ s − 1, B = 0 the reflection symmetry for perturbations U → −U, V → −V is broken, and hexagons emerge at the threshold of the instability.To clarify this point we perform weakly-nonlinear analysis of Eqs. ( 4),( 5) for s = s c − ǫ, and ǫ ≪ 1.At ǫ → 0, the variables U and V are related as in linear system: where |k| = k c and η = 2(s c − 1) + k 2 c /(bk 2 c − 2ω).Corresponding adjoint eigenvector is where η + = − 2(s c − 1) + k 2 c /bk 2 c .Substituting Eq. (10) into Eqs.(4), (5) and performing the orthogonalization, we obtain equations for the slowly-varying complex amplitudes A j , j = 1, 2, 3 (we assume only three waves with triangular symmetry, favored by quadratic nonlinearity) where the coefficients a 2 , a 3 are Eqs. ( 12) are well-studied (see [20]).There are three critical values of ǫ: ǫ A = −a 2 2 /40a 3 , ǫ R = a 2 2 /2a 3 , and ǫ B = 2a 2  2 /a 3 .The hexagons are stable for ǫ A < ǫ < ǫ B , and the stripes are stable for ǫ > ǫ R .Thus, near s = s c the model exhibit stable hexagons [21].Since we have two symmetric fixed points, both up-and down-hexagons co-exist.For smaller s stripes are stable, and for larger s, flat layers are stable, in agreement with observations [3,15], as well as direct numerical simulations of Eq.( 3) [22].The above analysis requires the values of ǫ A,B,R to be small.For parameters ω, b = O(1), this requirement is satisfied for ǫ A , but not for ǫ B,R .The estimates can be improved by substituting s = s c + ǫ instead of s = s c in (13).The resulting range of stable hexagons is plotted in the phase diagram (ω, γ) (see below Fig. 7).

III. INTERFACE SOLUTION
At s > 1, Eqs. ( 4),( 5) have an interface solution connecting two uniform states A = ± √ s − 1, B = 0. Let x be the coordinate perpendicular to the interface and y is parallel to the interface.The interface is the stationary solution to Eqs. ( 4),( 5): For b = 0 Eqs.( 14) have a solution of the form

Let us now consider the perturbed solution
where k is the transverse modulation wavenumber and λ(k) is the corresponding growth rate.For Ã, B we obtain linear equations where the matrix L is of the form In order to determine the spectrum of eigenvalues λ(k) we have solved Eq. ( 16) along with stationary Eqs. ( 14) numerically using numerical matching-shooting technique.We have found that for small ω the interface is stable with respect to transverse undulations, i.e. λ(k) ≤ 0. However, for the values of ω above certain critical value ω c (γ, b) the interface exhibits transverse instability: λ(k) > 0 in the band of wavenumbers |k| < k c , see Fig.This instability is confirmed by direct numerical simulations of Eq.( 3).An example of the evolution of slightly perturbed interface is shown in Fig. 4. Small perturbations grow to form a "decorated" interface.With time these decorations evolve slowly and eventually form a labyrinthine pattern.
The neutral curve for this instability can be determined as follows.Numerical analysis shows that at the threshold the most unstable wavenumber is k = 0 and we can expect that for k → 0 λ ∼ k 2 , see Fig. 2. Expanding Eqs.( 16) in power series of k 2 : in the zeroth order in k we obtain L(A (0) , B (0 ) = 0.The corresponding solution is the translation mode In the first order in k 2 we arrive at the linear inhomogeneous problem A bounded solution to Eq. ( 19) exists if the r.h.s. is orthogonal to the localized mode of the adjoint operator A + , B + .The adjoint operator L+ is of the form Since the operator L+ is not self-adjoint, the adjoint mode does not coincide with the translation mode, and, therefore, must be obtained numerically.
The orthogonality condition fixes the relation between λ and k: where σ is "surface tension" of the interface.The instability, corresponding to the negative surface tension of the interface, onsets at σ = 0.Although the interface solution A(x) is not localized (see Fig. 1), the integrals in Eq. ( 21) converge because the adjoint functions A + , B + are localized.The numerically-obtained adjoint eigen-mode is shown in Fig. 3.The neutral curve is shown in Fig. 7. Direct numerical simulations of Eq. (3) show that this instability leads to so called "decorated" interfaces (see Fig. 4,a), however At later stages the undulations grow and finally form labyrinthine patterns (Fig. 4, b-d).The emergence of the labyrinthine patterns as a result of interface instability contradicts experiments [4,5] where stable decorated interfaces are typically observed.Thus, the model Eq. ( 3) does not capture saturation of the interface instability and requires modification.We also checked that the coupling to the density field introduced in Ref. [13] also does not provide desired saturation.
Let us now discuss possible mechanisms of saturation of the transversal instability of the interface, which is not captured by Eq. ( 3).The reason for proliferation of the labyrinths is that local stabilization mechanisms cannot saturate the "negative surface tension instability" of the interface.Indeed, the behavior of small perturbation ζ of almost flat interface is described by the following linear equation: ζ t = −σ∆ζ, where σ is the surface tension.The linear growth must be counter-balanced by nonlinear terms.Due to translation invariance the local nonlinearity can be function of gradients of ζ only.In the lowest order the first term is |∇ζ| 3 , since the quadratic term does not saturate monotonic instability.Since for the modulation wavenumber k → 0 one has |∇ζ| 3 ∼ k 2 |ζ| 3 and ∆ζ ∼ k 2 ζ, the linear instability always overcomes local nonlinearity at small wavenumbers.As a result small perturbations grow and saturation happens due to non-local interactions as in a labyrinthine pattern which is stabilized by non-local repulsion of the fronts (see [23]).
As we see, in our model the interface instability does not result in stable decorated interfaces.One of the possibility for the saturation of this instability could be local convection induced by the interface.Indeed, due to shear motion of grains at the interface (flat regions moves in anti-phase) local convective flow can be created.The scale of this flow can be large then the interface thickness.One of the possible manifestation of the convection beside grain transport is reduction of the effective "viscosity" of flat layer, or increase in the sensitivity to external vibrations [24].As we can see from Fig. 7 increase in vibration magnitude γ suppresses the interface instability.As the instability develops, the length of the interface grows, which increases the convection and finally suppresses the instability.This feedback mechanism results in saturation of the interface instability and creates stable decorations.
Although the details of this coupling are not known yet, it can be implemented within the order-parameter approach in a phenomenological way.In Eq. 3 the forcing γψ * must be replaced by (γ 0 + w)ψ * , where w includes the effect of local convective flow and γ 0 = const.We take the simplest possible form of the coupling to the convective flow: where r 0 characterizes the typical scale of the convective flow and ǫ is the amplitude of the coupling, and S is the area of integration domain.Since |∇ψ(r ′ )| is nonzero only at the interface, this integral is proportional to the total length of the interface.With the growth of the interface length the value of "effective forcing" γ = γ 0 + w grows as well, and eventually leaves the instability domain.
We performed numerical simulations with Eq. ( 3) modified in this fashion.The convolution integral ( 22) was calculated using the fast Fourier transform.The results are presented in Fig. 5. Remarkably, for certain range of parameters the instability indeed saturates due to nonlocal term Eq. ( 22), and we find interfaces with stable decorations, see Fig. 5,a.

IV. LOCALIZED SOLUTIONS
In addition to the interface solution which exists for s > 1, Eq. ( 3) for even higher values of s possesses a variety of other nontrivial two-dimensional solutions including localized "super-oscillons".The typical solution and its corresponding localized linear mode are shown in Fig. 6.Although super-oscillon appears similar to the "conventional oscillon" discovered in Ref. [4], there is a significant difference.Asymptotic behavior of the super-oscillon is Ψ → A 0 = 0 for r → ∞ whereas the regular oscillon has Ψ → 0. As a result super-oscillons always have the same phase as their background.Oppositely-phased superoscillons must be separated by the interface connecting two opposite phases of the period-doubled flat layer whereas oppositely-phased oscillons can form stable bound states (dipoles, [4]).
Super-oscillons exist in a relatively narrow domain of parameters of Eq. ( 3).Decrease of γ leads to destabilization of the localized eigen-mode leading to subsequent nucleation of new "super-oscillons" on the periphery of the unstable super-oscillon, with the increase of γ super-oscillons disappear in a saddle-node bifurcation.
Interestingly, for certain values of γ 0 we also find an interesting structure consisting of a decorated interface and a chain of super-oscillons bound to the decorations, see Fig. 5,b.These objects can exist completely independently of the interface, Fig. 5 c, as one may expect for Eq. ( 3).Thus, we conclude that extended Eq. (3) qualitatively describes observed phenomenology.

V. PHASE DIAGRAM
The results of linear stability analysis can be summarized in the phase diagram shown in Fig. 7. Remarkably, the structure of the phase diagram Fig. 7 is qualitatively similar to that of the experiments for high frequencies of vibration, see Fig. 10 and Refs.[2][3][4][5]9,15].Increasing vibration amplitude γ leads to transition from a trivial state to stripes, hexagons, decorated interfaces, and finally, to stable interfaces.Transition from unstable to stable interfaces also occurs with decreasing ω (increasing vibration frequency f ), in agreement with Ref. [15,25].In experiments, at yet higher Γ, quarter-harmonic patterns appear, however these patterns are not described by our model.Our analysis also predicts co-existence of stripes and hexagons in a very narrow parameter range (see Sec. II).Super-oscillons were found in a narrow domain close to line 7 of Fig. 7.
Note that square patterns which are emerge at lower frequencies, are not included in the phase diagram.They are obtained within the full model Eqs.( 1),( 2) with additional density field ρ, see [13].

VI. SWIFT-HOHENBERG EQUATION
Near the line s = 1 (see Fig. 7) Eq. ( 3) can be simplified.In the vicinity of this line A ∼ (s − 1) 1/2 and B ∼ (s − 1) 3/2 ≪ A. In the leading order we can obtain from Eq. ( 5) B = b∇ 2 A/2, and Eq. ( 4) yields [26]: Rescaling the variables t → (s − 1)t, A → (s − 1) −1/2 A, x → (2(s − 1)/b 2 ) 1/4 , we arrive to the Swift-Hohenberg equation (SHE) where The description by SHE is valid if δ ∼ O(1) which implies additional condition for the validity of this approach ωb − 1 ≪ 1 at s → 1.Although this equation is simpler than Eq. ( 3), it captures many essential features of the original system dynamics, including existence and stability of stripes and hexagons in different parameter regions (see [27,28]), existence of the interface solutions, interface instability, super-oscillons and emergence of labyrinthine patterns [29].Indeed, simple analysis shows that the growth rate of the instability of the uniform state A = 1 as a function of the perturbation wavenumber is determined by the formula λ As in the original model, near the threshold of this instability, subcritical hexagonal patterns are preferred.Interface stability can also be analyzed more simply as the linearized operator corresponding to the model ( 24) is self-adjoint.The threshold value of δ is obtained from the following solvability condition which yields δ th = 1.011.It yields the equation for the neutral curve of the interface instability in the Swift-Hohenberg limit in the form: As it was shown recently, Swift-Hohenberg equation also possesses localized solutions (see [29,30]).These localized states are analogous to the super-oscillons of Eq. (3).

VII. SUBHARMONIC FORCING AND INTERFACE MOTION
Let us focus on the effect of additional subharmonic driving on the motion of the interface.In the original model Eq.(3) the interface does not move due to symmetry between left and right halves of the interface.However, if the plate oscillates with two frequencies, f and f /2, the symmetry between these two states connected by the interface, is broken, and interface moves.The velocity of interface motion depends on the relative phase of the sub-harmonic forcing with respect to the forcing at f .The effect of a small external subharmonic driving applied on the back-ground of primary harmonic driving with the frequency f can be described by the following equation (cf.Eq. ( 3)) where q characterizes the amplitude of the sub-harmonic pumping, and Φ determines its relative phase.
For small q, we look for moving interface solution in the form ψ = ψ 0 (x − vt) + qψ 1 (x − vt) + ... and v = O(q), or, alternatively, Solvability condition (compare with Eq. ( 21)) fixes the interface velocity as a function of amplitude and the phase of external subharmonic driving where α = const has the meaning of susceptibility of the interface for external forcing and Φ 0 = const is some phase shift.The explicit answer is possible to obtain for b = 0 when A + = ∂ x A 0 , and Φ = 0, π, which yields the interface velocity v = ∓ 3 2 qA −2 0 = ∓ 3 2 q(s−1) −1 .In general, A, B, A + , B + , and hence v, can be found numerically.The interface velocity as function of q, Φ is shown in Fig. 9. 5 vs q at Φ = 0. Inset: v vs Φ at q = 0.01.(•) -numerical results, (-)analytical expression (30).
Thus, from this analysis we conclude that additional subharmonic driving results in controlled motion of the interface.The velocity of the interface depends on the amplitude of the driving, and the direction is determined by the relative phase Φ.

VIII. EXPERIMENTAL RESULTS
We performed experiments with a thin layer of granular material subjected to periodic driving.Our experimental setup is similar to that of [5].We used 0.15mm diameter bronze or copper balls, and the layer thickness in our experiments was 10 particles.We performed experiments in a rectangular cell of 4×12 cm 2 .We varied the acceleration Γ and the frequency f of the primary driving signal as well as acceleration γ, frequency f 1 and the phase Φ of the secondary (additional) driving signal.The interface position and vertical accelerations were acquired using a highspeed video camera and accelerometers and further analyzed on a Pentium computer.To maintain and measure the proper acceleration at f and f /2 we employ the lock-in technique with our signal originating from an accelerometer attached to the bottom of the cell.This allows for the simultaneous real-time feed back and control of the device.To reduce the interstitial gas effects [31] we reduced the pressure to 2 mTorr.Our visual data is acquired using a high speed digital camera (Kodak SR-1000c), in addition to high speed recording this also allows for a synchronization between the patterns and the image acquisition.

A. Phase Diagram
The experimental phase diagram is shown in Fig. 10.It is similar to that of Ref. [3,4], however the transition from hexagons to interfaces is elaborated in more detail.The dashed line indicates the stability line for interfaces with respect to periodic undulations.To determine this stability limit, we kept the value of acceleration Γ fixed and increased the driving frequency f .For each value of the frequency we extracted the interface width by averaging up to 10 images of the interface.In order to find that amplitude of the periodic undulations l we subtracted from the obtained value the thickness of the flat interface l 0 .The resultant amplitude or undulations l as function of driving frequency f is shown in Fig. 11.As one sees from the Figure, the amplitude of undulations l decreases gradually approaching the critical value of the frequency.However, a small hysteresis at the transition point cannot be excluded because of large error bars.As it follows from the phase diagram, for small amplitude of the vertical acceleration Γ < 3.6 stripes are the only stable pattern (note that we focus on high-frequency patterns when squares do not exist, compare [3]).At slightly higher Γ = 3.6..3.7 stripes and hexagons coexist.For higher Γ = 3.7..4.0 the hexagons becomes stable.Due to the sub-harmonic character of motion both up and down hexagons may co-exist separated by a line phase defect, Fig. 12,a.There exists a narrow band, from Γ = 3.85 to Γ = 4.0, where localized states, or super-oscillons, appear both within the bulk of the material and also pinned to the the front between the interface and the bulk Fig. 12,b.For even higher values of the acceleration Γ > 4.0 the super-oscillons disappear (Fig. 12,c) leading to isolated interfaces with periodic decorations.

B. Controlled Motion of Interface
In the absence of additional driving the interface drifts toward the middle of the cell (see Fig. 12).We attribute this effect to the feedback between the oscillating granular layer and the plate vibrations due to the finite ratio of the mass of the granular material to the mass of the vibrating plate.Even in the absence of subharmonic drive, the vibrating cell can acquire subharmonic motion from the periodic impacts of the granular layer on the bottom plate at half the driving frequency.If the interface is located in the middle of the cell, the masses of material on both sides of the interface are equal, and due to the anti-phase character of the layer motion on both sides, an additional subharmonic driving is not produced.The displacement of the interface X from the center of the cell leads to a mass difference ∆m on opposite sides of the interface which in turn causes an additional driving proportional to ∆m.In a rectangular cell ∆m ∝ X.Our experiments show that the interface moves in such a way to decrease the subharmonic response, and the feedback provides an additional term −X/τ in the r.h.s. of Eq. ( 30), yielding The relaxation time constant τ depends on the mass ratio (this also holds for the circular cell if X is small compared to the cell radius).Thus, in the absence of an additional subharmonic drive (q = 0), the interface will eventually divide the cell into two regions of equal area (see Fig. 12).In order to verify the prediction of Eq. ( 31) we performed the following experiment.We positioned the interface off center by applying an additional subharmonic drive.Then we turned off the subharmonic drive and immediately measured the amplitude µ of the plate acceleration at the subharmonic frequency [32].The results are presented in Fig. 13.The subharmonic acceleration of the cell decreases exponentially as the interface propagates to the center of the cell.The measured relaxation time τ of the subharmonic acceleration increases with the mass ratio of the granular layer and of the cell with all other parameters fixed.The mass of the granular layer was varied by using two different cell sizes: circular, diameter 15.3 cm, and rectangular, 4 × 12 cm, while keeping the thickness of the layer unchanged.For these cells we found that the relaxation time τ in the rectangular cell is about four times greater than for the circular cell (see Fig. 13).This is consistent with the ratio of the total masses of granular material (52 grams in rectangular cell and 198 grams in circular one).In a separate experiment the mass of the cell was changed by attaching an additional weight of 250 grams to the moving shaft, which weights 2300 grams.This led to an increase of the corresponding relaxation times of 15-25 %.The relaxation time τ increases rapidly with Γ (see Fig. 13   When an additional subharmonic driving is applied, the interface is displaced from the middle of the cell.For small amplitude of the subharmonic driving γ, the stationary interface position is X = qατ sin(Φ − Φ 0 ), since the restoring force balances the external driving force.The equilibrium position X as function of Φ is shown in Fig. 14a.The solid line depicts the sinusoidal fit predicted by the theory.Because of the finite mass ratio effect, the amplitude of the measured plate acceleration µ at frequency f 1 also shows periodic behavior with Φ, (see Fig. 14b), enabling us to infer the interface displacement from the acceleration measurements.For even larger amplitude of subharmonic driving (more than 4-5 % of the primary driving) extended patterns (hexagons) re-emerge throughout the cell.
The velocity V 0 = qα which the interface would have in an infinite system, can be found from the measurements of the relaxation time τ and maximum displacement X m at a given amplitude of subharmonic acceleration γ, V 0 = X m /τ , see Eq. (31).We verified that in the rectangular cell the displacement X depends linearly on γ almost up to values at which the interface disappears at the short side wall of the cell (see Fig. 15, inset).Figure 15 shows the susceptibility α as a function of the amplitude of the primary acceleration Γ.The susceptibility decreases with Γ.The cusp-like features in the Γ-dependence of α (and τ , see inset to Fig. 13), are presumably related to the commensurability between the lateral size of the cell and the wavelength of the interface decorations.time τ and the "asymptotic" velocity V 0 .This was achieved by a small detuning ∆f of the additional frequency f 1 from the exact subharmonic frequency f /2, i.e. ∆f = f 1 − f /2 ≪ f .It is equivalent to the linear increase of phase shift Φ with the rate 2π∆f .This linear growth of the phase results in a periodic motion of the interface with frequency ∆f and amplitude X m = V 0 / τ −2 + (2π∆f ) 2 (see Eq (31)).The measurements of the "response functions" X m (∆f ) are presented in Fig. 16.From the dependence of X m on ∆f we can extract parameters V 0 , α and τ by a fit to the theoretical function.The measurements are in a very good agreement with previous independent measurements of relaxation time τ and susceptibility α.For comparison with the previous results, we indicate the values for τ and α, obtained from the response function measurements of Figs. 13 and 15 (stars).The measurements agree within 5 %.

IX. CONCLUSIONS
In this paper we studied the dynamics of thin vibrated granular layers in the framework of phenomenological model based on the Ginzburg-Landau type equation for the order parameter characterizing the amplitude of sub-harmonic sand vibrations.Although rigorous derivation of this model is not possible up to date, our approach gives a significant insight in the problem in hand and results in testable predictions.Our model is complementary to more elaborate fluid dynamics-like continuum descriptions or direct molecular dynamics simulations and allows us to carry out analytical calculations of various regimes and predict their stability.A challenging problem is to establish a quantitative relation between our order parameter model and other models and experiments.
We have shown in this paper that the generic parametric Ginzburg-Landau model Eq.( 3) captures not only patterns of vibrated sand near the primary bifurcation, but also large acceleration patterns such as hexagons, interfaces, and super-oscillons.The structure of the phase diagram Fig. 7 for high frequencies and amplitudes of vibrations is qualitatively similar to that of previous experiments Refs.[2,3,5,15], as well as the experiments reported here.Increasing vibration amplitude leads to transition from a trivial state to stripes, hexagons, decorated interfaces, and finally, to stable interfaces.Interface decorations were found to be a result of the transversal instability of a flat interface which is analogous to the negative surface tension.We found that Eq. ( 3) is not sufficient to describe the saturation of interface instability.We propose a possible non-local mechanism of saturation of this instability which takes into account the dependence of the overall length of the interface and magnitude of parametric forcing, which may occur due to a large scale convective flow induced by the sand motion near the interface.We described analytically the motion of the interface under the symmetry-breaking influence of the small subharmonic driving.
In our experimental study we found that on a qualitative level the theoretical phase diagram is similar to the experimental one.The experiments also confirmed existence of super-oscillons and their bound states with the interface as it was predicted by the theory.Further experimental study is necessary to elucidate the specific mechanism for saturation of the interface instability.
Stimulated by the theoretical prediction, we also performed experimental studies of interface motion under additional subharmonic driving.The experiment confirmed that the direction and magnitude of the interface displacement depend sensitively on the amplitude and the relative phase of the subharmonic driving.Moreover, we found that period-doubling motion of the flat layers produces subharmonic driving because of the finite ratio of the mass of the granular layer and the cell.This in turn leads to the restoring force driving the interface towards the middle of the cell.

Fig. 8 ,FIG. 8 .
Fig. 8, a-d shows the development of the interface instability within SHE with subsequent transition to labyrinthine patterns, similar to the dynamics of the full parametric equation (3).Again, to saturate the instability, an additional non-local mechanism is required.a

FIG. 11 .
FIG. 11.The amplitude of periodic undulations l as a function of driving frequency f .The acceleration amplitude Γ = 4.1

FIG. 12 .
FIG. 12. Representative high-acceleration patterns in the rectangular 4 × 12 cm cell, driving frequency f = 40 Hz: (a) Γ = 3.75 shows the onset of the interface with up and down hexagons, (b) Γ = 3.94 super-oscillons are present in the flat layer and are pinned to the front near the interface.(c) Γ = 4, isolated decorated interface.

FIG. 13 .
FIG.13.Amplitude of the subharmonic acceleration µ of the cell averaged over 4 measurements vs time for the interface propagating to the center of the cell for Γ = 3.97 and f = 40 Hz.Circles/squares correspond to the circular/rectangular cells, open/closed symbols correspond to light/heavy cells, respectively.Heavy cells differ from light ones by an additional weight of 250 g attached to the moving shaft.Solid lines show exponential fit µ ∼ exp(t/τ ) + const.Inset: τ vs Γ for light rectangular cell.