1. Introduction
Orbital shaking is a method to gently mix the liquid content of a container by its displacement at fixed container orientation along a circular trajectory and at a constant angular velocity. It is used in biological and chemical industrial applications, notably bacterial and cellular cultures (McDaniel & Bailey Reference McDaniel and Bailey1969; Wurm Reference Wurm2004), as an alternative to stirred tanks, where liquid agitation results from a rotating impeller or the rotation of a magnetic rod. In these cultivation protocols, cells are in suspension in an extracellular liquid medium, which serves as buffer for consumables from which they feed and for their secretions. The motion of the liquid prevents sedimentation and homogenizes the concentration of dissolved oxygen and nutrients and of secreted proteins and carbon dioxide. Because of the possible gas exchanges at the free surface, oxygen supply from the container bottom can possibly be circumvented, avoiding the formation of bubbles and thereby the damage that their collapse can exert on cells (Handa-Corrigan, Emery & Spier Reference Handa-Corrigan, Emery and Spier1989; Kretzmer & Schügerl Reference Kretzmer and Schügerl1991; Papoutsakis Reference Papoutsakis1991), sparking interest in the development of large-scale, in the hectolitre range, orbital-shaken bioreactors (Liu & Hong Reference Liu and Hong2001; Jesus et al. Reference Jesus, Girard, Bourgeois, Baumgartner, Jacko, Amstutz and Wurm2004; Muller et al. Reference Muller, Derouazi, van Tilborgh, Wulhfard, Hacker, Jordan and Wurm2007). It is therefore not a surprise that a significant body of research on gas exchange and mixing in these devices has emerged over the last two decades (Büchs et al. Reference Büchs, Maier, Milbradt and Zoels2000a,Reference Büchs, Maier, Milbradt and Zoelsb; Büchs Reference Büchs2001; Maier, Losen & Büchs Reference Maier, Losen and Büchs2004; Muller et al. Reference Muller, Girard, Hacker, Jordan and Wurm2005; Micheletti et al. Reference Micheletti, Barrett, Doig, Baganz, Levy, Woodley and Lye2006; Zhang et al. Reference Zhang2009; Tissot et al. Reference Tissot, Farhat, Hacker, Anderlei, Kühner, Comninellis and Wurm2010; Tan, Eberhard & Büchs Reference Tan, Eberhard and Büchs2011; Tissot et al. Reference Tissot, Oberbek, Reclari, Dreyer, Hacker, Baldi, Farhat and Wurm2011; Klöckner & Büchs Reference Klöckner and Büchs2012).
Since the shear stresses and, therefore, the mixing are proportional to the velocity gradients in the liquid phase, most of the gas exchange phenomena listed above are directly linked to the liquid motion, with the optimal working conditions essentially dictated by the wave pattern (Reclari Reference Reclari2013). For these reasons, at a more fundamental level, the hydrodynamics of these orbital shaking devices has received recent attention, from both experimental (Reclari et al. Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014; Bouvard, Herreman & Moisy Reference Bouvard, Herreman and Moisy2017; Moisy, Bouvard & Herreman Reference Moisy, Bouvard and Herreman2018) and theoretical (Reclari et al. Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014; Horstmann, Herreman & Weier Reference Horstmann, Herreman and Weier2020) perspectives, predominantly using linear potential flow models. These models are often complemented with effective viscous damping rates to incorporate the energy dissipation responsible for the phase shifts between wave and shaker, which was also seen to be sometimes responsible for damping-induced symmetry-breaking linear mechanisms resulting in linear spiral wave patterns (Horstmann et al. Reference Horstmann, Herreman and Weier2020, Reference Horstmann, Anders, Kelley and Weier2021). Previous studies mostly made use of classical existing theories for general linear sloshing dynamics, reviewed for instance in Ibrahim (Reference Ibrahim2005) and Faltinsen & Timokha (Reference Faltinsen and Timokha2009).
In order to refine the linear potential model and, specifically, to predict the occurrence of the super-harmonic wave dynamics observed experimentally (by super-harmonic, we mean here a wave of a certain frequency 
$\omega$ emerging from an excitation at 
$\varOmega =\omega /2$, with 
$\varOmega$ the driving angular frequency), Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) proposed an inviscid weakly nonlinear (WNL) analysis based on a second-order straightforward asymptotic expansion procedure, which was shown to be capable of capturing the observed resonance frequencies and of characterizing different multiple-crest wave patterns. Among these patterns, the super-harmonic double-crest (DC) wave dynamics is particularly relevant, as it appears to be the most stable and the one that displays the largest amplitude response. However, their analysis, as typical of straightforward asymptotic expansions, suffers from secular terms (Castaing Reference Castaing2005; Nayfeh Reference Nayfeh2008) and, therefore, it still fails in describing the correct nonlinear behaviour close to both harmonic and super-harmonic resonances.
With regards to the experiments of Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), Timokha & Raynovskyy (Reference Timokha and Raynovskyy2017) and Raynovskyy & Timokha (Reference Raynovskyy and Timokha2018a,Reference Raynovskyy and Timokhab) have applied the Narimanov–Moiseev multimodal sloshing theory (Narimanov Reference Narimanov1957; Moiseev Reference Moiseev1958; Dodge, Kana & Abramson Reference Dodge, Kana and Abramson1965; Faltinsen Reference Faltinsen1974; Narimanov, Dokuchaev & Lukovsky Reference Narimanov, Dokuchaev and Lukovsky1977; Lukovsky Reference Lukovsky1990). The theory is capable of accurately describing the nonlinear wave dynamics near the primary harmonic resonance, when no secondary resonances occur (Faltinsen, Rognebakke & Timokha Reference Faltinsen, Rognebakke and Timokha2005; Faltinsen, Lukovsky & Timokha Reference Faltinsen, Lukovsky and Timokha2016). Despite the fact that the experiments performed by Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) were made for non-dimensional fluid depths 
$H=h/R=1.04$ and 
$1$, which lie slightly beyond the applicability threshold of the multimodal theory (
$H_{th}$ should be 
$\gtrsim 1.05$ as stated by Raynovskyy & Timokha (Reference Raynovskyy and Timokha2020)) and imposed by the occurrence of secondary resonances, the authors found a quantitative good agreement with the experimental observations associated with hard-spring-type single-crest (SC) swirling.
In the spirit of the aforementioned multimodal theory but with regards to square-base basins, the resonant amplification of higher-order modes for forcing frequency in the vicinity of the primary resonance (secondary or internal resonances) was investigated by Faltinsen et al. (Reference Faltinsen, Rognebakke and Timokha2005), who formalized a so-called adaptive asymptotic modal approach capable of improving the agreement with earlier experiments. A thorough discussion on this regard is also outlined in Chapters 8 and 9 of Faltinsen & Timokha (Reference Faltinsen and Timokha2009), where the importance of the ratio of tank liquid depth to tank breadth for the occurrence of the internal resonance phenomenon is carefully discussed. Generally speaking, secondary resonance is a broader concept and it may occur even far from the primary resonance zone, as in the case of the DC swirling observed in Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). To the knowledge of the authors, the adaptive modal approach has never been extended to super-harmonic system responses of orbital-shaken circular cylindrical containers far from the primary resonance.
For these reasons, it appears that a quantitatively accurate model for the prediction of the diverse wave dynamics observed during the thorough experimental campaign carried out by Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) has not been provided yet.
The present work is precisely dedicated to the development of a WNL analysis based on the multiple-time-scale method, which will be seen to successfully capture nonlinear effects for the main additive harmonic resonances as well as the more subtle additive and multiplicative resonance governing the super-harmonic DC swirling. Amplitude equations are rigorously derived in an inviscid framework, which, once amended with an ad hoc damping term as the only tuning parameter, well match the experimental findings of Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). Lastly, the obtained amplitude equations for harmonic SC and super-harmonic DC waves are found to be compatible with two well-known one-degree-of-freedom systems: the Duffing and the Helmholtz–Duffing oscillators, respectively.
The paper is organized as follows. The flow configuration and governing equations are introduced in § 2. Section 3 is dedicated to a brief summary of the salient points of the asymptotic model proposed by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), the limitations of which motivated the present work. After tackling the more common case of harmonic SC waves in § 4.1, the WNL amplitude equation governing the super-harmonic DC wave dynamics is derived in § 4.2. Final comments and conclusions are outlined in § 5.
2. Flow configuration and governing equations: potential model
We consider a cylindrical container of diameter 
$D=2R$ filled to a depth 
$h$ with a liquid of density 
$\rho$. The air–liquid surface tension is denoted by 
$\gamma$. The orbital (circular) shaking motion (see sketch in figure 1) can be represented as the combination of two sinusoidal translations with a 
${\rm \pi} /2$ phase shift, thus leading to the following equations of motion for the container axis intersection with the 
$z=0$ plane, parametrized in cylindrical coordinates (
$r$, 
$\theta$):
\begin{equation} \dot{\boldsymbol{X}}_0 = \begin{cases} -\dfrac{d_s}{2}\varOmega_d\sin{\left(\varOmega_d t-\theta\right)}\,\boldsymbol{e}_r,\\ \dfrac{d_s}{2}\varOmega_d\cos{\left(\varOmega_d t-\theta\right)}\,\boldsymbol{e}_{\theta}. \end{cases} \end{equation}
In the classical potential flow limit, i.e. the flow is assumed to be inviscid, irrotational and incompressible, the motion is described in terms of free-surface deformation, 
$\eta$, and a potential velocity field, 
$\varPhi _{tot}$, which is typically separated into a container component, 
$\varPhi _c$, and a fluid component, 
$\varPhi$. Hence, the liquid motion within the moving container is governed by the Laplace equation,
\begin{equation} \Delta\varPhi=\frac{1}{r}\frac{\partial\varPhi}{\partial r}+\frac{\partial^2\varPhi}{\partial r^2}+\frac{1}{r^2}\frac{\partial^2\varPhi}{\partial \theta^2}+\frac{\partial^2\varPhi}{\partial z^2}=0, \end{equation}
subjected to the homogeneous no-penetration condition, 
$\boldsymbol {\nabla }\varPhi \boldsymbol {\cdot }\boldsymbol {n}=\boldsymbol {0}$, at the solid sidewall and bottom, and by the dynamic and kinematic free-surface boundary conditions at 
$z=\eta$ (see Ibrahim Reference Ibrahim2005),
\begin{gather} \frac{\partial\varPhi}{\partial t}+\frac{1}{2}\boldsymbol{\nabla}\varPhi\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi+\eta-\frac{\kappa\left(\eta\right)}{Bo}=rf\cos{\left(\varOmega t-\theta\right)}, \end{gather}
\begin{gather} \frac{\partial \eta}{\partial t}+\frac{\partial\varPhi}{\partial r}\frac{\partial \eta}{\partial r}+\frac{1}{r^2}\frac{\partial\varPhi}{\partial \theta}\frac{\partial \eta}{\partial \theta}-\frac{\partial\varPhi}{\partial z}=0, \end{gather}
which have been made non-dimensional by using the container's characteristic length 
$R$, the characteristic velocity 
$\sqrt {gR}$ and the time scale 
$\sqrt {R/g}$. In (2.3a), 
$\kappa (\eta )$ denotes the fully nonlinear curvature, while 
$Bo=\rho gR^2/\gamma$ is the Bond number. The non-dimensional driving amplitude and angular frequency read 
$f=d_s\varOmega _d^2/(2g)$ and 
$\varOmega =\varOmega _d/\sqrt {g/R}$, respectively. When surface tension is accounted for, an additional contact line boundary condition is required at 
$z=\eta$ and 
$r=1$, typically written as 
$\partial \eta /\partial r=\cot {\vartheta }$, where 
$\vartheta$ is the macroscopic contact angle. Under the classic free-end edge contact line assumption with 
$\vartheta ={\rm \pi} /2$ adopted here, the latter dynamic equation simply reduces to 
$\partial \eta /\partial r=0$. This means that the free surface at rest is flat and that a 
${\rm \pi} /2$ static contact angle is maintained when the contact line elevation changes dynamically.
Sketch of a cylindrical container of diameter 
$D=2R$ and filled to a depth 
$h$. The gravity acceleration is denoted by 
$g$. The Cartesian inertial reference frame is 
$O'\boldsymbol {e}'_x\boldsymbol {e}'_y\boldsymbol {e}'_z$, while 
$O\boldsymbol {e}_x\boldsymbol {e}_y\boldsymbol {e}_z$ is the Cartesian reference frame moving with the container. The origin of the moving cylindrical reference frame 
$(r,\theta,z)$ is placed at the container revolution axis and, specifically, at the unperturbed liquid height, 
$z=0$. The perturbed free surface and contact line elevation are denoted by 
$\eta$ and 
$\delta$, respectively. Parameter 
$d_s$ is the diameter of the circular shaking trajectory, characterized by a driving angular frequency 
$\varOmega _d$.
3. Linear solution and second-order straightforward asymptotic expansion
In order to enlighten the limitations of the expansion procedure developed by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), which motivates the formalization of the new theoretical framework proposed in the present paper, we briefly recall the salient points. Let us consider the following asymptotic expansion for the flow quantities:
\begin{gather} \varPhi=\varPhi_0+\epsilon\varPhi_1+\epsilon^2\varPhi_2+\text{O}\left(\epsilon^3\right), \end{gather}
\begin{gather} \eta=\eta_0+\epsilon\eta_1+\epsilon^2\eta_2+\text{O}\left(\epsilon^3\right), \end{gather}
together with the further assumption of small driving forcing amplitudes of 
$\text {O}(\epsilon )$, i.e. 
$f=\epsilon F$, with 
$\epsilon$ a small parameter 
$\epsilon \ll 1$ and the auxiliary variable 
$F$ of 
$\text {O}(1)$. Solution 
$\boldsymbol {q}_0=\{\varPhi _0,\eta _0\}^{\rm T}$ represents the rest state, which has a potential velocity field null everywhere, 
$\varPhi _0=0$, and a flat interface, 
$\eta _0=0$, as the contact angle is here assumed to be 
$\vartheta ={\rm \pi} /2$. Substituting the expansions above in (2.2)–(2.3b), a series of systems at the various orders in 
$\epsilon$ is obtained. At leading order, (2.2)–(2.3b) reduce to a forced linear system, whose matrix compact form reads
\begin{equation} \left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}\right)\boldsymbol{q}_1=\boldsymbol{\mathcal{F}}_1, \end{equation}
with 
$\boldsymbol {q}_1=\{\varPhi _1,\eta _1\}^{\rm T}$, 
$\boldsymbol {{\mathcal {F}}_1}=F\{0,r/2\}^{\rm T}{\rm e}^{\text {i}(\varOmega t-\theta )}+{\rm c.c.}=F\boldsymbol {\hat {\mathcal {F}}}_1^F{\rm e}^{\text {i}(\varOmega t-\theta )}+{\rm c.c.}$ and
\begin{equation} \boldsymbol{\mathsf{B}} = \begin{pmatrix} 0 & 0 \\ {\rm I}_{\eta} & 0 \end{pmatrix},\quad \boldsymbol{\mathsf{A}} = \begin{pmatrix} \Delta & 0 \\ 0 & -{\rm I}_{\eta}+\dfrac{1}{Bo}\dfrac{\partial\kappa}{\partial\eta} \end{pmatrix}, \end{equation}
where 
${\rm c.c.}$ stands for complex conjugate, 
$\partial \kappa /\partial \eta$ represents the first-order variation of the curvature associated with the small perturbation 
$\epsilon \eta _1$ and 
${\rm I}_{\eta }$ is the identity matrix associated with the interface 
$\eta$. Note that the kinematic condition does not explicitly appear in (3.3a,b), but it is enforced as a boundary condition at the interface (Viola, Brun & Gallaire Reference Viola, Brun and Gallaire2018). In the limit of zero external forcing, i.e. 
$F=0$, system (3.2) is a linear homogeneous problem which, by seeking for solutions having the following normal form:
\begin{equation} \hat{\boldsymbol{q}}_{{mn}}\left(r,z\right)\exp({\text{i}\left(\omega_{mn} t-m\theta\right)})+{\rm c.c.}, \end{equation}reduces to the classic generalized eigenvalue problem for inviscid capillary–gravity waves:
\begin{equation} \left(\text{i}\omega_{mn}\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\hat{\boldsymbol{q}}_{{mn}}=\boldsymbol{0}, \end{equation}
where indices 
$(m,n)$ represent the number of nodal circles and nodal diameters, respectively, with 
$m$ also commonly known as azimuthal wavenumber. Owing to the normal mode expansion, we note that the operator 
$\mathcal {A}$ depends on the azimuthal wavenumber, 
$m$, and, therefore, we denote it by 
$\mathcal {A}_m$. An exact analytical solution to (3.5) can be readily obtained via a Bessel–Fourier series representation leading to the well-known dispersion relation (Lamb Reference Lamb1993)
\begin{equation} \omega_{mn}^2=(k_{mn}+k_{mn}^3/Bo)\tanh{(k_{mn}H)}, \end{equation}
with 
$H=h/R$ and where the wavenumber 
$k_{mn}$ is given by the nth root of the first derivative of the mth-order Bessel function of the first kind satisfying 
${\rm J}'_{m}(k_{mn})=0$.
Despite the existence of this analytical solution, in this work we opt for a numerical scheme based on a discretization technique, where linear operators 
$\mathcal {B}$ and 
$\mathcal {A}_m$ are discretized in space by means of a Chebyshev pseudospectral collocation method with a two-dimensional mapping implemented in Matlab, which is analogous to that described by Viola et al. (Reference Viola, Brun and Gallaire2018). This numerical technique will enable us to avoid straightforward, but cumbersome calculations, otherwise required in the development of the rest of this work and, particularly, of § 4.2. One must note that when (3.5) is solved numerically as in the present case, additional boundary conditions need to be made explicit in order to regularize the problem on the revolution axis (
$r=0$), i.e.
\begin{gather} m=0:\quad \frac{\partial\hat{\eta}_{{mn}}}{\partial r}=\frac{\partial\hat{\varPhi}_{{mn}}}{\partial r}=0, \end{gather}
\begin{gather} |m|\ge1:\quad \hat{\eta}_{{mn}}=\hat{\varPhi}_{{mn}}=0. \end{gather}It is also useful to note that owing to the symmetries of the problem, system (3.5) is invariant under the transformation
\begin{equation} (\hat{\boldsymbol{q}}_{{mn}},+m,\text{i}\omega_{mn})\longrightarrow(\hat{\boldsymbol{q}}_{{mn}},-m,\text{i}\omega_{mn}). \end{equation} Convergence of the numerical solution was checked by computing the first 16 modes (
$m=0,2,3,4$ with 
$n=1,2,3,4$), whose corresponding natural frequency values, 
$\omega _{mn}$, matched the analytical ones given by (3.6) up to the fourth digit for a computational grid 
$N_r=N_z=60$, with 
$N_r$ and 
$N_z$ the number of radial and axial grid points, respectively.
Let us now reintroduce the forcing term on the right-hand side of (3.2). In contradistinction with the cases of unidirectional forcing (Miles Reference Miles1984a,Reference Milesb), for circular orbits, given the azimuthal periodicity of the associated forcing, the shaking at linear order is expected to excite non-axisymmetric modes only and, specifically, those with 
$m=1$. Therefore, the linear response to the external forcing can be sought as
\begin{equation} \boldsymbol{q}_{1}=F\hat{\boldsymbol{q}}_1^F \exp\left({\text{i}\left(\varOmega t-\theta\right)}\right)+{\rm c.c.}, \end{equation}
with 
$\hat {\boldsymbol {q}}_1^F$ being the solution of the following forced problem:
\begin{equation} \left(\text{i}\varOmega\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{1}\right)\hat{\boldsymbol{q}}_1^F=\boldsymbol{\hat{\mathcal{F}}}_1^F. \end{equation}
The response structure 
$\hat {\boldsymbol {q}}_1^F$ is here computed numerically, but, in practice, it is formally equivalent to that obtained analytically by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) by projecting the forcing term 
$\boldsymbol {\hat {\mathcal {F}}}_1$ onto the basis formed by the first-order Bessel functions of the first kind, except that surface tension is retained here because of the finite Bond number. Noting that 
$\epsilon F=f=d_s\varOmega ^2/(2g)$, in figure 2 the linear solution 
$\epsilon \boldsymbol {q}_1^F$ from (3.9) is shown (black solid lines) and compared with experimental measurements reported by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) in terms of maximum non-dimensional crest-to-trough contact line amplitudes, 
$\tilde {\delta }=\delta R/D$, with 
$\delta (\theta,t)=\eta (r=1,\theta,t)$. Measurements for different values of the non-dimensional shaking diameters, 
$\tilde {d}_s=d_s/D$, are shown. Blue and green markers in figure 2 correspond to highly nonlinear scenarios manifesting a free-surface breaking, which is therefore subsequently ignored. As discussed by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) and reproduced here, the linear solution describes well the SC wave dynamics for driving frequencies far enough from harmonic resonances and, particularly, for small 
$\tilde {d}_s$. However, as typical of undamped forced oscillators, the amplitude of the inviscid linear response to the external forcing is proportional to 
$1/|\varOmega ^2-\omega _{1n}^2|$ and therefore it naturally diverges close to 
$\omega _{1n}$, thus failing in predicting the close-to-resonance behaviour, e.g. for 
$\tilde {d}_s=0.02$ at 
$\varOmega \approx \omega _{11}$. Introduction of viscous dissipation would regularize the divergent behaviour at 
$\varOmega =\omega _{11}$; however, in the absence of any nonlinear restoring term, the hardening nonlinearity displayed in figure 2 cannot be retrieved.
Markers correspond to the experimentally measured maximum crest-to-trough contact line amplitude (non-dimensional), with 
$\tilde {\delta }=\delta R/D=\delta /2$, reported by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) for two container diameters, 
$D=0.144\,\text {m}$ and 
$D=0.287\,\text {m}$, a non-dimensional depth 
$\tilde {H}=h/D=0.52$ and five values of 
$\tilde {d}_s=d_s/D$, as a function of the non-dimensional shaking frequency 
$\varOmega$ normalized by the natural frequency of the first non-axisymmetric mode, 
$\omega _{11}=1.3286$ (
$m=1$), on the bottom 
$x$ axis and by that of first non-axisymmetric mode with 
$m=2$, 
$\omega _{21}=1.7475$, on the top 
$x$ axis (the frequency values correspond to 
$D=0.287\,\text {m}$). Colours denote different wave conditions. Black solid lines: linear potential model solution, from (3.9), computed by solving numerically equation (3.10). Red solid lines: WNL solution close to 
$\varOmega \approx \omega _{21}/2$, obtained by computing (3.14). Note that in order to ease the comparison with experiments, the non-dimensional 
$\delta$ was rescaled by a factor 
$R/D=1/2$, as the container diameter 
$D$, rather than the container radius 
$R$, was used by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) to make the equations non-dimensional.
Furthermore, in experiments multiple-crested waves were observed at fractions of the natural frequencies (red markers in figure 2), i.e. the system responses with a frequency which is n times (with n a positive integer) that of the external forcing. Here we refer to such conditions as super-harmonic dynamics (note that the terminology sub-harmonic was used by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) instead). Among these super-harmonics, the DC wave dynamics, occurring at a driving frequency 
$\varOmega \approx \omega _{21}/2$, was seen to be the most relevant (see figure 2), i.e. the most stable and the one displaying the largest deviation from the linear approximation. This specific multiple-crest dynamics, which is intrinsically nonlinear, is indeed favoured by the azimuthal symmetry of the external forcing. Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) tentatively described such a DC dynamics by pursuing the asymptotic analysis up to the second order in 
$\epsilon$, as in (3.1a) and (3.1b), in order to account for second-order system weak nonlinearities.
At the second order in 
$\epsilon$, one obtains the following forced linear system:
\begin{equation} (\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}})\boldsymbol{q}_2=\boldsymbol{\mathcal{F}}_2=F^2(\boldsymbol{\hat{\mathcal{F}}}_2^{FF} \exp({\text{i}(2\varOmega t-2\theta)})+{\rm c.c.})+F^2\boldsymbol{\hat{\mathcal{F}}}_2^{F\bar{F}}, \end{equation}
where 
$\boldsymbol {{\mathcal {F}}_2}$ gathers a series of terms produced by the first-order solution through the second-order system nonlinearities. For the sake of brevity, the explicit expression of these forcing terms is here omitted (see Appendix D for details). The overbar denotes the complex conjugate. Also note that amplitude 
$F$ is actually a real quantity and in the following the superscript 
${\bar {F}}$ is used only to indicate forcing terms produced by the combination of the direct and complex conjugate contributions of the first-order response to the external forcing. The right-hand side of (3.11) clearly shows how second-order terms naturally induce a super-harmonic response, whose spatial periodicity is 
$m=2$, hence precisely corresponding to the DC dynamics experimentally observed. The second forcing term on the right-hand side of (3.11) has 
$\omega =0$ and 
$m=0$, i.e. it is steady and axisymmetric. It originates in the leading-order contribution in time- and azimuthal-averaged flow, the so-called mean flow. Equation (3.11) was solved analytically by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) by retaining for convenience only two modes, namely those with 
$(m,n)=(2,1)$ and 
$(0,1)$, expected to be the relevant ones. The numerical scheme employed in this work allows us to effortlessly account for all the 
$(2,n)$ and 
$(0,n)$ modes simultaneously, as their contribution is directly encompassed in the spatial function 
$\hat {\boldsymbol {q}}_2^{FF}$ and 
$\hat {\boldsymbol {q}}_2^{F\bar {F}}$, appearing in the second-order solution,
\begin{equation} \boldsymbol{q}_{2}=(F^2\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(2\varOmega t-2\theta\right)})+{\rm c.c.})+F^2\hat{\boldsymbol{q}}_2^{F\bar{F}}, \end{equation}whose contributions are computed by solving the following systems:
\begin{equation} \left(\text{i}2\varOmega\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{2}\right)\hat{\boldsymbol{q}}_2^{FF}=\hat{\boldsymbol{\mathcal{F}}}_2^{FF},\quad -\boldsymbol{\mathsf{A}}_0\hat{\boldsymbol{q}}_2^{F\bar{F}}=\hat{\boldsymbol{\mathcal{F}}}_2^{F\bar{F}}. \end{equation}The total flow field, obtained through the asymptotic model, is then given by the sum of the first- and second-order solutions:
\begin{equation} \boldsymbol{q}=(f\hat{\boldsymbol{q}}_1^F \exp({\text{i}\left(\varOmega t-\theta\right)})+f^2\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(2\varOmega t-2\theta\right)})+{\rm c.c.})+f^2\hat{\boldsymbol{q}}_2^{F\bar{F}}, \end{equation}
where, in order to eliminate the implicit small parameter 
$\epsilon$, the amplitudes 
$\epsilon F$ and 
$\epsilon ^2 F^2$ are recast in terms of the physical amplitudes 
$f$ and 
$f^2$, respectively. The resulting prediction of the maximum crest-to-trough contact line amplitude 
$\delta (\theta,t)=\eta (r=1,\theta,t)$ is shown in figure 2 for driving frequencies close to 
$\varOmega /\omega _{21}\approx 0.5$ (see top 
$x$ axis) as red solid lines. Although this straightforward asymptotic expansion detects the emergence of the super-harmonic DC wave in that frequency window, it completely fails in capturing the correct nonlinear wave amplitude saturation displaying a hardening behaviour clearly visible in figure 2. Once again, the amplitude of the inviscid second-harmonic response is proportional to 
$1/|\varOmega ^2-\omega _{2n}^2/4|$ and the total solution tends to diverge close to the DC super-harmonic at 
$\omega _{21}/2$.
Such a symmetric and, in the absence of dissipation, close-to-resonance divergent behaviour is actually expected when performing straightforward asymptotic expansions as they typically suffer from secular (or resonating) terms that must be properly treated (see Castaing (Reference Castaing2005) and Nayfeh (Reference Nayfeh2008) among many other references).
4. Weakly nonlinear analysis via multiple-time-scale method
In order to overcome the aforementioned limitations of the straightforward asymptotic expansion procedure and thus to attempt to bridge the gap between theoretical predictions and experimental observations, we conduct in this section a WNL analysis based on the multiple-time-scale method. With the aim of deriving a WNL amplitude equation governing the DC dynamics, we first tackle the simpler problem of SC waves. In both cases we look for a third-order asymptotic solution of the system
\begin{equation} \boldsymbol{q}=\left\{\varPhi,\eta\right\}^{\rm T}=\epsilon\boldsymbol{q}_1+\epsilon^2\boldsymbol{q}_2+\epsilon^3\boldsymbol{q}_3+\text{O}(\epsilon^4), \end{equation}
where the zeroth-order solution, 
$\boldsymbol {q}_0=\boldsymbol {0}$, is omitted.
4.1. Single-crest dynamics
In § 3 the forcing amplitude 
$f$ was assumed of order 
$\epsilon$, thus leading to a linear first-order problem directly forced by external shaking, which produces a divergent response close to harmonic resonances. With regards to SC waves and specifically to the harmonic response at a driving frequency close to that of one of the non-axisymmetric modes, 
$\omega _{1n}$, we assume here a small forcing amplitude of order 
$\epsilon ^3$. This assumption is justified by the fact that close to resonance, 
$\varOmega \approx \omega _{1n}$, and in the absence of dissipation, even a small forcing will induce a large system response. The following analysis is therefore expected to hold for 
$\varOmega =\omega _{1n}+\lambda$, where 
$\lambda$ is a small detuning parameter assumed of order 
$\epsilon ^2$. Lastly, in the spirit of the multiple-scale technique, we introduce the slow time scale 
$T_2 = \epsilon ^2t$, with 
$t$ being the fast time scale at which the free surface oscillates with angular frequency 
$\approx \omega _{1n}$. Hence, the following scalings are assumed:
\begin{equation} f=\epsilon^3 F,\quad \lambda=\varOmega-\omega_{1n}=\epsilon^2\varLambda,\quad T_2=\epsilon^2 t, \end{equation}
with 
$F$ and 
$\varLambda$ of 
$\text {O}(1)$. We note that the forcing amplitude could be assumed of order 
$\epsilon ^2$ (as the other parameters); however, this complicates unnecessarily the second-order problem without modifying the final amplitude equation, even if the values of its coefficients end up being slightly, up to order 
$\epsilon$, different.
Although the asymptotic expansion is here pursued up to the third order in 
$\epsilon$, the procedure of the WNL analysis is essentially equivalent to that of the straightforward asymptotic model discussed in § 3. The major difference lies in the solution form of the leading-order problem that is now a homogeneous problem, as in (3.5). Given the azimuthal periodicity of the external forcing, among all possible natural eigenmodes we assume a leading-order solution as
\begin{equation} \boldsymbol{q}_1=A_1\left(T_2\right)\hat{\boldsymbol{q}}_{1}^{A_1} \exp({\text{i}\left(\omega_{1n} t-\theta\right)})+{\rm c.c.}, \end{equation}
where 
$\hat {\boldsymbol {q}}_1^{A_1}$ is the eigenmode (computed by solving (3.5)) associated with 
$(m,n)=(1,n)$ and 
$\omega _{1n}$ is the corresponding natural frequency (solution of (3.6)).
The complex amplitude 
$A_1$, a function of the slow time scale 
$T_2$ and still unknown at this stage of the expansion, describes the slow-time amplitude modulation of the oscillating wave 
$\hat {\boldsymbol {q}}_1^{A_1}$ and introduces a new arbitrariness in the problem, which must be fixed at a higher order. Eigensurface 
$\hat {\eta }_1^{A_1}$ and eigenpotential field 
$\hat {\varPhi }_1^{A_1}$, computed for 
$\omega _{1n}=\omega _{11}$, are shown in figures 3(a) and 3(b), respectively.
(a–c) Real part of the first-order 
$\hat {\eta }_1^{A_1}$ and second-order 
$\hat {\eta }_2^{A_1A_1}$ and 
$\hat {\eta }_2^{A_1\overline {A_1}}$ free-surface deformations computed for 
$\omega _{1n}=\omega _{11}$. (d–f) Imaginary part of the associated first-order 
$\hat {\varPhi }_1^{A_1}$ and second-order 
$\hat {\varPhi }_2^{A_1A_1}$ and 
$\hat {\varPhi }_2^{A_1\overline {A_1}}$ potential velocity field. Each response is denoted by its amplitude dependence, 
$\epsilon A_1$, 
$\epsilon ^2 A_1\bar {A}_1$ and 
$\epsilon ^2A_1A_1$. The first-order eigenmode is normalized with the amplitude and phase of the contact line (at 
$r=1$), such that the free surface 
$\hat {\eta }_1^{A_1}$ is purely real, whereas 
$\hat {\varPhi }_1^{A_1}$ is purely imaginary. Note that the second-order mean flow constantly induces an upside-down bell-like axisymmetric interface deformation pushing the free surface downward at the centre of the moving container. Calculations are performed for the case of figure 2, i.e. pure water with 
$\rho =1000\,\text {m}\,\text {m}^{-3}$, 
$\gamma =0.072\,\text {N}\,\text {m}^{-1}$, 
$D=0.287\,\text {m}$ and 
$\tilde {H}=h/D=0.52$, for which 
$Bo=2802.8$ and 
$\omega _{11}=1.3286$.
By pursuing the expansion to the second order, a linear system forced by the first-order solution and analogous to that of (3.11) is obtained (see Reclari (Reference Reclari2013) for the full expansion of the original nonlinear governing equation up to the second order). Nevertheless, the forcing terms on the right-hand side are here proportional to 
$A_1^2$ (super- or second-harmonic) and to 
$A_1\bar {A}_1$ (mean flow). Thus, we seek for a second-order solution of the form
\begin{equation} \boldsymbol{q}_2=A_1\bar{A}_1\hat{\boldsymbol{q}}_2^{A_1\bar{A}_1}+(A_1^2\hat{\boldsymbol{q}}_{2}^{A_1A_1} \exp({\text{i}\left(2\omega_{1n} t-2\theta\right)})+{\rm c.c.}), \end{equation}
with 
$\hat {\boldsymbol {q}}_2^{A_1\bar {A}_1}$ and 
$\hat {\boldsymbol {q}}_{2}^{A_1A_1}$ computed numerically and displayed in figures 3(b,d) and 3(c, f), respectively, in terms of second-order free-surface deformations and potential velocity fields evaluated for 
$\omega _{1n}=\omega _{11}$. From a numerical perspective, we note that the second-order responses can be straightforwardly computed as long as the pairs 
$(\omega,m)=(2\omega _{1n},2)$ and 
$(0,0)$ do not correspond to eigenvalues of (3.5), i.e. the second-order operators 
$(\text {i}2\omega _{1n}\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{2})$ and 
$-\boldsymbol{\mathsf{A}}_0$ are non-singular and hence invertible.
With regards to figure 3, it is interesting to note how the second-order mean flow potential velocity field is null everywhere. This can be rigorously proven by first noticing that the mean flow corresponds to a time- and azimuthal-averaged flow, i.e. 
$\partial /\partial t=\partial /\partial \theta =0$. Moreover, in the inviscid limit, free-surface elevation and potential field have a 
${\rm \pi} /2$ phase shift, meaning that the first-order eigenmode can be normalized such that the eigensurface is purely real, whereas the eigenpotential is purely imaginary. Under these conditions, the mean flow forcing term on the right-hand side of the kinematic equation cancels out, so that the associated Laplace equation appears to be constrained by homogeneous Neumann conditions on all the domain boundaries, thus prescribing a trivial constant potential field and therefore a null velocity field. In other words, the second-order mean flow system reduces to forced linear meniscus equation (resulting from (2.3a)) and its conditions at 
$r=0$ and 
$r=1$ (both 
$\partial \hat {\eta }_2^{A_1\bar {A}_1}/\partial r=0$), which prescribes a static mean interface deformation only. Such a result was expected since the second-order mean flow response represents the Eulerian mean flow, which, together with the so-called Stokes drift, contribute to the overall Lagrangian mean flow (see van den Bremer & Breivik (Reference van den Bremer and Breivik2018) for a thorough review).
While the Stokes drift is a pure kinematic concept, the Eulerian mean flow, often referred to as streaming flow (Bouvard et al. Reference Bouvard, Herreman and Moisy2017), is generally believed to be of viscous origin, although another appealing interpretation has been recently proposed (Faltinsen & Timokha Reference Faltinsen and Timokha2019). Sticking to the well-accepted viscous Eulerian mean flow generation mechanism, it is not a surprise that the absence of viscous boundary layers results in a vanishing Eulerian mean flow.
We now move forward to the 
$\epsilon ^3$-order problem, which is once again a linear problem forced by combinations of the first- and second-order solutions as well as by the slow time derivative of the leading-order solution and by the external forcing, which was assumed of order 
$\epsilon ^3$:
\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_3=\mathcal{F}_3 ={-}\partial_{T_2} A_1\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_1} \exp({\text{i}\left(\omega_{1n} t-\theta\right)}) + |A_1|^2A_1 \boldsymbol{\hat{\mathcal{F}}}_3^{A_1\bar{A}_1A_1} \exp({\text{i}\left(\omega_{1n} t-\theta\right)})\nonumber\\ &\quad +F\boldsymbol{\hat{\mathcal{F}}}_3^F \exp({\text{i}{\varLambda} T_2})\exp({\text{i}\left(\omega_{1n} t-\theta\right)}) +\text{NRT}+{\rm c.c.}, \end{align}
with 
$\boldsymbol {\hat {\mathcal {F}}}_3^F=\{0,r/2\}^{\rm T}$ and where 
$\text {NRT}$ stands for non-resonating terms, which are not relevant for further analysis. As standard in multiple-scale analysis, the indeterminacy introduced by the unknown amplitude 
$A_1$ is resolved by requiring that secular terms do not appear in the solution to (4.5). Secularity results from all resonant forcing terms in 
$\mathcal {F}_3$, i.e. all terms sharing the same frequency and wavenumber 
$(\omega _{1n},1)$ of 
$\boldsymbol {q}_1$, and in effect all terms explicitly written in (4.5). It follows that a compatibility condition must be enforced through the Fredholm alternative (Friedrichs Reference Friedrichs2012). Such a compatibility condition imposes the amplitude 
$B=\epsilon A_1 {\rm e}^{\text {i}{\lambda } t}$ to obey the following normal form:
\begin{equation} \frac{{\rm d}B}{{\rm d}t}={-}\text{i}{\lambda} B + \text{i}\,\mu_{{SC}} f + \text{i}\,\nu_{{SC}} |B|^2B, \end{equation}
where the physical time 
$t=T_2/\epsilon ^2$ has been reintroduced and where forcing amplitude and detuning parameter are recast in terms of their corresponding physical value, 
$f=\epsilon ^3 F$ and 
$\lambda =\epsilon ^2\varLambda$. Moreover, the small implicit parameter 
$\epsilon$ is eliminated by defining the total physical amplitude 
$A=\epsilon A_1$ (Bongarzone et al. Reference Bongarzone, Bertsch, Renaud and Gallaire2021a). The various normal form coefficients, which turn out to be real-valued quantities owing to the absence of dissipation, are computed as scalar products between the adjoint mode, 
$\hat {\boldsymbol {q}}_1^{A_1 {\dagger} }$, associated with 
$\hat {\boldsymbol {q}}_1^{A_1}$, and the third-order resonant forcing terms as follows:
\begin{gather} \text{i}\,\mu_{{SC}}=\frac{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\boldsymbol{\mathsf{B}}\hat{\mathcal{F}}_3^F\rangle}{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\mathcal{B}\hat{\boldsymbol{q}}_1^{A_1}\rangle}=\frac{\int_{z=0}^{}r\bar{\hat{\eta}}_1^{A_1 {\dagger}}/2\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_1 {\dagger}}\hat{\varPhi}_1^{A_1}+\hat{\varPhi}_1^{A_1 {\dagger}}\hat{\eta}_1^{A_1})\, r\text{d}r}, \end{gather}
\begin{gather} \text{i}\,\nu_{{SC}}=\frac{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\boldsymbol{\mathsf{B}}\boldsymbol{\hat{\mathcal{F}}}_3^{A_1\bar{A}_1A_1}\rangle}{\langle\hat{\boldsymbol{q}}_1^{A_1 {\dagger}},\mathcal{B}\hat{\boldsymbol{q}}_1^{A_1}\rangle}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_1 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{A_1\bar{A}_1A_1}+\hat{\varPhi}_1^{A_1 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{A_1\bar{A}_1A_1})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_1 {\dagger}}\hat{\varPhi}_1^{A_1}+\hat{\varPhi}_1^{A_1 {\dagger}}\hat{\eta}_1^{A_1})\, r\text{d}r}.\end{gather}
Here 
$\hat {\boldsymbol {q}}_1^{A_1 {\dagger} }=\bar {\hat {\boldsymbol {q}}}_1^{A_1}$, since the inviscid problem is self-adjoint with respect to the Hermitian scalar product 
$\langle \boldsymbol {a},\boldsymbol {b}\rangle =\int _{\varSigma }^{}\bar {\boldsymbol {a}}\boldsymbol {\cdot }\boldsymbol {b}\,\text {d}V$, with 
$\boldsymbol {a}$ and 
$\boldsymbol {b}$ two generic vectors (see Viola et al. (Reference Viola, Brun and Gallaire2018) for a thorough discussion and derivation of the adjoint problem). For the sake of brevity, the explicit expression of 
$\boldsymbol {\hat {\mathcal {F}}}_3^{A_1\bar {A}_1A_1}$ is omitted, as it only involves straightforward calculations, i.e. a Taylor expansion of nonlinear governing equations and boundary conditions (2.2)–(2.3b) around the rest state 
$\boldsymbol {q}_0=\boldsymbol {0}$. Here we simply denote with the subscripts ‘
${dyn}$’ and ‘
${kin}$’ the forcing components appearing in the dynamic and kinematic boundary conditions, respectively.
By turning to polar coordinates, 
$B=|B|{\rm e}^{\text {i}\varTheta }$, splitting the modulus and phase parts of (4.6) and looking for stationary solution, 
${\rm d}/{\rm d}t=0$ with 
$|B|\ne 0$, the following implicit relation is obtained:
\begin{equation} \tilde{d}_s\varOmega^2\mp \frac{\left({\lambda}-\nu_{{SC}}|B|^2\right)|B|}{\mu_{{SC}}}=0, \end{equation}or, in a more common polynomial form,
\begin{equation} P\left(|B|\right)=|B|^3-\frac{{\lambda}}{\nu_{{SC}}}|B|\pm\frac{\mu_{{SC}}\tilde{d}_s\varOmega^2}{\nu_{{SC}}}=0, \end{equation}
where 
$f=\tilde {d}_s\varOmega ^2$, 
$\lambda =(\varOmega -\omega _{1n})$ and the 
$\mp$ signs correspond to the phases 
$\varTheta =0$ and 
${\rm \pi}$, respectively. The two branches prescribed by (4.8) for 
$|B|$ as a function of 
$\varOmega$ at a fixed non-dimensional shaking diameter 
$\tilde {d}_s$ can be easily computed using the Matlab function fimplicit. After evaluating the stable and unstable stationary solutions for 
$|B|$ and 
$\varTheta$, the total SC wave solution is reconstructed as
\begin{equation} \boldsymbol{q}_{SC}=\left\{\varPhi,\eta\right\}^{\rm T}=\boldsymbol{q}_1+\boldsymbol{q}_2. \end{equation}4.1.1. Experiments versus WNL prediction: wave amplitude
In figure 4 the WNL prediction in terms of maximum crest-to-trough contact line amplitude, 
$\Delta \tilde {\delta }$, for SC waves is compared with two sets of experimental measurements and with the potential linear solution (3.9). In comparison with the linear theory presented in § 3, the agreement with experiments improves for different shaking diameters and for different harmonic resonances, e.g. those associated with modes 
$(m,n)=(1,1)$ and 
$(1,2)$ of figure 4(a). The hardening nonlinearity is correctly captured and the amplitude prediction matches well the measurements until the free surface eventually breaks and the wave regime leaves the WNL regime, hence suggesting the little relevance of dissipative effects attributable to viscosity in this regime.
Comparison of the WNL prediction for SC waves with experiments in terms of maximum crest-to-trough contact line amplitude (non-dimensional), 
$\Delta \tilde {\delta }$. (a) Black pentagons correspond to the experimental measurements presented in figure 4.30 of Reclari (Reference Reclari2013) (R13) for 
$\tilde {d}_s=d_s/D=0.01$, where the first two non-axisymmetric modes 
$(m,n)=(1,1)$ and 
$(1,2)$ are detected. Dotted black lines: solution of the linear potential model according to (3.9). Light blue and blue lines: WNL SC prediction (4.10). Unstable branches are represented as a dashed line. Violet solid lines: theoretical prediction by Raynovskyy & Timokha (Reference Raynovskyy and Timokha2018a) (R&T18) (see § 4.1.2 for further comments). (b) Same as (a) with the black filled circles corresponding to the measurements of Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) (R14) reported in figure 2 for 
$\tilde {d}_s=0.02$. The values of normal form coefficients 
$\mu _{SC}$ and 
$\nu _{SC}$ computed for 
$\varOmega \approx \omega _{11}$ and 
$\omega _{12}$ (see bottom and top 
$x$ axes) are given in table 1, together with the corresponding values of 
$\tilde {H}$, 
$D$ and natural frequencies 
$\omega _{mn}$ used in this calculation.
Values of the amplitude equation coefficients 
$\mu _{SC}$ and 
$\nu _{SC}$ used to produce figure 4.
However, one must note that in this WNL approach the driving frequency is essentially fixed around that of a unique non-axisymmetric natural mode, 
$\varOmega \approx \omega _{1n}$. Consequently, when performing the analysis for a mode 
$(1,n)$, the influence of all other modes is completely overlooked. In consequence, the accuracy of the asymptotic solution rapidly deteriorates moving away from harmonic resonances, when compared with the linear solution (3.9), which turns out to be more accurate. This is visible looking at the bottom stable branch in the multi-solution range of figure 4(b) or by looking at the driving frequency window 
$\varOmega \in [0.7\omega _{12},0.9\omega _{12}]$ in figure 4(a). In other words, the detuning parameter should be small in order for the present WNL analysis, based on a single-mode expansion, to hold. On this regard, as no other natural frequencies are encountered for 
$\varOmega <\omega _{11}$, an exception is made for the left-hand branch associated with the harmonic resonance of the first (or fundamental) non-axisymmetric mode, where an excellent agreement of the single-mode prediction, comparable to that of the linear solution, lasts until 
$\varOmega \approx 0$.
4.1.2. Comparison with the multimodal theory of Raynovskyy & Timokha (Reference Raynovskyy and Timokha2018a)
The violet solid curves reported in figure 4 correspond to the predictions associated with the 
$\omega _{11}$ SC swirling from the Narimanov–Moiseev multimodal sloshing theory employed by Raynovskyy & Timokha (Reference Raynovskyy and Timokha2018a) (R&T18) (only the stable branches are reported). Their curves have been here carefully reproduced by manually sampling those reported in figure 3 of R&T18 in the range of frequency available, i.e. 
$\varOmega /\omega _{11}\in [0.8,1.3]$. By looking at the upper stable branch, although an increasing discrepancy is observed for increasing wave amplitudes, one can see that both analyses are in fairly good agreement with experiments and with each other until wave breaking eventually occurs. Such a discrepancy could be tentatively attributed to the different definition of the detuning parameter employed in R&T18, i.e. 
$\epsilon ^2\varLambda _{R\&T18}=\omega _{1n}^2/\varOmega ^2-1$. On the other hand, by looking at the lower stable branch, one sees that, whereas the jump-up frequencies according to R&T18 and to the present model essentially coincide, the discrepancy between the two predictions increases at larger frequency, i.e. 
$\varOmega /\omega _{11}>1$, with that of R&T18 being closer to the linear potential model. One should also note that, in contradistinction with our analysis, that of R&T18 accounts for damping and predicts the jump-down transition visible in figure 4(a). This damping was essentially fitted from the experimental measurements and, specifically, from the jump-down frequency occurring at a larger frequency, once the wave-breaking regime is fully established, i.e. 
$\varOmega /\omega _{11}=1.27$ for 
$\tilde {d}_s=0.01$ and 
$\varOmega /\omega _{11}=1.45$ for 
$\tilde {d}_s=0.02$ (see figure 4). However, experiments suggest that the damping effect on the curves displayed in figure 4 would not be easy to observe, even for 
$\tilde {d}_s=0.01$. Indeed, the motion undergoes a SC wave breaking, thus entering in a fully nonlinear regime, where both our analysis and that of R&T18 loose predictive power. We therefore decided to discard damping while comparing our results with the close-to-harmonic resonance experiments from Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014).
4.1.3. The Duffing oscillator analogy
Mass–spring models are widely employed in several engineering fields, e.g. in aerospace engineering, for the description of close-to-resonance sloshing motions (Moiseev Reference Moiseev1958; Bauer Reference Bauer1966; Dodge Reference Dodge2000), where nonlinearities are of crucial importance. The most popular driven nonlinear mass–spring model is that developed by Duffing (Reference Duffing1918), who added a cubic nonlinear spring deformation (cubic term) to the classical driven harmonic oscillator:
\begin{equation} \ddot{x}+2\sigma\dot{x}+ x+c_3 x^3 = p\cos{\varOmega t}, \end{equation}
where 
$\sigma$ is the damping coefficient and where, depending on the sign of 
$c_3$, the resonance curve bends and the nonlinear resonance frequency shifts, i.e. it decreases for a softening spring (
$c_3 < 0$), whereas it increases for a hardening spring (
$c_3 > 0$), thus explaining the original observation of Duffing on vibration mechanism. Ockendon & Ockendon (Reference Ockendon and Ockendon1973) showed via asymptotic expansion of the potential flow solution in the neighbourhood of a harmonic resonance that for small external forcing amplitudes, sloshing in a two-dimensional rectangular container responds exactly as an undamped Duffing oscillator (with 
$\sigma =0$). In Appendix B, we briefly show that, as expected, the same holds for close-to-harmonic-resonance sloshing in orbital-shaken cylindrical containers, whose formal amplitude equation, starting from the full inviscid hydrodynamic system, was derived in § 4.1 (see (4.6)). Typically, when the Duffing equation is employed to model close-to-resonance responses in sloshing dynamics and experimental measurements are available, the nonlinear coefficient is often computed by fitting the experimental measurements. Recently, with regards to quasi-two-dimensional rectangular containers laterally excited, Bäuerlein & Avila (Reference Bäuerlein and Avila2021) have carried out careful quantitative comparisons between experiments and theoretical predictions from the damped Duffing equation, showing that their actual sloshing system is remarkably well described by the forced-damped Duffing oscillator. Nevertheless, for increasing wave amplitude responses, experiments deviate from the Duffing solution, which is not capable of predicting correctly the phase lag between driving and response, shown to be the key factor for an accurate estimation of the sloshing amplitude of the maximal nonlinear resonance (Cenedese & Haller Reference Cenedese and Haller2020; Bäuerlein & Avila Reference Bäuerlein and Avila2021). We note that, by analogy with the undamped Duffing equation, the WNL analysis formalized in § 4.1 exacerbates this aspect, since, owing to the lack of dissipation, it can only predict the classic phase-lag bounds, 
$0$ and 
${\rm \pi}$ (see Appendix A for further comments in this regard). Nevertheless, one should notice that this intrinsic limitation turns out to be unimportant in cases as those of figure 4, where for increasing amplitude a wave breaking eventually occurs and the WNL theory as well as the Duffing mechanical analogy no longer apply.
4.2. Double-crest dynamics
We now tackle the DC wave response. Its formalization is slightly more subtle, as it requires a new reordering of the small control parameter magnitudes as well as an unusual form of the leading-order problem, involving both a homogeneous and a particular solution. We recall that the DC dynamics in figure 2 occurs at a driving frequency 
$\varOmega \approx \omega _{21}/2$. Results at the end of this section are presented for mode 
$(2,1)$, for which experiments are available; however, for the sake of generality, we formalize the analysis for any mode 
$(2,n)$, i.e. 
$\varOmega =\omega _{2n}/2+\lambda$, where 
$\lambda$ is the small detuning parameter.
4.2.1. Formalism
To determine a suitable scaling for the forcing amplitude 
$f$ and detuning parameter 
$\lambda$ it is instructive to look at the experimental measurements shown in figure 2 for 
$\varOmega$ close to 
$\omega _{21}/2$. One can see that approaching 
$\varOmega \approx \omega _{21}/2$ from lower frequencies, the DC wave emerges on the top of a SC dynamics, with the latter being correctly described by the linear solution, which still behaves well as 
$\omega _{21}/2$ is far enough from the primary harmonic resonance occurring at 
$\omega _{11}$. It follows that the forcing amplitude 
$f$ and detuning 
$\lambda$ could be retained at order 
$\epsilon$ and the first-order problem takes the form (3.2), with 
$\boldsymbol {\mathcal {F}}_1=F\{0,r/2\}^{\rm T}{\rm e}^{\text {i}(\varOmega t-\theta )}+{\rm c.c.}=F\boldsymbol {\hat {\mathcal {F}}}_1^F{\rm e}^{\text {i}(\varOmega t-\theta )}+{\rm c.c.}$, with 
$f=\epsilon F$ and 
$\lambda =\epsilon \varLambda$.
Furthermore, in § 3 we have shown how, close enough to the super-harmonic resonance, the divergent behaviour is produced by a second-order resonating term, which breaks the straightforward expansion, as 
$\epsilon ^2$-order terms should not become larger than the 
$\epsilon$-order ones. In the following, this asymptotic breakdown is overcome by assuming that the leading-order solution is given by the sum of (i) a particular solution, given by the linear response to the external forcing computed by solving (3.10), and (ii) a homogeneous solution, represented by the natural mode 
$(m,n)=(2,n)$, obtained by solving the generalized eigenvalue problem (3.5), up to an amplitude to be determined at higher orders. The second-order resonating term will then require, in the spirit of multiple-time-scale analysis, an additional second-order solvability condition, complementing the third-order non-resonance condition already obtained in the SC wave WNL model. This suggests that two slow time scales exist, namely 
$T_1$ and 
$T_2$, with 
$T_1$ one 
$\epsilon$-order faster than 
$T_2$, hence implying that quadratic nonlinearities are stronger than cubic ones. To summarize, the fundamental scalings underpinning the WNL multiple-scale expansion for DC waves are the following:
\begin{equation} f=\epsilon F,\quad {\lambda=\varOmega-\omega_{2n}/2=\epsilon\varLambda},\quad T_1=\epsilon t,\quad T_2=\epsilon^2 t. \end{equation}The first-order solution reads
\begin{equation} \boldsymbol{q}_{1}=A_2\left(T_1,T_2\right)\hat{\boldsymbol{q}}_1^{A_2} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})+F\hat{\boldsymbol{q}}_1^F \exp({\text{i}\left(\left(\omega_{2n}/2\right) t-\theta\right)}) {{\rm e}^{\text{i}\varLambda T_1}}+{\rm c.c.}, \end{equation}
where the unknown slow-time amplitude modulation, 
$A_2$, is here a function of the two time scales 
$T_1$ and 
$T_2$, while the amplitude of the particular solution simply equals the forcing amplitude and 
$\hat {\boldsymbol {q}}_1^F$ is computed from (3.10) for 
$\varOmega =\omega _{2n}/2$.
The second-order linearized forced problem reads
\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_2={\boldsymbol{\mathcal{F}}_2}={\boldsymbol{\mathcal{F}}_2}^{ij}-\frac{\partial A_2}{\partial T_1}\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2}\exp({\text{i}\left(\omega_{2n} t-2\theta\right)}) \nonumber\\ &\qquad - \text{i}{\varLambda} F\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{{F}} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}) {\exp({\text{i}\varLambda T_1})}+{\rm c.c.} \end{align}
The first-order solution is made up of four different contributions of amplitude 
$A_2$, 
$\bar {A}_2$, 
$F$ and 
$\bar {F}$; therefore it generates 10 different second-order forcing terms, here denoted by 
$\boldsymbol {\mathcal {F}}_2^{ij}$, which exhibit a certain frequency and azimuthal periodicity 
$(\breve {\omega },\breve {m})$.
The additional two forcing terms stem from the time derivative of the first-order solution (4.13) with respect to the first-order slow time scale 
$T_1$. In order to interpret the last term in (4.14), it is worth first noting that, while the amplitude of the linear solution (3.9), computed at a generic driving frequency, grows with 
$\varOmega$ as 
$F/|\varOmega ^2-\omega _{11}^2|=\tilde {d}_s\varOmega ^2/|\varOmega ^2-\omega _{11}^2|\propto \varOmega ^2/|\varOmega ^2-\omega _{11}^2|$, in the WNL model for DC waves, the amplitude of the particular solution (4.13) is proportional to 
$F/|\omega _{21}^2/4-\omega _{11}^2|=\tilde {d}_s\varOmega ^2/|\omega _{21}^2/4-\omega _{11}^2|\propto \varOmega ^2$, since the driving frequency was frozen at 
$\varOmega =\omega _{21}/2+\lambda$, with the small detuning parameter 
$\lambda$ contributing to modify its phase, but not its amplitude. This leads to an increasing discrepancy between (3.9) and the leading-order particular solution (4.13) away from the super-harmonic resonance. The response to the forcing term proportional to 
$\varLambda F$ in (4.14) can be then interpreted as a second-order correction of the amplitude of the first-order particular solution accounting for a detuning from the exact resonance through 
$\varLambda F\propto \tilde {d}_s\varOmega ^2(\varOmega -\omega _{2n}/2)$ and contributing to improve the asymptotic approximation in a wider range of driving frequency in the neighbourhood of the super-harmonic frequency.
None of the forcing terms in (4.14) is resonant, as their oscillation frequency and azimuthal wavenumber differ from those of the leading-order homogeneous solution, except the term produced by the second harmonic of the leading-order particular solution, i.e. 
$\boldsymbol {\mathcal {F}}_2^{FF}=F^2\boldsymbol {\hat {\mathcal {F}}}_2^{FF} \exp ({\text {i}(\omega _{2n}-2\theta )}){\exp ({\text {i}2\varLambda T_1})}+{\rm c.c.}$ To avoid secular terms, a second-order compatibility condition is imposed, requiring that the following normal form is verified:
\begin{equation} \frac{\partial A_2}{\partial T_1}=\text{i}\,\mu_{{DC}} F^2{\exp({\text{i}2\varLambda T_1})}, \end{equation}
with 
$\mu _{{DC}}$ computed as before, i.e.
\begin{equation} \text{i}\,\mu_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{2_{{dyn}}}^{FF}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{2_{{kin}}}^{FF})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}. \end{equation}
Taken alone, the dynamics resulting from (4.15) is, however, of little relevance, since the solution, i.e. the frequency–response curve, would still diverge (symmetrically) to infinity for 
$\Delta =\varOmega -\omega _{2n}/2\rightarrow 0$ in the absence of any restoring term, i.e. the nonlinear mechanism responsible for the finite-amplitude saturation, which only comes into play at order 
$\epsilon ^3$. This means that the expansion must be pursued up to the following order, and thereby that we must solve for the second-order solution.
By substituting (4.15) in the forcing expression, (4.14) can be rearranged as follows:
\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_2= \boldsymbol{\mathcal{F}}_{2_{NRT}}^{i,j}+\boldsymbol{\mathcal{F}}_{2_{RT}}^{i,j}\nonumber\\ &\quad =\boldsymbol{\mathcal{F}}_{2_{NRT}}^{i,j}+F^2(\boldsymbol{\hat{\mathcal{F}}}_2^{FF}-\text{i}\,\mu_{{DC}} \boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2}) \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})\,{\exp({\text{i}2\varLambda T_1})}+{\rm c.c.}, \end{align}
where the subscripts ‘
${NRT}$’ and ‘
${RT}$’ denote non-resonating (whose frequencies and azimuthal periodicities are gathered in table 2) and resonating terms, respectively. Note that the term proportional to 
${\varLambda } F$ has been included in the non-resonating forcing terms, whereas the resonant term is written explicitly. The compatibility condition is now trivially satisfied, meaning that the new forcing term is orthogonal to the adjoint mode, 
$\hat {\boldsymbol {q}}_1^{A_2{\dagger} }=\bar {\hat {\boldsymbol {q}}}_1^{A_2}$, by construction and therefore, according to the Fredholm alternative, a non-trivial solution exists. Hence, we seek for a second-order solution having the following form:
\begin{align} \boldsymbol{q}_2&=A_2\bar{A}_2\hat{\boldsymbol{q}}_2^{A_2{\bar{A}_2}}+F^2\hat{\boldsymbol{q}}_2^{F\bar{F}}\nonumber\\ &\quad +A_2^2\hat{\boldsymbol{q}}_2^{A_2A_2} \exp({\text{i}\left(2\omega_{2n}t-4\theta\right)})+{\varLambda} F\hat{\boldsymbol{q}}_2^{{\varLambda} F} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}){\exp({\text{i}\varLambda T_1})}\nonumber\\ &\quad +A_2F\hat{\boldsymbol{q}}_2^{A_2F}\exp({\text{i}\left(\left(3\omega_{2n}/2\right)t-3\theta\right)}){\exp({\text{i}\varLambda T_1})}\nonumber\\ &\quad +A_2\bar{F}\hat{\boldsymbol{q}}_2^{A_2\bar{F}} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}){\exp({-\text{i}\varLambda T_1})}\nonumber\\ &\quad +F^2\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(\omega_{2n}t-2\theta\right)}){\exp({\text{i}2\varLambda T_1})}+{\rm c.c.} \end{align}
All non-resonant responses in (4.18) are handled similarly, i.e. they are computed in Matlab by performing a simple matrix inversion using standard LU solvers (as in §§ 3 and 4.1). As anticipated above, although the operator associated with the resonant forcing term, i.e. 
$(\text {i}\omega _{2n}\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_{2})$, is singular, the value of the normal form coefficient (4.16) ensures that a non-trivial solution for 
$\hat {\boldsymbol {q}}_2^{FF}$ exists. Diverse approaches can be followed to compute this response. Here such a response is computed by using the pseudoinverse matrix of the singular operator (Orchini, Rigas & Juniper Reference Orchini, Rigas and Juniper2016). Another possible approach is given in Appendix A of Meliga, Gallaire & Chomaz (Reference Meliga, Gallaire and Chomaz2012), where a two-step regularization procedure, involving an intermediate factious solution for 
$\hat {\boldsymbol {q}}_2^{FF}$, is employed. We also note that in (4.18), exactly as in (4.4), a second-order homogeneous solution has not been accounted for as its introduction would be irrelevant to the final solution.
First-order linear solutions and second-order non-resonating forcing terms gathered by their amplitude dependency and corresponding azimuthal and temporal periodicity 
$(\breve {m},\breve {\omega })$. Six terms have been omitted as they are the complex conjugates of 
$\epsilon A_2$, 
$\epsilon F$, 
$\epsilon ^2 A_2A_2$, 
$\epsilon ^2\varLambda F$, 
$\epsilon ^2A_2F$ and 
$\epsilon ^2 A_2\bar {F}$.
The first-order solutions together with all the non-resonating second-order responses are shown in the various panels of figure 5, where the two leading-order contributions, 
$\epsilon A_2$ and 
$\epsilon F$, corresponding to the DC and SC waves, respectively, can be identified. Moreover, we note that the second-order response proportional to 
$\epsilon ^2{\varLambda } F$ has a spatial structure similar to that of the leading-order response 
$\epsilon F$, as it represents the second-order correction to the latter caused by small frequency shifts of order 
$\epsilon$.
Real part of the first-order (a) 
$\hat {\eta }_1^{A_2}$ and (e) 
$\hat {\eta }_1^{F}$ and non-resonating second-order (b) 
$\hat {\eta }_2^{A_1\bar {A}_2}$, (c) 
$\hat {\eta }_2^{A_2A_2}$, (d) 
$\hat {\xi }_2^{\varLambda _1 F}$, ( f) 
$\hat {\eta }_2^{F\bar {F}}$, (g) 
$\hat {\eta }_2^{A_2F}$ and (h) 
$\hat {\eta }_2^{A_2\bar {F}}$ free-surface deformations computed for 
$\omega _{2n}=\omega _{21}$. The first-order eigenmode is normalized with the amplitude and phase of the contact line (at 
$r=1$), such that the free surface 
$\hat {\eta }_1^{A_2}$ is purely real, whereas 
$\hat {\varPhi }_1^{A_2}$ is purely imaginary. Note that the second-order mean flow 
$\hat {\eta }_2^{F\bar {F}}$ constantly induces an upside-down bell-like axisymmetric interface deformation pushing the free surface downward at the centre of the moving container, by analogy with the effect produced by 
$\hat {\eta }_2^{A_1\bar {A}_1}$ for SC waves, as the two responses are essentially equivalent up to a pre-factor. Here the mean flow 
$\hat {\eta }_2^{A_2\overline {A_2}}$ for DC waves pushes the interface upward at the wall (same as 
$\hat {\eta }_2^{F\bar {F}}$) and, at the same time, downward in an annular region at intermediate radial coordinates, without altering the free-surface elevation at the container revolutions axis.
Lastly, at third order in 
$\epsilon$, the problem reads
\begin{align} &\left(\partial_t\boldsymbol{\mathsf{B}}-\boldsymbol{\mathsf{A}}_m\right)\boldsymbol{q}_3=\boldsymbol{\mathcal{F}}_3\nonumber\\ &\quad ={-}\frac{\partial A_2}{\partial T_2 }\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_1^{A_2} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)}) - \text{i}2{\varLambda} F^2\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_2^{FF} \exp({\text{i}\left(\omega_{2n}-2\theta\right)}){\exp({\text{i}2\varLambda T_1})}\nonumber\\ &\qquad + |A_2|^2A_2 \boldsymbol{\hat{\mathcal{F}}}_3^{A_2\bar{A}_2A_2} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})+A_2F^2\boldsymbol{\hat{\mathcal{F}}}_3^{A_2F\bar{F}} \exp({\text{i}\left(\omega_{2n} t-2\theta\right)})\nonumber\\ &\qquad +{\varLambda} F^2\boldsymbol{\hat{\mathcal{F}}}_3^{\varLambda FF} \exp({\text{i}\left(\omega_{2n}{t}-2\theta\right)}){\exp({\text{i}2\varLambda T_1})}+\text{NRT}+{\rm c.c.}, \end{align}
where the first two forcing terms arise from the time derivative of the first-order solution with respect to the second-order slow time scale 
$T_2$ and from that of the second-order solution with respect to the first-order slow time scale 
$T_1$, respectively (see Appendix D for the full expression of 
$\boldsymbol {\mathcal {F}}_2$ and 
$\boldsymbol {\mathcal {F}}_3$). By noticing that the second and last forcing terms share the same amplitude dependence, i.e. 
${\varLambda } F^2$, they can be recast into a single forcing term, say 
${\varLambda } F^2\hat {\boldsymbol {\mathcal {F}}}_3^{{\varLambda } FF} \exp ({\text {i}(\omega _{2n}{t}-2\theta )}){\exp ({\text {i}2\varLambda T_1})}+{\rm c.c.}$
Once again, all terms explicitly written in (4.19) are resonant, as they share the same pair 
$(\omega _{2n},2)$ as the first-order homogeneous solution; hence a third-order compatibility condition, leading to the following normal form, must be enforced:
\begin{equation} \frac{\partial A_2}{\partial T_2}=\text{i}\zeta_{{DC}} {\varLambda} F^2{ \exp({\text{i}2\varLambda T_1})}+\text{i}\,\chi_{{DC}}A_2F^2+\text{i}\,\nu_{{DC}} |A_2|^2A_2, \end{equation}with
\begin{gather} \text{i}\,\zeta_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{{\varLambda} FF}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{{\varLambda} FF})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}, \end{gather}
\begin{gather} \text{i}\,\chi_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{A_2F\bar{F}}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{A_2F\bar{F}})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}, \end{gather}
\begin{gather} \text{i}\,\nu_{{DC}}=\frac{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{dyn}}}^{A_2\bar{A}_2A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\mathcal{F}}_{3_{{kin}}}^{A_2\bar{A}_2A_2})\, r\text{d}r}{\int_{z=0}^{}(\hat{\eta}_1^{A_2 {\dagger}}\hat{\varPhi}_1^{A_2}+\hat{\varPhi}_1^{A_2 {\dagger}}\hat{\eta}_1^{A_2})\, r\text{d}r}. \end{gather}
As a last step in the derivation of the final amplitude equation for the DC waves and in order to eliminate the implicit small parameter 
$\epsilon$, we unify (4.15) and (4.20) into a single equation recast in terms of the physical time 
$t=T_1/\epsilon =T_2/\epsilon ^2$, physical forcing control parameters 
$f=\epsilon F$, 
$\lambda =\epsilon \varLambda$, and total amplitude 
$A=\epsilon A_2$. This is achieved by summing (4.15) and (4.20) along with their respective weights 
$\epsilon ^2$ and 
$\epsilon ^3$, thus obtaining
\begin{equation} \frac{{\rm d}B}{{\rm d}t}={-}\text{i}\,(2{\lambda}-\chi_{{DC}} f^2) B + \text{i}\,(\zeta_{{DC}}{\lambda} + \mu_{{DC}}) f^2 + \text{i}\,\nu_{{DC}} |B|^2B, \end{equation}
where the change of variable 
$A=B{\rm e}^{\text {i}2 {\lambda } t}$ has been introduced for convenience. As in § 4.1, by turning to polar coordinates, 
$B=|B| {\rm e}^{\text {i}\varTheta }$, splitting the modulus and phase parts of (4.22) and looking for a stationary solution, 
${\rm d}/{\rm d}t=0$ with 
$|B|\ne 0$, the following implicit relation is obtained:
\begin{equation} \tilde{d}_s\varOmega^2-\sqrt{\left(2{\lambda}-\nu_{{DC}}|B|^2\right)|B|/\left[\chi_{{DC}}|B|\pm \left(\zeta_{{DC}}{\lambda}+\mu_{{DC}}\right)\right]}=0, \end{equation}
where only the real solutions corresponding to 
$f=\tilde {d}_s\varOmega ^2>0$ are retained, as the combinations 
$\tilde {d}_s\varOmega ^2<0$ are not physically meaningful.
Although two more terms appear in (4.22) and the dependence on the forcing amplitude is different with respect to the SC case, i.e. 
$f^2$ instead of 
$f$, thus leading to the square root in (4.23), amplitude equation (4.22) is reminiscent of that given in (4.6). Indeed, (4.22) contains essentially three contributions,
\begin{equation} \lambda \leftrightarrow \left(2{\lambda}-\chi_{{DC}}f^2\right),\quad \mu_{{SC}}f\leftrightarrow\left(\zeta_{{DC}}{\lambda}+\mu_{{DC}}\right)f^2,\quad \nu_{{SC}}\leftrightarrow\nu_{{DC}}, \end{equation}in order: a detuning term (forcing-amplitude-dependent), an additive (quadratic) forcing term (forcing-frequency-dependent) and the classic cubic restoring term. Hence, the same qualitative hardening or softening nonlinear behaviours as well as hysteresis, typical features of the Duffing equation, are expected under the hypotheses of the present analysis.
The total flow solution predicted by the WNL model for DC waves and reconstructed as
\begin{equation} \boldsymbol{q}_{DC}=\left\{\varPhi,\eta\right\}^{\rm T}=\boldsymbol{q}_1+\boldsymbol{q}_2 \end{equation}is compared in figures 6 and 7 with experiments from Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) (see also figure 2).
(a–d) The WNL prediction for DC waves versus experiments by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) (R14) (reported in figure 2) in terms of 
$\Delta \tilde {\delta }$ for a container diameter 
$D=0.144\,\text {m}$, 
$\tilde {H}=0.52$ and for various 
$\tilde {d}_s$. Dotted black lines: linear model (3.9). Red lines: straightforward asymptotic (3.14). Light-blue solid and dashed lines: stable and unstable branches, respectively, predicted by the WNL solution (4.25). The normal form coefficient values for this configurations are 
$\chi _{{DC}}=2.755$, 
$\zeta _{{DC}}=0.150$, 
$\mu _{{DC}}=0.129$ and 
$\nu _{{DC}}=10.018$. Violet lines: WNL solution (4.28) including the 
$\epsilon ^3$-order correction discussed in § 4.2.3.
Same as figure 6 (same colour convention) but for a container diameter 
$D=0.287\,\text {m}$ and 
$\tilde {H}=0.50$ (experiments are reported in figure 4.19 of Reclari (Reference Reclari2013) (R13)). Normal form coefficients: 
$\chi _{{DC}}=2.709$, 
$\zeta _{{DC}}=0.169$, 
$\mu _{{DC}}=0.134$ and 
$\nu _{{DC}}=9.885$.
4.2.2. Experiments versus WNL prediction: wave amplitude
In figures 6 and 7, the WNL prediction of DC waves according to (4.25) (light-blue solid and dashed lines) is quantitatively compared with the experimental measurements from Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) in terms of maximum non-dimensional crest-to-trough contact line amplitude 
$\Delta \tilde {\delta }$ for different values of the shaking diameters 
$\tilde {d}_s$ corresponding to those of figure 2 in the frequency window close to 
$\omega _{21}/2$.
The improvement gained through the formal WNL analysis, when compared with the linear and straightforward asymptotic models, is striking. The amplitude equation model correctly predicts the surge of the DC swirling and the resulting finite-amplitude saturation via hardening nonlinear mechanism, thus markedly narrowing the gap with experiments for all values of 
$\tilde {d}_s$ considered and for different container configurations.
Notwithstanding such an improvement, figure 6 highlights the main limitation of the present amplitude equation model for DC waves. Indeed, one notices that, while at larger shaking diameters, i.e. 
$\tilde {d}_s=0.13$ and 
$0.20$, a DC wave first emerges on the top of a SC dynamics and eventually a DC wave breaking occurs at larger frequencies, a jump-down transition from DC to SC takes place with increasing 
$\varOmega$ at lower shaking diameters, i.e. 
$\tilde {d}_s=0.07$ and 
$0.10$ for 
$D=0.144\,\text {m}$. This well-known hysteretic behaviour can be reasonably ascribed to the viscous dissipation of the system. For instance, at sufficiently small shaking diameters, e.g. 
$\tilde {d}_s\approx 0.02$ (see figure 2), the DC dynamics does not manifest at all, as the energy pumped into the system by the external forcing is likely not sufficient to dominate over the system viscous dissipation, whose effect also depends on the container diameter 
$D$. Indeed, figure 7 clearly shows that larger diameters, i.e. 
$D=0.287\,\text {m}$, generate less dissipation. It follows that for larger 
$D$, on increasing the driving frequency at a fixed shaking diameter, e.g. 
$\tilde {d}_s=0.10$, the free surface is more likely to undergo wave breaking, rather than a jump-down transition (see figures 6b and 7c). Obviously, the inviscid model employed here is not capable of predicting the so-called jump-down frequency. In Appendix A, a heuristic viscous damping model is introduced to tentatively overcome these limitations.
Finally, we note that, for frequency moderately far from the super-harmonic resonance, the agreement of the WNL model with experiments and with the linear solution, which behaves well away from 
$\varOmega \approx \omega _{11}$, progressively deteriorates. This is particularly evident on the lower stable branch and can be ascribed to the fact that the asymptotic model is essentially formalized for a fixed driving frequency, i.e. 
$\varOmega \approx \omega _{21}/2$. The second-order correction (proportional to 
$\varLambda F$ and discussed in § 4.2.1) to the leading-order particular solution seems to be sufficient to guarantee a fairly good agreement of the upper stable branch in a relatively wide range of frequency 
$\varOmega <\omega _{21}/2$. However, for the lower stable branch, i.e. 
$\varOmega >\omega _{21}/2$, the agreement with non-breaking-wave experiments sufficiently far from resonance is still relatively poor, despite the fact that these measurements essentially follow the linear prediction.
4.2.3. The 
$\epsilon ^3$-order correction to the leading-order SC particular solution
The last consideration in § 4.2.2 suggests that the leading-order SC solution with its second-order correction is only accurate in a limited range of frequency and higher-order terms should be accounted for in order to better retrieve the linear prediction far from resonance. In this regard, as the present asymptotic expansion is pursued up to the 
$\epsilon ^3$ order, we note that the third-order forcing contains a term, namely
\begin{equation} \boldsymbol{\mathcal{F}}_3^{\varLambda\varLambda F}={-}\text{i}\varLambda^2F\boldsymbol{\mathsf{B}}\hat{\boldsymbol{q}}_2^{\varLambda F} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)})\exp({\text{i}\varLambda T_1})+{\rm c.c.} , \end{equation}
generated by the derivative of the second-order solution with respect to the slow time scale 
$T_1$. This term is not resonating in 
$\exp ({\text {i}(\omega _{2n}-2\theta )})$ and therefore it can be gathered together with the other third-order non-resonating terms (NRT) in (4.19), which in asymptotic models are typically ignored, unless one aims to proceed to higher orders. Nevertheless, such a forcing term produces a system response:
\begin{equation} \boldsymbol{q}_3^{\varLambda\varLambda F}=\varLambda^2F\hat{\boldsymbol{q}}_3^{\varLambda \varLambda F} \exp({\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)})\exp({\text{i}\varLambda T_1})+{\rm c.c.} , \end{equation}
which precisely represents the third-order frequency correction of the leading-order particular solution and acts similarly to the second-order frequency correction 
$\boldsymbol {q}_2^{\varLambda F}=\varLambda F\hat {\boldsymbol {q}}_2^{\varLambda F}\exp ({\text {i}((\omega _{2n}/2)t-\theta )})\exp ({\text {i}\varLambda T_1})+{\rm c.c.}$ discussed in § 4.2.1.
A better intuition about why this third-order correction should improve the prediction farther away from the super-harmonic resonance is given in the following. As anticipated in § 4.2.1, the first-order particular solution simply represents the linear system response to the external forcing 
$\epsilon F\cos {(\varOmega t-\theta )}$. In general, the resulting amplitude is 
$\propto \epsilon F/(\varOmega ^2-\omega _{1n}^2)$, but in our asymptotic analysis the leading-order forcing frequency is frozen to 
$\omega _{2n}/2$, which implies that the first-order particular solution has an amplitude fixed to 
$\epsilon F/(\omega _{2n}^2 /4-\omega _{1n}^2)$, which does not account for the frequency detuning. If one replaces 
$\varOmega =\omega _{2n}/2+\epsilon \varLambda$ in 
$\epsilon F/(\varOmega ^2-\omega _{1n}^2)$ and takes its Taylor expansion, the 
$\epsilon ^2$-order term is 
$\propto \epsilon ^2 F\varLambda$, while the 
$\epsilon ^3$-order term is 
$\propto \epsilon ^3F\varLambda ^2$. It naturally follows that accounting for this third-order correction (4.27) leads to a more accurate description of the linear response to the external forcing and, therefore, it should give a better model prediction farther away from resonance, where the DC amplitude 
$|B|\approx 0$ and the non-breaking-wave experiments follow the linear theory.
To conclude, taking the total solution as
\begin{equation} \boldsymbol{q}_{DC}=\left\{\varPhi,\eta\right\}^{\rm T}=\boldsymbol{q}_1+\boldsymbol{q}_2+\boldsymbol{q}_3^{\varLambda \varLambda F} \end{equation}
is expected to leave the amplitude saturation prediction for the DC wave, 
$|B|$, unaltered, as it does not contribute to the amplitude equation solution, and, simultaneously, to better describe the SC swirling farther away from 
$\varOmega \approx \omega _{2n}/2$ (at least where no breaking waves occur).
The influence of the aforementioned 
$\epsilon ^3$-order corrections on the prediction for DC waves is shown as violet lines in figures 6 and 7, where it can be seen that accounting for the additional term allows one to eventually close the gap with the experiments even farther away from resonance, hence ensuring to the WNL model a wider frequency range of validity in all cases examined.
4.2.4. Experiments versus WNL prediction: free-surface reconstruction
In figure 8, the WNL model prediction for the DC waves is compared with the straightforward asymptotic prediction discussed in § 3 and the experimental measurements for DC waves from Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). The direct quantitative comparison is here outlined in terms of dimensionless and phase-averaged wave height measured at the sidewall, 
$\tilde {\delta }(\theta )$.
(a,c,e,g) Comparison of the dimensionless and phase-averaged wave height measured at the wall 
$\tilde {\delta }(\theta,{\rm \pi} /\varOmega )$ (black circles) with the straightforward asymptotic solution rebuilt via (3.14) (grey solid line) and the WNL solution for the DC wave (4.25). Panels correspond to 
$\tilde {H}=0.52$, 
$\tilde {d}_s=0.11$ and 
$D=0.144\,\text {m}$. The experimental measurements, here shown as black circles, are available in Reclari (Reference Reclari2013), except for those of (e), which are provided in Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014). Note that (c) the nonlinear prediction has a very large amplitude. (b,d, f,h) Corresponding three-dimensional free-surface deformation, 
$\eta (r,\theta,{\rm \pi} /\varOmega )$, reconstructed via (4.25). The transition from SC to DC swirling via hardening nonlinearity is clearly visible moving from top to bottom, i.e. for increasing frequency.
We observe that, if at 
$\varOmega /\omega _{21}=0.490$ both models match satisfactorily the experimental points, as soon as 
$\varOmega /\omega _{21}=0.5$ is approached, the straightforward asymptotic solution diverges due to the resonant (second-order) super-harmonic term, while the WNL solution predicts correctly the finite-amplitude saturation and the emergence of a DC wave on the top of a SC one. The WNL model for DC waves remains in fairly good agreement even at larger driving frequency, although the increasing phase asymmetry between the two local peaks at 
$\theta ={\rm \pi} /2$ and 
$3{\rm \pi} /2$ is not retrieved by the present inviscid asymptotic analysis, where secondary effects, e.g. the phase shift induced by viscous dissipation and influence of other higher modes, as well as stronger nonlinear effects for increasing wave amplitudes are overlooked.
For completeness, the three-dimensional free surface, 
$\eta (r,\theta,{\rm \pi} /\varOmega )$, is reconstructed through (4.25) and shown in figure 8(b,d,f,h), where, for increasing shaking frequencies, the nonlinear transition from nearly SC to DC swirling is evident.
4.2.5. The Helmholtz–Duffing oscillator analogy
While the Duffing equation is known to capture period-3 and period-1/3 dynamics arising from the cubic nonlinearity (Jordan & Smith Reference Jordan and Smith1999; Kalmár-Nagy & Balachandran Reference Kalmár-Nagy and Balachandran2011), as those observed by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) and occasionally by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), it cannot predict the period-halving (the system responds at a frequency which is twice that of the external forcing) dynamics associated with the super-harmonic resonance investigated in this paper. Therefore, in connection with § 4.1.3, here we aim to identify the simplest possible mechanical oscillator that could mimic, at least from a qualitative perspective, the period-1/2 dynamics studied in this work.
The WNL analysis as well as the straightforward asymptotic model highlighted the crucial role of quadratic nonlinearities emerging at second order and from which the DC dynamics stems. At the same time, the WNL model revealed that second-order terms only are not sufficient to capture all the dynamics features owing to the lack of restoring terms and, therefore, cubic nonlinearities must be retained. These considerations suggest that the DC dynamics could be tentatively described by a driven oscillator with both quadratic (asymmetric) and cubic (symmetric) nonlinear terms, i.e.
\begin{equation} \ddot{x}+2\sigma\dot{x}+ x+c_2 x^2 + c_3 x^3 = p\cos{\varOmega t}. \end{equation}Equation (4.29), also commonly known as the Helmholtz–Duffing equation, has wide applications in engineering problems, such as those related to beams, plates and shells subjected to an initial static curvature (Mirzabeigy, Yazdi & Nasehi Reference Mirzabeigy, Yazdi and Nasehi2014; Askari et al. Reference Askari, Saadatnia, Younesian, Yildirim and Kalami-Yazdi2011), whose governing equations are reduced to a second-order nonlinear ordinary equation with quadratic and cubic nonlinear terms (Ke, Yang & Kitipornchai Reference Ke, Yang and Kitipornchai2010; Alijani, Bakhtiari-Nejad & Amabili Reference Alijani, Bakhtiari-Nejad and Amabili2011; Fallah & Aghdam Reference Fallah and Aghdam2011).
Among the diverse asymptotic solutions of (4.29) in different limits (Rega Reference Rega1995; Benedettini & Rega Reference Benedettini and Rega1989; Kovacic & Brennan Reference Kovacic and Brennan2011), the most relevant to our work is that of Benedettini & Rega (Reference Benedettini and Rega1989). Within the context of planar nonlinear response of suspended elastic cables to an external excitation, they derived an amplitude equation which is concerned with the first or fundamental super-harmonic excitation, i.e. 
$\varOmega \approx 1/2$, of (4.29). Their WNL approach is detailed in Appendix C, with the additional assumption of vanishing damping 
$\sigma =0$. Assuming 
$2\varOmega =1+\lambda =1+\epsilon \hat {\lambda }$, small nonlinearities, 
$c_2=\epsilon \widehat {c_2}$ and 
$c_3=\epsilon ^2\hat {c}_3$, and introducing two slow time scales, one obtains
\begin{equation} {\rm d}D/{\rm d}t={-}\text{i}\,({2\lambda} +c_5f^2 )D +\text{i}\,\left(1-{\lambda}\right)c_2f^2/2-\text{i}\,4c_4|D|^2D, \end{equation}
with 
$C=D{\rm e}^{\text {i}{2\lambda } t}$ and with the auxiliary coefficients 
$c_4$ and 
$c_5$ (both functions of 
$c_2$ and 
$c_3$) defined in Benedettini & Rega (Reference Benedettini and Rega1989). By comparing term by term, the analogy with (4.22) is evident.
To conclude, although the DC dynamics examined in this paper is intrinsically related to the simultaneous interplay of multiple waves, thus making particularly challenging an accurate one-to-one quantitative comparison with a single-degree-of-freedom mechanical model, (4.30) seems to suggest that the actual inviscid sloshing dynamics in the DC regime may be, at least qualitatively, described by the undamped Helmholtz–Duffing equation (4.29) driven super-harmonically.
5. Conclusion
With regards to orbital-shaken cylindrical containers and, specifically, to the careful experimental campaign reported in Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014), a WNL analysis via a multiple-time-scale method was formalized in § 4 in order to investigate diverse features of the steady-state free-surface dynamics and, particularly, the DC wave dynamics pertaining at half the frequency of the first 
$m=2$ natural mode.
After having discussed the substantial limitations of the straightforward expansion procedure proposed by Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) and summarized in § 3, the WNL analysis was first formulated under the most common condition of pure harmonic resonance. Despite the inviscid assumption, the WNL analysis developed for the SC wave dynamics was shown to be in fairly good agreement with all the experimental measurements. In fact, the present model correctly describes the close-to-resonance hardening nonlinear behaviour experimentally observed. The agreement remains sufficiently accurate until the free surface eventually breaks and a transition to a fully nonlinear regime occurs.
It is well assessed in the literature that the close-to-harmonic-resonance sloshing dynamics can be modelled (from both qualitative and quantitative perspectives; Bäuerlein & Avila Reference Bäuerlein and Avila2021) by a single-degree-of-freedom system with a cubic nonlinearity and driven harmonically, i.e. by the famous Duffing oscillator, as rigorously proved for a two-dimensional rectangular container laterally excited (Ockendon & Ockendon Reference Ockendon and Ockendon1973). Without surprise, this was shown to hold for the case of orbital-shaken cylindrical containers as well.
The WNL analysis was then extended to the more complex case of a DC swirling. The overall agreement with experiments and, especially, the improvements with respect to the simple straightforward asymptotic model are remarkable in all cases considered, although the slight asymmetry observed in the reconstruction of the periodic free-surface dynamics at the sidewall was not retrieved in the present model.
To the knowledge of the authors, a formal amplitude equation describing the super-harmonic DC sloshing dynamics in orbital-shaken containers and coupled with a thorough experimental validation has not been reported in the literature yet, hence representing the most significant finding of this work.
Lastly, by analogy with the close-to-harmonic-resonance dynamics for SC waves, for which the Duffing oscillator represents a suitable mechanical analogy, a one-degree-of-freedom mechanical oscillator having both quadratic and cubic nonlinear terms, commonly referred to as a Helmholtz–Duffing oscillator, driven super-harmonically, was tentatively identified as the simplest possible mechanical system that could mimic, at least qualitatively, the super-harmonic DC sloshing dynamics investigated in this paper. The Helmholtz–Duffing equation was largely adopted in the last few decades within the context of structural analysis, i.e. beams, plates and shells subjected to an initial static curvature as well as suspended elastic cables (Nayfeh Reference Nayfeh1984; Benedettini & Rega Reference Benedettini and Rega1989), and it was here proposed as a direct mechanical analogy with the present orbital sloshing system.
The main limitation of the models derived in this work is intrinsic to the fundamental assumption of an inviscid fluid. This precludes correctly accounting for the jump-down transition experimentally observed for DC waves at low shaking amplitudes and, therefore, for an accurate estimation of the maximum amplitude response when such a transition occurs. Furthermore, in the absence of viscous boundary layers, the WNL time- and azimuthal-averaged mean flow reduces to a free-surface deformation only. This is in stark contrast with the existence of the so-called Eulerian mean flow (van den Bremer & Breivik Reference van den Bremer and Breivik2018), also known as viscous streaming flow, typically observed in experiments (Bouvard et al. Reference Bouvard, Herreman and Moisy2017). Therefore, the present work overlooks one of the essential points of interest in applications of orbital shaking. The mean flow, which contributes to an efficient mixing, is not captured.
The extension of the asymptotic models developed in this work to a viscous analysis is desirable, as it would enable one to predict quantitatively these secondary but fundamental effects for both cases of harmonic and super-harmonic resonances. However, it presently hinges on the subtle modelling of the moving contact line condition.
Acknowledgements
We acknowledge M. Farhat for fruitful discussions.
Funding
We acknowledge the Swiss National Science Foundation under grant 200021_178971.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Heuristic damping model: jump-down frequency and DC dynamics suppression at low driving amplitudes
In § 4.2.2 the WNL model for DC waves was compared with experimental measurements from Reclari (Reference Reclari2013) and Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014) in terms of non-dimensional maximum crest-to-trough contact line amplitude, 
$\Delta \tilde {\delta }$, for different non-dimensional shaking diameters, 
$\tilde {d}_s$, and container diameters, 
$D$ (see figures 6 and 7). We have observed that at larger shaking amplitudes, 
$\tilde {d}_s$, and for larger container diameter, 
$D$, a DC wave first emerges on the top of a SC wave at 
$\varOmega \approx \omega _{21}/2$ and eventually wave breaking occurs at larger frequencies. On the contrary, a jump-down transition from DC to SC then takes place with increasing 
$\varOmega$ at lower values of 
$\tilde {d}_s$ and/or for smaller 
$D$. The latter well-known hysteretic behaviour can be ascribed to the viscous dissipation of the system, obviously overlooked by the present inviscid analysis. In this appendix, viscous dissipation is tentatively reintroduced by employing a simple heuristic viscous damping model, as described in the following.
The viscous dissipation essentially arises at three locations: (i) at the solid tank boundary layers, i.e. bottom and sidewall, (ii) in the fluid bulk and (iii) at the free surface, the latter being typically negligible for ideal surface waves (in the absence of any form of contamination). A well-known formula for the prediction of the viscous damping coefficient of capillary–gravity waves in upright cylindrical containers was provided by Case & Parkinson (Reference Case and Parkinson1957) and Miles (Reference Miles1967). Such an estimation is computed according to the following formula:
\begin{align} \sigma=\frac{2k^{{2}}_{mn}}{Re}+\left(\frac{\omega_{mn}}{2Re}\right)^{1/2}\frac{k_{mn}}{\sinh{\left(2k_{mn}H\right)}}+\left(\frac{\omega_{mn}}{2Re}\right)^{1/2}\left[\frac{1}{2}\frac{1+\left(m/k_{mn}\right)}{1-\left(m/k_{mn}\right)}-\frac{k_{mn}H}{\sinh{\left(2k_{mn}H\right)}}\right], \end{align}
where the first term represents the bulk dissipation and the second and third terms are related to the dissipation occurring at the solid bottom and sidewall, respectively. In (A1), 
$H=h/R$ is the non-dimensional fluid depth, 
$k_{mn}$ is the non-dimensional wavenumber associated with mode 
$(m,n)$, 
$\omega _{mn}$ is the corresponding natural frequency obeying to the dispersion relation (3.6) and 
$Re=g^{1/2}R^{3/2}/\nu$ is the Reynolds number (
$\nu$ denotes the kinematic viscosity of the fluid). In § 4.2.2 an amplitude equation, governing the dynamics of a natural mode 
$(2,n)$ (which leads to the DC wave dynamics observed close to 
$\varOmega \approx \omega _{21}/2$), was derived. For mode 
$(2,1)$ in the conditions of figure 6, i.e. pure water with 
$\rho =1000\,\text {kg}\,\text {m}^{-3}$, 
$\gamma =0.072\,\text {N}\,\text {m}^{-1}$, 
$\nu =1\times 10^{-6}\,\text {m}^2\,\text {s}^{-1}$, 
$D=0.144\,\text {m}$ (for which the Bond number is 
$Bo=705.6$) and 
$H=1.04=2\tilde {H}$, the values 
$Re=60\,480$, 
$k_{21}=3.0542$ and 
$\omega _{21}=1.7561$ give a non-dimensional viscous damping coefficient 
$\sigma =0.0051$, mostly produced by the sidewall boundary layer. Typically, as in the present case and as supported by experimental (Cocciaro, Faetti & Festa Reference Cocciaro, Faetti and Festa1993) and numerical (Viola et al. Reference Viola, Brun and Gallaire2018) evidence, the viscous damping rate can be interpreted as a slow damping process over a faster time scale represented by the wave oscillation. Under this hypothesis, which translates in the assumption of a viscous damping coefficient of order 
$\epsilon ^2$ within the present WNL framework, the damping coefficient can be added a posteriori, i.e. in a phenomenological way, to the final inviscid amplitude equation from (4.22), leading to
\begin{equation} \frac{{\rm d}B}{{\rm d}t}={-}[\sigma + \text{i}(2{\lambda}-\chi_{{DC}} f^2)] B + \text{i}(\zeta_{{DC}}{\lambda}+\mu_{{DC}}) f^2 + \text{i}\nu_{{DC}} |B|^2B. \end{equation}The stationary form of (A2) can be rearranged in the following implicit form:
\begin{equation} (2{\lambda}-\nu_{{DC}}|B|^2-\chi_{{DC}}f^2)|B| \pm\sqrt{f^4(\zeta_{{DC}}{\lambda}+\mu_{{DC}})^2-(\sigma|B|)^2}=0, \end{equation}which can be solved using the Matlab function fimplicit. The effect of viscous dissipation on the DC regime is investigated in figure 9 for two representative values of the shaking diameter.
(a) Same case as for figure 6(a) with 
$\tilde {d}_s=0.07$. (b) Same as (a), but for 
$\tilde {d}_s=0.02$ (from figure 2), a value at which the DC dynamics does not manifest. Solid and dashed lines correspond to stable and unstable branches, respectively, computed via the WNL analysis in the inviscid case and for different values of damping coefficient, with 
$\sigma$ given by (A1). Markers correspond to the experimental points shown in figure 2 and extracted from Reclari et al. (Reference Reclari, Dreyer, Tissot, Obreschkow, Wurm and Farhat2014).
The case of figure 9(a) shows that the so-called jump-down frequency is somewhere between 
$\varOmega \in [0.675,0.685]$. The damping value produced by (A1) appears to be too small to match the experimental jump-down frequency; hence we tentatively added a pre-factor in order to fit the measurements. It turns out that a pre-factor of 1.35 is sufficient to provide a fairly good prediction of the jump-down frequency. We note that prediction (A1) does not involve any dissipation mechanism associated with the contact line, i.e. contact line hysteresis (Keulegan Reference Keulegan1959; Miles Reference Miles1967; Dussan Reference Dussan1979; Hocking Reference Hocking1987; Cocciaro et al. Reference Cocciaro, Faetti and Festa1993; Kidambi Reference Kidambi2009; Viola et al. Reference Viola, Brun and Gallaire2018; Viola & Gallaire Reference Viola and Gallaire2018; Bongarzone, Viola & Gallaire Reference Bongarzone, Viola and Gallaire2021b) or possible surface contamination (Henderson & Miles Reference Henderson and Miles1990, Reference Henderson and Miles1994). Indeed, depending on the configuration, contact line dynamics may rule the overall dissipation, with a measured damping coefficient up to 10–20 times larger (Benjamin & Ursell Reference Benjamin and Ursell1954; Hocking Reference Hocking1987; Kidambi Reference Kidambi2009) than that predicted by (A1). Comparisons of the theoretical damping coefficient value with that measured in moving contact line experiments, due to unavoidable sources of uncertainty in the meniscus dynamics, have always been mostly qualitative, rather than quantitative, requiring often the use of fitting parameters. For instance, in their predictive theory for single-mode Faraday experiments, Henderson & Miles (Reference Henderson and Miles1990) used an effective fluid viscosity three times larger than the actual one. Recently, Bäuerlein & Avila (Reference Bäuerlein and Avila2021) have measured the damping coefficient of the first anti-symmetric sloshing mode in a quasi-two-dimensional rectangular container, which was seen to be approximately 1.5 larger than that predicted by the theory (Faltinsen & Timokha Reference Faltinsen and Timokha2009).
Even in the absence of strong contact line dissipation, free-surface contamination may strongly contribute to the overall damping coefficient. Henderson & Miles (Reference Henderson and Miles1990) have considered a fully contaminated free surface, which can be modelled by a surface film that is free to move vertically but cannot stretch horizontally. They have also derived an analytical expression for the associated damping, which reads 
$(\omega _{mn}/2Re)^{1/2}k_{mn}\cosh ^2{k_{mn}H}/\sinh {2k_{mn}H}$. If such a contribution was accounted for in (A1), it would produce a value of 
$0.0109$, which is approximately twice the damping obtained without contamination (and would correspond to the green lines in figure 9(
$a$) for 
$2\sigma$). The need for a pre-factor of 1.35 in figure 9(a), which approximately corresponds to a fictitious fluid with a dynamic viscosity 1.8 times larger, is therefore not surprising when the damping is computed via (A1) and contact line dissipation as well as free-surface contamination are neglected.
We remark that the reasonings outlined in this appendix in order to elucidate the effect of viscosity are in fact only qualitative. Many aspects are ignored in the present inviscid analysis with phenomenological damping, two of which are commented upon in the following.
Prediction (A1) is only valid for free capillary–gravity waves, whereas dissipation rates of forced wave motions are generally more complex. In particular, the DC wave evolves out of a non-resonant forced SC swirling motion. A proper viscous WNL analysis would produce complex eigenmodes and responses (due to the phase shift owing to viscosity) and hence complex-valued normal form coefficients. For instance, among these coefficients, the imaginary part of 
$\nu _{{DC}}$ (or 
$\nu _{{SC}}$) multiplied by 
$|B|^2$ in (A2) could be interpreted as a sort of nonlinear damping (Douady Reference Douady1990), 
$(\sigma -\text {Im}[\nu ]|B|^2)$, whose contribution to the overall dissipation mechanisms is expected to increase at larger wave amplitudes, hence influencing the location of the jump-down frequency (note that 
$-\text {Im}[\nu ]$ can be 
$>0$). Moreover, the imaginary part of 
$\chi f^2$ in DC waves would contribute to the overall damping, 
$(\sigma +\text {Im}[\chi ]f^2-\text {Im}[\nu ]|B|^2)$, by introducing a further effect, proportional to 
$f^2$, that finds its origin in the fact that the DC wave emerges at second order owing to nonlinear square terms in the first-order non-resonant SC motion.
In contradistinction with the case of a pinned (or fixed) contact line, a formal viscous analysis undertaking the case of a moving contact line would require the introduction of a slip length model in order to regularize the well-known contact line stress singularity (Navier Reference Navier1823; Huh & Scriven Reference Huh and Scriven1971; Davis Reference Davis1974; Lauga, Brenner & Stone Reference Lauga, Brenner and Stone2007; Viola & Gallaire Reference Viola and Gallaire2018).
Most importantly, the inviscid WNL model is not capable of describing the continuous modulation of the phase lag between the external forcing and the wave amplitude response, which has been recently demonstrated by Bäuerlein & Avila (Reference Bäuerlein and Avila2021) (for unidirectional sloshing waves in a three-dimensional rectangular container) to be of crucial importance in the correct prediction of the jump-down frequency, otherwise often inaccurate, even when the considered damping coefficient is that measured experimentally. In principle, a formal viscous analysis, as briefly introduced above, is expected to correctly capture such a phase lag.
Another interesting case, which is worth commenting upon, is that shown in figure 9(b). At a shaking diameter 
$\tilde {d}_s=0.02$ (the lowest reported in figure 2), the DC dynamics was not observed at all. This is in conflict with the inviscid straightforward asymptotic analysis, which always prescribes a divergent behaviour close to the dominant super-harmonic, 
$\varOmega \approx \omega _{21}/2$, even for vanishing 
$\tilde {d}_s$. However, as soon as viscous dissipation is introduced, the energy pumped into the system is not sufficient to overcome dissipative effects and DC waves are essentially suppressed, with a system response that follows satisfactorily the linear solution (see figure 2) showing a SC dynamics ranging over the whole frequency window, 
$\varOmega /\omega _{11}\in [0,1]$, in agreement with experimental evidence.
Appendix B. Asymptotic harmonic solution of the undamped Duffing equation
By analogy with the WNL analysis for harmonic SC wave dynamics presented in § 4.1, we look for an asymptotic solution of the undamped Duffing equation
\begin{equation} \ddot{x}+x+ c_3 x^3 =p\cos{\varOmega t}, \end{equation}
having the form 
$x=x_0+\epsilon x_1$. Additionally, as standard in asymptotic solutions of the Duffing equation, we assume a small external forcing amplitude 
$p=\epsilon \hat {p}$ and detuning from the exact resonance, i.e. 
$\varOmega =1+\lambda =1+\epsilon \hat {\lambda }$, small nonlinearities through 
$c_3=\epsilon \hat {c}_3$ and the existence of a characteristic slow time scale 
$\hat {t}_1=\epsilon t$. Under these assumptions, the 
$\epsilon ^0$-order homogeneous solution simply reads
\begin{equation} x_0=C\left(\hat{t}_1\right){\rm e}^{\text{i} t}+{\rm c.c.}, \end{equation}
with 
$C(\hat {t}_1)$ to be determined at next order. At order 
$\epsilon$ one can readily verify that, in order to avoid secular terms, a solvability condition must be satisfied. Such a condition leads to the very classical amplitude equation
\begin{equation} {\rm d}D/{\rm d}t={-}\text{i}\,{\lambda} D + \text{i}\,\left({-}1/4\right) p+\text{i}\,\left(3c_3/2\right)|D|^2D, \end{equation}
where the change of variable 
$C=D{\rm e}^{\text {i}\lambda t}$ was introduced and each quantity was recast in terms of the corresponding physical value (to eliminate the implicit small parameter 
$\epsilon$).
By noticing that
\begin{equation} -1/4\ \leftrightarrow\ \mu_{{SC}},\quad 3c_3/2\ \leftrightarrow\ \nu_{{SC}}, \end{equation}one immediately recognizes that (4.6) has indeed the same structure of the formal amplitude (B3), thus suggesting that the continuous sloshing system and the one-degree-of-freedom Duffing system, under the specific conditions listed above, behave essentially in the same way.
Appendix C. Asymptotic super-harmonic solution of the undamped Helmholtz–Duffing equation
In this appendix, although with the additional assumption of vanishing damping, we briefly summarize the super-harmonic WNL solution of the Helmholtz–Duffing equation,
\begin{equation} \ddot{x}+x+c_2 x^2 + c_3 x^3= p\cos{\varOmega t}, \end{equation}derived by Benedettini & Rega (Reference Benedettini and Rega1989) and introduced in § 4.1.
We look for an asymptotic solution of the form 
$x=x_0+\epsilon x_1+\epsilon ^2 x_2$ to (C1) with 
$\sigma =0$ (undamped oscillator), 
$2\varOmega =1+\lambda =1+\epsilon \hat {\lambda }$ and with small nonlinearities through 
$c_2=\epsilon \hat {c}_2$ and 
$c_3=\epsilon ^2\hat {c}_3$ (with the cubic term one order smaller than the quadratic one). The existence of two slow time scales is hypothesized, 
$\hat {t}_1=\epsilon t$ and 
$\hat {t}_2=\epsilon ^2 t$. Under these assumptions, the solution of the 
$\epsilon ^0$-order forced linear problem reads
\begin{equation} x_1=C\left(\hat{t}_1,\hat{t}_2\right){\rm e}^{\text{i} t}+f {\exp({\text{i}\left(1/2\right)t})\exp({\text{i}(\hat{\lambda}/2) \hat{t}_1})}+{\rm c.c.}, \end{equation}
with 
$f=(2/3)p$ and 
$C(\hat {t}_1,\hat {t}_2)$ to be determined at next order. At orders 
$\epsilon$ and 
$\epsilon ^2$, resonating terms produced by the weak quadratic and cubic nonlinearities, respectively, arise, thus requiring the imposition of two solvability conditions prescribing that amplitude 
$C(\hat {t})$ must obey to the following normal forms:
\begin{gather} \boxed{\epsilon^1}\text{ :}\quad {\rm d}C/{\rm d}\hat{t}_1=\text{i}\,\left(c_2/2\right)f^2 {\exp({\text{i}\hat{\lambda} \hat{t}_1})}, \end{gather}
\begin{gather} \boxed{\epsilon^2}\text{ :}\quad {\rm d}C/{\rm d}\hat{t}_2={-}\text{i}\,{\hat{\lambda}}\left(c_2/4\right) {\exp({\text{i}\hat{\lambda} \hat{t}_1})}-\text{i}\,c_5f^2A-\text{i}\,4c_4|C|^2C, \end{gather}
with the full expression of the auxiliary coefficients 
$c_4$ and 
$c_5$ (both functions of 
$c_2$ and 
$c_3$) given in Benedettini & Rega (Reference Benedettini and Rega1989). Combining (C3a) and (C3b) into a single amplitude equation (by summing the two expressions by their respective weights, i.e. 
$\epsilon$ and 
$\epsilon ^2$, and reintroducing the physical quantities in order to eliminate the dependence on the implicit small parameter 
$\epsilon$), one obtains
\begin{equation} {\rm d}D/{\rm d}t={-}\text{i}\,({\lambda} +c_5f^2 )D +\text{i}\,(1-{\lambda/2})c_2f^2/2-\text{i}\,4c_4|D|^2D, \end{equation}
with 
$C=D{\rm e}^{\text {i}{\lambda } t}$. Note that the procedure used in the perturbation analysis above and outlined in Benedettini & Rega (Reference Benedettini and Rega1989) is in fact equivalent to that followed in Nayfeh (Reference Nayfeh1984) for treating the same second-order super-harmonic resonance in a more general case of a two-term excitation. By comparing the various terms of (C4) with those of (4.22), the analogy is evident, thus suggesting that the actual inviscid sloshing dynamics in the DC wave regime may be, at least qualitatively, described by the undamped Helmholtz–Duffing equation (4.29) driven super-harmonically.
Appendix D. The $\epsilon ^2$- and $\epsilon ^3$-order dynamic and kinematic equations
For completeness, in this appendix we provide the second-and third-order asymptotic expansions of the free-surface boundary conditions with regards to the most complex formulation presented in this paper, i.e. that related to the DC swirling. At second order, the dynamic and kinematic equations evaluated at 
$z=\eta _0=0$ read, respectively,
\begin{gather} {\frac{\partial\varPhi_2}{\partial t}\,{+}\,\eta_2\,{-}\,\frac{1}{Bo}\frac{\partial\kappa}{\partial\eta}\left(\eta_2\right)={-}\frac{\partial\varPhi_1}{\partial T_1}\,{-}\,\eta_1\frac{\partial}{\partial z}\frac{\partial\varPhi_1}{\partial t}\,{-}\,\frac{1}{2}\boldsymbol{\nabla}\varPhi_1\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi_1\,{+}\,\frac{1}{Bo}\frac{1}{2}\frac{\partial^2\kappa}{\partial \eta^2}\left(\eta_1\right)^2\,{=}\,\mathcal{F}_{2_{{dyn}}},} \end{gather}
\begin{gather} {\frac{\partial\eta_2}{\partial t}-\frac{\partial\varPhi_2}{\partial z}={-}\frac{\partial\eta_1}{\partial T_1}-\boldsymbol{\nabla}\varPhi_1\boldsymbol{\cdot}\boldsymbol{\nabla}\eta_1+\eta_1\frac{\partial}{\partial z}\frac{\partial\varPhi_1}{\partial z}=\mathcal{F}_{2_{{kin}}},} \end{gather}
where 
$\boldsymbol {\nabla }\eta =\{\partial \eta /\partial r,r^{-1}\partial \eta /\partial \theta,0\}^{\rm T}$, 
$\partial \kappa /\partial \eta$ represents the first-order variation of the curvature,
\begin{equation} {\frac{\partial\kappa}{\partial\eta}=\frac{\partial^2}{\partial r^2}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{1}{r^2}\frac{\partial^2}{\partial\theta^2},} \end{equation}
and it is applied to 
$\eta _2$, while 
$\partial ^2\kappa /\partial \eta ^2$ is its second-order variation applied to 
$(\eta _1)^2$. However, as the system is expanded around 
$z=\eta _0=0$, it turns out that 
$\partial ^2\kappa /\partial \eta ^2=0$.
By pursuing the expansion up to the third order in 
$\epsilon$, one obtains
\begin{align} &\frac{\partial\varPhi_3}{\partial t}+\eta_3-\frac{1}{Bo}\frac{\partial\kappa}{\partial\eta}\left(\eta_3\right)={-}\frac{\partial\varPhi_1}{\partial T_2}-\frac{\partial\varPhi_2}{\partial T_1}-\eta_1\frac{\partial}{\partial z}\frac{\partial\varPhi_1}{\partial T_1}\nonumber\\ &\quad -\eta_2\frac{\partial}{\partial z}\frac{\partial\varPhi_1}{\partial t}-\eta_1\frac{\partial}{\partial z}\frac{\partial\varPhi_2}{\partial t} -\frac{1}{2}\left(\eta_1\right)^2\frac{\partial^2}{\partial z^2}\frac{\partial\varPhi_1}{\partial t}-\boldsymbol{\nabla}\varPhi_1\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi_2-\eta_1\boldsymbol{\nabla}\varPhi_1\boldsymbol{\cdot}\frac{\partial}{\partial z}\boldsymbol{\nabla}\varPhi_1\nonumber\\ &\quad +\frac{1}{Bo}\frac{1}{2}\frac{\partial^2\kappa}{\partial\eta^2}\left(2\eta_1\eta_2\right)+\frac{1}{Bo}\frac{1}{6}\frac{\partial^3\kappa}{\partial\eta^3}\left(\eta_1\right)^3=\mathcal{F}_{3_{{dyn}}}, \end{align}
\begin{align} &\frac{\partial\eta_3}{\partial t}-\frac{\partial\varPhi_3}{\partial z}={-}\frac{\partial\eta_1}{\partial T_2}-\frac{\partial\eta_2}{\partial T_1}-\boldsymbol{\nabla}\varPhi_1\boldsymbol{\cdot}\boldsymbol{\nabla}\eta_2-\boldsymbol{\nabla}\varPhi_2\boldsymbol{\cdot}\boldsymbol{\nabla}\eta_1-\eta_1\boldsymbol{\nabla}\eta_1\boldsymbol{\cdot}\frac{\partial}{\partial z}\boldsymbol{\nabla}\varPhi_1 \nonumber\\ &\quad +\eta_2\frac{\partial}{\partial z}\frac{\partial\varPhi_1}{\partial z}+\eta_1\frac{\partial }{\partial z}\frac{\partial\varPhi_2}{\partial z}+\frac{1}{2}\left(\eta_1\right)^2\frac{\partial^2}{\partial z^2}\frac{\partial\varPhi_1}{\partial z}=\mathcal{F}_{3_{{kin}}}, \end{align}
with 
$\partial ^3\kappa /\partial \eta ^3(\ne 0)$ the third-order variation of the curvature, whose explicit expression is here omitted for the sake of brevity, applied to 
$(\eta _1)^3$.
Note that the second order in the straightforward asymptotic model is retrieved by retaining the 
$\epsilon ^2$-order system only and imposing 
$\partial /\partial T_1=\partial /\partial T_2=0$. On the other hand, the equations above reduce to the second- and third-order problem in the SC swirling formulation when 
$\partial /\partial T_1=0$ and the external forcing term, 
$rF\cos {(\varOmega t-\theta )}$, appears on the right-hand side of (D4). At each order in 
$\epsilon$, by substituting the previous-order solutions, it is possible to explicitly separate the various forcing contributions by their temporal and azimuthal periodicities, thus leading to expressions (3.11) for the straightforward model, (4.5) for the SC model and (4.17) and (4.19) for the DC one.
For the calculation of the amplitude equation coefficients at 
$\epsilon ^3$ order, only resonant terms matter. These terms, with their corresponding amplitudes, are proportional to 
$\exp ({\text {i}(\omega _{1n}t-\theta )})$ for SC waves and to 
$\exp ({\text {i}(\text {i}\omega _{2n}t-2\theta )})$ for DC waves (or, more generally, proportional to 
$\exp ({\text {i}(\omega _{mn}t-m\theta )})$). As an example, in the following we provide the expression of 
$\hat {\mathcal {F}}_{3,\text {kin}}^{A\bar {A}A}$ used in (4.7b) (with 
$A=A_1$) and in (4.21c) (with 
$A=A_2$) and which, together with 
$\hat {\mathcal {F}}_{3,{dyn}}^{A\bar {A}A}$, represents the most involved third-order resonant forcing term:
\begin{align} &\hat{\mathcal{F}}_{3,\text{kin}}^{A\bar{A}A}={-}\boldsymbol{\nabla}\hat{\varPhi}_1^A\boldsymbol{\cdot}\boldsymbol{\nabla}\hat{\eta}_2^{A\bar{A}} -\boldsymbol{\nabla}\overline{\hat{\varPhi}_1^A}\boldsymbol{\cdot}\boldsymbol{\nabla}\hat{\eta}_2^{AA}-\boldsymbol{\nabla}\hat{\varPhi}_2^{AA}\boldsymbol{\cdot}\boldsymbol{\nabla}\overline{\hat{\eta}_1^A}-\boldsymbol{\nabla}\hat{\varPhi}_2^{A\bar{A}}\boldsymbol{\cdot}\boldsymbol{\nabla}\hat{\eta}_1^A\nonumber\\ &\quad -\hat{\eta}_1^A\boldsymbol{\nabla}\hat{\eta}_1^A\boldsymbol{\cdot}\frac{\partial}{\partial z}\boldsymbol{\nabla}\overline{\hat{\varPhi}_1^A}-\hat{\eta}_1^A\boldsymbol{\nabla}\overline{\hat{\eta}_1^A}\boldsymbol{\cdot}\frac{\partial }{\partial z}\boldsymbol{\nabla}\hat{\varPhi}_1^A-\overline{\hat{\eta}_1^A}\boldsymbol{\nabla}\hat{\eta}_1^A\boldsymbol{\cdot}\frac{\partial}{\partial z}\boldsymbol{\nabla}\hat{\varPhi}_1^A\nonumber\\ &\quad +\hat{\eta}_2^{A\bar{A}}\frac{\partial}{\partial z}\frac{\partial\hat{\varPhi}_1^A}{\partial z}+\hat{\eta}_2^{AA}\frac{\partial}{\partial z}\frac{\partial \overline{\hat{\varPhi}_1^A}}{\partial z}+\hat{\eta}_1^A\frac{\partial}{\partial z}\frac{\partial\hat{\varPhi}_2^{A\bar{A}}}{\partial z}+\overline{\hat{\eta}_1^A}\frac{\partial}{\partial z}\frac{\partial\hat{\varPhi}_2^{AA}}{\partial z}\nonumber\\ &\quad +\left(\hat{\eta}_1^A\hat{\eta}_1^A\right)\frac{1}{2}\frac{\partial^2}{\partial z^2}\frac{\partial\overline{\hat{\varPhi}_1^A}}{\partial z}+\left(\hat{\eta}_1^A\overline{\hat{\eta}_1^A}\right)\frac{1}{2}\frac{\partial^2}{\partial z^2}\frac{\partial\hat{\varPhi}_1^A}{\partial z}+\left(\overline{\hat{\eta}_1^A}\hat{\eta}_1^A\right)\frac{1}{2}\frac{\partial^2}{\partial z^2}\frac{\partial\hat{\varPhi}_1^A}{\partial z}. \end{align}The expression of 
$\hat {\mathcal {F}}_{3,{dyn}}^{A\bar {A}A}$ (not reported here for the sake of brevity) is computed analogously. The extraction of resonant terms, especially those appearing at third order, involves tedious calculations due to several possible combinations of the previous-order solutions. Nevertheless, the procedure is straightforward and systematic, so that tools of symbolic calculus, e.g. the software Wolfram Mathematica, which was indeed used in this work, can be employed to ease such a procedure.
Open access
