# An amplitude equation modeling the single–to–double crest wave transition
in orbital shaken cylindrical containers

Alessandro Bongarzone1  Margherita Guido2,3  François Gallaire1
<EMAIL_ADDRESS>1Laboratory of Fluid Mechanics and Instabilities,
École Polytechnique Fédérale de Lausanne, Lausanne, CH-1015, Switzerland
2Swiss Plasma Center, École Polytechnique Fédérale de Lausanne, Lausanne,
CH-1015, Switzerland 3Numerical Algorithms and High-Performance Computing,
École Polytechnique Fédérale de Lausanne, Lausanne, CH-1015, Switzerland

###### Abstract

The container motion along a planar circular trajectory at a constant angular
velocity, i.e. orbital shaking, is of interest in several industrial
applications, e.g. for fermentation processes or in cultivation of stem cells,
where good mixing and efficient gas exchange are the main targets. Under these
external forcing conditions, the free surface typically exhibits a primary
steady state motion through a single–crest dynamics, whose wave amplitude, as
a function of the external forcing parameters, shows a Duffing-like behaviour.
However, previous experiments in lab-scale cylindrical containers have
unveiled that, owing to the excitation of super-harmonics, diverse dynamics
are observable in certain driving-frequency ranges. Among these super-
harmonics, the double-crest dynamics is particularly relevant, as it displays
a notably large amplitude response, that is strongly favored by the spatial
structure of the external forcing. In the inviscid limit and with regards to
circular cylindrical containers, we formalize here a weakly nonlinear analysis
via multiple timescale method of the full hydrodynamic sloshing system,
leading to an amplitude equation suitable to describe such a super-harmonic
dynamics and the resulting single–to–double crest wave transition. The weakly
nonlinear prediction is shown to be in fairly good agreement with previous
experiments described in the literature. Lastly, we discuss how an analogous
amplitude equation can be derived by solving asymptotically for the first
super-harmonic of the forced Helmholtz–Duffing equation with small
nonlinearities.

###### keywords:

##  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, 1969; Wurm, 2004), as an alternative to stirred tanks, where the
liquid agitation results from a rotating impeller or the rotation of magnetic
rod. In these cultivation protocols, cells are in suspension in the
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. Thanks to 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
damages that their collapse can exert on cells (Handa-Corrigan et al., 1989;
Kretzmer & Schügerl, 1991; Papoutsakis, 1991), sparking interest in the
development of large-scale, in the hectoliter range, orbital shaken
bioreactors (Liu & Hong, 2001; Jesus et al., 2004; Muller et al., 2007). It is
therefore not a surprise if a significant body of research on the gas exchange
and mixing in these devices has emerged over the last two decades (Büchs et
al., 2000a, b; Büchs, 2001; Maier et al., 2004; Muller et al., 2005;
Micheletti et al., 2006; Zhang et al., 2009; Tissot et al., 2010; Tan et al.,
2011; Tissot et al., 2011; Klöckner & Büchs, 2012).
At a more fundamental level, the hydrodynamics of these orbital shaking
devices has received recent attention, from both experimental (Reclari et al.,
2014; Bouvard et al., 2017; Moisy et al., 2018) and theoretical (Reclari et
al., 2014; Horstmann et al., 2020) 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., 2020, 2021). Previous
studies make mostly use of classical existing theories for general linear
sloshing dynamics, reviewed for instance in Ibrahim (2005) or Faltinsen &
Timokha (2009).
In order to refine the linear potential model and, specifically, to predict
the occurrence of the super-harmonic wave dynamics observed experimentally,
Reclari (2013) and Reclari et al. (2014) proposed an inviscid weakly nonlinear
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 wave dynamics is
particularly relevant, as it appears to be the most stable and the one which
displays the largest amplitude response. However, their analysis, as typical
of straightforward asymptotic expansions, suffers from secular terms
(Castaing, 2005; Nayfeh, 2008) and, therefore, it still fails in describing
the correct nonlinear behaviour close to both harmonic and super-harmonic
resonances.
This limitation was partially overcome by Timokha & Raynovskyy (2017) and
Raynovskyy & Timokha (2018), who have applied the Narimanov–Moiseev multimodal
sloshing theory (Narimanov, 1957; Moiseev, 1958; Dodge et al., 1965;
Faltinsen, 1974; Narimanov et al., 1977; Lukovsky, 1990) describing the
nonlinear wave dynamics near the first harmonic resonance, which was proven to
be of the hard-spring type, in qualitative agreement with the observations of
Reclari (2013) and Reclari et al. (2014). Nevertheless, they were not able to
quantitatively compare their predictions with the latter experimental
measurements because these were done for $H=h/R\leq 1.2$, e.g. $H=1$ and
$1.04$ (with $h$ the fluid depth and $R$ the container radius), where, as
stated by the authors themselves, the adopted nonlinear modal equations are
not applicable due to the secondary resonance phenomena (Faltinsen et al.,
2016).
In the spirit of the aforementioned multimodal theory but with regards to
square-base basins, the resonant amplification of higher order modes was
investigated by Faltinsen et al. (2005), who formalized a so-called adaptive
asymptotic modal approach capable to improve the agreements with earlier
experiments. Yet, to the knowledge of the authors, the adaptive modal approach
was never extended to orbital shaken circular cylindrical containers.
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 (2013) and Reclari et al. (2014)
has not been provided yet.
The present work is precisely dedicated to the development of a weakly
nonlinear analysis based on the multiple timescale 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 double–crest saturation. Amplitude equations are
rigorously derived in an inviscid framework, which once amended with an ad-hoc
damping term as only tuning parameter, well match the experimental findings of
Reclari (2013) and Reclari et al. (2014). Lastly, the obtained amplitude
equations for harmonic single–crest and super-harmonic double–crest waves are
found to be compatible with the two well known one-degree-of-freedom (1dof)
systems, the Duffing and the Helmholtz-Duffing oscillators, respectively.
The manuscript is organized as follows. The flow configuration and governing
equations are introduced in §2. §3 is dedicated to briefly summarize the
salient points of the asymptotic model proposed by Reclari et al. (2014),
whose limitations motivated the present work. After tackling the more common
case of harmonic single–crest wave in §4.1, the weakly nonlinear amplitude
equation governing the super-harmonic double–crest wave dynamics is derived in
§4.2. Final comments and conclusions are outlined in §5.

##  2 Flow configuration and governing equations: potential model

Figure 1: Sketch of a cylindrical container of diameter $D=2R$ and filled to a
depth $h$. The gravity acceleration is denoted by $g$.
$O^{\prime}\mathbf{e}^{\prime}_{x}\mathbf{e}^{\prime}_{y}\mathbf{e}^{\prime}_{z}$
is the Cartesian inertial reference frame, while
$O\mathbf{e}_{x}\mathbf{e}_{y}\mathbf{e}_{z}$ is the Cartesian reference frame
moving with the container. The origin of the moving cylindrical reference
frame $\left(r,\theta,z\right)$ 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. $d_{s}$ is the diameter of the circular shaking trajectory,
characterized by a driving angular frequency $\Omega_{d}$.

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
$\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$)

$\dot{\mathbf{X}}_{0}=\begin{cases}-\frac{d_{s}}{2}\Omega_{d}\sin{\left(\Omega_{d}t-\theta\right)}\,\mathbf{e}_{r}\\\
\frac{d_{s}}{2}\Omega_{d}\cos{\left(\Omega_{d}t-\theta\right)}\,\mathbf{e}_{\theta}\end{cases}.$
(1)

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,
$\Phi_{tot}$, which is typically separated into a container, $\Phi_{c}$, and a
fluid component, $\Phi$. Hence, the liquid motion within the moving container
is governed by the Laplace equation,

$\Delta\Phi=\frac{1}{r}\frac{\partial\Phi}{\partial
r}+\frac{\partial^{2}\Phi}{\partial
r^{2}}+\frac{1}{r^{2}}\frac{\partial^{2}\Phi}{\partial\theta^{2}}+\frac{\partial^{2}\Phi}{\partial
z^{2}}=0,$ (2)

subjected to the homogeneous no-penetration condition,
$\nabla\Phi\cdot\mathbf{n}=\mathbf{0}$, at the solid sidewall and bottom, and
by the dynamic and kinematic free surface boundary conditions at $z=\eta$ (see
Ibrahim (2005)),

$\frac{\partial\Phi}{\partial
t}+\frac{1}{2}\nabla\Phi\cdot\nabla\Phi+\eta-\frac{\kappa\left(\eta\right)}{Bo}=rf\cos{\left(\Omega
t-\theta\right)},$ (3a) $\frac{\partial\eta}{\partial
t}+\frac{\partial\Phi}{\partial r}\frac{\partial\eta}{\partial
r}+\frac{1}{r^{2}}\frac{\partial\Phi}{\partial\theta}\frac{\partial\eta}{\partial\theta}-\frac{\partial\Phi}{\partial
z}=0,$ (3b)

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 (3a), $\kappa\left(\eta\right)$ 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}\Omega_{d}^{2}/\left(2g\right)$ and, $\Omega=\Omega_{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=\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 $\pi/2$ static contact angle is maintained
when the contact line elevation changes dynamically.

##  3 Linear solution and second-order straightforward asymptotic expansion

In order to enlighten the limitations of the expansion procedure developed by
Reclari et al. (2014), 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,

$\Phi=\Phi_{0}+\epsilon\Phi_{1}+\epsilon^{2}\Phi_{2}+\text{O}\left(\epsilon^{3}\right),$
(4a)
$\eta=\eta_{0}+\epsilon\eta_{1}+\epsilon^{2}\eta_{2}+\text{O}\left(\epsilon^{3}\right),$
(4b)

together with the further assumption of small driving forcing amplitudes of
order $\text{O}\left(\epsilon\right)$, i.e. $f=\epsilon F$, with $\epsilon$ a
small parameter $\epsilon\ll 1$. Solution
$\mathbf{q}_{0}=\left(\Phi_{0},\eta_{0}\right)$ represents the rest state,
which has a potential velocity field null everywhere, $\Phi_{0}=0$, and a flat
interface, $\eta_{0}=0$, as the contact angle is here assumed to be
$\vartheta=\pi/2$. Substituting the expansions above in equations (2)-(3b), a
series of system at the various order in $\epsilon$ is obtained. At leading
order, equations (2)-(3b) reduce to a forced linear system, whose matrix
compact form reads,

$\left(\partial_{t}\mathcal{B}-\mathcal{A}\right)\mathbf{q}_{1}=\boldsymbol{\mathcal{F}}_{1},$
(5)

with $\mathbf{q}_{1}=\left\\{\Phi_{1},\eta_{1}\right\\}^{T}$,
$\boldsymbol{\mathcal{F}}_{1}=F\left\\{0,r/2\right\\}^{T}e^{\text{i}\left(\Omega
t-\theta\right)}+c.c.=F\boldsymbol{\hat{\mathcal{F}}}_{1}^{F}e^{\text{i}\left(\Omega
t-\theta\right)}+c.c.$ and

$\mathcal{B}=\begin{pmatrix}0&0\\\ I_{\eta}&0\\\
\end{pmatrix},\mathcal{A}=\begin{pmatrix}\Delta&0\\\
0&-I_{\eta}+\frac{1}{Bo}\frac{\partial\kappa\left(\eta\right)}{\partial\eta}\\\
\end{pmatrix},$ (6)

where $c.c.$ stands for complex conjugate,
$\partial\kappa\left(\eta\right)/\partial\eta$ represents the first order
variation of the curvature associated with the small perturbation
$\epsilon\eta_{1}$ and $I_{\eta}$ is the identity matrix associated with the
interface $\eta$. Note that the kinematic condition does not explicitly appear
in (6), but it is enforced as a boundary condition at the interface (Viola et
al., 2018). In the limit of zero external forcing, i.e. $F=0$, system (5) is a
linear homogeneous problem which, by seeking for solutions having the
following normal form

$\hat{\mathbf{q}}^{mn}\left(r,z\right)e^{\text{i}\left(\omega_{mn}t-m\theta\right)}+c.c.,$
(7)

reduces to the classic generalized eigenvalue problem for inviscid
capillary–gravity waves

$\left(\text{i}\omega_{mn}\mathcal{B}-\mathcal{A}_{m}\right)\hat{\mathbf{q}}^{mn}=\mathbf{0},$
(8)

where indices $\left(m,n\right)$ 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 notice that the operator
$\mathcal{A}_{m}$ is complex, since $\theta$ derivatives produce $-\text{i}m$
terms. An exact analytical solution to equation (8) can be readily obtained
via a Bessel-Fourier-series representation leading to the well-known
dispersion relation (Lamb, 1993)

$\omega_{mn}^{2}=\left(k_{mn}+k_{mn}^{3}/Bo\right)\tanh{\left(k_{mn}H\right)},$
(9)

with $H=h/R$ and where the wavenumbers $k_{mn}$ is given by the nth-root of
the first derivative of the mth-order Bessel function of the first kind
satisfying $J^{\prime}_{m}\left(k_{mn}\right)=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.
(2018). 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 section §4.2. One must note that when (8) 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.

$m=0:\ \ \ \ \ \frac{\partial\hat{\eta}^{mn}}{\partial
r}=\frac{\partial\hat{\Phi}^{mn}}{\partial r}=0,$ (10a) $|m|\geq 1:\ \ \ \ \ \
\ \ \ \ \hat{\eta}^{mn}=\hat{\Phi}^{mn}=0.$ (10b)

It is also useful to note that owing to the symmetries of the problem, system
(8) is invariant under the transformation

$\left(\hat{\mathbf{q}}^{mn},+m,\text{i}\omega_{mn}\right)\longrightarrow\left(\hat{\mathbf{q}}^{mn},-m,\text{i}\omega_{mn}\right).$
(11)

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 (9) 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 r.h.s. of equation (5). In
contradistinction with the cases of unidirectional forcing (Miles, 1984a, b),
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

$\mathbf{q}_{1}=F\hat{\mathbf{q}}_{1}^{F}e^{\text{i}\left(\Omega
t-\theta\right)}+c.c.,$ (12)

with $\hat{\mathbf{q}}_{1}^{F}$ being the solution of the following forced
problem

$\left(\text{i}\Omega\mathcal{B}-\mathcal{A}_{1}\right)\hat{\mathbf{q}}_{1}^{F}=\boldsymbol{\hat{\mathcal{F}}}_{1}^{F}.$
(13)

The response structure $\hat{\mathbf{q}}_{1}^{F}$ is here computed
numerically, but, in practice, it is formally equivalent to that obtained
analytically by Reclari et al. (2014) 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.

Figure 2: 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. (2014) 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 $\Omega$ 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}$). Colors denote different wave conditions. Black solid
lines: linear potential model solution, from (12), computed by solving
numerically equation (13). Red solid lines: weakly nonlinear solution close to
the $\Omega\approx\omega_{21}/2$, obtained by computing (17). 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. (2014) to make the
equations non-dimensional.

Noting that $\epsilon F=f=d_{s}\Omega^{2}/\left(2g\right)$, in figure 2 the
linear solution $\epsilon\mathbf{q}_{1}^{F}$ from (12) is shown (black solid
lines) and compared with experimental measurements reported by Reclari et al.
(2014) in terms of maximum non-dimensional crest-to-trough contact line
amplitudes, $\tilde{\delta}=\delta R/D$, with
$\delta\left(\theta,t\right)=\eta\left(r=1,\theta,t\right)$. 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 will be therefore ignored thereafter. As discussed by Reclari et al.
(2014) and reproduced here, the linear solution describes well the
single–crest (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 $\propto
1/\left(\omega_{1n}^{2}-\Omega^{2}\right)$ 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 $\Omega\approx\omega_{11}$.
Introduction of viscous dissipation would regularize the divergent behaviour
at $\Omega=\omega_{11}$, however, in absence of any nonlinear restoring term,
the hardening nonlinearity displayed in figure 2 cannot be retrieved.
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 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.
(2014) instead). Among these super-harmonics, the double–crest (DC) wave
dynamics, occurring at a driving frequency $\Omega\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 favored
by the azimuthal symmetry of the external forcing. Reclari et al. (2014)
tentatively described such a double–crest dynamics by pursuing the asymptotic
analysis up to the second order in $\epsilon$, as in equations (4a) and (4b),
in order to account for second order system weak nonlinearities.
At the second order in $\epsilon$, one obtains the following forced linear
system,

$\left(\partial_{t}\mathcal{B}-\mathcal{A}\right)\mathbf{q}_{2}=\boldsymbol{\mathcal{F}}_{2}=F^{2}\left(\boldsymbol{\hat{\mathcal{F}}}_{2}^{FF}e^{\text{i}\left(2\Omega
t-2\theta\right)}+c.c.\right)+F^{2}\boldsymbol{\hat{\mathcal{F}}}_{2}^{F\overline{F}},$
(14)

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 Ibrahim (2005) and Reclari (2013), among others, for a full
derivation up to the second and third order). The bar denotes the complex
conjugate. Also note that amplitude $F$ is actually a real quantity and in the
following the superscript ${}^{\overline{F}}$ will be 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 r.h.s.
of equation (14) clearly shows how second order terms naturally induce a
super-harmonic response, whose spatial periodicity is $m=2$, hence precisely
corresponding to the double–crest dynamics experimentally observed. The second
forcing term on the r.h.s. of (14) has $\omega=0$ and $m=0$, i.e. it is steady
and axisymmetric. It originates in the leading order contribution in the time
and azimuthal averaged flow, the so-called mean flow. Equation (14) was solved
analytically by Reclari et al. (2014) by retaining for convenience only two
modes, namely those with $\left(m,n\right)=\left(2,1\right)$ and
$\left(0,1\right)$, expected to be the relevant ones. The numerical scheme
employed in this work allows us to effortlessly account for all the
$\left(2,n\right)$ and $\left(0,n\right)$ modes simultaneously, as their
contribution will be directly encompassed in the spatial function
$\hat{\mathbf{q}}_{2}^{FF}$ and $\hat{\mathbf{q}}_{2}^{F\overline{F}}$,
appearing in the second order solution,

$\mathbf{q}_{2}=\left(F^{2}\hat{\mathbf{q}}_{2}^{FF}e^{\text{i}\left(2\Omega
t-2\theta\right)}+c.c.\right)+F^{2}\hat{\mathbf{q}}_{2}^{F\overline{F}},$ (15)

whose contributions are computed by solving the following systems

$\left(\text{i}2\Omega\mathcal{B}-\mathcal{A}_{2}\right)\hat{\mathbf{q}}_{2}^{FF}=\hat{\boldsymbol{\mathcal{F}}}_{2}^{FF},\
\ \ \ \ \
-\mathcal{A}_{0}\hat{\mathbf{q}}_{2}^{F\overline{F}}=\hat{\boldsymbol{\mathcal{F}}}_{2}^{F\overline{F}}$
(16)

The total flow field, obtained through the asymptotic model is then given by
the sum of the first and second order solutions,

$\mathbf{q}=\left(f\hat{\mathbf{q}}_{1}^{F}e^{\text{i}\left(\Omega
t-\theta\right)}+f^{2}\hat{\mathbf{q}}_{2}^{FF}e^{\text{i}\left(2\Omega
t-2\theta\right)}+c.c.\right)+f^{2}\hat{\mathbf{q}}_{2}^{F\overline{F}},$ (17)

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\left(\theta,t\right)=\eta\left(r=1,\theta,t\right)$ is shown in figure
2 for driving frequencies close to $\Omega/\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 double–crest 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 $\propto 1/\left(\omega_{2n}^{2}-4\Omega^{2}\right)$ and
the total solution tends to diverge close to the double–crest super-harmonic
at $\omega_{21}/2$.
Such a close-to-resonance divergent behaviour is actually expected
(particularly in absence of any form of dissipation) when performing
straightforward asymptotic expansions, as they typically suffer from secular
(or resonating) terms that must be properly treated (see Castaing (2005) and
Nayfeh (2008) among many other references).

##  4 Weakly nonlinear analysis via multiple timescale 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 weakly nonlinear analysis (WNL) based on the multiple timescale
method. With the aim to derive a weakly nonlinear amplitude equation governing
the double–crest dynamics (DC), we first tackle the simpler problem of
single–crest waves (SC). In both cases we look for a third order asymptotic
solution of the system

$\mathbf{q}=\left\\{\Phi,\eta\right\\}^{T}=\epsilon\mathbf{q}_{1}+\epsilon^{2}\mathbf{q}_{2}+\epsilon^{3}\mathbf{q}_{3}+\text{O}\left(\epsilon^{4}\right),$
(18)

where the zero order solution, $\mathbf{q}_{0}=\mathbf{0}$, is omitted.

###  4.1 Single–crest dynamics (SC)

In §3 the forcing amplitude $f$ was assumed of order $\epsilon$, thus leading
to a linear first order problem directly forced by the external shaking, which
produces a divergent response close to harmonic resonances. With regards to
single–crest 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,
$\Omega\approx\omega_{1n}$, and in absence of dissipation, even a small
forcing will induce a large system response. The following analysis is
therefore expected to hold for $\Omega=\omega_{1n}+\lambda$, where $\lambda$
is a small detuning parameter, here assumed of order $\epsilon^{2}$. Lastly,
in the spirit of the multiple scale technique, we introduce the slow time
scale $T_{2}=\epsilon^{2}t$, 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:

$f=\epsilon^{3}F,\ \ \ \ \ \lambda=\epsilon^{2}\Lambda,\ \ \ \ \
T_{2}=\epsilon^{2}t,$ (19)

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.
Although the asymptotic expansion is here pursued up to the third order in
$\epsilon$, the procedure of the weakly nonlinear analysis is essentially
equivalent to that of the straightforward asymptotic analysis discussed in §3.
The major difference lies in the solution form of the leading order problem
that is now a homogenous problem, as in equation (8). Given the azimuthal
periodicity of the external forcing, among all possible natural eigenmodes we
assume a leading order solution as

$\mathbf{q}_{1}=A_{1}\left(T_{2}\right)\hat{\mathbf{q}}_{1}^{A_{1}}e^{\text{i}\left(\omega_{1n}t-\theta\right)}+c.c.,$
(20)

where $\hat{\mathbf{q}}_{1}^{A_{1}}$ is the eigenmode (computed by solving
(8)) associated with $\left(m,n\right)=\left(1,n\right)$ and $\omega_{1n}$ is
the corresponding natural frequency (solution of (9)).
The complex amplitude $A_{1}$, 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{\mathbf{q}}_{1}^{A_{1}}$
and introduces a new arbitrariness in the problem, which must be fixed at a
higher order. Eigen-surface, $\hat{\eta}_{1}^{A_{1}}$, and eigen-potential
field, $\hat{\Phi}_{1}^{A_{1}}$, computed for $\omega_{1n}=\omega_{11}$, are
shown in figure 3(a) and (b), respectively.
By pursuing the expansion to the second order, a linear system forced by the
first order solution and analogous to that of equation (14) is obtained (see
Reclari (2013) for the full expansion of the original nonlinear governing
equation up to the second order). Nevertheless, the forcing terms on the
r.h.s. are here proportional to $A_{1}^{2}$ (super- or second-harmonic) and to
$A_{1}\overline{A}_{1}$ (mean flow). Thus, we seek for a second order solution
of the form

$\mathbf{q}_{2}=A_{1}\overline{A}_{1}\hat{\mathbf{q}}_{2}^{A_{1}\overline{A}_{1}}+\left(A_{1}^{2}\hat{\mathbf{q}}_{2}^{A_{1}A_{1}}e^{\text{i}\left(2\omega_{1n}t-2\theta\right)}+c.c.\right),$
(21)

Figure 3: (a), (b) and (c): real part of the first, $\hat{\eta}_{1}^{A_{1}}$,
and second order, $\hat{\eta}_{2}^{A_{1}A_{1}}$ and
$\hat{\eta}_{2}^{A_{1}\overline{A_{1}}}$, free surface deformations computed
for $\omega_{1n}=\omega_{11}$. (d), (e) and (f): imaginary part of the
associated first order, $\hat{\Phi}_{1}^{A_{1}}$, and second order,
$\hat{\Phi}_{2}^{A_{1}A_{1}}$ and $\hat{\Phi}_{2}^{A_{1}\overline{A_{1}}}$,
potential velocity field. Each response is denoted by its amplitude
dependence, $\epsilon A_{1}$, $\epsilon^{2}A_{1}\overline{A}_{1}$ and
$\epsilon^{2}A_{1}A_{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{\Phi}_{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 center 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}$, $D=0.287\,\text{m}$ and
$\tilde{H}=h/D=0.52$, for which $Bo=2\,802.8$ and $\omega_{11}=1.3286$.

with $\hat{\mathbf{q}}_{2}^{A_{1}\overline{A}_{1}}$ and
$\hat{\mathbf{q}}_{2}^{A_{1}A_{1}}$ computed numerically and displayed in
figure 3(c)-(d) and (e)-(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
$\left(\omega,m\right)=\left(2\omega_{1n},2\right)$ and $\left(0,0\right)$ do
not correspond to eigenvalues of (8), i.e. the second order operators
$\left(\text{i}2\omega_{1n}\mathcal{B}-\mathcal{A}_{2}\right)$ and
$-\mathcal{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 $\pi/2$ phase shift, meaning that the first order eigenmode can be
normalized such that the eigen-surface is purely real, whereas the eigen-
potential is purely imaginary. Under these conditions, the mean flow forcing
term on the r.h.s. 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 (3a)) and its conditions at $r=0$ and $r=1$ (both
$\partial\hat{\eta}_{2}^{A_{1}\overline{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 Bremer & Breivik (2018) 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., 2017), is
intrinsically a viscous concept. In other words, viscous boundary layers must
be reintroduced if one aims to account for it.
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}$,

$\displaystyle\left(\partial_{t}\mathcal{B}-\mathcal{A}_{m}\right)\mathbf{q}_{3}=\mathcal{F}_{3}=$
(22)
$\displaystyle=-\partial_{T_{2}}A_{1}\mathcal{B}\hat{\mathbf{q}}_{1}^{A_{1}}e^{\text{i}\left(\omega_{1n}t-\theta\right)}+|A_{1}|^{2}A_{1}\boldsymbol{\hat{\mathcal{F}}}_{3}^{A_{1}\overline{A}_{1}A_{1}}e^{\text{i}\left(\omega_{1n}t-\theta\right)}+F\boldsymbol{\hat{\mathcal{F}}}_{3}^{F}e^{\text{i}\Lambda
T_{2}}e^{\text{i}\left(\omega_{1n}t-\theta\right)}$
$\displaystyle+\text{N.R.T.}+c.c.,$

with $\boldsymbol{\hat{\mathcal{F}}}_{3}^{F}=\left\\{0,r/2\right\\}^{T}$ and
where N.R.T. 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 equation (22). Secularity
results from all resonant forcing terms in $\mathcal{F}_{3}$, i.e. all terms
sharing the same frequency and wavenumber $\left(\omega_{1n},1\right)$ of
$\mathbf{q}_{1}$, and in effect all terms explicitly written in (22). It
follows that a compatibility condition must be enforced through the Fredholm
alternative (Friedrichs, 2012). Such a compatibility condition imposes the
amplitude $B=\epsilon A_{1}e^{\text{i}\lambda t}$ to obey the following normal
form

$\frac{dB}{dt}=-\text{i}\lambda
B+\text{i}\,\mu_{{}_{SC}}f+\text{i}\,\nu_{{}_{SC}}|B|^{2}B,$ (23)

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}\Lambda$. Moreover, by considering that
$\mathbf{q}=\epsilon A_{1}\hat{\mathbf{q}}_{1}^{A_{1}}+\ldots$, the small
implicit parameter $\epsilon$ is eliminated by defining the total physical
amplitude $A=\epsilon A_{1}$. The subscript SC stands for single–crest (SC).
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{\mathbf{q}}_{1}^{A_{1}\dagger}$,
associated with $\hat{\mathbf{q}}_{1}^{A_{1}}$, and the third order resonant
forcing terms as follows

$\text{i}\,\mu_{{}_{SC}}=\frac{<\hat{\mathbf{q}}_{1}^{A_{1}\dagger},\mathcal{B}\hat{\mathcal{F}}_{3}^{F}>}{<\hat{\mathbf{q}}_{1}^{A_{1}\dagger},\mathcal{B}\hat{\mathbf{q}}_{1}^{A_{1}}>}=\frac{\int_{z=0}r\overline{\hat{\eta}}_{1}^{A_{1}\dagger}/2\,r\text{d}r}{\int_{z=0}\left(\hat{\eta}_{1}^{A_{1}\dagger}\hat{\Phi}_{1}^{A_{1}}+\hat{\Phi}_{1}^{A_{1}\dagger}\hat{\eta}_{1}^{A_{1}}\right)\,r\text{d}r},$
(24a)
$\text{i}\,\nu_{{}_{SC}}=\frac{<\hat{\mathbf{q}}_{1}^{A_{1}\dagger},\mathcal{B}\boldsymbol{\hat{\mathcal{F}}}_{3}^{A_{1}\overline{A}_{1}A_{1}}>}{<\hat{\mathbf{q}}_{1}^{A_{1}\dagger},\mathcal{B}\hat{\mathbf{q}}_{1}^{A_{1}}>}=\frac{\int_{z=0}\left(\hat{\eta}_{1}^{A_{1}\dagger}\hat{\mathcal{F}}_{3_{\text{dyn}}}^{A_{1}\overline{A}_{1}A_{1}}+\hat{\Phi}_{1}^{A_{1}\dagger}\hat{\mathcal{F}}_{3_{\text{kin}}}^{A_{1}\overline{A}_{1}A_{1}}\right)\,r\text{d}r}{\int_{z=0}\left(\hat{\eta}_{1}^{A_{1}\dagger}\hat{\Phi}_{1}^{A_{1}}+\hat{\Phi}_{1}^{A_{1}\dagger}\hat{\eta}_{1}^{A_{1}}\right)\,r\text{d}r}.$
(24b)

Here
$\hat{\mathbf{q}}_{1}^{A_{1}\dagger}=\overline{\hat{\mathbf{q}}}_{1}^{A_{1}}$,
since the inviscid problem is self–adjoint with respect to the Hermitian
scalar product
$<\mathbf{a},\mathbf{b}>=\int_{\Sigma}\overline{\mathbf{a}}\cdot\mathbf{b}\,\text{d}V$,
with $\mathbf{a}$ and $\mathbf{b}$ two generic vector (see Viola et al. (2018)
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}\overline{A}_{1}A_{1}}$ is omitted,
as it only involves straightforward calculations, i.e. a Taylor expansion of
nonlinear governing equations and boundary conditions (2)-(3b) around the rest
state $\mathbf{q}_{0}=\mathbf{0}$. Here we simply denote with the subscript
dyn and kin the forcing components appearing in the dynamic and kinematic
boundary condition, respectively.
By turning to polar coordinates, $B=|B|e^{\text{i}\Theta}$, splitting the
modulus and phase parts of (23) and looking for stationary solution, $d/dt=0$
with $|B|\neq 0$, the following implicit relation is obtained,

$\tilde{d}_{s}\Omega^{2}\mp\frac{\left(\lambda-\nu_{{}_{SC}}|B|^{2}\right)|B|}{\mu_{{}_{SC}}}=0,$
(25)

or, in a more common polynomial form,

$P\left(|B|\right)=|B|^{3}-\frac{\lambda}{\nu_{{}_{SC}}}|B|\pm\frac{\mu_{{}_{SC}}\tilde{d}_{s}\Omega^{2}}{\nu_{{}_{SC}}}=0,$
(26)

where $f=\tilde{d}_{s}\Omega^{2}$, $\lambda=\Omega-\omega_{1n}$ and the $\mp$
signs correspond to the phases $\Theta=0$ and $\pi$, respectively. The two
branches prescribed by (25) for $|B|$ as a function of $\Omega$ 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 $\Theta$, the total single–crest wave
solution is reconstructed as

$\mathbf{q}_{SC}=\left\\{\Phi,\eta\right\\}^{T}=\mathbf{q}_{1}+\mathbf{q}_{2}.$
(27)

####  4.1.1 Experiments vs weakly nonlinear prediction: wave amplitude

In figure 4(a) and (b) the weakly nonlinear (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 (12). In comparison to 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
$\left(m,n\right)=\left(1,1\right)$ and $\left(1,2\right)$ 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 weakly nonlinear regime, hence suggesting the little
relevance of dissipative effects attributable to viscosity in this regime.
However, one must note that in this weakly nonlinear approach the driving
frequency is essentially fixed around that of a unique non-axisymmetric
natural mode, $\Omega\approx\omega_{1n}$. Consequently, when performing the
analysis for a mode $\left(1,n\right)$, the influence of all other modes is
completely overlooked.

$\left(m,n\right)$ | $\tilde{H}=h/D$ |  $D$ $\left[\text{m}\right]$ | $\omega_{mn}$ | $\mu_{SC}$ | $\nu_{SC}$
---|---|---|---|---|---
$\left(1,1\right)$ | 0.50 | 0.287 |  $\omega_{11}=\,$1.3239 | 0.27665 | 1.5265
$\left(1,2\right)$ | 0.50 | 0.287 |  $\omega_{12}=\,$2.3206 | 0.04189 | 17.0246
$\left(1,1\right)$ | 0.52 | 0.287 |  $\omega_{11}=\,$1.3286 | 0.27761 | 1.4847
Table 1: Value of the amplitude equation coefficient $\mu_{SC}$ and $\nu_{SC}$
in the conditions of figure 4. Figure 4: Comparison of the WNL-MTS prediction
with experiments in terms of maximum crest-to-trough contact line amplitude
(non-dimensional), $\Delta\tilde{\delta}$. (a) Black filled pentagons
correspond to the experimental measurements presented by Reclari (2013) for
$\tilde{d}_{s}=d_{s}/D=0.01$, $\tilde{H}=h/D=0.5$ and $D=0.287\,\text{m}$,
where the first two non-axisymmetric mode with
$\left(m,n\right)=\left(1,1\right)$ and $\left(1,2\right)$ are detected. Solid
black lines: solution of the linear potential model according to (12). Colored
lines: WNL-MTS prediction (27). The unstable branch is represented as a dashed
line. The values of normal form coefficient $\mu_{SC}$ and $\nu_{SC}$ computed
for $\Omega\approx\omega_{11}$ and $\omega_{12}$ are given in table 1. (b)
Same as (a) with the black filled circles corresponding to the case of figure
2 with $\tilde{d}_{s}=0.02$, $\tilde{H}=0.52$ and $D=0.287\,\text{m}$ and for
mode $\left(1,1\right)$, $\Omega\approx\omega_{11}$.

In consequence, the accuracy of the asymptotic solution rapidly deteriorates
moving away from harmonic resonances, when compared to the linear solution
(12), 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
$\Omega\in\left[0.7\omega_{12},0.9\omega_{12}\right]$ in figure 4(a). In other
words, the detuning parameter should be small in order for the present weakly
nonlinear analysis, based on a single mode expansion, to hold. On this regard,
as no other natural frequencies are encountered for $\Omega<\omega_{11}$, an
exception is made for the left 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 $\Omega\approx 0$, i.e. there is no need to employ a
leading order multimodal expansion.

####  4.1.2 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, 1958; Bauer, 1966; Dodge, 2000), where nonlinearities are of
crucial importance. The most popular driven nonlinear mass-spring model is
that developed by Duffing (1918), who added a cubic nonlinear spring
deformation (cubic term) to the classical driven harmonic oscillator

$\ddot{x}+2\sigma\dot{x}+x+c_{3}x^{3}=p\cos{\Omega t},$ (28)

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 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 (1973)
showed via asymptotic expansion of the potential flow solution in the
neighborhood 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 equation (23)). 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 (2021)
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 to predict
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 (Bäuerlein & Avila, 2021; Cenedese & Haller, 2020). We
note that, by analogy with the undamped Duffing equation, the weakly nonlinear
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
$\pi$ (see Appendix A for further comments on this regard). Nevertheless, one
should notice that this intrinsic limitation turns to be unimportant in cases
as those of figure 4, where for increasing amplitude a wave breaking
eventually occurs and the weakly nonlinear theory as well as the Duffing
mechanical analogy no longer apply.

###  4.2 Double–crest dynamics (DC)

We now tackle the double–crest (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 homogenous and a particular solution. We remind that the
double–crest dynamics in figure 2 occurs at a driving frequency
$\Omega\approx\omega_{21}/2$. Results at the end of this section will be
presented for mode $\left(2,1\right)$, for which experiments are available,
however, for the sake of generality, we formalize the analysis for any mode
$\left(2,n\right)$, i.e. $\Omega=\omega_{2n}/2+\lambda$, where $\lambda$ is
the 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 $\Omega$ close to $\omega_{21}/2$. One can see that
approaching $\Omega\approx\omega_{21}/2$ from lower frequencies, the
double–crest wave emerges on the top of a single–crest dynamics, i.e a
single–to–double wave transition takes place, with the latter being correctly
described by the linear solution, which still behaves well as $\omega_{21}/2$
is far enough from the 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 (5),
with
$\boldsymbol{\mathcal{F}}_{1}=F\left\\{0,r/2\right\\}^{T}e^{\text{i}\left(\Omega
t-\theta\right)}+c.c.=F\boldsymbol{\hat{\mathcal{F}}}_{1}^{F}e^{\text{i}\left(\omega_{2n}/2t-\theta\right)}e^{\text{i}\epsilon\Lambda
t}+c.c.\,$, with $f=\epsilon F$ and $\lambda=\epsilon\Lambda$.
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 (13), and (ii) a
homogeneous solution, represented by the natural mode
$\left(m,n\right)=\left(2,n\right)$, obtained by solving the generalized
eigenvalue problem (8), up to an amplitude to be determined at higher orders.
The second order resonating term will then require, in the spirit of multiple
timescale analysis, an additional second order solvability condition,
complementing the third order non-resonance condition already obtained in the
single–crest wave weakly nonlinear 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
weakly nonlinear multiple scale expansion for double–crest waves are the
following:

$f=\epsilon F,\ \ \ \ \ \lambda=\epsilon\Lambda,\ \ \ \ \ T_{1}=\epsilon t,\ \
\ \ \ T_{2}=\epsilon^{2}t,$ (29)
$\mathbf{q}_{1}=A_{2}\left(T_{1},T_{2}\right)\hat{\mathbf{q}}_{1}^{A_{2}}e^{\text{i}\left(\omega_{2n}t-2\theta\right)}+F\hat{\mathbf{q}}_{1}^{F}e^{\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}e^{\text{i}\Lambda
T_{1}}+c.c.,$ (30)

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{\mathbf{q}}_{1}^{F}$ is computed from (13) for $\Omega=\omega_{2n}/2$.
The second order linearized forced problem reads

$\left(\partial_{t}\mathcal{B}-\mathcal{A}_{m}\right)\mathbf{q}_{2}=\mathcal{F}_{2}=\mathcal{F}_{2}^{i,j}-\frac{\partial
A_{2}}{\partial
T_{1}}\mathcal{B}\hat{\mathbf{q}}_{1}^{A_{2}}e^{\text{i}\left(\omega_{2n}t-2\theta\right)}-\text{i}\Lambda
F\mathcal{B}\hat{\mathbf{q}}_{1}^{\Lambda
F}e^{\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}e^{\text{i}\Lambda
T_{1}}+c.c.\,.$ (31)

The first order solution is made of four different contributions of amplitude
$A_{2}$, $\overline{A}_{2}$, $F$ and $\overline{F}$, therefore it generates 10
different second order forcing terms, here denoted by
$\boldsymbol{\mathcal{F}}_{2}^{i,j}$, which exhibit a certain frequency and
azimuthal periodicity, $\left(\breve{\omega},\breve{m}\right)$.

| $\epsilon A_{2}$ | $\epsilon F$ | $\epsilon^{2}A_{2}A_{2}$ | $\epsilon^{2}\Lambda F$ | $\epsilon^{2}A_{2}\overline{A}_{2}$ | $\epsilon^{2}F\overline{F}$ | $\epsilon^{2}A_{2}F$ | $\epsilon^{2}A_{2}\overline{F}$
---|---|---|---|---|---|---|---|---
$\breve{m}$ | 2 | 1 | 4 | 1 | 0 | 0 | 3 | 1
$\breve{\omega}$ | $\omega_{2n}$ |  $\omega_{2n}$/2 | $2\omega_{2n}$ | $\omega_{2n}/2$ | 0 | 0 | 3$\omega_{2n}$/2 |  $\omega_{2n}$/2
Table 2: First order linear solutions and second order non-resonating forcing
terms gathered by their amplitude dependency and corresponding azimuthal and
temporal periodicity $\left(\breve{m},\breve{\omega}\right)$. Six terms have
been omitted as they are the complex conjugates of $\epsilon A_{2}$, $\epsilon
F$, $\epsilon^{2}A_{2}A_{2}$, $\epsilon^{2}FF$, $\epsilon^{2}A_{2}F$ and
$\epsilon^{2}A_{2}\overline{F}$.

The additional two forcing terms stem from the time-derivative of the first
order solution (30) with respect to the first order slow time scale $T_{1}$.
In order to interpret the last term in (31), it is worth first noting that,
while the amplitude of the linear solution (12), computed at a generic driving
frequency, grows with $\Omega$ as
$F/\left(\omega_{11}^{2}-\Omega^{2}\right)=\tilde{d}_{s}\Omega^{2}/\left(\omega_{11}^{2}-\Omega^{2}\right)\propto\Omega^{2}/\left(\omega_{11}^{2}-\Omega^{2}\right)$,
in the weakly nonlinear model for double–crest waves, the amplitude of the
particular solution (30) is proportional to
$F/\left(\omega_{11}^{2}-\omega_{21}^{2}/4\right)=\tilde{d}_{s}\Omega^{2}/\left(\omega_{11}^{2}-\omega_{21}^{2}/4\right)\sim\Omega^{2}$,
since the driving frequency was frozen at $\Omega=\omega_{21}+\lambda$, with
the small detuning parameter, $\lambda$, contributing to modify its phase, but
not its amplitude. This leads to an increasing discrepancy between (12) and
the leading order particular solution (30) away from the super-harmonic
resonance. The response to the forcing term proportional to $\Lambda F$ in
(31) 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 $\Lambda
F\propto\tilde{d}_{s}\Omega^{2}\left(\Omega-\omega_{2n}/2\right)$ 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 (31) 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}e^{\text{i}\left(\omega_{2n}-2\theta\right)}e^{\text{i}2\Lambda
T_{1}}+c.c.\,$. To avoid secular terms, a second order compatibility condition
is imposed, requiring that the following normal form is verified

$\frac{\partial A_{2}}{\partial
T_{1}}=\text{i}\,\mu_{{}_{DC}}F^{2}e^{\text{i}2\Lambda T_{1}},$ (32)

with $\mu_{{}_{DC}}$ computed as before, i.e.

$\text{i}\,\mu_{{}_{DC}}=\frac{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{2_{\text{dyn}}}^{FF}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{2_{\text{kin}}}^{FF}\right)\,r\text{d}r}{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\Phi}_{1}^{A_{2}}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\eta}_{1}^{A_{2}}\right)\,r\text{d}r},$
(33)

Taken alone, the dynamics resulting from (32) is however of little relevance,
since the solution, i.e. the frequency-response curve, would still diverge
(symmetrically) to infinity for $\Lambda=\Omega-\omega_{2n}/2\rightarrow 0$ in
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 (32) in the forcing expression, equation (31) can be
rearranged as follows

$\displaystyle\left(\partial_{t}\mathcal{B}-\mathcal{A}_{m}\right)\mathbf{q}_{2}=\boldsymbol{\mathcal{F}}_{2_{NRT}}^{i,j}+\boldsymbol{\mathcal{F}}_{2_{RT}}^{i,j}=$
(34)
$\displaystyle=\boldsymbol{\mathcal{F}}_{2_{NRT}}^{i,j}+F^{2}\left(\boldsymbol{\hat{\mathcal{F}}}_{2}^{FF}-\text{i}\,\mu_{{}_{DC}}\mathcal{B}\hat{\mathbf{q}}_{1}^{A_{2}}\right)e^{\text{i}\left(\omega_{2n}t-2\theta\right)}e^{\text{i}2\Lambda
T_{1}}+c.c.\,,$

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 $\Lambda 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{\mathbf{q}}_{1}^{A_{2}\dagger}=\overline{\hat{\mathbf{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

$\displaystyle\mathbf{q}_{2}=A_{2}\overline{A}_{2}\hat{\mathbf{q}}_{2}^{A_{2}\overline{A}_{2}}+F^{2}\hat{\mathbf{q}}_{2}^{F\overline{F}}+$
$\displaystyle+A_{2}^{2}\hat{\mathbf{q}}_{2}^{A_{2}A_{2}}e^{\text{i}\left(2\omega_{2n}t-4\theta\right)}+\Lambda
F\hat{\mathbf{q}}_{2}^{\Lambda
F}e^{\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}e^{\text{i}\Lambda
T_{1}}+$
$\displaystyle+A_{2}F\hat{\mathbf{q}}_{2}^{A_{2}F}e^{\text{i}\left(\left(3\omega_{2n}/2\right)t-3\theta\right)}e^{\text{i}\Lambda
T_{1}}+A_{2}\overline{F}\hat{\mathbf{q}}_{2}^{A_{2}\overline{F}}e^{\text{i}\left(\left(\omega_{2n}/2\right)t-\theta\right)}e^{-\text{i}\Lambda
T_{1}}+$
$\displaystyle+F^{2}\hat{\mathbf{q}}_{2}^{FF}e^{\text{i}\left(\omega_{2n}t-2\theta\right)}e^{\text{i}2\Lambda
T_{1}}+c.c.\,.$ (35)

Figure 5: Real part of the first, (a) $\hat{\eta}_{1}^{A_{2}}$ and (e)
$\hat{\eta}_{1}^{F}$, and non-resonating second order, (b)
$\hat{\eta}_{2}^{A_{1}\overline{A}_{2}}$, (c) $\hat{\eta}_{2}^{A_{2}A_{2}}$,
(d) $\hat{\xi}_{2}^{\Lambda F}$, (f) $\hat{\eta}_{2}^{F\overline{F}}$, (g)
$\hat{\eta}_{2}^{A_{2}F}$ and (h) $\hat{\eta}_{2}^{A_{2}\overline{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{\Phi}_{1}^{A_{2}}$ is purely imaginary. Note that the second
order mean flow $\hat{\eta}_{2}^{F\overline{F}}$ constantly induces an upside
down bell-like an axisymmetric interface deformation pushing the free surface
downward at the center of the moving container, by analogy with the effect
produced by $\hat{\eta}_{2}^{A_{1}\overline{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 pushes the interface upward at
the wall (same as $\hat{\eta}_{2}^{F\overline{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.

All non-resonant responses in (4.2.1) 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.
$\left(\text{i}\omega_{2n}\mathcal{B}-\mathcal{A}_{2}\right)$, is singular,
the value of the normal form coefficient (33) ensures that a non-trivial
solution for $\hat{\mathbf{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 et al., 2016).
Another possible approach is given in Appendix A of Meliga et al. (2012),
where a two-step regularization procedure, involving an intermediate factious
solution for $\hat{\mathbf{q}}_{2}^{FF}$ is employed. We also note that in
(4.2.1), exactly as in (21), a second order homogeneous solution has not been
accounted for as its introduction would be irrelevant to the final solution.
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
double– and single–crest wave, respectively, can be identified. Moreover, we
note that the second order response proportional to $\epsilon^{2}\Lambda 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$.
Lastly, at third order in $\epsilon$, the problem reads

$\displaystyle\left(\partial_{t}\mathcal{B}-\mathcal{A}_{m}\right)\mathbf{q}_{3}=\boldsymbol{\mathcal{F}}_{3}=$
(36) $\displaystyle=-\frac{\partial A_{2}}{\partial
T_{2}}\mathcal{B}\hat{\mathbf{q}}_{1}^{A_{2}}e^{\text{i}\left(\omega_{2n}t-2\theta\right)}-\text{i}2\Lambda
F^{2}\mathcal{B}\hat{\mathbf{q}}_{2}^{FF}e^{\text{i}\left(\omega_{2n}-2\theta\right)}e^{\text{i}2\Lambda
T_{1}}+$
$\displaystyle+|A_{2}|^{2}A_{2}\boldsymbol{\hat{\mathcal{F}}}_{3}^{A_{2}\overline{A}_{2}A_{2}}e^{\text{i}\left(\omega_{2n}t-2\theta\right)}+A_{2}F^{2}\boldsymbol{\hat{\mathcal{F}}}_{3}^{A_{2}F\overline{F}}e^{\text{i}\left(\omega_{2n}t-2\theta\right)}+$
$\displaystyle+\Lambda F^{2}\boldsymbol{\hat{\mathcal{F}}}_{3}^{\Lambda
FF}e^{\text{i}\left(\omega_{2n}-2\theta\right)}e^{\text{i}2\Lambda
T_{1}}+\text{N.R.T.}+c.c.\,,$

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. By noticing that the second and last forcing
terms share the same amplitude dependence, i.e. $\Lambda F^{2}$, they can be
recast into a single forcing term, say $\Lambda
F^{2}\hat{\mathcal{F}}_{3}^{\Lambda
FF}e^{\text{i}\left(\omega_{2n}-2\theta\right)}e^{\text{i}2\Lambda
T_{1}}+c.c.\,$.
Once again, all terms explicitly written in (36) are resonant, as they share
the same pair $\left(\omega_{2n},2\right)$ than the first order homogeneous
solution, hence a third order compatibility condition, leading to the
following normal form, must be enforced

$\frac{\partial A_{2}}{\partial T_{2}}=\text{i}\zeta_{{}_{DC}}\Lambda
F^{2}e^{\text{i}2\Lambda
T_{1}}+\text{i}\,\chi_{{}_{DC}}A_{2}F^{2}+\text{i}\,\nu_{{}_{DC}}|A_{2}|^{2}A_{2},$
(37)

with

$\text{i}\,\zeta_{{}_{DC}}=\frac{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{3_{\text{dyn}}}^{\Lambda
FF}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{3_{\text{kin}}}^{\Lambda
FF}\right)\,r\text{d}r}{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\Phi}_{1}^{A_{2}}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\eta}_{1}^{A_{2}}\right)\,r\text{d}r},$
(38a)
$\text{i}\,\chi_{{}_{DC}}=\frac{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{3_{\text{dyn}}}^{A_{2}F\overline{F}}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{3_{\text{kin}}}^{A_{2}F\overline{F}}\right)\,r\text{d}r}{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\Phi}_{1}^{A_{2}}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\eta}_{1}^{A_{2}}\right)\,r\text{d}r},$
(38b)
$\text{i}\,\nu_{{}_{DC}}=\frac{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{3_{\text{dyn}}}^{A_{2}\overline{A}_{2}A_{2}}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\mathcal{F}}_{3_{\text{kin}}}^{A_{2}\overline{A}_{2}A_{2}}\right)\,r\text{d}r}{\int_{z=0}\left(\hat{\eta}_{1}^{A_{2}\dagger}\hat{\Phi}_{1}^{A_{2}}+\hat{\Phi}_{1}^{A_{2}\dagger}\hat{\eta}_{1}^{A_{2}}\right)\,r\text{d}r}.$
(38c)

As a last step in the derivation of the final amplitude equation for the
double–crest (DC) waves and in order to eliminate the implicit small parameter
$\epsilon$, we unify (32) and (37) 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$ and $\lambda=\epsilon\Lambda$, and total
amplitude, $A=\epsilon A_{2}$. This is achieved by summing (32) and (37) along
with their respective weights $\epsilon^{2}$ and $\epsilon^{3}$, thus
obtaining

$\frac{dB}{dt}=-\text{i}\,\left(2\lambda-\chi_{{}_{DC}}f^{2}\right)B+\text{i}\,\left(\zeta_{{}_{DC}}\lambda+\mu_{{}_{DC}}\right)f^{2}+\text{i}\,\nu_{{}_{DC}}|B|^{2}B,$
(39)

where the change of variable $A=Be^{\text{i}2\lambda t}$ has been introduced
for convenience. As in §4.1, by turning to polar coordinates,
$B=|B|e^{\text{i}\Theta}$, splitting the modulus and phase parts of (39) and
looking for stationary solution, $d/dt=0$ with $|B|\neq 0$, the following
implicit relation is obtained,

$\tilde{d}_{s}\Omega^{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,$
(40)

where only the real solutions corresponding to $f=\tilde{d}_{s}\Omega^{2}>0$
are retained, as the combinations $\tilde{d}_{s}\Omega^{2}<0$ are not
physically meaningful.
Although two more terms appear in (39) 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 (40), amplitude equation (39) is
reminiscent of that given in (23). Indeed, equation (39) contains essentially
three contributions,

$\lambda\leftrightarrow\left(2\lambda-\chi_{{}_{DC}}f^{2}\right),\ \ \
\mu_{{}_{SC}}f\leftrightarrow\left(\zeta_{{}_{DC}}\lambda+\mu_{{}_{DC}}\right)f^{2},\
\ \ \nu_{{}_{SC}}\leftrightarrow\nu_{{}_{DC}},$ (41)

in order, a detuning term (forcing amplitude dependent), an additive
(quadratic) forcing term (forcing frequency dependent) and the classic cubic
restoring term, respectively. 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 for DC waves and reconstructed as

$\mathbf{q}_{DC}=\left\\{\Phi,\eta\right\\}^{T}=\mathbf{q}_{1}+\mathbf{q}_{2},$
(42)

is compared in figures 6 and 7 with experiments from Reclari (2013) and
Reclari et al. (2014) (see also figure 2).

Figure 6: Weakly nonlinear (WNL) prediction for double–crest (DC) waves versus
experiments by Reclari et al. (2014) (reported in figure 2 of this paper) in
terms $\Delta\tilde{\delta}$ for the smallest container diameter
$D=0.144\,\text{m}$, for different shaking diameters and at a driving
frequency close to $\Omega\approx\omega_{21}/2$ ($\Omega/\omega_{11}=0.6576$).
Solid black lines: linear potential solution (12). Red solid lines:
straightforward asymptotic solution (17). Light blue solid and dashed lines:
stable and unstable branches, respectively, predicted by the WNL via (40). The
normal form coefficient values for this configurations are
$\chi_{{}_{DC}}=3.4572$, $\zeta_{{}_{DC}}=0.8708$, $\mu_{{}_{DC}}=0.1292$ and
$\nu_{{}_{DC}}=10.0181$. Figure 7: Same as figure 6 but for the largest
container diameter $D=0.287\,\text{m}$ and $\tilde{H}=0.50$ (see figure 4.19
of Reclari (2013)). The normal form coefficient values for this configurations
are $\chi_{{}_{DC}}=3.4718$, $\zeta_{{}_{DC}}=0.9059$, $\mu_{{}_{DC}}=0.1338$
and $\nu_{{}_{DC}}=9.8847$.

####  4.2.2 Experiments vs weakly nonlinear prediction: wave amplitude

In figures 6 and 7, the weakly nonlinear (WNL) prediction of double–crest (DC)
waves is quantitatively compared with the experimental measurements from
Reclari (2013) and Reclari et al. (2014) 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 transition from a single– to a
double–crest wave and the resulting finite amplitude saturation via hardening
nonlinear mechanism, thus remarkably narrowing the gap with experiments for
all the 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 a single–crest (SC) dynamics and eventually
a double–crest wave breaking occurs at larger frequencies, a jump-down
transition from DC to SC takes place by increasing $\Omega$ 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$, by increasing the
driving frequency at a fixed shaking diameter, e.g. $\tilde{d}_{s}=0.10$, the
free surface is more likely to undergo a wave breaking, rather than a jump-
down transition (see figures 6(b) and 7(c)). Obviously, the inviscid model
employed here is not capable to predict the so-called jump-down frequency. In
Appendix A, a heuristic viscous damping model is introduced to tentatively
overcome the beforehand mentioned 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 $\Omega\approx\omega_{11}$,
progressively deteriorates. This is particularly evident on the right stable
branch and can be ascribed to the fact that the asymptotic model is
essentially formalized for a fixed driving frequency, i.e.
$\Omega\approx\omega_{21}/2$, thus filtering out the existence of other
resonances. Nevertheless, owing to the assumption of detuning parameter of
order $\epsilon$, the second order correction to the leading order particular
solution guarantees a fairly good agreement in a relatively wide range of
frequency around $\omega_{21}/2$.

Figure 8: Left-panels: comparison of the dimensionless and phase-averaged wave
height measured at the wall, $\tilde{\delta}\left(\theta,\pi/\Omega\right)$,
(black circles) with the straightforward asymptotic solution rebuilt via
(3.14) (gray solid line) and the weakly nonlinear (WNL) solution for the
double–crest (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 (2013), except for panel
(c), which is provided in Reclari et al. (2014). Note that (b) the nonlinear
prediction has a very large amplitude. Right-panels: corresponding three-
dimensional free surface deformation, $\eta\left(r,\theta,\pi/\Omega\right)$,
reconstruct via (4.25). The single–to–double crest transition via hardening
nonlinearity is clearly visible moving from top to bottom, i.e. for increasing
frequency.

####  4.2.3 Experiments vs weakly nonlinear prediction: free surface
reconstruction

In figure 8, the weakly nonlinear (WNL) models for the double–crest waves (DC)
is compared versus the straightforward asymptotic prediction discussed in §3
and the experimental measurements for DC waves from Reclari (2013) and Reclari
et al. (2014). The direct quantitative comparison is here outlined in terms of
dimensionless and phase-averaged wave height measured at the sidewall,
$\tilde{\delta}\left(\theta\right)$.
We observe that, if at $\Omega/\omega_{21}=0.490$ both models match
satisfactorily the experimental points, as soon as $\Omega/\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 single–crest (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=\pi/2$ and
$3\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\left(r,\theta,\pi/\Omega\right)$, is reconstructed through (42) and
shown in the right-panels of figure 8, where, for increasing shaking
frequencies, the nonlinear transition from a nearly single–crest wave dynamics
to a double–crest wave dynamics is enlightened.

####  4.2.4 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, 1999; Kalmár-
Nagy & Balachandran, 2011), as those observed by Bäuerlein & Avila (2021) and
occasionally by Reclari et al. (2014), it cannot predict the period–halving
dynamics associated with the super-harmonic resonance investigated in this
paper. Therefore, in connection with §4.1.2, 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 weakly nonlinear analysis (WNL) as well as the straightforward asymptotic
model highlighted the crucial role of quadratic nonlinearities emerging at
second order and from which the double–crest (DC) dynamics stems. At the same
time, the WNL model enlightened 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.

$\ddot{x}+2\sigma\dot{x}+x+c_{2}x^{2}+c_{3}x^{3}=p\cos{\Omega t}.$ (43)

Equation (43), also commonly known as Helmholtz–Duffing equation, has wide
applications in engineering problems as those related to beams, plates and
shells subjected to an initial static curvature (Mirzabeigy et al., 2014;
Askari et al., 2011), whose governing equations are reconduced to a second-
order nonlinear ordinary equation with quadratic and cubic nonlinear terms (Ke
et al., 2010; Alijani et al., 2011; Fallah & Aghdam, 2011).
Among the diverse asymptotic solutions of (43) in different limits (Rega,
1995; Benedettini & Rega, 1989; Kovacic & Brennan, 2011), the most relevant to
our work is that of Benedettini & Rega (1989). Within the context of planar
nonlinear response of suspended elastic cables to an external excitation, they
derived an amplitude equation which concerns with the first or fundamental
super-harmonic excitation, i.e. $\Omega\approx 1/2$, of (43). Their weakly
nonlinear approach is detailed in Appendix C, with the additional assumption
of vanishing damping $\sigma=0$. Assuming
$2\Omega=1+\lambda=1+\epsilon\hat{\lambda}$, small nonlinearities,
$c_{2}=\epsilon\hat{c_{2}}$ and $c_{3}=\epsilon^{2}\hat{c}_{3}$, and
introducing two slow time scales, one obtains

$dD/dt=-\text{i}\,\left(\lambda+c_{5}f^{2}\right)D+\text{i}\,\left(1-\lambda/2\right)c_{2}f^{2}/2-\text{i}\,4c_{4}|D|^{2}D,$
(44)

with $C=De^{\text{i}\lambda t}$ and with the auxiliary coefficients $c_{4}$
and $c_{5}$ (both functions of $c_{2}$ and $c_{3}$) defined in Benedettini &
Rega (1989). By comparing term by term, the analogy with equation (39) 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 vis-à-vis quantitative comparison with a
single-degree-of-freedom mechanical model, equation (44) 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 (43)
driven super-harmonically.

##  5 Conclusion

With regards to orbital shaken cylindrical containers and, specifically, to
the careful experimental campaign reported in Reclari (2013) and Reclari et
al. (2014), a weakly nonlinear analysis (WNL) via multiple timescale method
was formalized in §4 in order to investigate diverse features of the steady
state free surface dynamics and, particularly, the double–crest (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 propose by Reclari et al. (2014) 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 single–crest (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 modeled (from both qualitative and quantitative
perspectives (Bäuerlein & Avila, 2021)) by a single degree of freedom (1dof)
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, 1973). 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 double–crest
wave dynamics and to the resulting single–to–double crest wave transition. 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 represent the suitable mechanical analogy, a
one-degree-of-freedom (1dof) mechanical oscillator having both quadratic and
cubic nonlinear terms, commonly referred to as Helmholtz-Duffing (HD)
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 HD
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, 1984;
Benedettini & Rega, 1989), and it was here proposed as 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 one to correctly
account 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
absence of viscous boundary layers, the weakly nonlinear time– and
azimuthal–averaged mean flow reduces to a free surface deformation only. This
is in stark contrast with existence of the so-called Eulerian mean flow
(Bremer & Breivik, 2018), also known as viscous streaming flow, typically
observed in experiments (Bouvard et al., 2017). 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 modeling of
the moving contact line condition.

##  Appendix A Heuristic damping model: jump–down frequency and DC dynamics
suppression at low driving amplitudes

In §4.2.2 the weakly nonlinear (WNL) model for double–crest (DC) waves was
compared with experimental measurements from Reclari (2013) and Reclari et al.
(2014) in terms 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
single–crest (SC) wave at $\Omega\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 by increasing $\Omega$ 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 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
(1957) and Miles (1967). Such an estimation is computed according to the
following formula

$\sigma=\frac{2k_{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],$
(45)

where the first term represents the bulk dissipation, whereas the second and
third terms are related to the dissipation occurring at the solid bottom and
sidewall, respectively. In equation (45), $H=h/R$ is the non-dimensional fluid
depth, $k_{mn}$ is the non-dimensional wavenumber associated with mode
$\left(m,n\right)$, $\omega_{mn}$ is the corresponding natural frequency
obeying to the dispersion relation (9) 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
$\left(2,n\right)$ (which leads the DC wave dynamics observed close to
$\Omega\approx\omega_{21}/2$), was derived. For mode $\left(2,1\right)$ 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}$,
$\nu=1\times 10^{-6}\,\text{m${}^{2}$}/\text{s}$, $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 et al., 1993) and numerical (Viola et al., 2018)
evidences, 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 (39), leading to

$\frac{dB}{dt}=-\left[\sigma+\text{i}\left(2\lambda-\chi_{{}_{DC}}f^{2}\right)\right]B+\text{i}\left(\zeta_{{}_{DC}}\lambda+\mu_{{}_{DC}}\right)f^{2}+\text{i}\nu_{{}_{DC}}|B|^{2}B.$
(46)

The stationary form of (46) can be rearranged in the following implicit form

$\left(2\lambda-\nu_{{}_{DC}}|B|^{2}-\chi_{{}_{DC}}f^{2}\right)|B|\pm\sqrt{f^{4}\left(\zeta_{{}_{DC}}\lambda+\mu_{{}_{DC}}\right)^{2}-\left(\sigma|B|\right)^{2}}=0,$
(47)

Figure 9: (a) Same case of figure 6-(a) with $\tilde{d}_{s}=0.07$. (b) same as
(a), but for $\tilde{d}_{s}=0.02$ (from figure 2), value at which the
double–crest dynamics does not manifest. Solid and dashed lines correspond to
stable and unstable branches, respectively, computed via the weakly nonlinear
analysis in the inviscid case and for different values of damping coefficient,
with $\sigma$ given by (45). Markers correspond to the experimental points
shown in figure 2 and extracted from Reclari et al. (2014).

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.
The case of figure 9(a) shows that the so-called jump-down frequency is
somewhere in between $\Omega\in\left[0.675,0.685\right]$. The damping value
produced by (45) 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
(45) does not involve any dissipation mechanism associated with the contact
line, i.e. contact line hysteresis (Miles, 1967; Cocciaro et al., 1993;
Dussan, 1979; Hocking, 1987; Keulegan, 1959; Kidambi, 2009; Viola et al.,
2018; Viola & Gallaire, 2018) or possible surface contamination (Henderson &
Miles, 1990, 1994). 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, 1954; Hocking, 1987; Kidambi,
2009) than that predicted by (45). Comparison 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 (1990) used an effective fluid
viscosity 3 times larger than the actual one. Recently, Bäuerlein & Avila
(2021) 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, 2009). 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
time larger, is therefore not surprising when the damping is computed via (45)
and contact line dissipation is 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 in the following.
Prediction (45) is only valid for free capillary–gravity waves, whereas
dissipation rates of forced wave motions are generally more complex. 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 (46), could be
interpreted as a sort of nonlinear damping (Douady, 1990),
$\left(\sigma+\text{Im}\left[\nu\right]|B|^{2}\right)$, whose contribution to
the overall dissipation mechanisms is expected to increase at larger wave
amplitudes, hence influencing the location of the jump-down frequency. 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 (Huh & Scriven, 1971; Davis, 1974; Lauga et
al., 2007; Navier, 1823; Viola & Gallaire, 2018).
Most importantly, the inviscid WNL model is not capable to describe 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 (2021) (for uni-directional sloshing waves in 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, that is worth to be commented, 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,
$\Omega\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 responses that follows satisfactorily the linear
solution (see figure 2) showing a single–crest dynamics ranging over the whole
frequency window, $\Omega/\omega_{11}\in\left[0,1\right]$, in agreement with
experimental evidences.

##  Appendix B Asymptotic harmonic solution of the undamped Duffing equation

By analogy with the weakly nonlinear analysis for harmonic single–crest wave
dynamics presented in §4.1, we look for an asymptotic solution of the undamped
Duffing equation

$\ddot{x}+x+c_{3}x^{3}=p\cos{\Omega t},$ (48)

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. $\Omega=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 solutions simply reads

$x_{0}=C\left(\hat{t}_{1}\right)e^{\text{i}t}+c.c.\,.$ (49)

with $C\left(\hat{t}_{1}\right)$ 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

$dD/dt=-\text{i}\,\lambda
D+\text{i}\,\left(-1/4\right)p+\text{i}\,\left(3c_{3}/2\right)|D|^{2}D,$ (50)

where the change of variable $C=De^{\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

$-1/4\ \leftrightarrow\ \mu_{{}_{SC}},\ \ \ \ \ 3c_{3}/2\ \leftrightarrow\
\nu_{{}_{SC}},$ (51)

one immediately recognizes that equation (23) has indeed the same structure of
the formal amplitude equation (50), thus suggesting that the continuous
sloshing system and the one-degree-of-freedom (1dof) 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 weakly nonlinear solution of
the Helmholtz–Duffing equation,

$\ddot{x}+x+c_{2}x^{2}+c_{3}x^{3}=p\cos{\Omega t},$ (52)

derived by Benedettini & Rega (1989) 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 equation (52) with $\sigma=0$ (undamped
oscillator), $2\Omega=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 a 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

$x_{1}=C\left(\hat{t}_{1},\hat{t}_{2}\right)e^{\text{i}t}+fe^{\text{i}\left(1/2\right)t}e^{\text{i}\left(\hat{\lambda}/2\right)\hat{t}_{1}}+c.c.,$
(53)

with $f=\left(2/3\right)p$ and $C\left(\hat{t}_{1},\hat{t}_{2}\right)$ 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\left(\hat{t}\right)$ must obey to the following normal forms

$\boxed{\epsilon^{1}}\text{ :}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \
dC/d\hat{t}_{1}=\text{i}\,\left(c_{2}/2\right)f^{2}e^{\text{i}\hat{\lambda}\hat{t}_{1}},$
(54a) $\boxed{\epsilon^{2}}\text{ :}\ \ \ \ \
dC/d\hat{t}_{2}=-\text{i}\,\hat{\lambda}\left(c_{2}/4\right)e^{\text{i}\hat{\lambda}\hat{t}_{1}}-\text{i}\,c_{5}f^{2}A-\text{i}\,4c_{4}|C|^{2}C,$
(54b)

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 (1989).
Combining (54a) and (54b) into a single amplitude equation (by summing the two
expression 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

$dD/dt=-\text{i}\,\left(\lambda+c_{5}f^{2}\right)D+\text{i}\,\left(1-\lambda/2\right)c_{2}f^{2}/2-\text{i}\,4c_{4}|D|^{2}D,$
(55)

with $C=De^{\text{i}\lambda t}$. Note that the procedure used in the
perturbation analysis above and outlined in Benedettini & Rega (1989) is in
fact equivalent to that followed in Nayfeh (1984) 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 (55) with those of (39), the
analogy is evident, thus suggesting that the actual inviscid sloshing dynamics
in the double–crest wave regime may be, at least qualitatively, described by
the undamped Helmholtz–Duffing equation (43) driven super-harmonically.

#### Acknowledgements

We acknowledge Mohamed 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.

## References

  * Alijani et al. (2011) Alijani, F., Bakhtiari-Nejad, F. & Amabili, M. 2011 Nonlinear vibrations of fgm rectangular plates in thermal environments. Nonlinear Dyn. 66 (3), 251–270.
  * Askari et al. (2011) Askari, H., Saadatnia, Z., Younesian, D., Yildirim, A. & Kalami-Yazdi, M. 2011 Approximate periodic solutions for the Helmholtz–Duffing equation. Comput. Math. Appl. 62 (10), 3894–3901.
  * Bauer (1966) Bauer, H. F. 1966 Nonlinear mechanical model for the description of propellant sloshing. AIAA 4 (9), 1662–1668.
  * Bäuerlein & Avila (2021) Bäuerlein, B. & Avila, K. 2021 Phase lag predicts nonlinear response maxima in liquid-sloshing experiments. J. Fluid Mech. 925.
  * Benedettini & Rega (1989) Benedettini, F. & Rega, G. 1989 Planar non-linear oscillations of elastic cables under superharmonic resonance conditions. J. Sound Vib. 132 (3), 353–366.
  * Benjamin & Ursell (1954) Benjamin, T. B. & Ursell, F. J. 1954 The stability of the plane free surface of a liquid in vertical periodic motion. Proc. R. Soc. A: Math. Phys. Eng. Sci. 225 (1163), 505–515.
  * Bouvard et al. (2017) Bouvard, J., Herreman, W. & Moisy, F. 2017 Mean mass transport in an orbitally shaken cylindrical container. Phys. Rev. Fluids 2 (8), 084801.
  * Bremer & Breivik (2018) Bremer, T. S. Van Den & Breivik, Ø. 2018 Stokes drift. Philos. Trans. Royal Soc. A PHILOS T R SOC A 376 (2111).
  * Büchs (2001) Büchs, J. 2001 Introduction to advantages and problems of shaken cultures. Biochem. Eng. J. 7 (2), 91–98.
  * Büchs et al. (2000a) Büchs, J., Maier, U., Milbradt, C. & Zoels, B. 2000a Power consumption in shaking flasks on rotary shaking machines: I. power consumption measurement in unbaffled flasks at low liquid viscosity. Biotechnol. Bioeng. 68 (6), 589–593.
  * Büchs et al. (2000b) Büchs, J., Maier, U., Milbradt, C. & Zoels, B. 2000b Power consumption in shaking flasks on rotary shaking machines: Ii. nondimensional description of specific power consumption and flow regimes in unbaffled flasks at elevated liquid viscosity. Biotechnol. Bioeng. 68 (6), 594–601.
  * Case & Parkinson (1957) Case, K. M. & Parkinson, W. C. 1957 Damping of surface waves in an incompressible liquid. J. Fluid Mech. 2 (2), 172–184.
  * Castaing (2005) Castaing, B. 2005 Hydrodynamics and nonlinear instabilities, , vol. 3. Cambridge University Press.
  * Cenedese & Haller (2020) Cenedese, M. & Haller, G. 2020 How do conservative backbone curves perturb into forced responses? a melnikov function analysis. Proc. R. Soc. A: Math. Phys. Eng. Sci. 476 (2234), 20190494\.
  * Cocciaro et al. (1993) Cocciaro, B., Faetti, S. & Festa, C. 1993 Experimental investigation of capillarity effects on surface gravity waves: non-wetting boundary conditions. J. Fluid Mech. 246, 43–66.
  * Davis (1974) Davis, S. H. 1974 On the motion of a fluid-fluid interface along a solid surface. J. Fluid Mech. 65 (1), 71–95.
  * Dodge (2000) Dodge, F. T. 2000 The new ”Dynamic Behavior of Liquids in Moving Containers”. Southwest Research Inst. San Antonio, TX.
  * Dodge et al. (1965) Dodge, F. T., Kana, D. D. & Abramson, H. N. 1965 Liquid surface oscillations in longitudinally excited rigid cylindrical containers. AIAA 3 (4), 685–695.
  * Douady (1990) Douady, S. 1990 Experimental study of the Faraday instability. J. Fluid Mech. 221, 383–409.
  * Duffing (1918) Duffing, G. 1918 Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedeutung. Vieweg.
  * Dussan (1979) Dussan, E. B. 1979 On the spreading of liquids on solid surfaces: static and dynamic contact lines. Annu. Rev. Fluid Mech. 11 (1), 371–400.
  * Fallah & Aghdam (2011) Fallah, A. & Aghdam, M. M. 2011 Nonlinear free vibration and post-buckling analysis of functionally graded beams on nonlinear elastic foundation. Eur. J. Mech. A. Solids 30 (4), 571–583.
  * Faltinsen (1974) Faltinsen, O. M. 1974 A nonlinear theory of sloshing in rectangular tanks. J. Sh. Res. 18 (04), 224–241.
  * Faltinsen et al. (2016) Faltinsen, O. M., Lukovsky, I. A. & Timokha, A. N. 2016 Resonant sloshing in an upright annular tank. J. Fluid Mech. 804, 608–645.
  * Faltinsen et al. (2005) Faltinsen, O. M., Rognebakke, O. F. & Timokha, A. N. 2005 Resonant three-dimensional nonlinear sloshing in a square-base basin. part 2. effect of higher modes. J. Fluid Mech. 523, 199–218.
  * Faltinsen & Timokha (2009) Faltinsen, O. M. & Timokha, A. N. 2009 Sloshing. Cambridge University Press.
  * Friedrichs (2012) Friedrichs, K. O. 2012 Spectral theory of operators in Hilbert space. Springer Science & Business Media.
  * Handa-Corrigan et al. (1989) Handa-Corrigan, A., Emery, A. N. & Spier, R. E. 1989 Effect of gas—liquid interfaces on the growth of suspended mammalian cells: mechanisms of cell damage by bubbles. Enzyme Microb. Technol. 11 (4), 230–235.
  * Henderson & Miles (1990) Henderson, D. M. & Miles, J. W. 1990 Single-mode faraday waves in small cylinders. J. Fluid Mech. 213, 95–109.
  * Henderson & Miles (1994) Henderson, D. M. & Miles, J. W. 1994 Surface-wave damping in a circular cylinder with a fixed contact line. J. Fluid Mech. 275, 285–299.
  * Hocking (1987) Hocking, L. M. 1987 The damping of capillary–gravity waves at a rigid boundary. J. Fluid Mech. 179, 253–266.
  * Horstmann et al. (2021) Horstmann, G. M., Anders, S., Kelley, D. H. & Weier, T. 2021 Formation of spiral waves in cylindrical containers under orbital excitation. J. Fluid Mech. 925.
  * Horstmann et al. (2020) Horstmann, G. M., Herreman, W. & Weier, T. 2020 Linear damped interfacial wave theory for an orbitally shaken upright circular cylinder. J. Fluid Mech. 891.
  * Huh & Scriven (1971) Huh, C. & Scriven, L. E. 1971 Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. J. Colloid Interface Sci. 35 (1), 85–101.
  * Ibrahim (2005) Ibrahim, R. A. 2005 Liquid sloshing dynamics: theory and applications. Cambridge University Press.
  * Jesus et al. (2004) Jesus, M. J. De, Girard, P., Bourgeois, M., Baumgartner, G., Jacko, B., Amstutz, H. & Wurm, F. M. 2004 Tubespin satellites: a fast track approach for process development with animal cells using shaking technology. Biochem. Eng. J. 17 (3), 217–223.
  * Jordan & Smith (1999) Jordan, D. W. & Smith, P. 1999 Nonlinear ordinary differential equations: an introduction to dynamical systems, , vol. 2. Oxford University Press, USA.
  * Kalmár-Nagy & Balachandran (2011) Kalmár-Nagy, T. & Balachandran, B. 2011 Forced harmonic vibration of a Duffing oscillator with linear viscous damping. The duffing equation: nonlinear oscillators and their behaviour pp. 139–174.
  * Ke et al. (2010) Ke, L.-L., Yang, J. & Kitipornchai, S. 2010 An analytical study on the nonlinear vibration of functionally graded beams. Meccanica 45 (6), 743–752.
  * Keulegan (1959) Keulegan, G. H. 1959 Energy dissipation in standing waves in rectangular basins. J. Fluid Mech. 6 (1), 33–50.
  * Kidambi (2009) Kidambi, R. 2009 Capillary damping of inviscid surface waves in a circular cylinder. J. Fluid Mech. 627, 323–340.
  * Klöckner & Büchs (2012) Klöckner, W. & Büchs, J. 2012 Advances in shaking technologies. Trends Biotechnol. 30 (6), 307–314.
  * Kovacic & Brennan (2011) Kovacic, I. & Brennan, M. J. 2011 The Duffing equation: nonlinear oscillators and their behaviour. John Wiley & Sons.
  * Kretzmer & Schügerl (1991) Kretzmer, G. & Schügerl, K. 1991 Response of mammalian cells to shear stress. Appl. Microbiol. Biotechnol. 34 (5), 613–616.
  * Lamb (1993) Lamb, H. 1993 Hydrodynamics. Cambridge university press.
  * Lauga et al. (2007) Lauga, Eric, Brenner, Michael & Stone, Howard 2007 Microfluidics: The No-Slip Boundary Condition, pp. 1219–1240. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Liu & Hong (2001) Liu, C.-M. & Hong, L.-N. 2001 Development of a shaking bioreactor system for animal cell cultures. Biochem. Eng. J. 7 (2), 121–125.
  * Lukovsky (1990) Lukovsky, I. A. 1990 Introduction to nonlinear dynamics of a solid body with a cavity including a liquid. Kiev: Naukova dumka (in Russian) .
  * Maier et al. (2004) Maier, U., Losen, M. & Büchs, J. 2004 Advances in understanding and modeling the gas–liquid mass transfer in shake flasks. Biochem. Eng. J. 17 (3), 155–167.
  * McDaniel & Bailey (1969) McDaniel, L. E. & Bailey, E. G. 1969 Effect of shaking speed and type of closure on shake flask cultures. Appl. Microbiol. 17 (2), 286–290.
  * Meliga et al. (2012) Meliga, P., Gallaire, F. & Chomaz, J. M. 2012 A weakly nonlinear mechanism for mode selection in swirling jets. J. Fluid Mech. 699, 216–262.
  * Micheletti et al. (2006) Micheletti, M., Barrett, T., Doig, S. D., Baganz, F., Levy, M. S., Woodley, J. M. & Lye, G. J. 2006 Fluid mixing in shaken bioreactors: Implications for scale-up predictions from microlitre-scale microbial and mammalian cell cultures. Chem. Eng. Sci. 61 (9), 2939–2949.
  * Miles (1967) Miles, J. W. 1967 Surface-wave damping in closed basins. Proc. R. Soc. A: Math. Phys. Eng. Sci. 297, 459–475.
  * Miles (1984a) Miles, J. W. 1984a Internally resonant surface waves in a circular cylinder. J. Fluid Mech. 149, 1–14.
  * Miles (1984b) Miles, J. W. 1984b Resonantly forced surface waves in a circular cylinder. J. Fluid Mech. 149, 15–31.
  * Mirzabeigy et al. (2014) Mirzabeigy, A., Yazdi, M. K. & Nasehi, M. H. 2014 Approximate periodic solution for the nonlinear Helmholtz-Duffing oscillator via analytical approaches. Int. J. Comput. Math. 2014.
  * Moiseev (1958) Moiseev, N. N. 1958 On the theory of nonlinear vibrations of a liquid of finite volume. J. Appl. Math. Mech. 22 (5), 860–872.
  * Moisy et al. (2018) Moisy, F., Bouvard, J. & Herreman, W. 2018 Counter-rotation in an orbitally shaken glass of beer. EPL 122 (3), 34002.
  * Muller et al. (2007) Muller, N., Derouazi, M., Tilborgh, F. Van, Wulhfard, S., Hacker, D. L., Jordan, M. & Wurm, F. M. 2007 Scalable transient gene expression in chinese hamster ovary cells in instrumented and non-instrumented cultivation systems. Biotechnol. Lett. 29 (5), 703–711.
  * Muller et al. (2005) Muller, N., Girard, P., Hacker, D. L., Jordan, M. & Wurm, F. M. 2005 Orbital shaker technology for the cultivation of mammalian cells in suspension. Biotechnol. Bioeng. 89 (4), 400–406.
  * Narimanov (1957) Narimanov, G. S. 1957 Movement of a tank partly filled by a fluid: the taking into account of non-smallness of amplitude. Prikl. Math. Mech. (in Russian) 21, 513–524.
  * Narimanov et al. (1977) Narimanov, G. S., Dokuchaev, L. V. & Lukovsky, I. A. 1977 Nonlinear dynamics of flying apparatus with liquid. moscow: Mashinostroenie. (in Russian) .
  * Navier (1823) Navier, C. L. M. H. 1823 Mémoire sur les lois du mouvement des fluides. Mém. Acad. R. des Sci. Inst. France 6 (1823), 389–440.
  * Nayfeh (1984) Nayfeh, A. H. 1984 Quenching of primary resonance by a superharmonic resonance. J. Sound Vib. 92 (3), 363–377.
  * Nayfeh (2008) Nayfeh, A. H. 2008 Perturbation Methods. Wiley.
  * Ockendon & Ockendon (1973) Ockendon, J. R. & Ockendon, H. 1973 Resonant surface waves. J. Fluid Mech. 59 (2), 397–413.
  * Orchini et al. (2016) Orchini, A., Rigas, G. & Juniper, M. P. 2016 Weakly nonlinear analysis of thermoacoustic bifurcations in the rijke tube. J. Fluid Mech. 805, 523–550.
  * Papoutsakis (1991) Papoutsakis, E. T. 1991 Fluid-mechanical damage of animal cells in bioreactors. Trends Biotechnol. 9 (1), 427–437.
  * Raynovskyy & Timokha (2018) Raynovskyy, I. & Timokha, A. N. 2018 Steady-state resonant sloshing in an upright cylindrical container performing a circular orbital motion. Math. Probl. Eng. .
  * Reclari (2013) Reclari, M. 2013 Hydrodynamics of orbital shaken bioreactors. Tech. Rep.. EPFL.
  * Reclari et al. (2014) Reclari, M., Dreyer, M., Tissot, S., Obreschkow, D., Wurm, F. M. & Farhat, M. 2014 Surface wave dynamics in orbital shaken cylindrical containers. Phys. Fluids 26 (5), 052104\.
  * Rega (1995) Rega, G. 1995 Bifurcation and chaos in the Helmholtz-Duffing oscillator. In Bifurc. Chaos, pp. 191–215. Springer.
  * Tan et al. (2011) Tan, R.-K., Eberhard, W. & Büchs, J. 2011 Measurement and characterization of mixing time in shake flasks. Chem. Eng. Sci. 66 (3), 440–447.
  * Timokha & Raynovskyy (2017) Timokha, A. N. & Raynovskyy, I. 2017 The damped sloshing in an upright circular tank due to an orbital forcing. Dopov. Nats. Akad. Mauk. Ukr. 10, 48–53.
  * Tissot et al. (2010) Tissot, S., Farhat, M., Hacker, D. L., Anderlei, T., Kühner, M., Comninellis, C. & Wurm, F. M. 2010 Determination of a scale-up factor from mixing time studies in orbitally shaken bioreactors. Biochem. Eng. J. 52 (2-3), 181–186.
  * Tissot et al. (2011) Tissot, S., Oberbek, A., Reclari, M., Dreyer, M., Hacker, D. L., Baldi, L., Farhat, M. & Wurm, F. M. 2011 Efficient and reproducible mammalian cell bioprocesses without probes and controllers? New Biotechnol. 28 (4), 382–390.
  * Viola et al. (2018) Viola, F., Brun, P.-T. & Gallaire, F. 2018 Capillary hysteresis in sloshing dynamics: a weakly nonlinear analysis. J. Fluid Mech. 837, 788–818.
  * Viola & Gallaire (2018) Viola, F. & Gallaire, F. 2018 Theoretical framework to analyze the combined effect of surface tension and viscosity on the damping rate of sloshing waves. Phys. Rev. Fluids 3 (9), 094801\.
  * Wurm (2004) Wurm, F. M. 2004 Production of recombinant protein therapeutics in cultivated mammalian cells. Nat. Biotechnol. 22 (11), 1393–1398.
  * Zhang et al. (2009) Zhang, X., Bürki, C., Stettler, M., Sanctis, D. De, Perrone, M., Discacciati, M., Parolini, N., DeJesus, M., Hacker, D. L., Quarteroni, A. & Wurm, F. M. 2009 Efficient oxygen transfer by surface aeration in shaken cylindrical containers for mammalian cell cultivation at volumetric scales up to 1000 l. Biochem. Eng. J. 45 (1), 41–47.