A revised gap-averaged Floquet analysis of Faraday waves in Hele-Shaw cells

Abstract Existing theoretical analyses of Faraday waves in Hele-Shaw cells rely on the Darcy approximation and assume a parabolic flow profile in the narrow direction. However, Darcy's model is known to be inaccurate when convective or unsteady inertial effects are important. In this work, we propose a gap-averaged Floquet theory accounting for inertial effects induced by the unsteady terms in the Navier–Stokes equations, a scenario that corresponds to a pulsatile flow where the fluid motion reduces to a two-dimensional oscillating Poiseuille flow, similarly to the Womersley flow in arteries. When gap-averaging the linearised Navier–Stokes equation, this results in a modified damping coefficient, which is a function of the ratio between the Stokes boundary layer thickness and the cell's gap, and whose complex value depends on the frequency of the wave response specific to each unstable parametric region. We first revisit the standard case of horizontally infinite rectangular Hele-Shaw cells by also accounting for a dynamic contact angle model. A comparison with existing experiments shows the predictive improvement brought by the present theory and points out how the standard gap-averaged model often underestimates the Faraday threshold. The analysis is then extended to the less conventional case of thin annuli. A series of dedicated experiments for this configuration highlights how Darcy's thin-gap approximation overlooks a frequency detuning that is essential to correctly predict the locations of the Faraday tongues in the frequency–amplitude parameter plane. These findings are well rationalised and captured by the present model.


Introduction
Recent Hele-Shaw cell experiments have enriched the knowledge of Faraday waves (Faraday 1831).Researchers have uncovered a new type of highly localized standing waves, referred to as oscillons, that are both steep and solitary-like in nature (Rajchenbach et al. 2011).These findings have spurred further experimentations with Hele-Shaw cells filled with one or more liquid layers, using a variety of fluids, ranging from silicone oil, and water-ethanol mixtures to pure ethanol (Li et al. 2018b).Through these experiments, new combined patterns produced by triadic interactions of oscillons were discovered by Li et al. (2014).Additionally, another new family of waves was observed in a cell filled solely with pure ethanol and at extremely shallow liquid depths (Li et al. 2015(Li et al. , 2016)).
All these findings contribute to the understanding of the wave behaviour in Hele-Shaw configurations and call for a reliable stability theory that can explain and predict the instability onset for the emergence of initial wave patterns.
Notwithstanding two-dimensional direct numerical simulations (Périnet et al. 2016;Ubal et al. 2003) have been able to qualitatively replicate standing wave patterns reminiscent of those observed in experiments (Li et al. 2014), these simulations overlook the impact of wall attenuation, hence resulting in a simplified model that cannot accurately predict the instability regions (Benjamin & Ursell 1954;Kumar & Tuckerman 1994) and is therefore not suitable for modelling Hele-Shaw flows.On the other hand, attempting to conduct threedimensional simulations of fluid motions in a Hele-Shaw cell poses a major challenge due to the high computational cost associated with the narrow dimension of the cell, which requires a smaller grid cell size to capture the shear dissipation accurately.Consequently, the cost of performing such simulations increases rapidly.
In order to tackle the challenges associated with resolving fluid dynamics within such systems, researchers have utilized Darcy's law as an approach to treating the confined fluid between two vertical walls.This approximation, also used in the context of porous medium, considers the fluid to be flowing through a porous medium, resulting in a steady parabolic flow in the short dimension.When gap-averaging the linearized Navier-Stokes equation, this approximation translates into a damping coefficient  that scales as 12/ 2 , with  the fluid kinematic viscosity and  the cell's gap-size, which represents the boundary layer dissipation at the lateral walls.However, Darcy's model is known to be inaccurate when convective and unsteady inertial effects are not negligible, such as in waves (Kalogirou et al. 2016).It is challenging to reintroduce convective terms consistently into the gap-averaged Hele-Shaw equations from a mathematical standpoint (Ruyer-Quil 2001;Plouraboué & Hinch 2002;Luchini & Charru 2010).
In their research, Li et al. (2019) applied the Kelvin-Helmholtz-Darcy theory proposed by Gondret & Rabaud (1997) to reintroduce advection and derive the nonlinear gap-averaged Navier-Stokes equations.These equations were then implemented in the open-source code Gerris developed by (Popinet 2003(Popinet , 2009) ) to simulate Faraday waves in a Hele-Shaw cell.Although this gap-averaged model was compared to several experiments and demonstrated fairly good agreement, it should be noted that the surface tension term remains twodimensional, as the out-of-plane interface shape is not directly taken into account.This simplified treatment neglects the contact line dynamics and may lead to miscalculations in certain situations.Advances in this direction were made by Li et al. (2018a), who found that the out-of-plane capillary forces associated with the meniscus curvature across the thin-gap direction should be retained in order to improve the description of the wave dynamics, as experimental evidence suggests.By employing a more sophisticated model, coming from molecular kinetics theory (Blake 1993;Hamraoui et al. 2000;Blake 2006) and similar to the macroscopic one introduced by Hocking (1987), to include the capillary contact line motion arising from the small scale of the gap-size between the two walls of a Hele-Shaw cell, they derived a novel dispersion relation, which indeed better predicts the observed instability onset.
However, discrepancies in the instability thresholds were still found.This mismatch was tentatively attributed to factors that are not accounted for in the gap-averaged model, such as the extra dissipation on the lateral walls in the elongated direction.Of course, a lab-scale experiment using a rectangular cell cannot entirely replace an infinite-length model, but if the container is sufficiently long, then this extra dissipation should be negligible.Other candidates were identified in the phenomenological contact line model or free surface contaminations.
If these factors can certainly be sources of discrepancies, we believe that a pure hydrody-namic effect could be at the origin of the discordance between theory and experiments in the first place.Despite the use of the Darcy approximation is well-assessed in the literature, the choice of a steady Poiseuille flow profile as an ansatz to build the gap-averaged model appears in fundamental contrast with the unsteady nature of oscillatory Hele-Shaw flows, such as Faraday waves.At low enough oscillation frequencies or for sufficiently viscous fluids, the thickness of the oscillating Stokes boundary layer becomes comparable to the cell gap: the Stokes layers over the lateral solid faces of the cell merge and eventually invade the entire fluid bulk.In such scenarios, the Poiseuille profile gives an adequate flow description, but this pre-requisite is rarely met in the above-cited experimental campaigns.It appears, thus, very natural to ask oneself whether a more appropriate description of the oscillating boundary layer impacts the prediction of stability boundaries.This study is precisely devoted to answering this question by proposing a revised gap-averaged Floquet analysis, based on the classical Womersley-like solution for the pulsating flow in a channel (Womersley 1955;San & Staples 2012).
Following the approach taken by Viola et al. (2017), we examine the impact of inertial effects on the instability threshold of Faraday waves in Hele-Shaw cells, with a focus on the unsteady term of the Navier-Stokes equations.This scenario corresponds to a pulsatile flow where the fluid's motion reduces to a two-dimensional oscillating channel flow, which seems better suited than the steady Poiseuille profile to investigate the stability properties of the system.When gap-averaging the linearized Navier-Stokes equation, this results in a modified damping coefficient becoming a function of the ratio between the Stokes boundary layer thickness and the cell's gap, and whose complex value will depend on the frequency of the wave response specific to each unstable parametric region.
First, we consider the case of horizontally infinite rectangular Hele-Shaw cells by also accounting for the same dynamic contact angle model employed by Li et al. (2019), so as to quantify the predictive improvement brought by the present theory.A vis-à-vis comparison with experiments by Li et al. (2019) points out how the standard Darcy model often underestimates the Faraday threshold, whereas the present theory can explain and close the gap with these experiments.
The analysis is then extended to the case of thin annuli.This less common configuration has been already used to investigate oscillatory phase modulation of parametrically forced surface waves (Douady et al. 1989) and drift instability of cellular patterns (Fauve et al. 1991).For our interest, an annular cell is convenient as it naturally filters out the extra dissipation that could take place on the lateral boundary layer in the elongated direction, hence allowing us to reduce the sources of extra uncontrolled dissipation and perform a cleaner comparison with experiments.Our homemade experiments for this configuration highlight how Darcy's theory overlooks a frequency detuning that is essential to correctly predict the locations of the Faraday's tongues in the frequency spectrum.These findings are well rationalized and captured by the present model.
The paper is organized as follows.In §2 we revisit the classical case of horizontally infinite rectangular Hele-Shaw cells.The present model is compared with theoretical predictions from the standard Darcy theory and with existing experiments.The case of thin annuli is then considered.The model for the latter unusual configuration is formulated in §3 and compared with homemade experiments in §4.Conclusions are outlined in §5.

Horizontally infinite Hele-Shaw cells
Let us begin by considering the case of a horizontally infinite Hele-Shaw cell of width  filled to a depth ℎ with an incompressible fluid of density , dynamic viscosity  (kinematic viscosity  = /) and liquid-air surface tension  (see also sketch in figure 1(a)).The vessel undergoes a vertical sinusoidal oscillation of amplitude  and angular frequency Ω.In a frame of reference which moves with the oscillating container, the free liquid interface is flat and stationary for small forcing amplitudes, and the oscillation is equivalent to a temporally modulated gravitational acceleration,  ( ′ ) =  − Ω 2 cos Ω ′ .The equation of motion for the fluid bulk are Linearizing about the rest state U ′ = 0 and  ′ ( ′ ,  ′ ) = − ( ′ )  ′ , the equations for the perturbation velocity, u ′ ( ′ ,  ′ ,  ′ ,  ′ ) = { ′ ,  ′ ,  ′ }  , and pressure,  ′ ( ′ ,  ′ ,  ′ ,  ′ ), fields, associated with a certain perturbation's wavelength  ∼  −1 (, wavenumber), read Assuming that  ≪ 1, then the velocity along the narrow  ′ -dimension  ′ ≪  ′ ,  ′ and, by employing the Hele-Shaw approximation as in, for instance, Viola et al. (2017), one can simplify the linearized Navier-Stokes equations as follows: Focus on Fluids articles must not exceed this page length Equations (2.3a)-(2.3b)are made dimensionless using  −1 for the directions  and , and  for .The forcing amplitude and frequency provide a scale Ω for the in-plane -velocity components, whereas the continuity equation imposes the transverse component  ′ to scale as Ω ≪ Ω, due to the strong confinement in the -direction ( ≪ 1).With these choices, dimensionless spatial scales, velocity components and pressure write: (2.4) The first two equations in (2.3b) in a non-dimensional form are

Floquet analysis of the gap-averaged equations
Given its periodic nature, the stability of the base flow, represented by a time-periodic modulation of the hydrostatic pressure, can be investigated via Floquet analysis.We therefore introduce the following Floquet ansatz (Kumar & Tuckerman 1994) where   is the real part of the non-dimensional Floquet exponent and represents the growth rate of the perturbation.We have rewritten ( + /Ω) =  to better explicit the parametric nature of the oscillation frequency of the wave response.In the following, we will focus on the condition for marginal stability (boundaries of the Faraday's tongues), which require the growth rate   = 0.In addition, values of  = 0 and Ω/2 correspond, respectively, to harmonic and sub-harmonic parametric resonances (Kumar & Tuckerman 1994).This implies that  is a parameter whose value is either , for harmonics, or  + 1/2, for subharmonics, with  an integer  = 0, 1, 2, . ... The gap-averaged velocity along the -direction satisfies a Darcy-like equation, (2.9) In order to obtain a governing equation for the pressure p , we average the continuity equation and we impose the impermeability condition for the span-wise velocity,  = 0 at  = ±1/2, Since < ũ >= i (  /  ) ∇ p , the pressure field p must obey the Laplace equation (2.11) It is now useful to expand each Fourier component p (, ) in the infinite -direction as sin  such that the -average implies, p = p () sin , (2.12a) which admits the solution form p =  1 cosh  +  2 sinh . (2.14) The presence of a solid bottom imposes that ŵ = 0 and, therefore, that  p / = 0, at a non-dimensional fluid depth  = −ℎ, hence giving (2.15) Let us now invoke the linearized kinematic boundary condition (2.16) Note that free surface elevation,  ′ ( ′ ,  ′ ,  ′ ), has been rescaled by the forcing amplitude , i.e.  ′ / = , and represents the projection of the bottom of the transverse concave meniscus on the -plane of figure 1(a).Moreover, by recalling the Floquet ansatzs (2.6a)-(2.6b)(with   = 0), here specified for the interface, we get an equation for each Fourier component , (2.17) Expanding η in the -direction as sin  and averaging in , i.e. < η >= η sin , leads to (2.18) Lastly, we consider the linearized dynamic condition (or linearized normal stress), evaluated at  ′ =  ′ and where the term associated with the curvature of the free surface appears, In (2.19),  ′ / ′ represents the first-order variation of the curvature associated with the small perturbation  ′ .Capillary force in the -direction is only important at large enough wavenumbers, although the associated term can be retained in the analysis in order to retrieve the dispersion relation for capillary-gravity waves (Li et al. 2019).On the other hand, the small gap of Hele-Shaw cells is such that surface tension effects in the narrow -direction are strongly exacerbated.In general, the curvature can be divided into two parts (Saffman & Taylor 1958;Chuoke et al. 1959): where the first term indicates the principal radius of curvature and the second term represents the out-of-plane curvature of the meniscus (see figure 1(a)).A common treatment of Hele-Shaw cells assumes the out-of-plane interface shape to be semicircular (Saffman & Taylor 1958;McLean & Saffman 1981;Park & Homsy 1984;Afkhami & Renardy 2013).
Nevertheless, laboratory observations have unveiled that liquid oscillations in Hele-Shaw cells experience an up-and-down driving force with  constantly changing (Jiang et al. 2004), hence giving rise to a dynamic contact angle.Here, as in Li et al. (2019), we use the following model (Hamraoui et al. 2000) to evaluate the cosine of the dynamic contact angle  as where  =  ′ / is the Capillary number defined using the vertical contact line velocity  ′ =  ′ / ′ .The friction coefficient , sometimes referred to as mobility parameter  (Xia & Steen 2018), can be interpreted in the framework of molecular kinetics theory (O.V. Voinov 1976;Hocking 1987;Blake 1993Blake , 2006;;Johansson & Hess 2018), but here, in the same spirit of Li et al. (2019), we simply view this coefficient as a constant phenomenological fitting parameter that defines the energy dissipation rate per unit length of the contact line.By combining equations (2.20)-(2.21)and taking their first-order curvature variation applied to the small perturbation, one can express After turning to non-dimensional quantities using the scaling in (2.4), equations (2.19) reads where the viscous stress term has been eliminated, as it is negligible compared to the others.With introduction of the Floquet ansatz (2.6b)-(2.17)and by recalling the -expansion of the interface and pressure as sin , the averaged normal stress equations become where the decomposition cos Ω ′ =  iΩ ′ +  −iΩ ′ /2 =  i +  −i /2 has also been used.Equations (2.15) and (2.18) are finally used to express the dynamic equation as a function of the non-dimensional averaged interface only, (2.25) with the auxiliary variables  = Ω 2 / and Γ =  2 /, such that (1 + Γ)  tanh  ℎ =  2 0 , the well-known dispersion relation for capillary-gravity waves (Lamb 1993).
As in the present form the interpretation of coefficient   does not appear straightforward, it is useful to define the damping coefficients where   is used to help rewriting 1 These auxiliary definitions allows one to express (2.25) as (2.27) or, equivalently, (2.28) Subscripts  and   in (2.26a) denote, respectively, the boundary layers and contact line contributions to the total damping coefficient   .
At the end of this long mathematical derivation, the main result is the modified damping coefficient   .Since the boundary layer contribution,   depends on the th Fourier component, the overall damping,   , is mode dependent and its value is different for each specific th parametric resonant tongue considered.This is in stark contrast with the standard Darcy approximation, where   is the same for each resonance and amounts to 12/ 2 .In our model, the case of  = 0 with  = 0 constitutes a peculiar case, as   =  0 = 0 and  0 → +∞.In such a situation,  0 () tends to the steady Poiseuille profile, so that we take  0 = 12.
Similarly to Kumar & Tuckerman (1994), equation (2.28) is rewritten as with The non-dimensional amplitude of the external forcing,  = Ω 2 / appears linearly, therefore (2.29) can be considered to be a generalized eigenvalue problem with eigenvalues  and eigenvectors whose components are the real and imaginary parts of η .See Kumar & Tuckerman (1994) for the structure of matrices A and B. 2/Ω, whose values are specified by the filled circles in (a) with matching colors.The Poiseuille profile is also reported for completeness.In drawing these figures we let the oscillation frequency of the wave, Ω, free to assume any value, but we recall that the parameter  can only assume discrete values, and so do  and  ().
For one frequency forcing we use a truncation number  = 10, which produces 2 ( + 1) × 2 ( + 1) = 22×22 matrices.Eigen-problem (2.31) is then solved in Matlab using the built-in function eigs, by asking for the eigenvalue (or few eigenvalues) with the smallest real part.
Figure 3 shows the results of this procedure for one of the configurations considered by Li et al. (2019) and neglecting the dissipation associated with the contact line motion, i.e.  = 0.In each panel, associated with a fixed forcing frequency, the black regions correspond to the unstable Faraday tongues computed using   = 12/ 2 as given by Darcy's approximation, whereas the red regions are the unstable tongues computed with the modified   =   / 2 .At a forcing frequency 4 Hz the first sub-harmonic tongues computed using the two models essentially overlap.Yet, successive resonances display an increasing departure from Darcy's model due to the newly introduced complex coefficient   .Particularly, the real part of   is responsible for the higher onset acceleration, while the imaginary part is expected to act as a detuning term, which shifts the resonant wavenumbers .

Asymptotic approximations
The main result of this analysis consists in the derivation of the modified damping coefficient   =  , + i , associated with each parametric resonance.Aiming at better elucidating how this modified complex damping influences the stability properties of the system, we would like to derive in this section an asymptotic approximation, valid in the limit of small forcing amplitudes, damping and detuning, of the first sub-harmonic (SH1) and harmonic (H1) Faraday tongues.
Unfortunately, the dependence of   on the parametric resonance considered and, more specifically, on the th Fourier component, does not allow one to convert (2.27), expressed in a discrete frequency domain, back into the continuous temporal domain.By keeping this in mind, we can still imagine fixing the value of   to that corresponding to the parametric resonance of interest, e.g. 0 (with  = 0 and  0 Ω = Ω/2) for SH1 or  1 (with  = 1 and  1 Ω = Ω) for H1.By considering then that for the SH1 and H1 tongues, the system responds in time as exp (iΩ/2) and exp (iΩ), respectively, we can recast, for these two whereas red regions are the unstable tongues computed with the present modified   =   / 2 .For this example, we consider ethanol 99.7% (see table 1) in a Hele-Shaw cell of gap size  = 2 mm filled to a depth ℎ = 60 mm. denotes the non-dimensional forcing acceleration,  = Ω 2 /, with dimensional forcing amplitude  and angular frequency Ω.For plotting, we define a small scale-separation parameter  =  /2 and we arbitrarily set its maximum acceptable value to 0.2.Contact line dissipation is not included, i.e.  =    = 0. SH stands for sub-harmonic, whereas  stands for harmonic.specific cases, equations (2.27) into a damped Mathieu equation (Benjamin & Ursell 1954;Kumar & Tuckerman 1994;Müller et al. 1997) (2.32) with either σ =  0 (SH1) or σ =  1 (H1) and where one can recognize that − (  Ω) 2 η ↔  2 η/ ′2 and i (  Ω) η ↔  η/ ′ .Asymptotic approximations can be then computed by expanding asymptotically the interface as η = η0 +  η1 +  2 η2 + . .., with  a small parameter ≪ 1.

First sub-harmonic tongue
As anticipated above, when looking at the first or fundamental sub-harmonic tongue (SH1), one should take σ →  0 (with  0 Ω = Ω/2), which is assumed small of order .The forcing amplitude  is assumed of order  as well.Furthermore, a small detuning ∼ , such that Ω = 2 0 + , is also considered, and, in the spirit of the multiple timescale analysis, a slow time scale  =  ′ (Nayfeh 2008) is introduced.At leading order, the solution reads η0 =  ()  i 0  ′ + .., with ..denoting the complex conjugate part.At the second order in , the imposition of a solvability condition necessary to avoid secular terms prescribes the amplitude  () =  ()  −i/2 to obey the following amplitude equation (2.33) Turning to polar coordinates, i.e.  = || iΦ , keeping in mind that  0 =  0, + i 0, and looking for stationary solutions with || ≠ 0 (we skip the straightforward mathematical steps), one ends up with the following approximation for the marginal stability boundaries associated with the first sub-harmonic Faraday tongue whose onset acceleration value, min  1  , amounts to Note that the final approximation on the right-hand-side of (2.35)only holds if  ℎ ≫ 1, so that tanh  ℎ ≈ 1 (deep water regime).Given that  0, > 12 and  0, > 0 always, the asymptotic approximation (2.35), in its range of validity, suggests that Darcy's model underestimates the sub-harmonic stability threshold.Moreover, from (2.34), the critical wavenumber , associated with min  1 , would correspond to that prescribed by the Darcy approximation but at an effective forcing frequency Ω +  0, = 2 0 instead of at Ω = 2 0 .This explains why the modified tongues appear shifted towards higher wavenumbers.These observations are well visible in figure 4.

First harmonic tongue
By analogy with §2.2.1, an analytical approximation of the first harmonic tongue (H1) can be provided.In the same spirit of Rajchenbach & Clamond (2015), we adapt the asymptotic  et al. (2019).The value of the friction parameter  for ethanol-70% is fitted from the experimental measurements reported in §4, but lies well within the range of values used by Li et al. (2019) and agrees with the linear trend displayed in figure 5 of Hamraoui et al. (2000).
scaling such that  is still of order , but  =  2 , σ =  1 ∼  2 (with  1 Ω = Ω) and Ω =  0 +  2 .Pursuing the expansion up to  2 -order, with η0 =  ()  i 0  ′ + ..and  () =  ()  −i , will provide the amplitude equation (2.36) The approximation for the marginal stability boundaries derived from (2.36) takes the form with a minimum onset acceleration, min and where, as before, the final approximation on the right-hand side is only valid in the deep water regime.Similarly to the sub-harmonic case, the critical wavenumber  corresponds to that prescribed by the Darcy approximation but at an effective forcing frequency Ω +  1, /2 =  0 instead of at Ω =  0 and the onset acceleration is larger than that predicted from the Darcy approximation (as  1, > 12).

Comparison with experiments by Li et al. (2019)
Results presented so far were produced by assuming the absence of contact line dissipation, i.e. coefficient  was set to  = 0, so that    = 0.In this section, we reintroduce such a dissipative contribution and we compare our theoretical predictions with a set of experimental measurements reported by Li et al. (2019), using the values they have proposed for M.This comparison, shown in figure 5, is outlined in terms of non-dimensional minimum onset acceleration, min  = min  1 , versus driving frequency.These authors performed experiments in two different Hele-Shaw cells of length  = 300 mm, fluid depth ℎ = 60 mm and gap-size  = 2 mm or  = 5 mm.Two fluids, whose properties are reported in table 1, were used: ethanol 99.7% and ethanol 50%.The empty squares in figure 5 are computed via Floquet stability analysis (2.31) using the Darcy approximation for   = 12/ 2 and correspond to the theoretical prediction by Li et al. (2019), while the colored triangles are computed using the corrected   =   / 2 .Although the trend is approximately the same, the Darcy approximation underestimates the onset acceleration with respect to the present model, which overall compares better with the experimental measurements (blackfilled circles).Some disagreement still exists, especially at smaller cell gaps, i.e.  = 2 mm, where surface tension effects are even larger.This is likely attributable to an imperfect phenomenological contact line model (Bongarzone et al. 2021(Bongarzone et al. , 2022)), whose definition falls beyond the scope of this work.Yet, this comparison shows how the modifications introduced by the present model contribute to closing the gap between theoretical Faraday onset estimates and these experiments.

The case of thin annuli
We now consider the case of a thin annular container, whose nominal radius is  and the actual inner and outer radii are −/2 and +/2, respectively (see the sketch in figure 1(b)).In the limit of / ≪ 1, the wall curvature is negligible and the annular container can be considered a Hele-Shaw cell.The following change of variable for the radial coordinate,  ′ =  +  ′ =  (1 +  ′ /) with  ′ ∈ [−/2, /2], will be useful in the rest of the analysis.As in §2, we first linearize around the rest state.Successively, we introduce the following non-dimensional quantities, It follows that, at leading order,  =  (1 + /) ∼  −→ 1/ = 1/( (1 + /)) ∼ 1/ but /  = (/) /  ∼ / ≫ 1.With this scaling and introducing the Floquet ansatzs (2.6a)-(2.6b),one obtains the following simplified governing equations, which are fully equivalent to those for the case of conventional rectangular cells if the transformation  →  is introduced.Averaging the continuity equation with the imposition of the no-penetration condition at  = ∓1/2,  (∓1/2), eventually leads to (3.3) identically to (2.11).Expanding p in the azimuthal direction as p = p sin , with  the azimuthal wavenumber, provides and the no-penetration condition at the solid bottom located at  = −ℎ/, ŵ =   p = 0, prescribes p =  1 (cosh  + tanh ℎ/ sinh ) .
(3.5)Although so far the rectangular and the annular cases are indistinguishable from each other, here it is crucial to observe that the axisymmetric container geometry translates into a periodicity condition according to which sin (−) = sin () −→ sin  = 0, (3.6) and that always imposes the azimuthal wavenumber to be an integer.In other words, in contradistinction with the case of §2, where the absence of lateral wall ideally allows for any wavenumber , here we have  = 0, 1, 2, 3, . . .∈ N.
By repeating the calculations outlined in §2, one ends up with the very same equation (2.28) (and subsequent (2.29)-(2.31)),but where  0 obeys to the quantized dispersion relation with Γ =  2 / 2 .In this context, a representation of Faraday's tongues in the forcing frequency-amplitude plane appears most natural, as each parametric tongue will correspond to a fixed wavenumber .Consequently, instead of fixing Ω and varying the wavenumber, here we solve (2.31) by fixing  and varying Ω.

Floquet analysis and asymptotic approximation
The results from this procedure are reported in figure 6, where, as in figure 3, the black regions correspond to the unstable tongues obtained according to the standard gap-averaged Darcy model, while the red ones are computed using the corrected   =   / 2 .The Faraday threshold is represented in terms of forcing (panels (a) and (b)) and forcing amplitude (panels (c) and (d)).In figure 6(a)-(c) no contact line model is included, whereas in (b)-(d) a mobility parameter  = 0.0485 is accounted for.This specific value for  will be used in the next section when comparing the theory with dedicated experiments.The regions with the lowest thresholds in each panel are sub-harmonic tongues associated with modes from  = 1 to 14.For the case with  ≠ 0, the instability onset acceleration associated with each wavenumber  appears to follow a linear trend, as already reported in figure 5.
In general, the present model gives a higher instability threshold, consistent with the results reported in the previous section.However, the tongues are here shifted to the left.This apparent opposite correction is a natural consequence of the different representations: varying wavenumber at a fixed forcing frequency (as in figure 3) versus varying forcing frequency at a fixed wavenumber (figure 6).
The asymptotic approximation for the sub-harmonic onset acceleration, adapted to this case from (2.34) in §2.2.1, helps us indeed in rationalizing the influence of the modified complex damping coefficient.
The apparent opposite shift shown in figure 6 in comparison to that displayed in figure 3, is clarified by the asymptotic relation (3.8) and, particularly by the term Ω+ 0, 2 0 − 1 .In §2, the analysis is based on a fixed forcing frequency, while the wavenumber  and, hence, the natural frequency  0 , are let free to vary.The first sub-harmonic Faraday tongue occurs when Ω +  0, ≈ 2 0 .Since Ω is fixed and  0, > 0, Ω +  0, > Ω such that  0 and therefore  have to increase in order to satisfy the relation.On the other hand, if the wavenumber  and, hence,  0 are fixed as in this section, then 2 0 −  0, < 2 0 and the forcing frequency around which the sub-harmonic resonance is centered, decreases of a contribution  0, , which introduces a frequency detuning responsible for the negative frequency shift displayed in figure 6.

Discussion on the system's spatial quantization
A first aspect that needs to be better discussed is the frequency-dependence of the damping coefficient   associated with each Faraday's tongue.In the case of horizontally infinite cells, the most natural description for investigating the system's stability properties is in the (,  ) plane for a fixed forcing angular frequency Ω (Kumar & Tuckerman 1994).According to our model, the oscillating system's response occurring within each tongue is characterized by a Stokes boundary layer thickness   = √︁ 2/(Ω + )/.For instance, let us consider sub-harmonic resonances with  = Ω/2.As Ω is fixed (see any sub-panel of figure 3), each unstable region sees a constant   (with  = 0, 1, 2, . ..) and hence a constant damping   .
On the other hand, in the case of quantized wavenumber as for the annular cell of §3, the most suitable description is in the driving frequency-driving amplitude plane at fixed wavenumber  (see figure 6) (Batson et al. 2013).In this description, each sub-harmonic ( = Ω/2) or harmonic ( = Ω) th tongue associated with a wavenumber , sees a   , and thus a   , changing with Ω along the tongue itself.
Furthermore, it is important to realize that in a real lab-scale experiment, the horizontal size of rectangular cells is never actually infinite.It follows that if the analysis of §2 is restrained to horizontally finite cells of overall length , then one must impose the no-penetration condition for  ′ at  ′ = ±/2, which would set the admissible wavenumbers to  = / only, with  = 0, 1, 2, . . .∈ N, thus completing the analogy with the annular configuration.
In such a case however, the solution form (2.9) prevents the no-slip condition for the in-plane -velocity components to be imposed (Viola et al. 2017).This always translates into an underestimation of the overall damping of the system in standard Hele-Shaw cells, although the sidewall contribution is expected to be negligible for sufficiently long cells.
On the other hand, the case of a thin annulus, by naturally filtering out this extra dissipation owing to the periodicity condition, offers a prototype configuration that can allow one to better quantify the correction introduced by the present gap-averaged model when compared to dedicated experiments, as outlined in the next section.

Setup
The experimental apparatus, shown in figure 7, is very simple.We used a Plexiglas annular container of height 100 mm, nominal radius  = 44 mm and gap-size  = 7 mm, which is then filled to a depth ℎ = 65 mm with ethanol 70% (see table 1 for the fluid properties).An air conditioning system helps in maintaining the temperature of the room at around 22 • .The container is mounted on a loudspeaker VISATON TIW 360 8Ω placed on a flat table and connected to a wave generator TEKTRONIX AFG 1022, whose output signal is amplified using a wideband amplifier THURKBY THANDER WA301.The motion of the free surface is recorded with a digital camera NIKON D850 coupled with a 60mm f/2.8D lens and operated in slow motion mode, allowing for an acquisition frequency of 120 frames per second.A LED panel placed behind the apparatus provides back illumination of the fluid interface for better optimal contrast.The wave generator imposes a sinusoidal alternating voltage,  = (  /2) cos (Ω ′ ), with Ω the angular frequency and    the full peak-to-peak voltage.The response of the loudspeaker to this input translates into a vertical harmonic motion of the container,  cos (Ω ′ ), whose amplitude,  [mm], is measured with a chromatic confocal displacement sensor STI CCS PRIMA/CLS-MG20.This optical pen, which is placed around 2 cm (within the admissible working range of 2.5 cm) above the container and points at the top flat surface of the outer container's wall, can detect the time-varying distance between the fixed sensor and the oscillating container's surface with a sampling rate in the order of kHz (i) First, we need to ensure that the loudspeaker's output translates into a vertical container's displacement following a sinusoidal time signal.To this end, the optical sensor is used to measure the container motion at different driving frequencies.These time signals are then fitted with a sinusoidal law. Figure 8 shows how below a forcing frequency of 8 Hz, the loudspeaker's output begins to depart from a sinusoidal signal.This check imposes a first lower bound on the explorable frequency range.
(ii) In addition, as Faraday waves only appear above a threshold amplitude, it is convenient to measure a priori the maximal vertical displacement  achievable.The loudspeaker response curve is reported in the bottom part of figure 8.A superposition of this curve with the predicted Faraday's tongues immediately identifies the experimental frequency range within which the maximal achievable  is larger than the predicted Faraday threshold so that standing waves are expected to emerge in our experiments.Assuming the herein proposed gap-averaged model (red regions) to give a good prediction of the actual instability onset, the experimental range explored in the next section is limited to approximately ∈ [10.2, 15.6] Hz.

Procedure
Given the constraints discussed in §4.2, experiments have been carried out in a frequency range between 10.2 Hz and 15.6 Hz with a frequency step of 0.1 Hz.For each fixed forcing frequency, the Faraday threshold is determined as follows: the forcing amplitude  is set to the maximal value achievable by the loudspeaker, so as to quickly trigger the emergence of the unstable Faraday wave.The amplitude is then progressively decreased until the wave disappears and the surface becomes flat again.
More precisely, a first quick pass across the threshold is made to determine an estimate (b) (a) of the sought amplitude.A second pass is then made by starting again from the maximum amplitude and decreasing it.When we approach the value determined during the first pass, we perform finer amplitude decrements, and we wait several minutes between each amplitude change to ensure that the wave stably persists.We eventually identify two values: the last amplitude where the instabilities were present (see figure 9(a)) and the first one where the surface becomes flat again (see figure 9(b)).Two more runs following an identical procedure are then performed to verify the values previously found.Lastly, an average between the smallest unstable amplitude and the largest stable one gives us the desired threshold.
Once the threshold amplitude value is found for the considered frequency, the output of the wave generator is switched off, the frequency is changed, and the steps presented above are implemented again for the new frequency.In this way we always start from a stable configuration, hence limiting the possibility of nonlinear interaction between different modes.
For each forcing frequency, the two limiting amplitude values, identified as described above, are used to define the error bars reported in figure 10.Those error bars must also account for the optical pen's measurement error (0.1 m), as well as the non-uniformity of the output signal.By looking at the measured average, minimum, and maximum amplitude values in the temporal output signal, it is noteworthy that the average value typically deviates from the minimum and maximum by around 10 m.Consequently, we incorporate in the error bars this additional 10 m of uncertainty in the value of .The uncertainty in the frequency of the output signal is not included in the definition of the error bars, as it is extremely small, on the order of 0.001 Hz.

Instability onset and wave patterns
The experimentally detected threshold at each measured frequency is reported in figure 10 in terms of forcing acceleration  and amplitude .Once again, the black unstable regions are calculated according to the standard gap-averaged model with   = 12/ 2 , whereas red regions are the unstable tongues computed using the modified damping   =   / 2 .Both scenarios include contact line dissipation    = (2/) (/) tanh (ℎ/), with a value of  equal to 0.0485 for ethanol 70%.Although, at first, this value has been simply selected in order to fit well our experimental measurements, it is in perfect agreement with the linear relation linking  to the liquid's surface tension reported in figure 5 of Hamraoui et al. (2000) and used by Li et al. (2019)

(see table 1).
As figure 10 strikingly shows, the present theoretical thresholds match well our experimental measurements.On the contrary, the poor description of the oscillating boundary layer in the classical Darcy model translates into a lack of dissipation.The arbitrary choice of a higher fitting parameter  value, e.g. ≈ 0.09 would increase contact line dissipation and compensate for the underestimated Stokes boundary layer one, hence bringing these predictions much closer to experiments; however, such a value would lie well beyond the typical values reported in the literature.Furthermore, the real damping coefficient   = 12/ 2 given by the Darcy theory does not account for the frequency detuning displayed by experiments.This frequency shift is instead well captured by the imaginary part of the new damping   =   / 2 (  =  , + i , ).
Within the experimental frequency range considered, five different standing waves, corresponding to  = 5, 6, 7, 8 and 9, have emerged.The identification of the wavenumber  has been simply performed by visual inspection of the free surface patterns reported in figure 11.Indeed, by looking at two time snapshots separated by a forcing period , it is possible to count the various wave peaks along the azimuthal direction.
When looking at figure 10, it is worth commenting that on the left sides of the marginal stability boundaries associated with modes  = 5 and 6 we still have a little discrepancy between experiments and the model.Particularly, the experimental thresholds are slightly lower than the predicted ones.A possible explanation can be given by noticing that our experimental protocol is agnostic to the possibility of subcritical bifurcations and hysteresis, while such behaviour has been predicted by Douady (1990).
As a last comment, one has to keep in mind that the Hele-Shaw approximation remains good only if the wavelength, /2 does not become too small, i.e. comparable to the cell's gap, .In other words, one must check that the ratio /2 is of the order of the small  separation-of-scale parameter, .For the largest wavenumber observed in our experiments,  = 9, the ratio /2 amounts to 0.23, which is not exactly small.Yet, the Hele-Shaw approximation is seen to remain fairly good.

Contact angle variation and thin film deposition
Before concluding, it is worth commenting on why the use of dynamic contact angle model (2.21) is justifiable and seen to give good estimates of the Faraday thresholds.
Existing lab experiments have revealed that liquid oscillations in Hele-Shaw cells constantly experience an up-and-down driving force with an apparent contact angle  constantly changing (Jiang et al. 2004).Our experiments are consistent with such evidence.In figure 12 we report seven snapshots, (i)-(vii), covering one oscillation period, , for the container motion.These snapshots illustrate a zoom of the dynamic meniscus profile and show how the macroscopic contact angle changes in time during the second half of the advancing cycle (i)-(v) and the first half of the receding cycle (vi)-(x), hence highlighting the importance of the out-of-plane meniscus curvature variations.Thus, on the basis of our observations, it seems appropriate to introduce in the theory a contact angle model so as to justify this associated additional dissipation, which would be neglected by assuming  = 0.The model   used in this study, and already implemented by Li et al. (2019), is very simple; it assumes the cosine of the dynamic contact angle to linearly depend on the contact line speed through the capillary number  (Hamraoui et al. 2000).Accounting for such a model is shown, both in Li et al. (2019) and in this study, to supplement the theoretical predictions by a sufficient extra dissipation suitable to match experimental measurements.This dissipation eventually reduces to a simple damping coefficient    as it is of linear nature.A unique constant value of the mobility parameter  is sufficient to fit all our experimental measurements at once, suggesting that the meniscus dynamics is not significantly affected by the evolution of the wave in the azimuthal direction, i.e. by the wavenumber, and  can be seen as an intrinsic property of the liquid-substrate interface.
Several studies have discussed the dependence of the system's dissipation on the substrate material (Huh & Scriven 1971;Dussan 1979;Cocciaro et al. 1993;Ting & Perlin 1995;Eral et al. 2013;Viola et al. 2018;Viola & Gallaire 2018;Xia & Steen 2018).These authors, among others, have unveiled and rationalized interesting features such as solid-like friction induced by contact angle hysteresis.This strongly nonlinear contact line behaviour does not seem to be present in our experiments.This can be tentatively explained by looking at figure 13.These snapshots illustrate how the contact line constantly flows over a wetted substrate, due to the presence of a stable thin film deposited and alimented at each oscillation cycle.This feature has been also recently described by Dollet et al. (2020), who showed that the relaxation dynamics of liquid oscillation in a U-shaped tube filled with ethanol, due to the presence of a similar thin film, obey an exponential law that can be well-fitted by introducing a simple linear damping, as done in this work.

Conclusions
Previous theoretical analyses for Faraday waves in Hele-Shaw cells have so far relied on the Darcy approximation, which is based on the parabolic flow profile assumption in the narrow direction and that translates into a real-valued damping coefficient   = 12/ 2 , with  the fluid kinematic viscosity and  the cell's gap-size, that englobes the dissipation originated from the Stokes boundary layers over the two lateral walls.However, Darcy's model is known to be inaccurate whenever inertia is not negligible, e.g. in unsteady flows such as oscillating (i) (ii) (iii) Figure 13: These three snapshots correspond to snapshots (ii), (iii) and (iv) of figure 12 and show, using a different light contrast, how the contact line constantly moves over a wetted substrate due to the presence of a stable thin film deposited and alimented at each cycle.See also supplementary movies 7 at: LINK.

Figure 1 :
Figure 1: (a) Sketch of Faraday waves in a rectangular Hele-Shaw cell of width  and length  filled to a depth ℎ with a liquid.Here  denotes the gap size of the Hele-Shaw cell,  is the wavelength of a certain wave, such that / ≪ 1, and  is the dynamic contact angle of the liquid on the lateral walls.The vessel undergoes a vertical sinusoidal oscillation of amplitude  and angular frequency Ω.The free surface elevation is denoted by  ′ ( ′ ).(b) Same as (a), but in an annular Hele-Shaw cell with internal and external radii, respectively,  − /2 and  + /2.Here, / ≪ 1 and the free surface elevation is a function of the azimuthal coordinate  ′ , i.e.  ′ ( ′ ).
Figure 2: (a) Real and imaginary parts of the complex auxiliary coefficient  =   + i  versus twice the non-dimensional Stokes boundary layer thickness .The horizontal black dotted line indicates the constant value 12 given by the Darcy approximation.(b) Normalized profile  () (Womersley profile) for different  =  −1 √︁2/Ω, whose values are specified by the filled circles in (a) with matching colors.The Poiseuille profile is also reported for completeness.In drawing these figures we let the oscillation frequency of the wave, Ω, free to assume any value, but we recall that the parameter  can only assume discrete values, and so do  and  ().

Figure 3 :
Figure3: Faraday tongues computed via Floquet analysis at different fixed driving frequencies (reported on the top of each panel).Black regions correspond to the unstable Faraday tongues computed using   = 12/ 2 as in the standard Darcy approximation, whereas red regions are the unstable tongues computed with the present modified   =   / 2 .For this example, we consider ethanol 99.7% (see table 1) in a Hele-Shaw cell of gap size  = 2 mm filled to a depth ℎ = 60 mm. denotes the non-dimensional forcing acceleration,  = Ω 2 /, with dimensional forcing amplitude  and angular frequency Ω.For plotting, we define a small scale-separation parameter  =  /2 and we arbitrarily set its maximum acceptable value to 0.2.Contact line dissipation is not included, i.e.  =    = 0. SH stands for sub-harmonic, whereas  stands for harmonic.

Figure 4 :
Figure 4: First sub-harmonic and harmonic Faraday tongues at a driving frquency 1/ = 18 Hz for the same configuration of figure 3. Black and red regions show unstable tongues computed via Floquet analysis by using, respectively,   = 12/ 2 and the modified   =  1 / 2 from the present model.Dashed and solid light-blue lines correspond to the asymptotic approximations according to (2.34) and (2.37).

Figure 5 :
Figure 5: Sub-harmonic instability onset, min  , versus forcing frequency.Comparison between theoretical data (empty squares: standard Darcy model,   = 12/ 2 ; colored triangles: present model,   =   / 2 ) and experimental measurements by Li et al. (2019).The values of the mobility parameter  here employed are reported in the figure.

Figure 6 :
Figure 6: Faraday tongues computed via Floquet analysis (2.31) at different fixed azimuthal wavenumber  and varying the driving frequency.(a)-(b) Faraday thresholds in terms of forcing acceleration  = Ω 2 /; (c)-(d) Threshold in terms of forcing amplitude .Black regions correspond to the unstable Faraday tongues computed using   = 12/ 2 , whereas red regions are the unstable tongues computed with the present modified   =   / 2 .The fluid parameters used here correspond to those given in table 1 for ethanol 70%.The gap-size is set to  = 7 mm, the fluid depth to ℎ = 65 mm and the nominal radius to  = 44 mm.Contact line dissipation is included in (b) and (d) by accounting for a mobility coefficient  = 0.0485.The regions with the lowest thresholds in each panel are sub-harmonic tongues associated with modes from  = 1 to 14.

Figure 7 :
Figure 7: Photo of the experimental setup

Figure 8 :
Figure 8: Top: vertical container displacement  versus time at different forcing frequencies.The black curves are the measured signal, while the green dash-dotted curves are sinusoidal fitting.Below a forcing frequency of 8 Hz, the loudspeaker's output begins to depart from a sinusoidal signal.Bottom: same as in figure 6(d): sub-harmonic Faraday tongues computed by accounting for contact line dissipation with a mobility parameter  = 0.0485.The light blue curve here superposed corresponds to the maximal vertical displacement  achievable with our setup.With this constraint, Faraday waves are expected to be observable only in the frequency range highlighted in blue.

Figure 9 :
Figure 9: Free surface shape at a forcing frequency 1/ = 11.7 Hz and corresponding to: (a) the lowest forcing amplitude value,  = 0.4693 mm, for which the  = 6 standing wave is present (the figure shows a temporal snapshot); (b) the largest forcing amplitude value,  = 0.4158 mm, for which the surface becomes flat and stable again.Despite the small forcing amplitude variation, the change in amplitude is large enough to allow for a visual inspection of the instability threshold with sufficient accuracy.

Figure 10 :
Figure 10: Experiments (empty circles) are compared to the theoretically predicted sub-harmonic Faraday threshold computed via Floquet analysis (2.31) for different fixed azimuthal wavenumber  and according to the standard (black region) and revised (red regions) gap-averaged models.The shaded band around the instability onset indicates the error bar for the threshold amplitudes at each measured driving frequency.The tongues are computed by including contact line dissipation with a value of  equal to 0.0485 as in figures 6(b)-(d) and 8.As explained in §4.3, the vertical error bars indicate the amplitude range between the smallest measured forcing amplitude at which the instability was detected and the largest one at which the surface remains stable and flat.These two limiting values are successively corrected by accounting for the ptical pen's measurement error and the non-uniformity of the output signal of the loudspeaker.

Figure 12 :
Figure 12: Zoom of the meniscus dynamics recorded at a driving frequency 11.6 Hz and amplitude  = 1.2 mm for  = 6.Seven snapshots, (i)-(vii), covering one oscillation period, , for the container motion are illustrated.These snapshots show how the meniscus profile and the macroscopic contact angle change in time during the second half of the advancing cycle and the first half of the receding cycle, hence highlighting the importance of the out-of-plane curvature or capillary effects. .See also supplementary movie 6 at: LINK.
, is a rescaled Stokes boundary layer thickness specific to the th Fourier component.The function   () is displayed in figure 2(b), which depicts how a decrease in the value of   starting from large values corresponds to a progressive transition from a fully developed flow profile to a plug flow connected to thin boundary layers.

Table 1 :
Characteristic fluid parameters for the three ethanol-water mixtures considered in this study.Data for the pure ethanol and ethanol-water mixture (50%) are taken from Li