2.1. Experimental set-up

Figure 1(a) shows the experimental set-up we used to excite and visualize Faraday waves. A circular cylinder of radius $R=45\, \mathrm {mm}$ and height $H=5.6\, \mathrm {mm}$ was assembled from a circular plate and two ring-shaped plates (see figure 1b) manufactured from laser-cut transparent acrylic sheets. The three acrylic plates were tightly bonded together in order that the liquid depth $H$ was equal to the thickness $h_0$ of the middle plate. The upper plate of a larger inner radius prevented spillage of oscillating liquid. In each test case, a random-dot pattern printed on white paper was inserted beneath the lower plate (thickness $h_b=2.85 \, \mathrm {mm}$) to enable the displacement field to be determined using a DIC algorithm. The cylinder was mounted on an electromechanical shaker (ESS-050) capable of vibrating sinusoidally in the vertical direction for a prescribed driving frequency in the range $\varOmega _0/2{\rm \pi} =5\sim 10\,000\, \mathrm {Hz}$. The shaker was driven by an acceleration generator (Amber) and a power amplifier (PA-1200). The stability of the shaker oscillation was guaranteed by a closed-loop control system, consisting of a function generator, amplifier, shaker and acceleration sensor. A prescribed oscillating amplitude $\xi$ and driving angular frequency $\varOmega _0$ determined the time-varying acceleration $A_0 \cos {\varOmega _0 t}$ with $A_0=\xi \varOmega _0^2$. In the experiments a stable acceleration amplitude $A_0$ was guaranteed by sustained feedback between the acceleration sensor and shaker as they oscillated together.

Figure 1. Experimental set-up: (a) schematic layout; (b) photograph of liquid container composed of three acrylic plates with a random-dot pattern on paper between the base of the lowest plate and the top panel of the shaker; (c) transverse geometry of the liquid-filled cylindrical container.

The liquid surface pattern was determined from data obtained using a visualization system comprising a set of LED lamps and a SpeedCam MacroVis EoSens (Germany) camera with Tokina atx-i 100 mm F2.8 FF MACRO (Japan) lens (see figure 1a). A ring-shaped light source placed at a suitably raised elevation produced uniform white light. The high-speed camera was positioned so that it pointed vertically downward at a height $h_c$ $= 1.2\ \mathrm {m}$ above the sheet containing the random-dot pattern to capture a square area covering the overall surface domain. High-resolution images ($1328\times 1330\ {\mathrm {pixel}}^2$) were collected to ensure accurate surface reconstruction. In practice, a frame rate of 500 fps and an exposure time of 2 ms were found to guarantee proper capture of the wave patterns throughout the driving cycle. Before parametric oscillations were generated, a reference image was obtained of the flat surface of the liquid in almost still conditions except for a slight disturbance caused by low-level noise emanating from the standby mechanism of the shaker. The trigger for camera capture was delayed until the liquid surface pattern had begun to evolve steadily with periodicity. The elapsed time from shaking onset to stable pattern formation, herein called the growth time $\tau$, was different for each parametric vibration case. As with the previous regulation of cross-waves (Moisy et al. Reference Moisy, Michon, Rabaud and Sultan2012), $\tau$ was found to have a roughly inverse-proportional relationship with acceleration. A sufficiently long onset-trigger time of 4 min ($>\tau$) was set to capture stable surface patterns.

Previous research has demonstrated that simultaneous superposition of harmonic meniscus waves and subharmonic Faraday waves would occur when the container is either over- or under-filled (Shao et al. Reference Shao, Wilson, Saylor and Bostwick2021b). Hence, great care had to be taken to ensure that meniscus waves would not form. As illustrated in figure 1(c), the cylindrical container was filled with ultra-pure water to the inside brim of the middle ring plate. Subsequent addition or removal of liquid was undertaken until the liquid surface was almost flat, such that the contact line was pinned at the edge and the static contact angle was ${\rm \pi} /2$, suppressing both meniscus waves and dynamic contact-line effects. All experiments were carried out at room temperature ($25\pm 0.5 \, ^\circ \mathrm {C}$) so that the properties of pure water remained constant as follows: water density $\rho = 997\, \mathrm {kg}\,\mathrm {m}^{-3}$, dynamic viscosity $\mu =10^{-3}\,\mathrm {Pa}\,\mathrm {s}$ and surface tension $\sigma =72\, \mathrm {mN}\,\mathrm {m}^{-1}$. Hence, the Bond number $Bo=\rho g R^2/\sigma \approx$275, where $g$ is the gravitational constant $g=9.8\, \mathrm {m}\,\mathrm {s}^{-2}$. Due to the gradual evaporation of water, a meniscus initially appeared along the interior wall about 15 minutes after filling and so pure water was added until the meniscus disappeared. Additionally, to confirm the reproducibility of our experimental observations in an open environment (similar to Henderson & Miles Reference Henderson and Miles1991), at about 90 min intervals we emptied and cleaned the container, and changed the liquid to avoid the apparent drop in surface tension should surface pollution be present. Hence, it could reasonably be assumed that the surface tension $\sigma$, although not measured, was approximately the same as that of pure water at room temperature.

Stable Faraday waves exist in a limited acceleration window spanning from the critical Faraday acceleration ($A_F$) to chaotic acceleration ($A_c$). Furthermore, in order for the FS-SS method (introduced in § 2.2) to be valid, it was necessary to restrict the deformed liquid surface to have weak slopes, thus limiting the range of experimental acceleration values. This was because high-sloping waves, such as large standing gravity waves, alter the light intensity distribution through refraction causing visible light and dark streaks and modal patterns to be superimposed on the images. Moreover, the high curvature of the liquid surface can give rise to ray crossings and light reflection spots in the recorded images. Upper limits were placed on the driving parameters, amplitude $\xi$ and frequency $\varOmega _0/2{\rm \pi}$, to prevent contamination of the images by uneven intensity distribution, ray crossings and light spots.

2.2. Surface reconstruction

The FS-SS method (figure 2a) developed by Moisy et al. (Reference Moisy, Rabaud and Salsac2009) was used to determine the liquid free-surface topography. The FS-SS method relies on knowledge of the refracted light and the displacement field of the random-dot pattern to reconstruct the surface gradient field. Figure 2(b–e) illustrates the application of the FS-SS method. From refracted images of the flat and deformed interfaces shown in figure 2(b,c), the displacement vector field $\delta \boldsymbol {r}$ (figure 2d) was computed using the augmented Lagrangian digital image correlation (AL-DIC) algorithm proposed by Yang & Bhattacharya (Reference Yang and Bhattacharya2019). Interrogation windows of $16\times 16\ \text {pixel}^2$ with a subset size of $8\times 8\ \text {pixel}^2$ were used in the computation of the AL-DIC algorithm, resulting in a spatial resolution of $0.68\, \mathrm {mm}$. By multiplying the pixel scale by the conversion factor determined for each set of experimental photographs, the spatial length in pixels could convert to physical scale ($\mathrm {mm}$).

Figure 2. (a) Side elevation schematic of light refraction in a Faraday wave whereby the local free-surface slope $i$ results in a corresponding refracted pattern. The reference pattern (b) and deformed pattern (c), captured experimentally for $\varOmega _0/2{\rm \pi} =18.5\, \mathrm {Hz}$ and $\xi =0.15\, \mathrm {mm}$, are used to compute the deformation field (d), which further leads to the 2-D overview of the reconstructed free-surface elevation (e).

Due to the high elevation of the camera above the free surface and the weak surface slope encountered in the present experiments, the surface gradient $\boldsymbol {\nabla } \eta$ was linearly related to the vector field $\delta \boldsymbol {r}$ as

(2.1)\begin{equation} \boldsymbol{\nabla} \eta={-}\frac{\delta \boldsymbol{r}}{(1-n_a/n_l)h_p}, \end{equation}

where $n_a$ and $n_l$ denote the optical indices of air and liquid, respectively, and $h_p$ is the effective surface-pattern distance including overall optical correction, given by

(2.2)\begin{equation} h_p=h_0+\frac{n_l}{n_b}h_b, \end{equation}

where $n_b$ is the optical index of the bottom acrylic plane. The derivation of relation (2.1) was based on three mandatory approximations required by the FS-SS method: a paraxial approximation, a weak-slope approximation and a weak-amplitude approximation. In our experiments, the pattern-camera distance $h_c$ was substantially larger than the field size $R$, thus satisfying the paraxial approximation. The wave slope $i$ was approximately determined as the ratio of wave amplitude $|\eta |$ to wavelength $\lambda$. In the present tests the experimental parameters (including the tiny driving amplitude $\xi$ and acceleration $\epsilon$) were selected to guarantee that the weak-slope approximation was satisfied. Parametric excitation of micro-amplitude waves and a sufficiently thick base plate ensured that the weak-amplitude approximation held. Due to $h_c\gg R$ and $\xi \ll R$, the reference patterns, after marginal alteration, had minimal effect on the calculated displacement field, and so it was feasible to select the reference pattern captured immediately before vibration commenced. The modest optical depth and weak surface slope helped prevent ray crossing (Moisy et al. Reference Moisy, Rabaud and Salsac2009).

Some additional, unavoidable effects arising from the relative motion between the visual plane and random-dot pattern were discerned when the FS-SS method resolved the surface interface in the vibrating system (Damiano et al. Reference Damiano, Brun, Harris, Galeano-Rios and Bush2016). Conventionally, the background pattern has been assumed to remain stationary in the derivation of the FS-SS method and in later implementations (Moisy et al. Reference Moisy, Rabaud and Salsac2009). However, in the present experiments the shaking amplitude varied from 0.1 mm to 0.25 mm, and scaled patterns were recorded as the background pattern vibrated. Although zoom and translation effects were indistinguishable by the naked eye, the error concealed in the deformation field revealed itself on the reconstructed surface. Detailed analysis of the vibrating patterns is given in the following evaluation of a typical test case where $\xi =0.15\, \mathrm {mm}$ and ${\varOmega _0/2{\rm \pi} =6\, \mathrm{Hz}}$.

Figure 3 shows the deformation field $\delta \boldsymbol {r}$ corresponding to a specific frame, and its decomposition into primary components by removing the mean slope translation and applying 2-D Butterworth filters. Theoretically, low-acceleration vibration below the Faraday threshold should not excite either a harmonic edge wave or a subharmonic Faraday wave; hence, the ideal reconstructed surface is flat ($\delta \boldsymbol {r}=\boldsymbol {0}$). In practice, however, concave/convex patterns emerged due to spurious interference of the deformation field. Moreover, the pinned boundary (static contact angle $\alpha ={\rm \pi} /2$) condition could not be perfectly satisfied because of a combination of slow evaporation of water and machining errors of the container, which led to some uniform ripples appearing in the reconstructed surface even when the shaker was in standby mode. Consequently, the raw deformation field was divided into the following parts:

(2.3)\begin{equation} \delta \boldsymbol{r} = \delta \boldsymbol{r}^{*}-\delta \boldsymbol{r}_v-\delta\boldsymbol{r}_c-\delta \boldsymbol{r}_o. \end{equation}

Here the actual deformation field $\delta \boldsymbol {r}$ was extracted from the raw deformation data $\delta \boldsymbol {r}^*$ by implementing three filters to eliminate the vibration-induced deformation field $\delta \boldsymbol {r}_v$, the high-frequency deformation field $\delta \boldsymbol {r}_c$ produced by machine standby and the (small) low-frequency noise field $\delta \boldsymbol {r}_o$.

Figure 3. Decomposition of deformation fields in the $\hat{x}$ direction (a) and in the $\hat{y}$ direction (b) in a low-acceleration experiment and corresponding reconstructed results $\eta$ (c) obtained by applying the FS-SS method (unit: mm). Here, $\hat{x}$ and $\hat{y}$ are two unit vectors parallel to the row direction and column direction of the image pixel, respectively. After multi-step filtering, the raw deformation field $\delta \boldsymbol {r}^*$ (see (a i–b i)) is decomposed into the vibration-induced deformation field $\delta \boldsymbol {r}_v$, low-frequency deformation field $\delta \boldsymbol {r}_o$ and high-frequency deformation field $\delta \boldsymbol {r}_c$, which are shown sequentially in (a ii–b ii) to (a iv–b iv).

In our experiments, the vibration amplitude was of the order of the wave topography, so the zoom effect of visual pictures could not be ignored. Besides, the divergence between the central axes of the camera and the cylindrical vessel caused the overall movement of the images, referred to as the translation effect. The zoom effect contributed to a scaled deformation whereas the translation effect resulted in an equivalent deformation in the visualization system. Both effects were directly related to relative motion between the camera and the background. For simplicity, the superimposed deformation was expressed in linear form as

(2.4)\begin{equation} \delta \boldsymbol{r}_v = \left(c_x\, x+d_x\right)\hat{\boldsymbol{x}}+\left(c_y\, y+d_y\right)\hat{\boldsymbol{y}}, \end{equation}

where $c_x$ and $c_y$ are zoom effect factors, and $d_x$ and $d_y$ are translation effect factors. Hence, $\delta \boldsymbol {r}_v$ was subtracted by filtering out the best-fit plane (2.4) from the raw deformation field. We call this a translation/zoom filter.

Two further effects were uncovered by selectively applying low-pass and high-pass Butterworth filters. An irregular deformation field $\delta \boldsymbol {r}_o$ was linked to ambient fluctuation arising from inhomogeneity of light intensity, etc. A ripple-shaped deformation field $\delta \boldsymbol {r}_c$ originated from the shaking standby mode as confirmed by matching the electric current frequency and the measured wavelength $\lambda _c$ (see figure 3). The ripples are harmonic resonant patterns, and the corresponding wave frequency of measured wavelength $\lambda _c$ ($\approx$5.78 mm) was roughly 51 Hz according to the capillary-gravity wave dispersion relation

(2.5)\begin{equation} \omega^2 = \left(g k+\frac{\sigma}{\rho}k^3\right)\tanh{k H}, \end{equation}

where $k$ is calculated as $2{\rm \pi} / \lambda _c$. Owing to the standby oscillation of the shaker, whose frequency is determined by the electric current ($50\sim 60$ Hz), the hypothesis that natural shaking of the vibrator would result in harmonic ripples was found to be valid. Furthermore, no ripple-shaped pattern occurred in the reconstructed surface when the shaker was not operating.

In the tests, the vibration-induced effect $\delta \boldsymbol {r}_v$ dominated the other two effects $\delta \boldsymbol {r}_c$ and $\delta \boldsymbol {r}_o$ in the reconstruction of resonant waves. The free-surface deformation field $\delta \boldsymbol {r}$ of Faraday waves was sufficiently large to mask the aforementioned spurious vibration-induced deformations and unavoidable multi-frequency deformations, which were hard to detect by the naked eye. Moreover, given that multiple frequency components coexist in resonant waves, it was not necessary to remove the relatively minor error deformations at a specific frequency.

In short, precise adjustment of the vibrating system and post-processing of the raw deformation field $\delta \boldsymbol {r}^*$ helped ensure the reconstructed results were acceptably accurate. In our systems, a translation/zoom filter was applied to acquire the post-processed deformation field $\delta \boldsymbol {r}^*-\delta \boldsymbol {r}_v$ and a Butterworth band-pass filter was implemented to identify monochromatic deformation.

The inverse gradient operator $\nabla ^{-1}$ was used to reconstruct the surface topography $\eta$ from the free-surface gradient field (see figure 2e). For a circular domain $D$ comprising $Q$ discrete points in figure 4(a), different schemes were utilised in the expression for height gradient at given points. For those points whose surrounding grid points were not entirely located in the domain $D$, a pinned boundary condition was applied to one or two components of the gradient, and so these points were artificially defined as near-boundary points. The remaining points inside $D$ were classified as internal points. Third-order four-point differences were used to approximate the $\hat {\boldsymbol {x}}$-component $\nabla _1 \eta$ and $\hat {\boldsymbol {y}}$-component $\nabla _2 \eta$ of the gradient $\boldsymbol {\nabla } \eta$ at internal points, with all four grid points involved in each computation determined according to the quadrant in which the target point lay. Taking the $\hat {\boldsymbol {x}}$ component of the gradient as an example, $\nabla _1 \eta$ at grid points in each of the four quadrants $D_i$ ($i=1,2,3,4$; see figure 4b) is given by

(2.6)\begin{equation} \nabla_1 \eta(P_{i,j})=\left\{ \begin{array}{@{}ll@{}} \displaystyle (\eta_{i-2,j}-6 \eta_{i-1,j}+3 \eta_{i,j}+2\eta_{i+1,j})/6s, & P_{i,j} \in D_1 \cup D_4,\\ \displaystyle ({-}2 \eta_{i-1,j}-3 \eta_{i,j}+6\eta_{i+1,j}-\eta_{i+2,j})/6s, & P_{i,j} \in D_2 \cup D_3, \\ \end{array}\right. \end{equation}

where $s$, defined as element size, is the spatial distance between each adjacent grid point. The $\hat {\boldsymbol {y}}$-component $\nabla _2 \eta$ of an internal point is obtained in a similar manner.

Figure 4. (a) Classified grid points near the circular boundary (circle, near-boundary point and square, internal point). (b) Four computational sectors where different differential formats are applied.

The gradient $\boldsymbol {\nabla } \eta$ components at a point close to the boundary were related to nearby interior and boundary grid points. We found that the central difference format used by Moisy et al. (Reference Moisy, Rabaud and Salsac2009) was not sufficient to reconstruct the surface with a curved boundary, leading to reconstruction without convergence. Hence, three-point difference formats of higher order were used to approximate the gradient $\boldsymbol {\nabla } \eta$ at point $P_{i,j}$ as

(2.7)\begin{align} \boldsymbol{\nabla} \eta_{i,j}&=\left[\frac{s^2-(\Delta x)^2}{s\Delta x (\Delta x+s)}\eta_{i,j}+\frac{\Delta x}{s (\Delta x+s)}\eta_{i+1,j}\right]\hat{\boldsymbol{x}}\nonumber\\ &\quad +\left[\frac{-\Delta y}{s (\Delta y+s)}\eta_{i,j-1}+\frac{(\Delta y)^2-s^2}{s \Delta y (\Delta y+s)}\eta_{i,j}\right]\hat{\boldsymbol{y}}, \end{align}

where $\Delta x$, $\Delta y$ are the distances from $P_{i,j}$ to the corresponding boundary points, respectively. Due to the pinned boundary of the wave surface, a homogeneous Dirichlet boundary condition ($\eta =0$) was applied at these points on $\partial D$ such that the corresponding gradient was related to two points. The gradient of a near-boundary point was similarly quadrant determined, with both components in a four-point format, as was the case with the $\hat {\boldsymbol {x}}$ component of $\boldsymbol {\nabla } \eta (P_{i+1,j})$. Occasionally, second-order centred differences, or forward or backward Euler differences, were used to approximate the gradient in cases where there were insufficient grid points to construct a higher-order difference scheme. Over-determined linear systems with $2Q\times Q$ elements were generated, and a least-squares solution was readily obtained using the ‘$\backslash$’ operator in Matlab.

The accuracy of the surface reconstruction was evaluated in part by considering volume conservation in terms of a height error defined as the ratio of a discrete integral of volume variation to the cross-sectional area of the cylindrical container,

(2.8)\begin{equation} h_e = \frac{1}{{\rm \pi} R^2} \sum_{i=1}^{Q} \eta_i s^2, \end{equation}

in which $s$ is the element size, which is related to the subset size in the AL-DIC computation. The height error was sensitive to the initial location of reference images, and so repeatability experiments were utilised to improve results with the smallest $h_e$. After confirming all reconstructed surfaces of experimental modes, the final result, $h_e \ll \xi$, was valid for every subharmonic case, indicating the remarkable accuracy and applicability of the FS-SS method in our system. The main reconstruction error arose from two sources: inevitably inaccurate identification of the circular boundary because of the discrete length scale (pixel); and possibly poor plane fitting $\delta \boldsymbol {r}_v$ for relatively steep-sloped waves by the translation/zoom filter.

4.1. Spatiotemporal surface fitting

We hypothesise that the surface function of a certain mode $(\nu,\zeta )$ excited above the Faraday threshold may be approximated by a series of $\nu$-basis modes $(m,n)_\nu$ (3.1), such that the spatial decomposition is given by

(4.1)\begin{equation} \eta(r,\theta,t)=\sum_{m=1}^{M}\sum_{n=1}^{N}a_{m,n}(t) \cos{m\nu\theta} \,\mathrm{J}_{m\nu}(\delta_{m\nu,n}r), \end{equation}

where the total modal number $MN$ is determined by the calculated limits $M$ and $N$ of two multipliers $m$ and $n$, respectively. Normally, $N$ is set to an integer value much greater than $\zeta$ in order to obtain a more accurate expression for radial structures in the deformed surface. Additionally, an adequately large value of $M$ ($M\geq 2$) should be identified because of the existence of circumferential structures with multiple azimuths, particularly for small local features that occur during the transition period. We commence by fitting the reconstructed surface, generated by the FS-SS method, with a multi-mode function (4.1) and then analyse the modal information by examining the fitted coefficients $a_{m,n}$. For example, we evaluate the observed mode $(3,3)$ in figure 8 to indicate the accuracy of the 3-basis modal expression (4.1) using the spatial surface-fitting procedure.

We establish a corresponding equation for every grid point in the circular domain to generate a linear system containing $Q$ equations and then solve this over-determined problem with $MN$ unknowns. A similar method ‘$\backslash$’, which has been employed in the calculation of the inverse gradient field, is performed to acquire the least-squares solution of the fitted surface. We set $M=6$ and $N=15$ (a total of 90 $\nu$-basis modes) to obtain surface-fitting results for the investigated modes $(3,3)$. It can be seen in figure 11 that the spatial expression in the form of (4.1) provides a satisfactory fit to the reconstructed surface despite the small, complicated displacement of the modal transition. For other explored modes, the results of spatial fitting with $\nu$-basis expressions are found to be qualitatively of high accuracy, confirming that resonant modes $(\nu,\zeta )$ excited above the Faraday threshold can be linearly decomposed into corresponding $\nu$-basis modes $(m,n)_\nu$.

Figure 11. Comparison between (a) reconstructed patterns (see figure 8e) and (b) fitted patterns using 3-basis modal expressions. A common colour bar is given for all patterns where $|\eta |_{max}$ denotes the maximal modulus of surface elevation calculated from the reconstructions of each pattern in the first row (a). The close agreement demonstrates the validity of the spatial surface-fitting scheme used herein. Complete surface-fitting patterns over a wave period can be seen in the supplementary animation.

We now examine the spatiotemporal characteristics of the $\nu$-basis modes and perform further analysis on the time-varying coefficients $a_{m,n}(t)$ to reveal the magnitude and periodicity of $(m,n)_\nu$. Many $\nu$-basis modes are applied in the surface-fitting computation, and intuitive speculation concerning the dominant modes of the observed $(\nu,\zeta )$ mode suggests that the largest two are $(1,\zeta )_\nu$ and $(1,\zeta -1)_\nu$ due to constraint by the cylindrical boundary on the spatial structures of $(\nu,\zeta )$. Furthermore, it is useful to discern other relatively dominant wave components in such mode superposition, and so we conduct a comparative analysis of the $a_{m,n}(t)$ coefficients. Figure 12 displays the magnitude $\gamma _{m,n}$ plots of $\nu$-basis modes for three investigated modes, where $\gamma _{m,n}$ denotes the amplitude of time varying $a_{m,n}(t)$,

(4.2)\begin{equation} \gamma_{m,n}=\max_t{\left|a_{m,n}(t)\right|}. \end{equation}

The magnitude of $\gamma _{m,n}$ determines the amplitude of the $\nu$-basis mode $(m,n)_\nu$ for a given experimental mode $(\nu,\zeta )$, which also demonstrates the corresponding degree of dominance of the $\nu$-basis mode. In figure 12 the discrete scattered curves of magnitude $\gamma _{m,n}$ for different $m$ contain multiple peaks, which allow us to deduce the dominant modes. These $\nu$-basis modes at peaks associated with each $m$ are referred to as peak modes. The magnitudes $\gamma _{1,\zeta }$ and $\gamma _{1,\zeta -1}$ are found to be substantially greater than the other coefficients, in accordance with our earlier hypothesis. We suggest that the non-peak modes modulate the radial structures (wavenumbers) of peak modes to satisfy the dispersion relation. For this reason, our paper focuses solely on these peak modes, which are related to surface instability or wave interactions.

Figure 12. Magnitudes $\gamma _{m,n}$ of fitting coefficients $a_{m,n}$ for three experimental modes (a–c) with indicated different vibration driving parameters displayed as a function of azimuthal multiplier $m$ and radial multiplier $n$. Major peaks are labelled according to the dominant $(m,n)_\nu$ modes. Magnitudes $\gamma _{m,n}$ associated with given $m$ are indicated by the coloured bands. Results are shown for (a) $(2,3)$ $\xi = 0.2$ mm, $\varOmega _0/2{\rm \pi} = 15.5$ Hz; (b) $(4,2)$ $\xi = 0.175$ mm, $\varOmega _0/2{\rm \pi} = 15.2$ Hz; (c) $(3,3)$ $\xi = 0.15$ mm, $\varOmega _0/2{\rm \pi} = 18.5$ Hz.

For a given $\nu$-basis mode with $m$ and $n$ non-equivalent to any of the peak modes, the magnitude $\gamma _{m,n}$ is negligible, suggesting that wave energy is primarily stored in peak modes. Hence, the spatial structures of the explored mode $(\nu,\zeta )$ mainly depend on those of peak modes, and a larger number of peak modes leads to a more deformed modal pattern throughout the wave period. Apart from the peak modes, the magnitudes of the $\nu$-basis modes $(m,n)_\nu$ rapidly reduce with progressively larger $m$ and $n$. Moreover, $\gamma _{m,n}$ disappears when $m$ and $n$ reach critical integer values, implying weak or no excitation of the corresponding modes in the experiments. As a result, the upper limits of $m$ and $n$ can be ascertained in order to obtain a sufficiently precise fitting surface. However, not all peak modes can be regarded as dominant modes triggered by surface instability, which should be further classified according to frequency distribution, and so a detailed analysis follows concerning the temporal fitting results for $a_{m,n}(t)$.

We consider several peak modes with $m\leq 3$ for three cases of interest, which are labelled $(m,n)_\nu$ in figure 12. Further analysis evaluates the periodicity relation for these peak modes, based on frequency decomposition of the time varying $a_{m,n}(t)$. The observed Faraday modes are subharmonic resonance, excited at frequency $\omega _0$ equal to half the driving angular frequency such that $\omega _0=\varOmega _0 / 2$. The temporal decomposition is expressed as the following sum of multi-frequency cosine functions:

(4.3)\begin{equation} a_{m,n}(t) = \sum_{l=0}^{\infty} \chi_l \cos{l \omega_0 t}. \end{equation}

Here the fitting coefficients $\chi _l$ are different for various $a_{m,n}(t)$. Due to the temporal difference between the selected frames and expression (4.3), we adopt a more general form to fit the time-varying coefficient $a_{m,n}(t)$ in practical calculations, which reads as

(4.4)\begin{equation} a_{m,n}(t) = \sum_{l=0}^{\infty} \chi_l \cos{l \omega_0 t}+\sum_{l=1}^{\infty} \hat{\chi}_l \sin{l \omega_0 t}, \end{equation}

where $\hat {\chi }_l$ are coefficients corresponding to the additional sine terms ${\sin {l \omega _0 t}}$ (${l=1,2,3,\ldots}$) in the expansion. For computational purposes, the upper limit of $l$ is truncated to $L$.

We now let $L=10$, and perform linear fitting using the Matlab function lsqcurvefit() of the approximate expression (4.4) to discrete $a_{m,n}(t)$ values for the three experimental modes. Figure 13(a,c,e) shows the fitting results of indicated dominant modes over an evolving period $T_F$ ($=2{\rm \pi} /\omega _0$). It can be seen that the fitted curves of $a_{m,n}(t)$ are in relatively close agreement with the calculated evolution. Close-fit results are also obtained for the remaining modes in figure 12, confirming that (4.4) (or (4.3)) precisely replicates the periodicity of $\nu$-basis modes $(m,n)_\nu$. We then focus on the fitting coefficients $\chi _l$ and $\widehat {\chi _l}$. The right pseudo-colour illustrations in figure 13 show the normalized magnitudes $\hat {\gamma }_{m,n,l}$ of frequency components for the indicated modes $(m,n)_\nu$, each of which is defined as

(4.5)\begin{equation} \hat{\gamma}_{m,n,l}=\left\{ \begin{array}{@{}ll@{}} \displaystyle\dfrac{\left|\chi_0\right|}{\gamma_{m,n}}, & l=0, \\ \displaystyle\dfrac{\sqrt{\chi_l^2+\hat{\chi}_l^2}}{\gamma_{m,n}}, & l>0, \end{array} \right. \end{equation}

where $\chi _0$, $\chi _l$ and $\hat {\chi }_l$ ($l=1,2,\ldots,L$) are the fitting coefficients of $a_{m,n}$. For simplicity, we define the component with angular frequency $l\omega _0$ of $a_{m,n}(t)$ as the corresponding $l\omega _0$ component of the $\nu$-basis mode $(m,n)_\nu$, and peak modes with dominant $l\omega _0$ component as $l\omega _0$-dominant modes. As shown in figure 13(b,d, f), the $\omega _0$ and $2\omega _0$ components dominate, so these modes are purely subharmonic or harmonic. The most dominant $(1,\zeta )_\nu$ modes of the three different $(\nu,\zeta )$ modes oscillate at an angular frequency $\omega _0$. The dominant components of peak modes $(2,5)_2$, $(2,5)_3$ and $(2,3)_4$ are $2\omega _0$ components. On this basis, peak modes with dominant $\omega _0$ or $2\omega _0$ components are denoted as dominant modes, with $(1,\zeta )_\nu ^{sh}$ invariably the most dominant mode in the modal decomposition of $(\nu,\zeta )$. For these $\nu$-basis modes $(m,n)_\nu$ with constant $m$, no more than two dominant modes usually exist, and one is subharmonic and the other is harmonic if two dominant modes emerge with the same $m$, such as $(2,5)^{sh}$ and $(2,2)^{h}$ in figure 13(a,b). Meanwhile, the subharmonic dominant mode $(3,1)_3$ is indeed observed, as inferred in the discussion for figure 10.

Figure 13. Temporal fitting of dominant modes for three experimental cases of Faraday waves: plots (a,c,e) show the resulting $a_{m,n}(t)$ time series; and (b,d, f) pseudo-colour diagrams, where $\omega /\omega _0$ equals $l$, show the normalized magnitude $\hat {\gamma }_{m,n,l}$ of multiple frequency components, demonstrating the frequency distribution of selected dominant $\nu$-basis modes $(m,n)_\nu$. Deeper-coloured blocks with larger $\hat {\gamma }_{m,n,l}$ represent dominant frequency components of given dominant modes. Results are shown for (a,b) $(2,3)$ $\xi = 0.2$ mm, $\varOmega _0/2{\rm \pi} = 15.5$ Hz; (c,d) $(4,2)$ $\xi = 0.175$ mm, $\varOmega _0/2{\rm \pi} = 15.2$ Hz; (e, f) $(3,3)$ $\xi = 0.15$ mm, $\varOmega _0/2{\rm \pi} = 18.5$ Hz.

4.2. Peak modes

It is significant to focus on the spatiotemporal characteristics of peak modes because they strongly correlate with multi-azimuth structures appearing in the fast wave evolution. The pseudo-colour illustrations in figure 14(a–c) depict the logarithmic magnitude $\log _{10}{\gamma _{m,n}}$ of $MN$ $\nu$-basis modes in a $m$–$n$ plane for three investigated cases. Deeper red blocks marked with circles represent peak modes of relatively large magnitude in figure 12; these dominant modes are marked by blue pentagrams. In the $m$–$n$ plane we observe that the peak modes appear to be regulated in a seemingly uniform arrangement within a constrained range of $m$ and $n$ values. We connect these peak modes with solid lines and dashed lines in two directions in figure 14(a–c). For simplicity, we heuristically refer to the sequences of uniformly arranged peak modes on solid lines and dashed lines as peak chains and frequency levels, respectively.

Figure 14. Pseudo-colour diagrams that indicate the magnitude $\log _{10}{\gamma _{m,n}}$ of the $\nu$-basis modes for three investigated modes (a–c), with corresponding peak modes in figure 12 marked by hollow circles. These dominant modes are also marked by blue pentagrams. Peak modes are uniformly arranged in the $m$–$n$ plane and connected by solid lines and dashed lines in two directions, which refer to peak chains and frequency levels, respectively. Four peak chains of mode $(3,3)$ are identified in (d). The pseudo-colour illustrations of $\hat {\gamma }_{m,n,l}$ demonstrate the frequency distribution of peak modes $(m,n)_\nu$ on a certain peak chain (e) and four specific frequency levels. Colour scales for $\log _{10}{\gamma _{m,n}}$ and $\hat {\gamma }_{m,n,l}$ are given on the right. Results are shown for (a) $(2,3)$ $\xi = 0.2$ mm, $\varOmega _0/2{\rm \pi} = 15.5$ Hz; (b) $(4,2)$ $\xi = 0.175$ mm, $\varOmega _0/2{\rm \pi} = 15.2$ Hz; (c) $(3,3)$ $\xi = 0.15$ mm, $\varOmega _0/2{\rm \pi} = 18.5$ Hz; (e) $(1,3)_3 \to (2,5)_3 \to (3,7)_3 \to (4,9)_3$.

The peak chains starting at $(1,\zeta )_\nu$ for modes $(2,3)$, $(3,3)$ and $(4,2)$ are given by

(4.6)\begin{equation} \left. \begin{aligned} (1,3)_2 & \rightarrow (2,5)_2 \rightarrow (3,7)_2 \rightarrow (4,9)_2 \rightarrow \cdots, \\ (1,3)_3 & \rightarrow (2,5)_3 \rightarrow (3,7)_3 \rightarrow (4,9)_3 \rightarrow \cdots, \\ (1,2)_4 & \rightarrow (2,3)_4 \rightarrow (3,4)_4 \rightarrow (4,5)_4 \rightarrow (5,6)_4 \rightarrow \cdots, \end{aligned} \right\} \end{equation}

where the right arrows specify the primary transmission direction of energy and spatiotemporal information as discussed in the following sections. Interestingly, the indices $m$ and $n$ of a particular peak chain are in the order of an arithmetic progression, i.e. the next peak mode after $(m,n)_\nu$ is always $(m+1,n+2)_\nu$ for modes $(2,3)$ and $(3,3)$, whereas that for $(4,2)$ is $(m+1,n+1)_\nu$, which is related to the first two peak modes in any given peak chain. Hence, we suggest that the uniform arrangement of peak modes is arithmetic-progression-like. Moreover, more than one uniformly arranged peak chain can be discerned for any arbitrary $(\nu,\zeta )$ mode. For example, four peak chains for the $(3,3)$ mode are indicated in figure 14(d), where the peak modes are arranged according to the solid lines and dashed lines in figure 14(c). A similar analysis is next carried out for the frequency characteristics of peak modes on solid lines (peak chains) and dashed lines (frequency levels) for $(3,3)$. As shown in figure 14(e), the frequency distribution for these modes in a certain peak chain demonstrates that the dominant frequency components $l\omega _0$ (deep blocks) are uniformly arranged with increasing $l$ along the peak chains, which is also arithmetic-progression-like. For example, the $l$th mode in the peak chains starting at $(1,\zeta )_\nu$ is $l\omega _0$ dominant. A more universal rule is that a mode in a certain chain is $l\omega _0$ dominant whereas the next mode is $(l+1)\omega _0$ dominant. Figure 14( f) shows the frequency distribution of these peak modes at four frequency levels. One can observe that those modes at the same frequency level possess the same dominant frequency ($l\omega _0$) provided that the frequency levels are sorted in a $1,2,\ldots,l,\ldots$ sequence corresponding to the increasing dominant frequency of modes along the peak chain. White blocks represent a lack of corresponding frequency components in the peak modes, whose uniform arrangements in figure 14(e, f) reveal that these peak modes on peak chains and frequency levels are either subharmonic (only odd $l$) or harmonic (only even $l$). However, an exception exists for $(2,2)$ whereby $(2,2)_2$ and $(2,5)_2$ have impure frequency distributions, as evident in figure 13(a,b), due to co-interference. This invariably occurs when two dominant modes occur with the same azimuthal structures excited by surface instability.

Hence, according to the frequency distributions, we define the modal decomposition of $(\nu,\zeta )$ by

(4.7)\begin{equation} (\nu,\zeta) = \sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\left[ (m,n)_\nu^{sh} + (m,n)_\nu^{h}\right], \end{equation}

where the quantities $m$, $n$ and $\nu$ indicate spatial characteristics, and the superscripts $sh$ and $h$ represent temporal characteristics. Accordingly, the temporal expression (4.3) is rewritten as the following combination of subharmonic and harmonic components:

(4.8)\begin{equation} a_{m,n}(t) = \underbrace{\sum_{l=0}^{\infty} \chi_l^n \cos{2 l \omega_0 t}}_{\text{harmonic}}+\underbrace{\sum_{l=0}^{\infty} \hat{\chi}_l^n \cos{(2 l+1)\omega_0 t}}_{\text{subharmonic}}. \end{equation}

4.3. Modal interaction with energy transfer

We next pay attention to the mechanism of arithmetic-progression-like arrangements of peak modes in the $m$–$n$ plane (see figure 14a–c) and their frequency distributions (see figure 14e, f). We recall that the spatiotemporal decomposition of mode $(\nu,\zeta )$ is in the form of $\nu$-basis modes $(m,n)_\nu$, and the minimum spatiotemporal component of $(m,n)_\nu$ is expressed as $\cos {l\omega _0 t}\cos {m\nu \theta }\,\mathrm {J}_{m\nu }(\delta _{m\nu,n}r)$. To describe the generation of the arithmetic combination of mode indices, we provide a qualitative derivation based on the simple multiplication of two minimum spatiotemporal components, which represents the interaction of the $l_1\omega _0$ component of mode $(m_1,n_1)_\nu$ and $l_2\omega _0$ component of mode $(m_2,n_2)_\nu$,

(4.9) $$\begin{gather} \cos{l_1 \omega_0 t}\cos{m_1 \nu \theta}\,\mathrm{J}_{m_1\nu}(\delta_{m_1\nu,n_1}r)\, \cos{l_2 \omega_0 t}\cos{m_2 \nu \theta}\,\mathrm{J}_{m_2\nu}(\delta_{m_2\nu,n_2}r) \nonumber\\ = \underbrace{\cos{m_1 \nu \theta}\cos{m_2 \nu \theta}\,\mathrm{J}_{m_1\nu}(\delta_{m_1\nu,n_1}r) \mathrm{J}_{m_2\nu}(\delta_{m_2\nu,n_2}r)}_{\text{(I)}} \,\underbrace{\cos{l_1 \omega_0 t} \cos{l_2 \omega_0 t}}_{\text{(II)}}, \end{gather}$$

where $m_1$, $m_2$, $n_1$, $n_2$, $l_1$ and $l_2$ are positive integers; and terms (I) and (II) are spatially and temporally related products. Self-interaction is approximated by the overall operation (4.9) with $m_1=m_2$, $n_1=n_2$ and $l_1=l_2$. Due to the pinned boundary condition, the zero points of the two Bessel functions in term (I) at $r=1$ are coincident, and so the number of nodes of the Bessel function product in the radial interval ($0< r\leq 1$) is $n_1+n_2-1$. Moreover, the azimuthal terms $\cos {(m_1-m_2)\nu \theta }$ and $\cos {(m_1+m_2) \nu \theta }$ are derived from the multiplication of two cosine functions. Hence, we obtain a criterion that the interaction of $(m_1,n_1)_\nu$ and $(m_2,n_2)_\nu$, which both lie on a certain peak chain, is necessary for a new mode $(m_1+m_2,n_1+n_2-1)_\nu$ to be created belonging to the same peak chain. It can be verified that peak chains in figure 14(a–c) satisfy this criterion, and thus, the rule of arithmetic-progression-like arrangement of peak modes in the $m$–$n$ plane follows accordingly. Similarly, for term (II), two frequency components $\cos {|l_1-l_2| \omega _0 t}$ and $\cos {|l_1+l_2| \omega _0 t}$ are generated from the interaction between the $l_1\omega _0$ component and $l_2\omega _0$ component. The dominant frequency of these peak modes along the peak chain (4.6) also increases progressively as shown in figure 14(e), and is consistent with the simple nonlinear interaction of modes on peak chains. Here, using the simple multiplication (4.9), we provide typical examples about $(4,9)_3$ to display mode interaction on peak chains. The generation of $(4,9)_3$ originates from (a) the interaction of $(1,3)_3$ and $(3,7)_3$ and (b) the self-interaction of $(2,5)_3$, which both induce three frequency components ($0$, $2\omega _0$ and $4\omega _0$) as evident in figure 14(e). Self-interaction is usually a relatively weak process in wave interaction. Similar nonlinear interactions appear between arbitrary, two or more, peak modes, and so other peak modes of larger $m$ and $n$ with specific frequency distributions in peak chains emerge accordingly, which result in the uniform arrangements of peak modes and dominant frequency components displayed in figure 14. Note that peak chains (4.6) labelled with right-directed arrows indicate that the two preceding modes can (self-)interact to generate a third mode, which we refer to as resonance involving three standing modes. Although nonlinear wave interactions in the present experimental system are far more complex than represented by simple multiplication (4.9), such schemes nevertheless enable us to understand the spatiotemporal structures of harmonics induced by wave interaction.

Furthermore, although the aforementioned results are derived from the surface fitting of symmetric Faraday waves at the final steady state in experiments where the driving parameters much exceed the instability threshold, the chain-like dynamics of superposed standing wave components, identified as $\nu$-basis modes, enable us to backtrack the pattern formation. Rajchenbach et al. (Reference Rajchenbach, Clamond and Leroux2013) suggested a mechanism based on three-wave resonance to explain the triggering of parametrically forced surface waves with certain symmetries that were observed in the experiments. However, such a nonlinear mechanism can only describe the initial response and does not account for the presence of high-order harmonics in the final periodic state. Using modal decomposition, we gain insight into the mechanism that enhances symmetries triggered during the initial response period as they evolve to reach steady states. Here, analogous resonant coupling between three standing waves is used to describe the observed spatiotemporal characteristics of high-order harmonics, as introduced in (4.9). The interaction between two primary modes, $(m_1,n_1)_\nu$ and $(m_2,n_2)_\nu$, can result in the emergence of a secondary mode $(m_1+m_2,n_1+n_2-1)_\nu$, allowing for continuous energy supply. By extending the mechanism in Rajchenbach et al. (Reference Rajchenbach, Clamond and Leroux2013) and incorporating nonlinear wave interactions, we provide a comprehensive explanation of overall pattern formation in our experiments. At first, using energy supplied from the forcing vibration of the shaker, surface instability is triggered resulting in the emergence of dominant modes, which are either subharmonic ($\omega _0$ dominant) or harmonic ($2\omega _0$ dominant). Their symmetries are selected by the contact angle between two waves involved in three-wave resonance. Once these dominant modes are established, they become locked and subsequently grow. Persistent nonlinear wave interactions, including self-interactions, take place with sustained energy transfer, and peak chains and frequency levels with specific spatiotemporal characteristics are then uniquely determined. Consequently, newly formed wave components with complex spatiotemporal structures originate from the coupled interaction of several dominant waves and simultaneously extract a small amount of energy. After a certain number of forcing cycles, the pattern evolution process becomes steady, and a balance is achieved among energy supply, transfer and dissipation. In this way, energy transfer in each investigated wave $(\nu,\zeta )$ primarily occurs along the peak chain (4.6), and the arithmetic-progression-like arrangement is determined by dominant modes and wave interactions. Hence, the wave components involved can be classified as being of two types: dominant modes determined by the effect of parametric instability at the air–liquid interface; and harmonics induced by wave interactions. The generation of dominant modes is still not fully understood, and therefore, it is crucial to provide theoretical verification regarding the selection of these dominant modes.

5.1. Theory of inviscid Faraday waves in a brimful cylinder

Consider the vertical forced oscillation of liquid in a brimful cylindrical container of radius $R$ and height $H$ shown in figure 9(b). The container is filled to the brim with a low-viscosity liquid (i.e. water) of density $\rho$ and surface tension $\sigma$ such that the contact line of the liquid surface and the lateral wall remains pinned during the course of continuous oscillation. Surface patterns and interface deformation $\eta (r,\theta,t)$, restored by gravity $g$ and surface tension $\sigma$, are generated under the parametric driving force expressed by dimensionless acceleration amplitude $\epsilon$ and dimensionless angular frequency $\varOmega$. Here, linear spatial dimensions are scaled by $R$, time by $\sqrt {\rho R^3/\sigma }$, velocity potential by $\sqrt {\sigma R/\rho }$ and acceleration by $g$. Dimensionless governing equations are used in the linear stability analysis of the free-surface disturbance. Assuming the liquid is inviscid and irrotational, the velocity potential $\phi$ satisfies Laplace's equation

(5.1)\begin{equation} \nabla^2 \phi=0. \end{equation}

No-penetration conditions at the bottom and lateral walls of the container are given by

(5.2a,b)\begin{equation} \left.\frac{\partial \phi}{\partial z}\right|_{z={-}h}=0,\quad \left.\frac{\partial \phi}{\partial r}\right|_{r=1}=0, \end{equation}

where dimensionless height $h=H/R$.

At the free surface, assumed to be at $z=0$, the kinematic condition,

(5.3)\begin{equation} \frac{\partial \phi}{\partial z}=\frac{\partial \eta}{\partial t}, \end{equation}

and the dynamic condition,

(5.4)\begin{equation} \frac{\partial\phi}{\partial t}+Bo(1-\epsilon\cos{\varOmega t})\eta-\left(\frac{\partial^2\eta}{\partial r^2}+\frac{1}{r}\frac{\partial\eta}{\partial r}+\frac{1}{r^2}\frac{\partial^2\eta}{\partial \theta^2}\right)=0, \end{equation}

both hold under linearisation for small deformation $|\eta |\ll 1$, where $Bo=\rho g R^2/\sigma$, $\epsilon =A_0/g$ and $\varOmega =\varOmega _0\sqrt {\rho R^3/\sigma }$.

The pinned contact-line condition and the volume conservation condition in a cylindrical domain are given by

(5.5)\begin{equation} \left.\eta\right|_{r=1}=0 \end{equation}

and

(5.6)\begin{equation} \int_0^{2{\rm \pi}} \int_{0}^1 \eta(r,\theta,t) r\, \mathrm{d}r\,\mathrm{d}\theta=0. \end{equation}

In our experiments we explored modes excited above the Faraday threshold where wave resonance is ubiquitous. Although the modal patterns remained symmetric and periodic, small deformations occurred (mainly in the transition period) compared with pure modes at critical conditions. Although the azimuthal mode number $\nu$ of a certain mode can be identified intuitively, $\nu$ is not in itself sufficient to express transitional patterns where azimuthal nodes are typically multiples of $\nu$. We recall that components of multi-azimuth modes were introduced to indicate the spatial structure of Faraday waves in § 3. Inspired by the forms of $\eta$ and $\phi$ given by Kidambi (Reference Kidambi2013) who derived a remarkably applicable expression for the onset of pure mode oscillation, we utilise the following pair of functions formed from eigenmodes $\cos {m\nu \theta }\,\mathrm {J}_{m\nu }$ with multiple azimuths in a cylindrical domain,

(5.7)\begin{equation} \eta(r,\theta,t)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}a_{m,n}(t) \cos{m\nu\theta}\, \mathrm{J}_{m\nu}(\delta_{m\nu,n}r) \end{equation}

and

(5.8)\begin{equation} \phi(r,\theta,z,t)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}b_{m,n}(t) \cos{m\nu\theta} \frac{\cosh{k_{m\nu,n}(z+h)}}{\sinh{k_{m\nu,n}h}}\,\mathrm{J}_{m\nu}(k_{m\nu,n}r), \end{equation}

where $\delta _{m\nu,n}$ is the $n\mathrm {th}$ root of the Bessel function $\mathrm {J}_{m\nu }$, and $k_{m\nu,n}$ is the $n\mathrm {th}$ root of function $\mathrm {J}_{m\nu }^{\prime }$. Equations (5.7) and (5.8) are two improved functions comprising complex components $\cos {m\nu \theta }$ coupled with radial parameter $r$ that satisfy governing equations ((5.1)–(5.2)) and (5.5). The truncated form of (5.7) is given by (4.1), which characterizes the spatial structures of the modes investigated in our experiments. It should be noted that the aforementioned pair of functions only ensure sufficient validity provided the following constraints are met.

The volume conservation condition (5.6) holds naturally for non-axisymmetric modes ($\nu \neq 0$). For axisymmetric modes ($\nu =0$), where multi-azimuth structures are not observed, volume conservation requires

(5.9)\begin{equation} \sum_{n=1}^{\infty} a_n \int_0^1 \mathrm{J}_{0}(\delta_{0,n}r)r\,\mathrm{d}r=0. \end{equation}

Similar to the derivation procedure by Kidambi (Reference Kidambi2013), we obtain coupled Mathieu equation systems for each $m$ by introducing ((5.7)–(5.8)) into ((5.3)–(5.4)) and solving coefficients $b_{m,n}$ from coefficients $a_{m,n}$. With truncated limits $M$ and $N$ respectively for $m$ and $n$ ($MN$ $\nu$-basis modes), the $m$-determined Mathieu equations are given by

(5.10)\begin{equation} \frac{\mathrm{d}^2\boldsymbol{a}_m}{\mathrm{d}T^2}+({{\boldsymbol{\mathsf{P}}}}-2{{\boldsymbol{\mathsf{Q}}}}\cos{2 T})\boldsymbol{a}_m=0, \end{equation}

where $T=\varOmega t/2$, $\boldsymbol {a}_m$ is the amplitude vector composed of $a_{m,n}\ (n=1,2,\ldots,N)$ and two non-diagonal functional matrices ${{\boldsymbol{\mathsf{P}}}}$ and ${{\boldsymbol{\mathsf{Q}}}}$, whose detailed elements are shown in Appendix A, are related to the driving parameters, Bond number and given azimuthal mode number $\nu$. Therefore, $M$ numerical systems comprising coupled Mathieu equations are obtained. Note that we have included multi-azimuth wave terms $\cos {m\nu \theta }\, \mathrm {J}_{m\nu }$ into the expressions for the free-surface elevation and velocity potential to describe steady-state pattern evolution, and thus, different coupled Mathieu equations are generated according to each azimuthal multiplier $m$. Herein, the stability analysis of parametric resonance above the Faraday threshold involves coupled analysis of corresponding stability results from $M$ numerical systems.

The Mathieu equation (5.10) has been widely applied to the mechanics of nonlinear oscillations (McLachlan Reference McLachlan1951; Nayfeh & Mook Reference Nayfeh and Mook1995). Floquet theory is typically used to identify the solutions and perform a stability analysis of the Mathieu equation. We use Hill's infinite determinant method (Hill Reference Hill1886) to determine the form of the series solution and the stability regions of the Mathieu equation (5.10) in the $\varOmega$–$\epsilon$ plane. In previous studies of Faraday waves, Nayfeh & Mook (Reference Nayfeh and Mook1995) and Kidambi (Reference Kidambi2013) utilised the infinite determinant method when conducting stability analysis of the $m=1$ mode; here we extend their work to determine solutions for larger-$m$ modes.

According to Floquet theory, we use a solution of the $m$th system (5.10) in a truncated form of Fourier expansion with $2L+1$ unknown coefficient vectors with $N$ components, $\boldsymbol {\alpha }_l$ ($l=0,1,\ldots,L$) and $\boldsymbol {\beta }_l$ ($l=1,2,\ldots,L$), such that

(5.11)\begin{equation} \boldsymbol{a}_m(T)={\rm e}^{\varrho_m T}\left(\sum_{l=0}^{L}\boldsymbol{\alpha}_l {\rm e}^{\mathrm{i}2 l T}+\sum_{l=1}^{L}\boldsymbol{\beta}_l {\rm e}^{-\mathrm{i}2 l T}\right), \end{equation}

in which $\varrho _m$ is the characteristic exponent related to the azimuthal multiplier $m$. The vector $\boldsymbol {a}_m$ has period $2{\rm \pi}$ for subharmonic resonance and ${\rm \pi}$ for harmonic resonance, and $\varrho _m=\mathrm {i}$ (or 0) is selected to obtain the subharmonic (or harmonic) solution (Kidambi Reference Kidambi2013). As previously mentioned, (4.8) provides a highly accurate description of the temporal decomposition of the $\nu$-basis mode $(m,n)_\nu$. In practice, the selection of $\varrho _m$ is closely linked to the periodicity of each mode mentioned in § 4.1, and so we set $\varrho _m=\mathrm {i}$ for subharmonic components $(m,n)_\nu ^{sh}$ (or 0 for harmonic components $(m,n)_\nu ^{h}$) to obtain subharmonic (or harmonic) resonance regions. Equation (4.8) combines the real form of (5.11) for $\varrho _m=\mathrm {i}$ with that for $\varrho _m=0$. In this way, the final subharmonic solution of the explored modes, which is $\omega _0$ dominant in our experiments, comprises the sum of subharmonic and harmonic components.

Substituting (5.11) into (5.10) and taking orthogonality of the Fourier series $\mathrm {e}^{\mathrm {i}l 2 T}$ ($l=0,\pm 1,\pm 2,\ldots$) into consideration, we obtain the $m$-determined equation system, whose detailed generation is in Appendix B, given by

(5.12)\begin{equation} {{\boldsymbol{\mathsf{G}}}}\boldsymbol{v}=\epsilon{{\boldsymbol{\mathsf{H}}}}\boldsymbol{v}, \end{equation}

where ${{\boldsymbol{\mathsf{G}}}}$ and ${{\boldsymbol{\mathsf{H}}}}$, both $(2L+3)N\times (2L+1)N$, are two functional matrices originating from ${{\boldsymbol{\mathsf{P}}}}$ and ${{\boldsymbol{\mathsf{Q}}}}/\epsilon$, respectively, and $\boldsymbol {v}$ is the solution vector given by

(5.13a,b)\begin{equation} \boldsymbol{v}=[\hat{\boldsymbol{v}}_1;\cdots;\hat{\boldsymbol{v}}_n;\cdots\hat{\boldsymbol{v}}_N], \quad \hat{\boldsymbol{v}}_n = [\alpha_0^n,\ldots,\alpha_L^n,\beta_1^n,\ldots,\beta_L^n]^{\rm T}. \end{equation}

Note that the additional $2N$-row elements in matrices ${{\boldsymbol{\mathsf{G}}}}$ and ${{\boldsymbol{\mathsf{H}}}}$ are relevant to the terms $\mathrm {e}^{\pm \mathrm {i}2(L+1)T}$ arising from the expansion of $\cos {2T}{\rm e}^{\pm \mathrm {i}2LT}$ and so the corresponding rows of ${{\boldsymbol{\mathsf{G}}}}$ are zero. We adjust these zero rows to the bottom to obtain the following sub-matrices:

(5.14a,b)\begin{equation} {{\boldsymbol{\mathsf{G}}}}= \begin{bmatrix} \hat{{{\boldsymbol{\mathsf{G}}}}}_1 \\ {{\boldsymbol{\mathsf{0}}}} \end{bmatrix},\quad {{\boldsymbol{\mathsf{H}}}}= \begin{bmatrix} \hat{{{\boldsymbol{\mathsf{H}}}}}_1 \\ \hat{{{\boldsymbol{\mathsf{H}}}}}_2 \end{bmatrix}. \end{equation}

Here $\hat {{{\boldsymbol{\mathsf{G}}}}}_1$ and $\hat {{{\boldsymbol{\mathsf{H}}}}}_1$ are $(2L+1)N$-tuple sparse matrices, and $\hat {{{\boldsymbol{\mathsf{H}}}}}_2$ is the coefficient matrix of the $2N$ equations related to terms $\mathrm {e}^{\pm \mathrm {i}2(L+1)T}$.

Since the matrices ${{\boldsymbol{\mathsf{G}}}}$ and ${{\boldsymbol{\mathsf{H}}}}$ are not square, the system (5.12) is over-determined and so is difficult to solve. Hence, (5.12) needs to be altered so that it is solvable as an eigenvalue problem. To achieve this, two approximate methods are introduced. The first involves pre-multiplying (5.12) by a $(2L+1)N\times (2L+3)N$ matrix ${{\boldsymbol{\mathsf{H}}}}^{\dagger}$ and ${{\boldsymbol{\mathsf{H}}}}^{\dagger} {{\boldsymbol{\mathsf{H}}}}={{\boldsymbol{\mathsf{I}}}}$ (identity matrix) so that

(5.15a,b)\begin{equation} {{\boldsymbol{\mathsf{C}}}}\boldsymbol{v}=\epsilon\boldsymbol{v},\quad {{\boldsymbol{\mathsf{C}}}}={{\boldsymbol{\mathsf{H}}}}^{\dagger} {{\boldsymbol{\mathsf{G}}}}. \end{equation}

Then a least-square solution to the generalized inverse matrix ${{\boldsymbol{\mathsf{H}}}}^{\dagger}$ can be easily obtained with Matlab operator ‘/’, i.e. ${{\boldsymbol{\mathsf{H}}}}^{\dagger} ={{\boldsymbol{\mathsf{I}}}}/{{\boldsymbol{\mathsf{H}}}}$. The second method involves neglecting equations involving $\mathrm {e}^{\mathrm {i}2(L+1)T}$ and $\mathrm {e}^{-\mathrm {i}2(L+1)T}$ so that a straightforward generalized eigenvalue problem is generated as

(5.16)\begin{equation} \hat{{{\boldsymbol{\mathsf{G}}}}}_1\boldsymbol{v}=\epsilon\hat{{{\boldsymbol{\mathsf{H}}}}}_1\boldsymbol{v}. \end{equation}

In practice, the predictions of eigenvalue $\epsilon$ by the two methods are relatively close. The accuracy of the first method mainly depends on the inverse calculation of ${{\boldsymbol{\mathsf{H}}}}^{\dagger}$, and the method may not work well when $N$ is very small. Although the second method incurs some error, the highly sparse structure of the bottom blocks of ${{\boldsymbol{\mathsf{G}}}}$ and ${{\boldsymbol{\mathsf{H}}}}$ indicate that there is little impact on the eigenvalues (see also Kidambi Reference Kidambi2013). Note that the stability analysis in this paper is based on the eigenvalue problem (5.16) for these non-axisymmetric modes ($\nu \neq 0$), but for the axisymmetric modes, additional equations (5.9) should be taken into account.

In the numerical procedure, for given $\varOmega$ and constant experimental parameters, we generate the matrices ${{\boldsymbol{\mathsf{P}}}}$, ${{\boldsymbol{\mathsf{Q}}}}$, ${{\boldsymbol{\mathsf{G}}}}$ and ${{\boldsymbol{\mathsf{H}}}}$ sequentially, and solve (5.16) using the Matlab eig() function to obtain the $\epsilon$ eigenvalues. Then, stability regions for given $\nu$ and $m$ are reproduced in the $\varOmega$–$\epsilon$ plane by real eigenvalues ($\mathrm {Im}(\epsilon )=0$) as they intercept the complex space of $\epsilon$. It should be noted that, for a certain $\nu$, different $m$ lead to different ${{\boldsymbol{\mathsf{G}}}}$, ${{\boldsymbol{\mathsf{H}}}}$ resulting in $2M$ stability diagrams corresponding to $\nu$-basis modes $(m,{\cdot })_\nu ^{sh}$ and $(m,{\cdot })_\nu ^{h}$, where $m=1,2,\ldots,M$. Hence, joint analysis of $2M$ stability diagrams is necessary to confirm the presence of the dominant modes in figure 12.

5.2. Dominant modes triggered by surface instability

To connect the theoretical predictions with the experimental data, we consider stability analysis of the excited modes in figure 12. Recall that the surface expression (5.7) comprises the sum of $\nu$-basis modes $(m,n)_\nu$ where the upper limits of the modal indices are set to $m=M$ and $n=N$. The surface-fitting results in figure 12 indicate that: $N$ should be sufficiently large to obtain converged results for each $m$; and the combination of $\nu$-basis modes with a restricted $m$ should provide a fair representation of the spatial structure of the modes of interest. To this end, we first compute the stable and unstable regions, known as stability tongues, with increasing $N$ $(=1,2,3)$ and constant $L$ $(=5)$ for subharmonic $(1,{\cdot })_3^{sh}$ and harmonic $(1,{\cdot })_3^{h}$ modes. Figure 15(a,b) depict the bottom boundary of unstable tongues for subharmonic cases and harmonic cases, respectively. Unstable tongues in the dimensionless $\varOmega$–$\epsilon$ plane, determined from the eigenvalue equation (5.16), intersect to form fundamental resonance tongues (FRTs) whose bottoms correspond to $\epsilon =0$ due to the neglect of viscous damping in the analysis. Combination resonances, which normally appear on the FRTs when $Bo$ ($\approx 275$ herein) is finite, have two effects: the parametric range of unstable resonances increases; and the oscillatory evolution of corresponding modes is non-periodic (Kidambi Reference Kidambi2013). The FRTs emerge at fundamental frequencies (the bottom of tongues) where pure Faraday modes $(m,n)_\nu$ are excited. The CRTs arising from the combination frequencies densely overlap with FRTs, thus contributing to larger unstable regions. Additionally, the narrower harmonic tongues indicate that subharmonic resonances are more readily excited in vibration experiments. Both FRTs and CRTs increase with radial multiplier $n$, and the number of FRTs calculated using $N$ modes is always equal to $N$. Given that these FRTs are closely correlated with $\nu$-basis modes $(m,n)_\nu$, we artificially label the FRTs (from left to right) as $(1,n)_3^{sh}$ or $(1,n)_3^{h}$ ($n=1,2,3$). For other cases with varying $m$ and $\nu$, FRTs in the $\varOmega$–$\epsilon$ plane can be named similarly to those in the selected cases. Kidambi (Reference Kidambi2013) calculated CRTs by means of mapping at a period and found the corresponding oscillatory states were non-periodic. Combination resonances at the boundaries of the CRTs may be regarded as wave couplings among different modes of the same symmetry and slowly evolving amplitude, unlike the pattern competition in different symmetries discussed by Ciliberto & Gollub (Reference Ciliberto and Gollub1984, Reference Ciliberto and Gollub1985a). It is challenging to observe such quasi-periodic or chaotic combination resonances in experiments due to the narrow frequency bandwidths of CRTs and the requirement for careful selection of driving parameters to ensure other symmetric modes are absent. An experimental investigation of non-periodic CRTs is beyond the scope of the present paper, and so we restrict our discussion of modal verification to comparison of periodic observations with corresponding FRTs.

Figure 15. Predicted resonance tongues of harmonic modes (a) and subharmonic modes (b) obtained for $N=1,2,3$, $m=1$, $\nu =3$ and $L=5$. The fundamental resonance tongues (FRTs) are labelled sequentially according to corresponding $\nu$-basis modes. The insets show CRTs that emerge as $N$ increases, and which usually intersect with FRTs. Only the bottom boundaries of the resonance tongues are displayed here.

Subsequently, larger truncation limits were incorporated to obtain converged calculations. In analysing the emergence of dominant modes through spatiotemporal surface fitting, we primarily focus on the narrow domain of $0\leq \epsilon \leq 0.4$ because our experimental condition $(\epsilon < 0.4)$ is slightly above the Faraday threshold. This restricted-$\epsilon$ range can embrace small CRTs while avoiding the perplexing intersection of resonance tongues at large $\epsilon$, thus improving stability analysis for verifying the foregoing dominant $\nu$-basis modes. In the case of small $Bo$ = 275 (as investigated in the experiments), convergence necessitates numerous modes (large $N$) (Kidambi Reference Kidambi2013). As demonstrated in figure 16, converged FRTs are achieved with $N=15$ and $L=5$, signifying the parameter selection for fitting spatiotemporal structures of the experimental modes considered herein. Herein, the phenomenological criterion of a close positional relationship between experimental observations and stability tongues in the $\varOmega$–$\epsilon$ plane is used to verify the existence of dominant modes. In figure 16 the symbols denoting the driving parameters of four experimental modes $(3,\zeta )$ $(\zeta =2,3,4,5)$ are situated reasonably near the left boundary of the corresponding subharmonic FRTs $(1,\zeta )_3^{sh}$ $(\zeta =2,3,4,5)$, indicating that these most dominant $\nu$-basis modes $(1,\zeta )_3^{sh}$ (as described earlier in the analysis of surface-fitting results) can be possibly excited by the indicated driving parameters.

Figure 16. Converged stability plot for four subharmonic dominant tongues $(1,{\cdot })_3$ in a local window of small $\epsilon$. The four corresponding experimental modes are in one-to-one correspondence with these FRTs.

By utilising the close positional relationship between an investigated mode $(\nu,\zeta )$ and the (sub)harmonic FRTs $(m,n)_\nu$ in the $\varOmega$–$\epsilon$ plane, we are able to confirm qualitatively the potential existence of such dominant modes, denoted in figure 12. Figure 17 compares stability plots obtained for the driving parameters $(\varOmega,\epsilon )$ and the corresponding FRTs calculated by varying $m$ and $\varrho _m$ for three investigated modes $(2,3)$, $(3,3)$ and $(4,2)$. Several $\nu$-basis modes can be identified based on their close positional relationships, which invariably agree with the spatiotemporal surface-fitting results. For modal interaction $(\nu,\zeta )$, the principal subharmonic mode is typically $(1,\zeta )_\nu ^{sh}$ and the principal harmonic mode is $(2,\hat {\zeta })_\nu ^{h}$, where $\hat {\zeta }$ is determined by the dispersion relation. The existence of secondary subharmonic modes $(2,2)_2^{sh}$ and $(3,1)_3^{sh}$ is also confirmed, indicating the potential emergence of more than one subharmonic mode. These identified wave components are considered to be dominant modes triggered by surface instability, which are coincident with $\omega _0$-dominant or $2\omega _0$-dominant peak modes. Meanwhile, other incompatible modes, such as $(1,2)^{sh}_3$, $(1,4)^{sh}_3$ and $(3,2)^{sh}_3$ for the $(3,3)$ case, appear to exhibit no or only a weak response in the cases studied. Additional wave components with relatively weak amplitudes at the $\omega _0$-dominant or $2\omega _0$-dominant frequency levels can be identified from stability analyses for different $m$ and $\varrho _m$. Hence, it is both significant and convenient to classify qualitatively these dominant peak modes through linear stability analysis, enabling us to achieve a deep understanding of the mechanism behind the full-life wave response, as explained in § 4.3.

Figure 17. Stability plots for subharmonic tongues (a,c,e) and harmonic tongues $(2,{\cdot })_\nu ^{h}$ (b,d, f) corresponding to three experimental modes $(2,3)$, $(3,3)$ and $(4,2)$. Black triangles denote the indicated parameters of the three investigated modes. Only overlapping CRTs associated with these FRTs representing tongues $(m,n)_\nu ^{sh}$ and $(m,n)_\nu ^{h}$ are displayed here. The close positional relationships between the experimental observations and stability tongues indicate the dominant modes excited by surface instability. Results are shown for (a,b) $(2,3)$ $\xi = 0.2$ mm, $\varOmega _0/2{\rm \pi} = 15.5$ Hz; (c,d) $(4,2)$ $\xi = 0.175$ mm, $\varOmega _0/2{\rm \pi} = 15.2$ Hz; (e, f) $(3,3)$ $\xi = 0.15$ mm, $\varOmega _0/2{\rm \pi} = 18.5$ Hz.

In addition, validation tests of $m=1$ were performed by comparisons with experimental observations by Shao et al. (Reference Shao, Wilson, Saylor and Bostwick2021b) and numerical results by Kidambi (Reference Kidambi2013). Our calculated results exhibit reasonable agreement with those of Shao et al. (Reference Shao, Wilson, Saylor and Bostwick2021b) and reproducibility with those of Kidambi (Reference Kidambi2013). However, the close positional relationships in our analyses present a slight rightward shift in the discrepancies of instability tongues, which are strongly influenced by nonlinear effects present in the experiments. Based on high-resolution surface reconstructions, averaged estimates of local surface steepness can readily be made to evaluate the effect of nonlinearity. For the surface field in figure 8(b), the maximum radial steepness $\Delta \eta /\Delta r|_{\theta =0}\approx 0.125$ and the maximum azimuthal steepness $\Delta \eta /(r\Delta \theta )|_{r=0.357}\approx 0.123$ can be calculated using the ratio of elevation differences ($\Delta \eta$) to transversal spatial differences ($\Delta r$ or $r\Delta \theta$) of neighbouring peaks and troughs. The maximum local steepness is small but not insignificant, indicating that the nonlinearity involved is not negligible but instead relatively weak. Additionally, natural frequencies for the pinned contact-line cases are invariably larger than those with freely sliding contact-line dynamics (Bostwick & Steen Reference Bostwick and Steen2015; Wilson et al. Reference Wilson, Shao, Saylor and Bostwick2022). This suggests that the frequency difference between the experimental case and theoretical prediction may be due to the non-ideal contact-line condition in our experiments, and the loss of the pinned end being exacerbated by the strong forcing. Meanwhile, decoupling analyses for different-$m$ cases cannot explain existing nonlinear wave behaviour, which is typically described by cubic or higher-order terms in amplitude equations (Meron Reference Meron1987). The foregoing highlights the need for a coupled, nonlinear analysis of subharmonic and harmonic solutions for comparison against experimental resonances excited above the Faraday threshold. However, the foregoing positional relationships are sufficient for qualitative confirmation of these dominant modes, thus, further supporting the validity and rationality of modal decomposition.