SS symmetry
ArticleThe Stability of a Hydrodynamic Bravais Lattice
Miles M. P. Couchman 1 , Davis J. Evans 2
1

2*

and John W. M. Bush 2, *Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Cambridge CB3 0WA, UK; mc2277@cam.ac.uk
Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA;
djevans@mit.edu
Correspondence: bush@math.mit.edu

Abstract: We present the results of a theoretical investigation of the stability and collective vibrationsof a two-dimensional hydrodynamic lattice comprised of millimetric droplets bouncing on the surface
of a vibrating liquid bath. We derive the linearized equations of motion describing the dynamics of ageneric Bravais lattice, as encompasses all possible tilings of parallelograms in an infinite plane-fillingarray. Focusing on square and triangular lattice geometries, we demonstrate that for relatively low
driving accelerations of the bath, only a subset of inter-drop spacings exist for which stable lattices
may be achieved. The range of stable spacings is prescribed by the structure of the underlying
wavefield. As the driving acceleration is increased progressively, the initially stationary lattices
destabilize into coherent oscillatory motion. Our analysis yields both the instability threshold and the
wavevector and polarization of the most unstable vibrational mode. The non-Markovian nature of
the droplet dynamics renders the stability analysis of the hydrodynamic lattice more rich and subtle
than that of its solid state counterpart.
Keywords: bouncing droplets; Faraday waves; lattice instability; normal-mode analysis; phonons

Citation: Couchman, M.M.P.;
Evans, D.J.; Bush, J.W.M. The
Stability of a Hydrodynamic Bravais
Lattice. Symmetry 2022, 1, 0.
https://doi.org/
Academic Editor: Andrzej
Herczyński
Received: 20 June 2022
Accepted: 21 July 2022
Published: 26 July 2022
Publisher’s Note: MDPI stays neutral
with regard to jurisdictional claims inpublished maps and institutional affiliations.

Copyright: © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms andconditions of the Creative Commons
Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

1. Introduction
The collective vibrations of atoms within a crystal lattice, referred to as phonons,
are of fundamental interest in materials science and solid-state physics, and govern bulk
properties of the crystal such as the heat capacity, thermal and electrical conductivity,and elastic modulus [1–4]. Complementary insight into lattice dynamics is often gained
through consideration of macroscopic analog systems, including mass-and-spring networks [5,6] and more complex metamaterials [7,8]. Numerous hydrodynamic analogues ofcrystals have been explored, from Bragg’s packed bubbles at an air–liquid interface [9] tomore recent studies of vibrations in lattices of colloidal particles [10–14] and microfluidic
droplets [15–18]. We here consider theoretically the dynamics of a two-dimensional latticecomprised of millimetric droplets that bounce on the surface of a vibrating liquid bath andare coupled through an underlying Faraday wavefield.
A fluid bath vibrating vertically with acceleration γ sin(2π f t) will destabilize into
a subharmonic field of standing Faraday waves, characterized by period TF = 2/ f andwavelength λ F (as is related to ω F = π f through the standard water-wave dispersion
relation), when the vibrational acceleration γ exceeds a critical value γF known as the Faraday threshold [19,20]. Below γF , but above the bouncing threshold γB < γF , a millimetricdroplet may bounce indefinitely on the bath’s surface, generating a spatially extended,
temporally decaying wavefield at each impact [21,22]. As γ is gradually increased beyond
γB , the droplet’s bouncing period increases until eventually becoming twice that of thevibrational driving, thus achieving resonance with the Faraday wavefield and prompting asubstantial increase of the system’s wave energy. In this period-doubled bouncing regime,
droplets may bounce in either low- or high-energy states, as distinguished by bouncing
amplitude and denoted by (2, 1)1 and (2, 1)2 , respectively [23]. Such period-doubled drops

Symmetry 2022, 1, 0. https://doi.org/10.3390/sym1010000

https://www.mdpi.com/journal/symmetry

Symmetry 2022, 1, 0

2 of 20

may bounce either in phase or out of phase with respect to each other. For γ > γW , thewalking threshold, a single droplet destabilizes from stationary bouncing into steady horizontal motion, propelled through a resonant interaction with its own wavefield [24,25]. As
γ approaches γF , the persistence time of the waves increases and the droplet’s dynamics
become more strongly influenced by its past trajectory, and the system has heightened‘path memory’ [26,27]. The resulting non-Markovian nature of the droplet dynamics gives
rise to numerous features reminiscent of quantum systems; consequently, this system has
provided the basis for the burgeoning field of hydrodynamic quantum analogues [28–31].
Multiple droplets may organize into static and dynamic bound states by virtue of
their shared wavefield. Specifically, the static bound states of droplet pairs [32], free rings[33], and radially confined rings (modelling a one-dimensional periodic lattice) [34] haveinter-drop spacings related to the Faraday wavelength λ F . As the vibrational acceleration isincreased, droplet pairs destabilize into either in-line oscillations, or orbital or promenading(side by side) motion, depending on the droplet size and inter-drop distance [32]. Radiallyconfined rings destabilize into out-of-phase angular oscillations or propagating soliton-likewaves [34], and free rings exhibit additional radial vibrational modes [33]. The stability andresonant oscillations of forced chains with a free end have also been considered theoretically[35]. Pairs of identical walking drops may form dynamic bound states consisting of either
orbital [24,25,36,37] or promenading motion [38,39].
A variety of two-dimensional bound structures have been studied, including therotational [40] and translational instabilities of droplet aggregates [41]. Eddi et al. [42]constructed eight of the eleven possible Archimedean tilings of the plane, some of which
required tuning the relative bouncing phase of neighbouring droplets. Eddi et al. [43]considered square and triangular lattices and observed the emergence of coherent modes ofoscillation, a hydrodynamic analog of phonons, beyond a critical vibrational acceleration.
Only one lattice spacing was considered for each geometry and a full characterization of
the emergent oscillations as a function of the lattice spacing was not undertaken. Edgeeffects were seen to influence the observed lattice dynamics owing to their finite size. To
rationalize their experimental observations, Eddi et al. [43] proposed a phenomenological,
one-dimensional model in which each droplet was connected to its nearest neighbours via
an effective spring force proportional to the wave amplitude. Notably, their model did notinclude an explicit waveform or the influence of the system memory, which precluded aquantitative characterization of the lattice stability.
A detailed linear stability analysis of an infinite one-dimensional droplet lattice was
performed by Thomson et al. [44], using the stroboscopic model of Oza et al. [45], in order
to rationalize the observed dynamics of a periodic droplet chain [34]. In this configuration,
only certain lattice spacings were found to remain stable below the instability threshold ofa single drop γW , and the lattice subsequently destabilized via either super- or sub-critical
Hopf bifurcations as the driving acceleration was increased. This linear stability analysis
was then extended to investigate both weakly nonlinear oscillations and solitary waves [46].
We here build upon the work of Eddi et al. [43] and Thomson et al. [44] by developinga theoretical framework for studying the stability and dynamics of two-dimensional droplet
lattices, which one expects to exhibit a richer set of instabilities than their one-dimensionalcounterparts. Specifically, we consider the Bravais lattice, a theoretical construct usedin solid-state physics to describe regular crystalline structures. The defining feature of
the Bravais lattice is its discrete translational symmetry, which allows its lattice points tobe expressed as integer multiples of two basis vectors. The Bravais lattice thus appears
to be identical from each constituent lattice point [4]. In two dimensions, five possible
geometries satisfy this required symmetry, specifically, square, triangular, rectangular,
centred rectangular, and oblique lattices.
In Section 2, we use the theoretical model of Couchman et al. [32] to derive thedispersion relation governing the stability of a generic Bravais lattice. In Section 3, we focus
on the square and triangular geometries, which admit analytical simplifications due to theirrotational symmetry. For these geometries, we predict the stability threshold and mode

Symmetry 2022, 1, 0

3 of 20of vibrational instability that emerge as the bath’s vibrational acceleration is increased
progressively. Our theoretical predictions are compared to the experimental results of
Eddi et al. [43]. In Section 4, we summarize our results and propose future avenues ofinvestigation.
2. Normal Mode Analysis
In this section, we characterize the linear stability of a Bravais droplet lattice. In
Section 2.1, we review the variable-phase stroboscopic model of Couchman et al. [32] forthe dynamics of multiple interacting droplets, which provides the basis of our analysis.
Definitions of relevant variables and parameters are provided in Table 1. In Section 2.2, wethen characterize the base state of the Bravais lattice and, in Section 2.3, perturb the basestate to derive the dispersion relation governing the lattice’s stability and normal modesof vibration.
Table 1. Definitions of relevant variables and parameters.

SymbolDefinition

Fluidρ, σ, νef , TF = 2/ f , λ F , k F = 2π/λ Fγ, γW , γ∗ , γF

Density, surface tension, effective kinematic viscosity [47]Bath driving frequency, Faraday period, wavelength, wavenumber
Peak driving acceleration of bath, walking threshold of single drop, lattice instability
threshold, Faraday threshold

Trajectory equationx = ( x, y), t, g
R, m = 4πρR3 /3h( x, t)f (r )

Td = 1/ νe k2Fd
Me = T (1−Tγ/γ
F)qF
3νe TF mgk FA=
2π 3k2 σ +ρge2

Horizontal position, time, gravitational acceleration
Droplet radius, mass
Wave amplitude strobed at bouncing period TF
Wave kernel
Wave decay timescale [47]
Memory parameter
Wave-amplitude coefficientF

2π f ρν k
α = 2ν (1+2e2 ) , e = 3k2 σ+eρgF
eqFζ = k2 TFαMeF

D = 0.17mg
κ = Tm
FD
β=
S, C

q

Spatio-temporal damping coefficient, viscosity induced wavenumber correction [48]Non-dimensional spatial-damping coefficientρR
σ + 6πRµ a , µ a

mgATF k2F
D

Latticea, bdmn = ma + nb
k, ξ, ωk

Horizontal drag coefficient [47], air viscosity
Non-dimensional droplet mass
Non-dimensional wave-force coefficient
Impact phase parameters (see Appendix A) [32]Basis vectors defining geometry of Bravais lattice
Horizontal position of droplets in base lattice, (m, n) ∈ Z
Wave vector, polarization vector, complex-valued frequency of vibrational mode

2.1. Review of Theoretical Model
Throughout this work, we choose physical parameters for the drop and bath liquids
corresponding to silicone oil with kinematic viscosity ν = 20 cSt, density ρ = 949 kg m−3 ,and surface tension σ = 20.6 × 10−3 N m−1 , which have been widely used in experimental
and theoretical studies of droplet–droplet interactions [32,33,37,39]. In a deep bath vibrating
vertically at f = 80 Hz, such a fluid will be characterized by a Faraday wavelength
λ F ≈ 4.75 mm and Faraday threshold γF ≈ 4.25g, where g denotes the gravitational
acceleration. We assume that all droplets have radius R = 0.36 mm and bounce in phasewith each other at the Faraday period TF in the higher-energy (2, 1)2 mode, as assumed inSymmetry 2022, 1, 0

4 of 20

the prior analysis of free rings [33]. An example of a triangular droplet lattice constructedin the laboratory is shown in Figure 1a.

5 mm

Figure 1. (a) An oblique view of a triangular lattice of millimetric droplets bouncing on the surface ofa vertically vibrating liquid bath. The droplets are coupled by a shared wavefield, as is visualized byilluminating the bath with a striped pattern of coloured light [49]. We note that beneath each drop, thedrop’s reflection from the bath surface is visible. (b) A Bravais lattice is generated via enumerationof the points dmn = ma + nb, (m, n) ∈ Z, with point (m, n) = (2, 1) highlighted for illustration. The√
basis vectors describing the triangular lattice shown here are a = a(1, 0) and b = a(−1/2, 3/2),
where a has units of length and sets the lattice spacing. We study the response of the lattice toperturbations δx away from its base state.

The trajectory equation for walking droplets was developed by Moláček and Bush [47].
The ‘stroboscopic’ trajectory equation of Oza et al. [45] is a simplification thereof, so called
because it effectively eliminates the vertical droplet motion from consideration by averaging
the drop dynamics over a bouncing period, and so describes the horizontal dynamics visible
in the laboratory when the system is strobed at the Faraday period [49]. We here employ thevariable-phase stroboscopic model of Couchman et al. [32], an extension of the stroboscopic
model that accounts for variations in the phase of impact between the droplet and bath.
Consideration of such variations has been found to be necessary in order to rationalize
the observed stability of multi-drop systems [48], including bound droplet pairs [32] andrings [33].
The variable-phase stroboscopic model predicts that the horizontal positions, xmn , ofinteracting drops of mass m in the lattice evolve according tom ẍmn + D ẋmn = −mgCmn ∇h( xmn , t),(1)where overdots denote time derivatives, and the wavefield h is described byh( x, t) = A ∑pq

with wave kernel

Z t−∞

S pq f (k F | x − x pq (s)|)e−(t−s)/(TF Me ) ds,−2
f (r ) = J0 (r ) 1 + (K1 (ζr )ζr − 1)e−r
,

(2)

(3)

where J0 and K1 denote Bessel functions of the first kind and modified second kind,respectively.
In Equation (1), drop motion is driven by a wave force, proportional to the local
gradient of the underlying wavefield h, and resisted by a linear drag force with coefficientD. The wavefield, defined in Equation (2), is modelled as the superposition of waves ofspatial form f (r ) generated by each droplet along its past trajectory, as summed over all of
the droplets in the lattice. The wave kernel f (r ) is based on the wave model of Moláček
and Bush [47] but more accurately captures the experimentally observed far-field decayof the wavefield [50], and so the long-range interactions between droplets [32,48,50,51].
For r < 1, f (r ) is well approximated by J0 (r ), but then decays more rapidly than J0 (r )Symmetry 2022, 1, 0

5 of 20

for r > 1, as prescribed by the spatial damping factor ζ which is defined in terms of
fluid parameters in Table 1. The memory timescale TF Me ∼ (1 − γ/γF )−1 appearing inEquation (2) characterizes the temporal decay of the waves, with higher values of thememory parameter Me indicating waves that decay more slowly and thus have a greater
influence on the lattice’s evolution.
The phase factors S and C capture the phase shift between the resonant oscillations of
the drop and wavefield, which may vary as a function of the system parameters. Such aphase shift influences both the wave amplitude generated by each droplet at impact (ascaptured by S ) and the horizontal wave force imparted to the droplet (as captured by C ).
Couchman et al. [32] determined the dependence of these phase factors on the droplet
radius, local wave height, and vibrational acceleration. They also found that, for drops inthe (2, 1)2 bouncing mode [52] commonly used in experiments and treated here, neglecting
variations in these phase factors may lead to the prediction that bound states destabilize
below the walking threshold of a single droplet γW , a prediction at odds with experimental
observations of drop–drop interactions [32–34,37,39]. Further details concerning the phaseparameters S and C , and their explicit functional forms, may be found in Appendix A.
By introducing the following non-dimensional variables x̄ = k F x, h̄ = h/A, t̄ = t/TF
for horizontal position, wave height, and time, respectively, the governing Equations (1)–(2)take the following non-dimensional form:κ x̄¨ mn + x̄˙ mn = − βCmn ∇h̄( x̄mn , t̄),h̄( x̄, t̄) = ∑

Z t̄

pq

−∞

S pq f (| x̄ − x̄ pq (s̄)|)e−(t̄−s̄)/Me ds̄.

(4)
(5)

Expressions for the non-dimensional mass κ and waveforce coefficient β are given in Table 1.
For the remainder of the paper, we drop the overbars denoting non-dimensional variables
for the sake of notational simplicity.
2.2. Base State of the Bravais Lattice
The unperturbed Bravais lattice has droplets located at horizontal positions
dmn = ma + nb,(6)

where m and n are integers and a and b are the basis vectors that define the lattice geometry,commonly referred to as ‘primitive vectors’ (see Figure 1b). We first demonstrate that ageneric Bravais lattice is a stationary solution to Equation (4). Substituting xmn (t) = dmn
into Equation (5) yields the wavefield for the Bravais lattice,h( x, t) = Me ∑ S(γ, h(d pq )) f (| x − d pq |).

(7)

pq

We note that the lattice symmetry ensures that all drops encounter the same wave heighth(d pq ) = h0 , so their phase factors are likewise identical, S(h0 ) = S0 . Evaluating Equation (7) at the lattice points yields an implicit expression for the local wave amplitude,
h0 = Me S(γ, h0 ) ∑ f ( d pq ).

(8)

pq

Making use of the functional form for S presented in Equation (A2), one may solve theimplicit Equation (8) numerically to obtain h0 . Having obtained h0 , the phase parameters
S0 = S(γ, h0 ) and C0 = C(γ, h0 ) may be computed using Equations (A2) and (A3), respectively. In Section 3, it will be shown that the values S0 and C0 depend on the initialgeometry of the Bravais lattice, which in turn strongly influence the lattice stability.

Symmetry 2022, 1, 0

6 of 20

Having solved for the base state wavefield, it is immediately evident that xmn (t) = dmn
satisfies Equation (4) when ∇h(dmn , t) = 0, which may be written explicitly as
d pq∑ f 0 (|d pq |) |d pq | = 0.

(9)

pq

Physically, Equation (9) signifies that the net waveforce on each droplet must vanish, or
equivalently that the slope of the wavefield beneath each droplet is zero. Noting thatd(− p,−q) = −d( p,q) , the pairwise terms ( p, q) and (− p, −q) cancel in Equation (9), leaving
only the ( p, q) = (0, 0) term which vanishes because the wave kernel is even and so
f 0 (0) = 0 (see Equation (3)). It is thus apparent that, by virtue of its translational symmetry,
any Bravais lattice is a stationary solution of the trajectory equation.
2.3. Dispersion Relation of the Perturbed Lattice
We may now assess the stability of the Bravais lattice by studying the growth of small
perturbations to the base state, xmn (t) = dmn + εδxmn (t), where ε  1. Substituting thisperturbation into Equation (4), using the equilibrium condition (9), and expanding to firstorder in ε, yields the following linearized equations of motion governing the perturbations
δxmn (t),¨ mn + δx
˙ mn = − ϕ ∑
κ δx
pq

Z t
−∞



H f d pq · δxmn (t) − δx pq (s) e−(t−s)/Me ds,

(10)

where
ϕ = β S0 C0 ,(11)

and H f denotes the Hessian matrix of the radially symmetric wave kernel f (Equation (3)),∂2 f
2H f ( x) =  ∂∂x2 f
∂x∂y



∂2 f
∂x∂y 
.
∂2 f
∂y2 x



(12)

The elements of H f may be expressed explicitly as
∂2 f
( x)
∂x2
∂2 f
( x)
∂y2
∂2 f
( x)
∂x∂y

=
=
=

f 00 (r ) 2
f 0 (r ) 2
x
+
y ,
r2
r3
f 00 (r ) 2
f 0 (r ) 2
y
+
x ,
r2
r3
 00

f (r )
f 0 (r )
−
xy,
r2
r3

(13)
(14)
(15)

where x = ( x, y) and r = | x|.
In deriving Equation (10), we note that variations in the phase parameters around thebase values S0 , C0 are of O(ε2 ) and so can be neglected. To see this, note that S(γ, h(dmn + eδxmn )) = S γ, h(dmn ) + ε∇h(dmn ) · δxmn + O ε2 ,(16)
but ∇h(dmn ) = 0 in the base state (Equation (9)). Therefore, the impact phase parameters
only enter our linear analysis through their base values S0 , C0 , which influence the resulting
stability of the lattice through the parameter ϕ (Equation (11)) that scales the waveforce inEquation (10).
To derive the dispersion relation governing the linear stability of the lattice, we
substitute the following normal mode perturbation into Equation (10)δxmn = ξeik·dmn +ωk t + c.c.,(17)Symmetry 2022, 1, 0

7 of 20

where c.c. denotes the complex conjugate of the preceding exponential term. Equation (17)represents a plane wave characterized by wavevector k, polarization vector ξ, and complex frequency ωk . The real and imaginary components of ωk represent the mode’s
growth rate and oscillation frequency, respectively. Making use of the following inteRt
gral, −∞ eωs e−(t−s)/Me ds = Me eωt /(1 + Me ω ), we thus obtain the dispersion relation"|

κωk2 + ωkI + ϕMe H

(0)

H(k)
−
1 + Me ω k{z
Dk ( ωk )!#ξ = 0,(18)}

where I denotes the identity matrix andH(k) = 2 ∑ H f (d pq ) cos(k · d pq )

(19)

pq

may be recognized as the elementwise discrete cosine transform of the Hessian matrix
over the lattice. For k = 0, H(0) is the Hessian matrix of the net wavefield at the origin, as
characterizes the wavefield’s local curvature.
Non-trivial solutions (ξ 6= 0) to the dispersion relation (18) exist when
det Dk (ωk ) = 0,(20)which takes the form of a sextic polynomial in the complex-valued frequency ωk and may
be solved numerically. If a wavevector k exists such that any root ωk∗ of Equation (20) has
a positive real part, then the associated mode will grow and the lattice will destabilize. If
no such k exist, the lattice remains stable. For an initially stable lattice, as the vibrational
forcing γ is increased, the instability threshold γ∗ is reached when there exists at least onek∗ such that Re(ωk∗ ) ≥ 0. The polarization vector ξ ∗ associated with the unstable mode k∗
may then be found by solving Dk (ωk∗ )ξ ∗ = 0 (Equation (18)). We note that the polarization
vector ξ must be real-valued as Dk is a symmetric matrix.
It is noteworthy that, owing to the discrete translational symmetry of the Bravais lattice,not all wavevectors k produce physically distinguishable oscillations when substituted
into the normal mode ansatz (17). For example, consider
 a square lattice characterized by
dmn = a(m x̂ + nŷ) and two wavevectors k = k x , k y and k0 = k x + 2π/a, k y + 2π/a .
Substituting k0 into Equation (17) yields
δxmn = ξe2πi(m+n) eik·dmn +ωk t = ξeik·dmn +ωk t ,(21)

an oscillation characterized by wavevector k. Therefore, k0 and k result in the same physicaloscillation. It thus suffices to consider only wavevectors k in the lattice’s so-called Brillouinzone, defined as the smallest set of k required to describe all distinguishable vibrationsof the discrete lattice [4]. The Brillouin zones for the square and triangular geometries
considered in Section 3 are illustrated in Figure 2, and the procedure for generating theBrillouin zone for a generic Bravais lattice may be found in standard reference texts onsolid-state physics [4].

Symmetry 2022, 1, 0

8 of 20

Figure 2. The extent of the Brillouin zone in wavevector space (k x , k y ) for (a) square and (b) triangular
Bravais lattices with lattice spacing a [4]. Wavevectors k0 outside the Brillouin zone are paired withan equivalent k within the Brillouin zone such that k0 and k yield the same physical vibration once
substituted into the normal mode ansatz (17), a result of the discrete translational symmetry of thelattice. In characterizing lattice stability, it thus suffices to consider only modes with wavevectors kwithin the Brillouin zone.

3. Results
3.1. Square and Triangular Lattices
We first consider the stability of square and triangular lattices, as are characterized bya single inter-drop spacing a with droplets at the base state positions(square)dmn(triangle)dmn= a(m x̂ + nŷ),√ = a m x̂ + (n/2) − x̂ + 3ŷ ,

(22)
(23)

respectively. As is demonstrated in Appendix B, the four- and six-fold symmetry of thesquare and triangular geometries, respectively, result in the Hessian matrix for the net
wavefield H(0) (Equation (19), with k = 0) reducing to a scalar multiple of the identity
matrix I,
H(0) = η0 I,(24)
where η0 is the unique (repeated) eigenvalue of H(0) . Geometrically, Equation (24) signifies that the curvature of the local wavefield beneath each droplet in the square andtriangular base lattices is isotropic, having the same value in all horizontal directions.
Using Equation (24), the general dispersion relation (18) reduces to the simpler eigenvalue
problem,1
H(k) ξ =(25)(1 + Me ωk ) κωk2 + ωk + ϕη0 Me ξ,
ϕMewhich admits non-trivial solutions when ξ is an eigenvector of H(k) with corresponding
eigenvalue ηk . When such is the case, one may equate ηk with the bracketed term on theright-hand side of Equation (25), yielding




1
κ
P1 (ωk ) ≡ κωk3 + 1 +ωk2 ++ Me ϕη0 ωk + ϕ(η0 − ηk ) = 0,
(26)
Me
Me
which may be solved for the three ωk associated with each of the two eigenvalues ηk . Since
κ > 0 and P1 (ωk ) → ∞ as ωk → ∞, the intermediate value theorem guarantees a positive
real root whenever P1 (0) < 0, i.e., ϕ(η0 − ηk ) < 0. We thus deduce that a lattice is always
unstable if there are any wavevectors k such that either eigenvalue ηk of H(k) is greaterthan η0 .
Based on Equation (26), we now have a simple procedure for assessing whether asquare or triangular lattice is unstable at the initial vibrational acceleration γ = 0.7γF ,which is below the walking threshold of a single drop γW and corresponds approximately
to the lowest γ at which the droplets bounce in a period-doubled mode [52]. Specifically,

Symmetry 2022, 1, 0

9 of 20

we compute the eigenvalues η0 and ηk of H(k) (Equation (19)) numerically, noting that
there are two ηk for each wave vector k in the Brillouin zone, and plot their dependence
on the lattice spacing a for both the square (Figure 3a) and triangular (Figure 4a) latticegeometries. Intervals of a where ηk > η0 for all k are shaded in red, indicating that thelattice is already unstable at γ = 0.7γF . We observe that for both the square and triangular
geometries, there are only discrete intervals of a for which the lattice is initially stable,
which roughly correspond to geometries in which drops bounce in minima of the local
wave amplitude (see Figures 3e and 4e). We note that this behaviour is in accord withprevious studies reporting that bound states are most stable when each drop bounces ina minimum of the net wavefield produced by its neighbours [32,33]. Thomson et al. [44]found similar swaths of initial instability as a function of the lattice spacing for infinite
one-dimensional droplet chains at low vibrational accelerations, which they referred to as‘geometric instabilities’.
For the set of lattices that are initially stable at γ = 0.7γF , we may then determine theinstability threshold γ∗ > 0.7γF at which destabilization occurs. At γ∗ , the real part of ωkwill vanish and the imaginary part will correspond to the frequency of the destabilizing
oscillation. Substituting ωk = iΩ (Ω ∈ R) into Equation (26), separating the complex
polynomial into real and imaginary parts, and eliminating Ω yields

Pk ( Me ) ≡ ϕη0 Me3 + ϕκηk Me2 + Me + κ = 0.

(27)

Noting that κ > 0, ϕ > 0 and η0 < 0, we observe that Pk ( Me ) = 0 has at least one
positive real root, the minimum of which corresponds to the memory Me∗ at which thelattice destabilizes. Because Qk ( Me ) ≡ Pk ( Me ) − CMe2 < Pk ( Me ) for every non-zero
Me provided C > 0, we infer that Qk must have a root at a lower Me than does Pk .
Furthermore, if there is a wavevector k0 such that ηk0 < ηk , then since Pk ( Me ) − Pk0 ( Me ) =ϕκMe2 (ηk − ηk0 ) > 0, the smallest real positive root of Pk0 ( Me ) must be smaller than that ofPk ( Me ). We thus conclude that the wave vector k that goes unstable at the lowest memory
value Me , corresponds to that k with the minimum eigenvalue ηk .
We may thus gain additional insight from the curves ηk plotted in Figures 3a and 4a forthe square and triangular geometries, respectively. Namely, in the initially stable intervals
of a, the lowermost curve ηk corresponds to the destabilizing vibrational wavemode kthat will emerge at the instability threshold γ∗ . In Figures 3 and 4, we plot the instability
thresholds γ∗ (panel b), and the magnitude and direction of the wavevector characterizing
this destabilizing mode (panels c and d, respectively) as a function of the lattice spacing a.
In all cases, we find that the polarization vector ξ is parallel to the wave vector k, signifyingthat the lattices always destabilize into longitudinal, as opposed to transverse, oscillations.
For both the square and triangular geometries, the instability threshold of the lattice,γ∗ , is almost always greater than the walking threshold γW of a single drop. The lattices
are most stable when the constituent drops bounce in the deepest minima of the wavefield
produced by their neighbours, as was the case for droplet pairs and rings [32,33]. This effectis a direct result of variations in the impact phase, with the product S0 C0 (see Equations(10) and (11)) decreasing with decreasing local wave amplitude (see Figures 3e and 4e) andthus reducing the horizontal waveforce exerted on each droplet at impact. Conversely,
we note that in their theoretical analysis of a one-dimensional lattice, Thomson et al. [44]predicted destabilization below γW , despite the fact that experimentally such lattices were
found to remain stable above γW [34]. This mismatch followed from their assumption ofa constant impact phase, and highlights the importance of accounting for variations in
the vertical dynamics by using a variable-impact-phase model when considering droplet–
droplet interactions. In Figures 3b and 4b we observe that γ∗ approaches γW in the limitof large a, consistent with the droplets in the lattice becoming effectively uncoupled fromtheir neighbours at sufficiently large distances.

Symmetry 2022, 1, 0

10 of 20

Figure 3. The square lattice. (a) The eigenvalues ηk of the matrix H(k) (19), governing the latticestability, are plotted as a function of the inter-drop spacing a, normalized by the Faraday wavelength
λ F . In intervals of a where η0 is greater than all ηk , the lattice is initially stable at γ/γF = 0.7.
Otherwise, the lattice is already unstable as denoted by red shading. Curves ηk for two common
modes of instability are highlighted. There are two eigenvalues ηk for a single mode k, denotedby dashed and solid lines of the same colour. We note that the two magenta curves are virtuallyindistinguishable. For a given a, the most negative ηk corresponds to the destabilizing wave mode.
(b) The instability threshold γ∗ , normalized by the walking threshold for a single drop γW , for
initially stable spacings a. The corresponding wavevector magnitude |k|, normalized by π/a (seeFigure 2a), and wave angle with respect to the x-axis, are shown for the destabilizing mode in panels(c,d), respectively. Green and magenta points indicate spacings a where one of the two wave modes
highlighted in panel (a) are found to be destabilizing. These modes are further illustrated in Figure 5.
(e) The black curve indicates the dependence on lattice spacing of the local wave amplitude h0
beneath each droplet in the base lattice, normalized by the drop radius R. The orange curve denotes
the corresponding product of impact phases S0 C0 that influences the lattice stability through ϕ inEquation (10).

In Figures 3c,d, we highlight two dominant modes of vibration for the square lattice,as are further illustrated in Figure 5. The first mode (green) corresponds to out-of-phaseoscillations of neighbouring lattice planes along either the x̂ or ŷ directions, which areSymmetry 2022, 1, 0

11 of 20

equivalent given the four-fold rotational symmetry of the square lattice. We note that Eddiet al. [43] observed a superposition of such oscillations in the x̂ and ŷ directions, resulting in
the appearance of each droplet exhibiting roughly circular orbits around its base point (seeFigure 4a of [43]). In magenta, we highlight a separate mode characterized√ by out-of-phaseoscillations along 45 degree planes in the lattice with a wavelength of 2a, corresponding
to the diagonal of a square cell. Apart from these two readily identifiable modes, the scaled
wavevector a|k|/π (Figure 3c) varies approximately linearly with the scaled lattice spacinga/λ F , thus representing an oscillation wavelength that is no longer solely governed by a,but is now also influenced by the Faraday wavelength λ F that characterizes the underlying
wavefield.

Figure 4. The triangular lattice. Quantities characterizing the stability of a triangular lattice areplotted as a function of the lattice spacing a, in the same manner as for the square lattice in Figure 3.
The dominant mode of instability is highlighted in magenta in panels (c,d), and is illustrated in
Figure 6.

In Figures 4c,d, we highlight the dominant mode of vibration for the triangular lattice(magenta) which, as illustrated in Figure 6, corresponds to out-of-phase oscillations ofneighbouring lattice planes along the π/6 direction (or equivalently, the π/6 + n(π/3),n = (1, 2, . . . , 5) directions given the six-fold rotational symmetry of the lattice). This mode

Symmetry 2022, 1, 0

12 of 20is consistent with that reported by Eddi et al. [43] (see Figure 4b of [43]). We note thatin their experiments, Eddi et al. [43] used different drop sizes and fluid parameters thanconsidered here, so we cannot present a quantitative comparison between their results andour theoretical predictions.
a)

b)

c)

d)

Figure 5. Common vibrational modes of the square lattice corresponding to longitudinal waves alonghigh-symmetry lattice directions. (a) An illustration of the k = (π/a, 0) wave mode highlighted bygreen markers in Figure 3, corresponding to the out-of-phase oscillations of neighbouring planes
along the x̂ direction with wavelength λ = 2a. (b) The maximum real part of the eigenvalues ηk areshown for a lattice spacing a/λ F = 3.40, with the most unstable modes k = (π/a, 0) and (0, π/a)marked by crosses. (c,d) Analogous plots for the wave mode k = (π/a, π/a), highlighted by magenta
markers in Figure 3, corresponding to the out-of-phase oscillations of neighbouring planes along the√
45 degree diagonal, with wavelength λ = 2a.

Symmetry 2022, 1, 0

13 of 20

a)

b)

√

Figure 6. Common vibrational mode of the triangular lattice. (a) An illustration of the k = √2π ( 23 , 12 )3awave mode highlighted by magenta markers in Figure 4, corresponding to the out-of-phase oscilla√tions of neighbouring planes along a line 30 degrees to the horizontal, with wavelength λ = 3a.
(b) The maximum real part of the eigenvalues ηk for the lattice spacing a/λ F = 2.5. The Brillouinzone boundaryis traced in black (see Figure 2b), and crosses indicate the most unstable modes:√k = √2π ( 23 , 12 ) and π/3 rotations thereof.
3a3.2. Geometric Instabilities in the Low-Memory Limit
We have demonstrated that, for the relatively low driving acceleration γ = 0.7γF ,square and triangular lattices are already unstable for certain values of the lattice spacing a.
Following Thomson et al. [44], we refer to such lattices as being ‘geometrically unstable’.
We proceed by demonstrating a method for predicting whether a generic Bravais lattice isgeometrically unstable on the basis of the shape of the lattice-induced wavefield. In the
low-memory limit Me  1, the dispersion relation (18) reduces to the eigenvalue problem




1  2
(0)
(k)
H −H
ξ= −
κωk + ωk ξ,
(28)
ϕMe


which admits non-trivial solutions when ξ is an eigenvector of the matrix H(0) − H(k)with corresponding eigenvalue αk . When such is the case, one may equate αk with thebracketed term on the right-hand side of Equation (28), yielding

P2 (ωk ) ≡ κωk2 + ωk + ϕMe αk = 0,(29)which may be solved to obtain the two ωk associated with each of the two eigenvalues αk .
We note that P2 (ωk ) → ∞ as ωk → ∞. In the case αk < 0, we also have P2 (0) < 0, and so
the intermediate value theorem guarantees the existence of a positive real root. Therefore,
αk < 0 is a sufficient condition for the geometric instability of a generic Bravais lattice;
while we have derived this result in the low-memory limit, one expects that increasing thememory will tend to promote lattice destabilization. Thus, a lattice that is geometricallyunstable in the low-memory limit should remain so at higher memories.
In Figure 7, we use the above criterion to determine the regions of initial stabilityof a rectangular lattice, as is described by two lattice spacings a and b. We note that the
diagonal line a = b in Figure 7 yields the intervals of stability shown in Figure 3b. As was
the case for the square and triangular geometries, the pockets of stability are correlatedwith minima in the local wave amplitude. Having identified the regions of initial stability,
the dispersion relation (18) could then be used to characterize the instability threshold andmost unstable vibrational mode for each lattice.

Symmetry 2022, 1, 0

14 of 20

0.11
0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03

Figure 7. The rectangular lattice is parameterised by two lattice spacings, a and b. Regions ofinitial stability at low memory are shaded in green. All other lattices are geometrically unstable.
The colourmap corresponds to the local wave amplitude beneath each drop h0 , normalized by thedrop radius R. Stable regions are roughly correlated with minima in h0 . The diagonal line a = b
corresponds to the square lattice, where the regions of initial stability correspond to those presented
in Figure 3.

4. Discussion
We have used the variable-phase stroboscopic model of Couchman et al. [32] to consider the linear stability of a group of droplets arranged in a Bravais lattice. All Bravais
lattices are stationary solutions to the stroboscopic model, since the lattice’s discrete translational symmetry ensures that the net waveforce on each droplet vanishes. By considering
the response of the lattice to normal mode perturbations, our analysis yielded an implicit
dispersion relation relating the vibrational mode’s wavevector k to its complex frequency
ωk , as captures both the mode’s growth rate and oscillation frequency. Particular attention
was given to the square and triangular Bravais geometries, for which the rotational symmetry of the lattice allowed our general dispersion relation to be reduced to a form similar tothat considered by Thomson et al. [44] in their investigation of a one-dimensional droplet
chain. A distinctive feature of our analysis is the inclusion of variations in the drop’s
vertical motion, as are required to capture the stabilizing influence of droplet–droplet interactions that allow lattices to remain stable above the walking threshold of the individual
constituent drops. This stabilizing phenomenon has also been reported for a variety ofother bound states [32–34,37,39].

Symmetry 2022, 1, 0

15 of 20

While all lattices are stationary solutions to the trajectory equation, not all are stable
at low vibrational accelerations; indeed, the majority are geometrically unstable [44]. For
the square and triangular lattices, we deduced a criterion for geometric instability solely
in terms of the local wave curvature beneath each droplet. This criterion predicts discreteintervals of stability in the lattice spacing. These stable regions correspond roughly togeometries that minimize the local wave amplitude beneath each droplet, a feature also
reported in bound droplet pairs [32], rings [33], and one-dimensional lattices [44]. By
considering the low-memory limit of our dispersion relation, we deduced a similar criterion
for assessing the geometric instability of a generic Bravais lattice, which we used to identify
the initial regions of stability for the rectangular lattice.
Consistent with the experimental observations reported by Eddi et al. [43], increasing the memory causes an initially stable lattice to destabilize into phonon-like motions
characterized by coherent small-amplitude oscillations. The most unstable modes werelongitudinal waves in all cases and were usually aligned along high-symmetry directionsof the lattice; shear modes only destabilize at higher vibrational forcings. We numericallycomputed the instability threshold of the most unstable modes, and found that they typically arise above the walking threshold γW for a single droplet, as highlights the stability
imparted by neighbouring droplets. Furthermore, we found local maxima for the instability
threshold at the centre of the initially stable regions, corresponding to minima in the local
wave amplitude, a result not captured in the theoretical modelling of the one-dimensionallattice [44] where variable phase factors were neglected.
The hydrodynamic lattice exhibits certain features that are distinct from canonicalmodels of phonons in crystalline lattices, such as those analysed by Blackman [53] andMontroll [54]. In these models, a discrete lattice of masses is connected to their nearest andnext-nearest neighbours via a linear spring force. If we were to remove the damping termfrom our dispersion relation (18), and take the low-memory limit considered in Section 3.2,
we would recover a dispersion relation of the same form as that of Blackman and Montroll,for which all modes are neutrally stable. For droplet lattices accessible in the laboratory, thedamping term in the trajectory Equation (1) dominates the inertia term in the low-memory
limit, resulting in overdamped oscillations. In our system, neutrally stable oscillations only
arise when the stabilizing effects of damping precisely balance the destabilizing effects ofmemory.
In the low-memory limit, the potential landscape is a sum of wave kernels centredat each droplet of the lattice. This potential induces a linear spring force with the local
curvature playing the role of the spring constant. Since the wave kernel is oscillatory, anindividual droplet may contribute a negative curvature (analogously, a negative spring
constant) to this sum, depending on the inter-droplet distance. This effect is not seen in a
canonical mass-spring lattice, as the springs are only perturbed about their equilibrium
lengths. The oscillatory wave kernel of the hydrodynamic lattice thus gives rise to muchricher dynamics than arise in generic crystal lattices, including the discrete windows ofinitial stability.
Our linear theory for the hydrodynamic lattice cannot predict the amplitude of theemergent oscillations. A weakly nonlinear extension of our theory, which parallels thatdeveloped for a one-dimensional hydrodynamic lattice [46], could yield insight into themore complex oscillations reported by Eddi et al. [43]. The potential for further theoretical
explorations abound. For example, one might consider how the introduction of defects,
such as holes or an additional droplet, might modify the resulting lattice dynamics [55].
In solid-state physics, more complex lattice geometries arise, for example, in ioniclattices comprised of more than one type of atom. In such cases, the ‘unit cell’ constitutes
the smallest non-repeating subcomponent of the lattice, and the structure may be defined interms of a Bravais lattice of such unit cells. As is well established in the phonon literature,having more than one element per unit cell is a prerequisite for optical modes and band gaps
[4,56], frequency ranges in which no vibrational states arise. In our study, we restricted ourattention to the case where there is a single droplet per unit cell, and focused on relatively

Symmetry 2022, 1, 0

16 of 20

simple geometries. In the future, one might extend this framework to describe lattices tiled
by unit cells containing multiple droplets, thereby describing arbitrary lattices in the plane.
We note that for such arrangements, the equilibrium condition ∇h(dmn ) = 0 (Equation (9))
would not necessarily be satisfied for all arbitrary unit cells. Nevertheless, an extended
theoretical framework might enable an investigation of the Archimedian tilings that Eddiet al. [42] were unable to access in the laboratory, such as the truncated hexagonal and greatrhombi-trihexagonal lattices.
Author Contributions: Conceptualization, M.M.P.C. and J.W.M.B.; methodology, software validation,
analysis, data curation, and visualization, M.M.P.C. and D.J.E.; writing—original draft preparation,
M.M.P.C. and D.J.E.; writing—review and editing, M.M.P.C., D.J.E., and J.W.M.B.; resources, supervision, administration, and funding acquisition, J.W.M.B. All authors have read and agreed to the
published version of the manuscript.
Funding: J.W.M.B. gratefully acknowledges the financial support of the NSF through grant CMMI2154151.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.
Conflicts of Interest: The authors declare no conflict of interest.

Appendix A. Impact-Phase Parameters
The strobing of the trajectory equation of Moláček and Bush [47] at the vertical bouncing period TF introduces two impact-phase parameters,

R t+TFS=t

FN (t0 ) sin(π f t0 )dt0

R t+TFt

FN (t0 )dt0

R t+TF, C=t

FN (t0 ) cos(π f t0 )dt0,
R t+TF0 ) dt0
F
t
(
N
t

(A1)

that represent averages of the sine and cosine of the bath’s phase of oscillation over theduration of droplet impact, weighted by the vertical contact force FN exerted on the dropby the bath [47]. The impact phase directly scales both the wave amplitude generated ateach impact (see Equation (2)) and the horizontal wave force imparted to the drop by thebath (see Equation (1)), as captured by S and C , respectively. While the phase parameters
are often combined into a constant fitting parameter sinΦ, modulations in a drop’s impact
phase have been shown to be critically important in accurately capturing the stability ofbound droplet states [32,33,37,39]. Couchman et al. [32] thus derived functional forms
for the dependencies of S and C on the bath’s vibrational acceleration γ, the local wave
amplitude beneath each drop h0 = h( x0 , t), and the drop radius R. To ensure that ourtheoretical predictions may provide the best possible comparison with future experimental
studies, we here include these variable-phase parameters in our analysis, noting that the
product SC may be set to a constant for a simpler description of the system.
The functional dependencies of S and C on the bath’s vibrational acceleration γ, localwave amplitude h0 , and droplet radius R, are presented in Table 2 of Couchman et al. [32].
To obtain the explicit results presented in Section 3, we here focus our theoretical analysis
on drops of radius R = 0.36 mm bouncing in a (2, 1)2 mode, as are typical parameters forexperimental studies [32,33,37,39]. In this case, the impact phase functions take the formS(γ, h0 ) = 1 − 1.32 exp{−3.52(γ/g − 5.73h0 /R − 2)},C(γ, h0 ) = 1.98 exp{−2.37(γ/g − 5.86h0 /R − 2)},

(A2)
(A3)

where the gravitational acceleration g and drop radius R are used to non-dimensionalize γ
and h0 , respectively.

Symmetry 2022, 1, 0

17 of 20

Appendix B. Rotational Symmetries of the Square and Triangular Lattices
We here derive the simplification H(0) = η0 I (Equation (24)) which holds for thesquare and triangular lattices, due to their respective four- and six-fold symmetries (seeFigure A1). We adopt the notationc pq ≡ x pq / d pq , s pq ≡ y pq / d pq ,(A4)

where d pq = x pq , y pq , and c pq and s pq thus represent the cosine and sine of the angle
between the position vector d pq and the x-axis. Equation (19), using k = 0, thus yields

H

(0)



"
#
"
#

 00
f 0 (|d pq |)c2pq
c pq s pq
s2pq−c pq s pq .
= 2 ∑ f (|d pq |)+c pq s pq
s2pq−c pq s pq
c2pq
|d pq |
pq 

|
{z
}
|
{z
}
A1

(A5)

A2

For both the square and triangular geometries, we now demonstrate that the right-handside of Equation (A5) is proportional to the identity matrix I.

Figure A1. Square (a) and triangular (b) lattices exhibit four- and six-fold symmetry, respectively,
and may be generated by rotating the bolded points around the origin in increments of π/2 and π/3.

Square Lattice
Consider performing the sums in Equation (A5) by grouping four points at a time,
selected by rotating an initial point ( p, q) in the first quadrant of Figure A1a through anglesof θ = lπ/2, l ∈ [0, 3]. Importantly, each of these four points maintains the same distance|d pq | from the origin. Using the identities
3

3

l =0
3

l =0

∑ sin2 (θ0 + lπ/2) = ∑ cos2 (θ0 + lπ/2) = 2,

(A6)

∑ sin(θ0 + lπ/2) cos(θ0 + lπ/2) = 0,

(A7)

l =0

yieldsSymmetry 2022, 1, 0

18 of 20

3

3  2
cl
A1 =
c
l sl
l =0
l =0

∑

∑

3

3 

l =0

l =0

∑ A2 = ∑


0
=
,
2
 

−cl sl
2 0
=
.
0 2
c2l

cl sl
s2l

s2l
−cl sl





2
0

(A8)
(A9)

Thus, the right-hand side of Equation (A5) is proportional to the identity matrix I.
Triangular Lattice
Performing an analogous procedure as for the square lattice, we now evaluate thesums in Equation (A5) by grouping six points at a time, selected by rotating an initial point( p, q) in the bolded sixth of Figure A1b through angles of θ = lπ/3, l ∈ [0, 5]. Using theidentities
5

5

l =0
5

l =0

∑ sin2 (θ0 + lπ/3) = ∑ cos2 (θ0 + lπ/3) = 3,

(A10)

∑ sin(θ0 + lπ/3) cos(θ0 + lπ/3) = 0,

(A11)

l =0

yields
5

5  2
cl
cl sl
l =0

∑ A1 = ∑

l =0
5

5 

∑ A2 = ∑

l =0

l =0

s2l
−cl sl


0
,
3
 

−cl sl
3 0
=
.
c2l
0 3

cl sl
s2l





=

3
0

(A12)
(A13)

Thus, once again, the right-hand side of Equation (A5) is proportional to the identity
matrix I.
References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.

Krumhansl, J. Lattice vibrations in solids. J. Appl. Phys. 1962, 33, 307–319.
Cochran, W. Lattice vibrations. Rep. Prog. Phys. 1963, 26, 1.
Dean, P. Atomic vibrations in solids. IMA J. Appl. Math. 1967, 3, 98–165.
Kittel, C. Introduction to Solid State Physics; Wiley: Hoboken, NJ, USA, 1996.
Kane, C.; Lubensky, T. Topological boundary modes in isostatic lattices. Nat. Phys. 2014, 10, 39–45.
Lubensky, T.; Kane, C.; Mao, X.; Souslov, A.; Sun, K. Phonons and elasticity in critically coordinated lattices. Rep. Prog. Phys.
2015, 78, 073901.
Nash, L.M.; Kleckner, D.; Read, A.; Vitelli, V.; Turner, A.M.; Irvine, W.T. Topological mechanics of gyroscopic metamaterials. Proc.
Natl. Acad. Sci. USA 2015, 112, 14495–14500.
Süsstrunk, R.; Huber, S.D. Classification of topological phonons in linear mechanical metamaterials. Proc. Natl. Acad. Sci. USA
2016, 113, E4767–E4775.
Bragg, W.L.; Nye, J.F. A dynamical model of a crystal structure. Proc. R. Soc. A 1947, 190, 474–481.
Penciu, R.; Kriegs, H.; Petekidis, G.; Fytas, G.; Economou, E. Phonons in colloidal systems. J. Chem. Phys. 2003, 118, 5224–5240.
Keim, P.; Maret, G.; Herz, U.; von Grünberg, H.H. Harmonic lattice behavior of two-dimensional colloidal crystals. Phys. Rev.
Lett. 2004, 92, 215504.
Kaya, D.; Green, N.; Maloney, C.; Islam, M. Normal modes and density of states of disordered colloidal solids. Science 2010,
329, 656–658.
Chen, K.; Still, T.; Schoenholz, S.; Aptowicz, K.B.; Schindler, M.; Maggs, A.; Liu, A.J.; Yodh, A. Phonons in two-dimensional soft
colloidal crystals. Phys. Rev. E 2013, 88, 022315.
Ghosh, A.; Mari, R.; Chikkadi, V.; Schall, P.; Maggs, A.; Bonn, D. Low-energy modes and Debye behavior in a colloidal crystal.
Phys. A Stat. Mech. Its Appl. 2011, 390, 3061–3068.
Beatus, T.; Tlusty, T.; Bar-Ziv, R. Phonons in a one-dimensional microfluidic crystal. Nat. Phys. 2006, 2, 743–748.
Schiller, U.D.; Fleury, J.B.; Seemann, R.; Gompper, G. Collective waves in dense and confined microfluidic droplet arrays. Soft
Matter 2015, 11, 5850–5861.

Symmetry 2022, 1, 0

17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.

19 of 20

Tsang, A.C.H.; Shelley, M.J.; Kanso, E. Activity-induced instability of phonons in 1D microfluidic crystals. Soft Matter 2018,14, 945–950.
Ono-dit-Biot, J.C.; Soulard, P.; Barkley, S.; Weeks, E.R.; Salez, T.; Raphaël, E.; Dalnoki-Veress, K. Mechanical properties of 2D
aggregates of oil droplets as model mono-crystals. Soft Matter 2021, 17, 1194–1201.
Benjamin, T.B.; Ursell, F. The stability of the plane free surface of a liquid in vertical periodic motion. Proc. R. Soc. A 1954,
225, 505–515.
Miles, J.; Henderson, D. Parametrically forced surface waves. Annu. Rev. Fluid Mech. 1990, 22, 143–165.
Walker, J. The amateur scientist. Drops of liquid can be made to float on liquid. What enables them to do so? Sci. Am. 1978,
238, 151.
Couder, Y.; Fort, E.; Gautier, C.H.; Boudaoud, A. From bouncing to floating: noncoalescence of drops on a fluid bath. Phys. Rev.
Lett. 2005, 94, 177801.
Moláček, J.; Bush, J.W.M. Drops bouncing on a vibrating bath. J. Fluid Mech. 2013, 727, 582–611.
Couder, Y.; Protiere, S.; Fort, E.; Boudaoud, A. Dynamical phenomena: walking and orbiting droplets. Nature 2005, 437, 208.
Protière, S.; Boudaoud, A.; Couder, Y. Particle–wave association on a fluid interface. J. Fluid Mech. 2006, 554, 85–108.
Eddi, A.; Sultan, E.; Moukhtar, J.; Fort, E.; Rossi, M.; Couder, Y. Information stored in Faraday waves: the origin of a path memory.
J. Fluid Mech. 2011, 674, 433–463.
Fort, E.; Eddi, A.; Boudaoud, A.; Moukhtar, J.; Couder, Y. Path-memory induced quantization of classical orbits. Proc. Natl. Acad.
Sci. USA 2010, 107, 17515–17520.
Bush, J.W.M. Pilot-wave hydrodynamics. Annu. Rev. Fluid Mech. 2015, 47, 269–292.
Bush, J.W.M.; Couder, Y.; Gilet, T.; Milewski, P.A.; Nachbin, A. Introduction to focus issue on hydrodynamic quantum analogs.
Chaos Interdiscip. J. Nonlinear Sci. 2018, 28, 096001.
Bush, J.W.M.; Oza, A.U. Hydrodynamic quantum analogs. Rep. Prog. Phys. 2020, 84, 017001.
Rahman, A.; Blackmore, D. Walking droplets through the lens of dynamical systems. Mod. Phys. Lett. B 2020, 34, 2030009.
Couchman, M.M.P.; Turton, S.E.; Bush, J.W.M. Bouncing phase variations in pilot-wave hydrodynamics and the stability ofdroplet pairs. J. Fluid Mech. 2019, 871, 212–243.
Couchman, M.M.P.; Bush, J.W.M. Free rings of bouncing droplets: stability and dynamics. J. Fluid Mech. 2020, 903.
Thomson, S.J.; Couchman, M.M.P.; Bush, J.W.M. Collective vibrations of confined levitating droplets. Phys. Rev. Fluids 2020,
5, 083601.
Barnes, L.; Pucci, G.; Oza, A.U. Resonant interactions in bouncing droplet chains. C. R. Mécanique 2020, 348, 573–589.
Protière, S.; Bohn, S.; Couder, Y. Exotic orbits of two interacting wave sources. Phys. Rev. E 2008, 78, 036204.
Oza, A.U.; Siéfert, E.; Harris, D.M.; Moláček, J.; Bush, J.W.M. Orbiting pairs of walking droplets: dynamics and stability. Phys.
Rev. Fluids 2017, 2, 053601.
Borghesi, C.; Moukhtar, J.; Labousse, M.; Eddi, A.; Fort, E.; Couder, Y. Interaction of two walkers: wave-mediated energy andforce. Phys. Rev. E 2014, 90, 063017.
Arbelaiz, J.; Oza, A.U.; Bush, J.W.M. Promenading pairs of walking droplets: dynamics and stability. Phys. Rev. Fluids 2018,
3, 013604.
Lieber, S.I.; Hendershott, M.C.; Pattanaporkratana, A.; Maclennan, J.E. Self-organization of bouncing oil drops: two-dimensionallattices and spinning clusters. Phys. Rev. E 2007, 75, 056308.
Eddi, A.; Terwagne, D.; Fort, E.; Couder, Y. Wave propelled ratchets and drifting rafts. Europhys. Lett. 2008, 82, 44001.
Eddi, A.; Decelle, A.; Fort, E.; Couder, Y. Archimedean lattices in the bound states of wave interacting particles. Europhys. Lett.
2009, 87, 56002.
Eddi, A.; Boudaoud, A.; Couder, Y. Oscillating instability in bouncing droplet crystals. Europhys. Lett. 2011, 94, 20004.
Thomson, S.J.; Durey, M.; Rosales, R.R. Collective vibrations of a hydrodynamic active lattice. Proc. R. Soc. A 2020, 476, 20200155.
Oza, A.U.; Rosales, R.R.; Bush, J.W.M. A trajectory equation for walking droplets: hydrodynamic pilot-wave theory. J. Fluid Mech.
2013, 737, 552–570.
Thomson, S.J.; Durey, M.; Rosales, R.R. Discrete and periodic complex Ginzburg-Landau equation for a hydrodynamic active
lattice. Phys. Rev. E 2021, 103, 062215.
Moláček, J.; Bush, J.W.M. Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory. J. Fluid Mech. 2013,
727, 612–647.
Turton, S.E.; Couchman, M.M.P.; Bush, J.W.M. A review of the theoretical modeling of walking droplets: towards a generalized
pilot-wave framework. Chaos Interdiscip. J. Nonlinear Sci. 2018, 28, 096111.
Harris, D.M.; Quintela, J.; Prost, V.; Brun, P.T.; Bush, J.W.M. Visualization of hydrodynamic pilot-wave phenomena. J. Vis. 2017,20, 13–15.
Damiano, A.P.; Brun, P.T.; Harris, D.M.; Galeano-Rios, C.A.; Bush, J.W.M. Surface topography measurements of the bouncing
droplet experiment. Exp. Fluids 2016, 57, 163.
Tadrist, L.; Shim, J.B.; Gilet, T.; Schlagheck, P. Faraday instability and subthreshold Faraday waves: surface waves emitted bywalkers. J. Fluid Mech. 2018, 848, 906–945.
Wind-Willassen, Ø.; Moláček, J.; Harris, D.M.; Bush, J.W.M. Exotic states of bouncing and walking droplets. Phys. Fluids 2013,
25, 082002.

Symmetry 2022, 1, 0

53.
54.
55.
56.

20 of 20

Blackman, M. Contributions to the theory of the specific heat of crystals II—On the vibrational spectrum of cubical lattices and its
application to the specific heat of crystals. Proc. R. Soc. A 1935, 148, 384–406.
Montroll, E.W. Dynamics of a square lattice I. Frequency spectrum. J. Chem. Phys. 1947, 15, 575–591.
Bell, R. The dynamics of disordered lattices. Rep. Prog. Phys. 1972, 35, 1315–1409.
Kosevich, A.M. The Crystal Lattice: Phonons, Solitons, Dislocations, Superlattices; John Wiley & Sons: Hoboken, NJ, USA, 2006.