# Statistical self-organization of an assembly of interacting walking drops in a confining potential

Adrien Hélias, Matthieu Labousse

► To cite this version:

Adrien Hélias, Matthieu Labousse. Statistical self-organization of an assembly of interacting walking drops in a confining potential. European Physical Journal E: Soft matter and biological physics, 2023, 46 (4), pp.29. 10.1140/epje/s10189-023-00288-5. hal-04070960

HAL Id: hal-04070960

https://hal.science/hal-04070960

Submitted on 20 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistical self-organization of an assembly of interacting walking drops in a confining potential

Adrien Hélias

Gulliver, UMR CNRS 7083, ESPCI Paris, Université PSL, 75005 Paris, France and

Current address: Department of Physics and Astronomy,
Western University, 1151 Richmond St,

N6A 3K7 London, Canada

ahelias@uwo.ca

Matthieu Labousse

Gulliver, UMR CNRS 7083, ESPCI Paris,

Université PSL, 75005 Paris, France

matthieu.labousse@espci.fr

A drop bouncing on a vertically-vibrated surface may self-propel forward by standing waves and travels along a fluid interface. This system called *walking drop* forms a non-quantum wave-particle association at the macroscopic scale. The dynamics of one particle has triggered many investigations and has resulted in spectacular experimental results in the last decade. We investigate numerically the dynamics of an assembly of walkers, i.e. a large number of walking drops evolving on a unbounded fluid interface in the presence of a confining potential acting on the particles. We show that even if the individual trajectories are erratic, the system presents well-defined ordered internal structure that remains invariant to parameter variations such as the number of drops, the memory time and the bath radius. We rationalize such non-stationary self-organization in terms of the symmetry of the waves and show that oscillatory pair potentials form a wavy collective state of active matter.

## I. INTRODUCTION

In a seminal article, Y. Couder and E. Fort [1] have shown that a submillimetric drop may bounce on a vertically-vibrated fluid interface. Above a certain fluid acceleration threshold the drop may self-propel and be guided by standing waves [2–13] which are the footprint of Faraday waves [14–16]. The system has triggered a flurry of thought-provoking experiments mimicking effects which were previously thought to be peculiar to the quantum scale [2, 3]. Specifically, the dynamics of one drop has been investigated in many configurations, like moving through a slit and double slits [17–19], in cavities [20–24], with a Coriolis force [25–27], in a harmonic potential [13, 26, 28–30], exhibiting Friedel-like oscillations [31], enabling statistical projections [23]. In the last decade, an important effort of research has led to a better understanding of both the emergence of wave-like statistics and the mechanisms responsible for classical quantization on the macroscopic scale [2, 3].

In a parallel stream of research, it has been shown, numerically and experimentally, that complex dynamical behaviors such as memory-induced diffusion, and run and tumbling dynamics [32–40] appear for a single droplet system, in which the degrees of freedom in the wave field play the role of a tailored thermal bath [40]. Generalized pilot-wave models [41] have shown that these behaviours are observed in parameter regimes beyond those accessible in the laboratory [36].

The dynamics of two drops has also been investigated and promenade modes [5, 13, 42–45], as well as

quantized orbiting states have been observed and rationalized [1, 13, 45–47]. In addition, correlated motions of two walkers trapped in a coupled set of cavities has led to analogs of superradiance [48, 49] and the violation of static Bell's inequality in a static test [50]. As for the many drops' dynamics, crystalline structures have been obtained in 2D [51–53], and the collective dynamics in toroidal channel [54, 55] as well as the spatially periodic potentials [56] have been investigated in the linear and non-linear regime [57–61]. Recently, drops trapped in circular cavities distributed over a two-dimensional lattice have shown the spontaneous emergence of long distance motion synchronisation in analogy with ferromagnetism and antiferromagnetism [62]. In all these investigations [51–55, 57–62], and because of the complexity of preventing the coalescence during drop collisions, the drops are carefully initially placed in either particular positions or specific relative distances to each other. As a result of these constraints, and although it would be a remarkable playground for active matter and statistical physics, a *active gas phase of walkers*, i.e. many interacting walking drops confined in a single large domain (see Fig. 1) is experimentally unexplored. A numerical investigation is a necessary first step to trigger and motivate further experimental realisations. In this article, we numerically investigate the dynamics of a large number of walking drops evolving on a unbounded fluid interface in the presence of a confining wall potential acting on the particles.

The article is organised as follows. In Sec. II, we introduce the numerical model. Then, in Sec. III, we characterize and rationalize the properties of the dynam-
ical phases and propose a mechanism to predict the main properties of this statistical quantized self-organization. Finally, we conclude in Sec. IV.

## II. NUMERICAL MODEL AND METHODS

We perform a numerical simulation of many walking drops in two dimensions. Several models [11, 12, 25, 41, 63], varying in their strategies to solve the fluid equation, the impact conditions, and the horizontal dynamics of the drops, have been proposed. They are in fairly good qualitative agreement as they share the same main essential features (see the review [3] for the comparison and relative merits of all these models). Here, we aim to simulate an assembly of walkers which have long-range interactions. As the connectivity matrix of particles is fully coupled, the computational cost becomes a critical aspect. As a result, we implement the walking drops dynamics by considering a numerically optimized [32, 40] version of a discrete time evolution [25], which can be summarized as follows.

First, we model the wave created by a single impact on an unbounded fluid interface by a Bessel function of order 0, $J_0$, centered at the point of impact. The standing wave frequency $f_F$ is directly linked to the vibration frequency of the bath $f_0$ by the relation $f_F = f_0/2$ [8, 9, 14–16, 64–68] corresponding to a quasi-monochromatic wave at a wavelength $\lambda_F = 6.1$ mm. As the drops impact the bath at constant relative phase every period $T_F = 1/f_F = 40$ ms, the wave field is updated periodically. Thus, a single impact at the surface position, $r_0$, gives rise to a surface standing wave field with a spatial structure proportional to $J_0(k_F \|r - r_0\|)$, where $r$ is a two-dimensional vector, and $k_F = 2\pi/\lambda_F$ the Faraday wave vector. Experimentally, taking into account the fluid bath viscosity is known to lead to a spatial decay more pronounced than a Bessel function [12, 41, 69]. However considering a simplified Bessel wave field is a widely used approximation which has provided most of the theoretical predictions in the field. As a consequence, we keep this simplified Bessel approximation for the wave field expression [3, 70].

Secondly, the persistence time of these waves is defined by the memory time, $Me \times T_F$, with Me the memory parameter. The larger the value, the longer the waves remain on the surface and influence the dynamics of the drops. Experimentally, the memory time is tuned by varying the acceleration amplitude of the bath relative to the critical Faraday acceleration threshold. In the numerical model, the memory is a scalar parameter Me which we can vary at will. Me ~ 10 is considered as a low memory regime, and Me ~ 50 as a moderately small memory regime [40]. We only investigate Me values ranging from 10 to 50. In practice, this memory range

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematics of the system. (a) Example of the system with $N = 3$ drops, $Me = 30$, $R_{max} = 3\lambda_F$, seen from a top view. The coloured surface illustrates the wave field (in arbitrary units) created by the bouncing drops dotted in grey. The waves are not restrained by the outer wall, only the drops are. (b) A sketch of the confining force per unit mass: it is equivalent to a step of force of 1 N/kg acting beyond a distance to the center larger than $R_{max}$.</figcaption></figure>

corresponds to an acceleration amplitude of the fluid bath $\gamma_m/\gamma_F = (1 - 1/Me) \in [0.9, 0.98]$, with $\gamma_F$ the fluid-dependent Faraday acceleration threshold.

Thirdly, we consider the system as sketched in Fig. 1. N walking drops are confined by a very steep external potential acting only beyond a radius $R_{max}$. The wave field is considered as not influenced by the presence of this confining potential and as a consequence, the waves are generated as if they were feeling an unbounded fluid interface. We determine the evolution of the position of the i-th drop ($i \in \{1, ..., N\}$) from time $t_p$ and at position $r_i(t_p)$ to time $t_{p+1} = t_p + T_F$, i.e. one bounce further, from the relation $r_i(t_{p+1}) = r_i(t_p) + v(t_p)T_F$ and by computing in parallel the evolution of the velocity of the i-th drop, $v_i$, as

$$v_i(t_{p+1}) = v_i(t_p) + T_F \left( -\eta v_i(t_p) + F_{\text{wave},i}(t_p) + F_{\text{wall},i} + F_{\text{drops}\to\text{drop } i} \right). \quad (1)$$

wherein $\eta$ is the dissipation coefficient is set to 4.72 s⁻¹ [26, 33], which corresponds to the situation of a silicone oil with a kinematic viscosity of 20 cSt. It takes into account both the friction due to the lubrication layer between the drop and the bath, and the transfer of momentum to the fluid bath. $F_{wave,i}$ is the force per unit mass propelling the i-th drop and is proportional to the local value of the two-dimensional gradient of the wave
169 field. Specifically, it writes

$$ \mathbf{F}_{\text{wave},i}(t_p) = -\frac{CT_F}{k_F} \nabla \mathbf{r}_i(t_p) \sum_{j=1}^{N} \left( \sum_{k=p-3\text{Me}}^{p} J_0(k_F \|\mathbf{r}_i(t_p) - \mathbf{r}_j(t_k)\|) e^{-\frac{(p-k)}{\text{Me}}} \right) \quad (2) $$

170

The first sum accounts for the linear superposition of the wave generated by the drops $j$ ($j \in \{1...N\}$). The second sum originates from the memory kernel which indicates that the amplitude of the waves generated by the $j$-th drop at time $t_k < t_p$ is exponentially damped. $C = 1.1$ m.s$^{-3}$ is a coupling coefficient [26, 33] expected from a drop of diameter 0.9 mm. Rigorously, the sum over $k$ should start from $k = 0$ but thanks to the exponential term, we start the summation from $p - 3\text{Me}$, neglecting the other residual terms (5%). This truncation enables us to reduce significantly the computational costs, without observable difference in the statistical convergence of the pair correlation function introduced later. $F_{\text{wall}}$ is the force per unit mass accounting for the presence of the outer wall. It is phenomenologically modelled by a step of force at a radius $R_{\text{max}}$. It is only acting for a radius larger than $R_{\text{max}}$ as sketched in Fig. 1b. It writes

188

$$ \mathbf{F}_{\text{wall}} = -F_0 e_r H (r_i(t_p) - R_{\text{max}}) \quad (3) $$

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

<figure><img src="image_3.png" /><figcaption>FIG. 2. Screenshots of 30 walking drops evolving after $10^3$ time steps for a memory parameter $Me = 30$, with different bath sizes (a) $3\lambda_F$ and (b) $5\lambda_F$. The grey circles represent the drops with diameter 0.9 mm. The blue lines link a drop to its $Me$ past positions to illustrate the magnitude of the memory parameter. Two movies showing the drops moving during a few tens of periods are available in Supplemental Materials.</figcaption></figure>

5 drops, $10^5$ Faraday periods are considered to decrease the measurements uncertainties. We discard the first $10^3$ Faraday periods to avoid any transient effects.

### III. RESULTS : DYNAMICAL SELF-ORGANIZATION

We investigate the dynamics by changing the number of drops $N$, the domain radius $R_{\text{max}}$, and the memory parameter $Me$. Figures 2a and 2b illustrate two characteristic snapshots of the positions of 30 drops confined after $10^3$ time steps in a domain of radius $R_{\text{max}} = 3\lambda_F$ and $5\lambda_F$, respectively (see also Supplemental Movie
1 and 2). The magnitude of the memory parameter is indicated by a blue tail which visually links a drop to its Me past positions. The mean density of drops, $\rho = N/(\pi R_{\text{max}}^2)$ is low in both Fig. 2a and b, specifically $\rho_{\text{drop}} \simeq 1.1 \lambda_F^{-2}$ and $0.4 \lambda_F^{-2}$, respectively, and we observe erratic trajectories, regardless of the bath size. We remark that sometimes, the drops walk together by forming small groups during a few tens of periods but we did not observe any reproducible patterns nor any obviously-identifiable long-standing stable emerging structure. Very small excursions outside $R_{\text{max}}$ can be observed in Figures 2 and are dictated by the precise details of the modelling of the wall repulsion. We also observe that the drops are overlaying and crossing each other sometimes. In this low drop density regime, we check that the dynamics is statistically not influenced by this effect by comparing with the case when the force $F_{\text{drops}} \to \text{drop i}$ is switched on, an investigation which we shall discuss further.

As we observe no apparent structure by looking solely at the dynamics, we proceed by investigating the statistical properties on the N drops system for a simulation time $10^4 T_F$. We first compute the probability density function (PDF), $F(r)$. We define the radial probability function, $f(r)$, by integrating $F$ over all possible angular components, specifically, $f(r) = \int_0^{2\pi} F(r, \theta) d\theta$. By construction, we get the normalization $\iint f(r)rdr = \iint F(r, \theta)rdrd\theta = 1$. As a consequence, Figure 3a shows the evolution of $rf(r)$, with the length expressed in $\lambda_F$ units, for $N = 5$ and 30, a confining radius $R_{\text{max}} = 3\lambda_F$ and $5\lambda_F$, and for a common memory parameter Me = 30. For all parameter configurations, we observe a linear dependency of $rf(r)$ with $r$ up to a radius $R^* \simeq R_{\text{max}} - \lambda_F$, where $rf(r)$ is maximum. For $r < R^*$, the slopes are independent of the drop density. Then, in this low density regime, we find that the slopes scale as $1/R_{\text{max}}^2$. For distances from the center $r > R^*$, the drops feel the influence of the boundary and $rf(r)$ sharply decreases in the vicinity of $R_{\text{max}}$. For radial distance to the center $r < R^*$, a linear behavior of $rf(r)$ means that $f$ is constant with $r$. So, in this region, the distribution of drops is radially homogenous. We note that this radial distribution is different from the case of a walker in a corral where both the walker motion and the wave modes are influenced by the boundary of the domain [20]. Additionally, we did not find any correlation between position and speed which is another difference with the corral system [20].

Then, we analyse the angular structure of the system by computing the angular probability distribution $h(\theta)$, specifically $h(\theta) = \int F(r, \theta)rdr$. Figure 3b shows that the drops angular distribution remains isotropic for all drop densities and all values of $R_{\text{max}}$ investigated.

Figure 3c shows the 2D probability density function established for a simulation time $10^4 T_F$, 30

<figure><img src="image_4.png" /><figcaption>FIG. 3. Statistical properties of the dynamics (a) Radial, (b) angular probability density function (PDF). (c) two-dimensional PDF for simulation time $10^4 T_F$ and the relative absolute value of the weight modes $|a_n|/|a_0|$ of the mean wave field $\bar{h}$. (d) speed PDF. Common color code: green $N = 30$ and red $N = 5$. The squares represent the situation $R_{\text{max}} = 3\lambda_F$ and the triangles $R_{\text{max}} = 5\lambda_F$. Me = 30 for all configurations.</figcaption></figure>
droplets and Me = 30. It does not have any particular structure and is homogenous up to some fluctuations which we assume to be due to the spatial coarse-graining. Then, we compute the convolution between the PDF and the Bessel wave kernel which gives the mean wave field, $\bar{h}$, according to the mean wave field theorem of Durey et al. [34]. We decompose $\bar{h}$ into a wave polar basis ($J_n(2\pi r/\lambda_F)e^{in\theta})_{n\in\mathbb{Z}}$ and calculate the weight modes $a_n$. The mean wave field is mainly dominated by the axisymmetric mode ($n=0$) which takes the spatial form of the Bessel function $J_0$. This is expected as the distribution of drops is found to be homogenous. So, because of this angular invariance, the only mode which contributes is the axisymmetric one. It is also a difference with [40] where the mean wave field of a single particle, in asymptotically large memory limit, is not only a function of the density but requires the knowledge of higher-order temporal correlation functions. The comparison of the mean wave field with that obtained in [71] with walking drops in a central circular well with a Bessel potential, is interesting. In both cases, the spatial structure of the mean wave field is expected to be the same as the trajectories are statistically invariant by rotation. However, the prefactors must be very different as they encode the particular shape of the radial distributions which strongly differ between the two.

Finally, we compute the speed probability distribution in Fig. 3d, from which we measure a mean speed $\bar{v} = 1.40, 1.59, 1.38, 1.49$ cm/s, respectively, for $(N, R_{max}) = (5, 3\lambda_F), (30, 3\lambda_F), (5, 5\lambda_F), (30, 5\lambda_F)$, respectively. The mean speeds are in the range expected experimentally. We observe that increasing the radius tends to decrease the mean drops' speed. Finally, we note that the systems with larger numbers of drops have more pronounced high-speed excursions. We investigate in more details this point in Figure 4 by plotting the speed density for increasing number of drops, up to 300.

We rationalize the evolution of the speed PDF by taking inspirations from the active statistical theory developed in [40]. In the large memory regime (Me $\gg 10^2$) investigated by the authors, the motion of one single drop is equivalent to a self-propelled particle coupled to a white noise thermal bath. Following [40], we fit, in Fig. 4a, the numerical speed probability density functions with $h(v) = \alpha v \exp(-\beta((v-v_0)^2))$, wherein $v_0$ originates from a constrain on the speed arising from the self-propulsion of the drops and $\beta$ the effective temperature of the system. $\alpha$ is a normalisation factor such that $\int h(v)dv = 1$. Note that in the theoretical asymptotic limit $v_0 \to 0$, $h$ identifies to a two-dimensional Maxwell-Boltzmann distribution as one would expect from a classical ideal gas in two dimensions. We observe some systematic quantitative differences between the numerical results and the fitting model especially with the tail of the distribution, suggesting strong remaining

correlations.

To reveal these correlations, we measure in Fig. 4b and 4c, the evolution of the mean speed and of its standard deviation as a function of the number of drops along with their fits with 95% confidence intervals. We observe a small and steady increase of the mean speed with the number of drops. It is a difference with the case of the promenade mode in which two correlated drops move slower than the case of one single drop. We conclude that we do not reach a statistical limit by increasing the number of drops, an important difference with an ideal gas and with the behavior of one single drop in the large memory regime [40]. We note that the increase of the mean speed should obviously reach an upper limit in experiments as nonlinear wave effects are expected to start playing a saturation role.

The speed fluctuations also increase with the number of drops (see Fig. 4c) and we rationalize its evolution as follows. We expect that the squared speed fluctuation of a particle $i$, $\sigma_i^2$, should be the sum of two contributions: one intrinsic, denoted $\sigma_{i,0}^2$, and one due to the presence of all the other drops $j \neq i$, denoted $\Sigma_i^2$:

$$ \sigma_i^2 = \sigma_{i,0}^2 + \Sigma_i^2. \quad (4) $$

In principle, $\Sigma_i^2$ is a function of $i$ the number of drops, $N$, and $ii$ all possible degrees of freedom of all the other drops $j \neq i$. All drops being indistinguishable, we assume that $\Sigma_i^2$ can be factorized as $\Sigma_i^2 = \Sigma_{m,i}^2 F(N)$ where $F$ is a function independent on $i$ and only dependent on the number of drops, while $\Sigma_{m,i}^2$ is a function which depends, in all generality, on all possible degrees of freedom of all the other drops $j \neq i$. This factorization ansatz is equivalent to a meanfield approach, which we expect to hold for large $N$. Additionally, the number of drops is a scale-free parameter, so $F(N)$ is expected to be an algebraic function of $N$, specifically $F(N) = N^\nu$. $\nu$ is an exponent to be determined from the fitting onto the numerical results. By combining Eq. 4 with the factorization ansatz, and then averaging on both time and on all drops $i=1...N$, we infer that the speed standard deviation $\delta v$ takes the form

$$ \delta v^2 = \delta v_0^2 + \delta v_m^2 N^\nu \quad (5) $$

with $\delta v_0$ and $\delta v_m$ two constants to be determined from the numerical results. Fig. 4c indicates a good agreement with the proposed fit (Eq. 5).

To investigate these internal correlations in more details, we aim at spotting the signature of any internal structure. To do so, we compute the drops pair correlation function

$$ g(\mathbf{r}) = \left\langle \sum_{i,j, i \neq j} \delta(\mathbf{r} - (\mathbf{r}_i - \mathbf{r}_j)) \right\rangle \quad (6) $$
<figure><img src="image_5.png" /><figcaption>FIG. 4. Analysis of speed probability density functions. (a) Speed PDF for different numbers of drops, Me = 30 and R<sub>max</sub> = 5λ<sub>F</sub>. The dashed lines represent the fits by h(v) = αv exp(−β(v − v<sub>0</sub>)<sup>2</sup>) of the numerical speed probability density functions (plain transparent lines). (b) Evolution of the mean speed as a function of the number of drops. N. The fitting functions used is of the form aN<sup>b</sup> + c with a = 0.0161, b = 0.7201, c = 0.3190 and a determination coefficient R<sup>2</sup> = 0.9987. (c) Evolution of the speed standard deviation as a function of the number of drops N. The fitting functions used is of the form of Eq. (5) with ν = 0.681, δv<sub>0</sub><sup>2</sup> = 1.36 10<sup>−2</sup> cm<sup>2</sup>/s<sup>2</sup>, δv<sub>m</sub><sup>2</sup> = 0.262 cm<sup>2</sup>/s<sup>2</sup> with a determination coefficient, R<sup>2</sup> = 0.9996.</figcaption></figure>

where ⟨...⟩ denotes an average over time. As there is no preferred orientation, we assume g(r) = g(r). The corresponding results are presented sequentially in Fig. 5, each panel corresponding to the variations of one single parameter. Fig. 5a indicates that the pair correlation function decreases with the radial distance, up to twice R<sub>max</sub>, with the presence of well-pronounced equally-spaced peaks being associated with preferred inter-drops distances. The number of discernible peaks increases with the radius of the domain. Peaks apart, the global decreasing trend is expected as a situation, for which the drops would be homogeneously distributed.

in a circular domain of radius R<sub>max</sub>, exhibit the same decreasing behavior. Fig. 5b indicates that varying memory in the range Me ∈ [10, 50] does not change the pair correlation function profile, and moderately alter the magnitude of the peak-to-peak amplitude. Finally, Fig. 5c shows that varying the number of drops only affects the signal-to-noise ratio, which is very low for 5 drops. We observe on Figure 5 local maxima on the curves. It means that drops tend to locate statistically from one to each other at a preferred set of quantized distances. Some of them are well pronounced with a well-measurable position and some of them are emerging which make their position difficult to estimate accurately. Among the peaks which are clearly formed, we measure, for example at R<sub>max</sub> = 7λ<sub>F</sub>, the position of maxima r*/λ<sub>F</sub> = 1.09, 2.07, 3.08, 4.06, 5.04, 6.06, 6.97, 7.91 (± 0.04).

This quantized distance originates from the force due to the wave field (Eq. 2), which we rationalize now by considering the dominant terms in Eq. 2. We introduce the notations (see also Fig. 6): $\mathbf{r}_{ij}(t_p) = \mathbf{r}_i(t_p) - \mathbf{r}_j(t_p)$ is a vector pointing from the drop j to the drop i, and $\mathbf{\Delta}_j(t_p, t_k) = \mathbf{r}_j(t_k) - \mathbf{r}_j(t_p)$ linking the current position of the drop j to one of its position in the past. Additionally, we denote the associated distances $r_{i,j}(t_p) = \|\mathbf{r}_{ij}(t_p)\|$, $\mathbf{\Delta}_j(t_p, t_k) = \|\mathbf{\Delta}_j(t_p, t_k)\|$, $J_n$ the Bessel function of first kind of order n and $\theta_{i,j}(t_p, t_k)$ the relative angle between $\mathbf{r}_{ij}(t_p)$ and $\mathbf{\Delta}_j(t_p, t_k)$. Note that the symbol i in the complex exponential in Eq. 7 denotes the imaginary number and differs from the italicized subscript i. We identify the dominant terms in Eq. 2 by proceeding as follows. As varying the memory leaves the position of maxima unchanged, we consider the regime of large Me which means that we neglect the exponential decay hereafter in Eq. 2. Then, we note that the wave term 'J<sub>0</sub>' in Eq. 2 can be written

$$ 
\begin{align} 
J_0(k_F \|\mathbf{r}_i(t_p) - \mathbf{r}_j(t_k)\|) 
&= J_0 (k_F \|(\mathbf{r}_i(t_p) - \mathbf{r}_j(t_p)) - (\mathbf{r}_j(t_k) - \mathbf{r}_j(t_p))\|) \\ 
&= J_0 (k_F \|\mathbf{r}_{ij}(t_p) - \mathbf{\Delta}_j(t_p, t_k)\|) \\ 
&= \sum_{n=-\infty}^{+\infty} J_n(k_F \mathbf{r}_{ij}(t_p)) J_n(k_F \mathbf{\Delta}_j(t_p, t_k)) e^{i n \theta_{i,j}(t_p, t_k)} 
\end{align} 
\tag{7} $$

In the first equality, we only added and subtracted $\mathbf{r}_j(t_p)$. To get the second equality, we use the notation defined above. Then, we use the Graf's addition theorem [72] to obtain the last equality. Then, to get the total wave field as in Eq. 2, we still have to perform a double summation over both the past positions (indexed by k) and droplets pairs ((indexed by i and j)). Because of both this double summation and the erratic form of the drops' individual dynamics, the terms in Eq. 7 containing a complex phase rapidly cancel out and so give rise to only small contributions to the wave field. As a consequence, we expect the term n = 0 in Eq. 7, namely $U_{0,j}(r_{ij}) = J_0(k_F r_{ij}(t_p)) J_0(k_F \mathbf{\Delta}_j(t_p, t_k))$, to be the leading term. As the wave-induced force, $F_{\text{wave},i}$ in Eq. 2 derives from the gradient of the wave field, we expect the equilibrium solution to be given by the extrema of
$U_{0,j}(r_{ij})$, apart from stability aspects. This can be found by looking for the zeros of $r \mapsto J_1(2\pi r/\lambda_F)$ which are

$$r^*/\lambda_F \approx 0, 0.63, 1.13, 1.63, 2.13, 2.63, 3.13, \dots \quad (8)$$

The even zeros correspond to those found numerically up to the uncertainty measurement (see vertical dashed lines in Fig. 5b). We recall that the situation is unstable and does not correspond to a crystalline structure. We conjecture that the odd zeros correspond to pair modes that are linearly unstable similarly to that observed in [45]. The peaks in the pair correlation function are the signature of a dynamically evolving internal structure and we have established its statistical wave origin.

Finally, we numerically predict that the pair correlation function exhibits a maximum at the origin. Experimentally, it is obviously impossible as short distance interactions must be taken into account. We present briefly in this last paragraph the influence of a short-distance elastic repulsion interaction. In Fig. 7, we observe that the systems with repulsion have pair correlation function $g(r)$ vanishing when $r$ tends to 0, contrary to the system in the absence of repulsion. The stronger the repulsion, the more efficient the convergence towards 0 at $r = 0$, is. We also observe that including short distance interactions let the positions of the peaks, outside the origin, unchanged while their aspects are mainly unaltered.

## IV. CONCLUSION AND DISCUSSION

We have numerically investigated the dynamics of a large number of walking drops that are confined by an external potential well while the surface Faraday waves are created as if they were in an unbounded fluid domain. The drops interact between each other by the wave they collectively generate and evolve in a complex chaotic dynamics which explores the whole accessible domain. A statistical analysis of the dynamics shows that for a simulation time smaller than $10^4 T_F$ the drops probability density distribution is homogeneous up to a radial distance approximately located at the edge $R_{max} - \lambda_F$. A fine analysis of the pairwise correlation functions shows that while being dynamic, time-evolving and presenting many indications of a good mixing in the phase space, the system adopts, on average, preferred inter-drops distances whose origin has been rationalized by analysing the internal symmetry of the waves. Thus, this numerical investigation sheds light on a statistical many-body wave self-organization in an apparent erratic dynamics.

It is relatively surprising that the internal order evidenced in the pair correlation function do not lead to any identifiable signature in the position probability density function as reported in corral experiments [20]. However, we remain in a parameter regime with a low

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 5. Evolution of the pair correlation function $g$ for various (a) confining potential radii $R_{max}$, (b) memory parameters Me, and (c) numbers of drops $N$. The dashed lines in (a) represent the theoretical situation for which all the drops where homogeneously distributed in a circular domain of radius $R_{max}$. Vertical black lines in (b) represent the theoretical predictions for the position of the stable maxima (even solutions in Eq. 8)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 6. Schematics to illustrate the notations used in Eq. 7. The drops (dark grey circles) *i* or *j* (*i*,*j* = 1,...,N) follow their path marked in blue. Shaded grey circles indicate past positions.</figcaption></figure>

time, even if already reasonably long, is currently
restricted by computational costs. Developing a more
parallelizable numerical scheme, by taking inspiration
from the field of molecular dynamics, would give access
to both larger system sizes and longer simulation
times for which the diluted phases of walkers may not
necessarily remain homogeneous. For example, in these
regimes, can we observe an ordering of the active gas
phase or a condensate? The exploration of these still
hypothetical regimes is the natural next step.

V. DATA AVAILABILITY STATEMENT

The C++ numerical code to adapt to various platforms with a Readme file is available at the link repository https://mycore.core-cloud.net/index.php/s/qlFPCxcMHDEjQr6 .
Supplementary Movie 1: movie associated to Fig. 2a.
Supplementary Movie 2: movie associated to Fig. 2b.

VI. AUTHORS CONTRIBUTIONS

Conceptualization: ML. Methodology: AH and ML.
Investigation: AH and ML. Visualization: AH. Project
administration: ML. Supervision: ML. Writing - original
draft: ML. Writing - review and editing: AH, ML.

density of drops and low memory, while the simulation

[1] Y. Couder, S. Protière, E. Fort, and A. Boudaoud, Na-574
ture 437, 208 (2005).

[2] J. W. M. Bush, Annual Review of Fluid Mechanics 47, 576
269 (2015).

[3] J. W. M. Bush and A. U. Oza, Report on Progress in Physics 84, 017001 (2021).

[4] J. Walker, Scientific American 238, 151 (1978).

[5] S. Protière, A. Boudaoud, and Y. Couder, Journal of Fluid Mechanics 554, 85 (2006).

[6] Y. Couder, E. Fort, C. Gautier, and A. Boudaoud, Physical Review Letters 94, 177801 (2005).

[7] N. Vandewalle, D. Terwagne, K. Mulleners, T. Gilet, and S. Dorbolo, Physical Review Letters 100 (2008).

[8] A. Eddi, E. Sultan, J. Moukhtar, E. Fort, M. Rossi, and Y. Couder, Journal of Fluid Mechanics 674, 433 (2011).

[9] J. Moláček and J. W. M. Bush, Journal of Fluid Mechanics 599, 727, 612 (2013).

[10] J. Moláček and J. W. M. Bush, Journal of Fluid Mechanics 727, 582 (2013).

[11] A. U. Oza, R. R. Rosales, and J. W. M. Bush, Journal of Fluid Mechanics 737, 552 (2013).

[12] P. A. Milewski, C. A. Galeano-Rios, A. Nachbin, and J. W. M. Bush, Journal of Fluid Mechanics 778, 361 (2015).

[13] M. Durey and P. A. Milewski, Journal of Fluid Mechanics 598, 821, 296 (2017).

[14] M. Faraday, Philosophical Transactions of the Royal Society of London 121, 299 (1831).

[15] J. Miles and D. Henderson, Annual Review of Fluid Mechanics 22, 143 (1990).

[16] T. B. Benjamin and F. Ursell, Proceedings of the Royal Society London A 225, 505 (1954).

[17] Y. Couder and E. Fort, Phys. Rev. Lett. 97, 1 (2006).

[18] G. Pucci, D. M. Harris, L. M. Faria, and J. W. M. Bush, Journal of Fluid Mechanics 835, 11361156 (2018).

[19] C. Ellegaard and M. T. Levinsen, Phys. Rev. E 102, 023115 (2020).

[20] D. M. Harris, J. Moukhtar, E. Fort, Y. Couder, and J. W. M. Bush, Physical Review E 88, 011001 (2013).

[21] T. Gilet, Physical Review E 90, 052917 (2014).

[22] T. Gilet, Physical Review E 93 (2016), 10.1103/Phys-RevE.93.042202.

[23] P. J. Sáenz, T. Cristea-Platon, and J. W. M. Bush, Nature Physics 14, 315 (2018).

[24] M. Durey, P. A. Milewski, and Z. Wang, Journal of Fluid Mechanics 891, A3 (2020).

[25] E. Fort, A. Eddi, A. Boudaoud, J. Moukhtar, and Y. Couder, Proceedings of the National Academy of Sciences 107, 17515 (2010).

[26] M. Labousse, A. U. Oza, S. Perrard, and J. W. M. Bush, Physical Review E 93 (2016).

[27] A. U. Oza, D. M. Harris, R. R. Rosales, and J. W. M. Bush, Journal of Fluid Mechanics 744, 404 (2014).

[28] S. Perrard, M. Labousse, M. Miskin, E. Fort, and Y. Couder, Nature Communications 5, 3219 (2014).

[29] S. Perrard, M. Labousse, E. Fort, and Y. Couder, Physical Review Letters 113, 104101 (2014).
<figure><img src="image_10.png" /><figcaption>FIG. 7. Evolution of pairwise density function *g* for the system *N* = 30, *Me* = 30, *R*<sub>max</sub> = 5*λ*<sub>*F*</sub> for various repulsion strengths. Color code: from light grey to black, the spring constant per unit mass is *K* = 500, 1000, 1500, 2000, 2500 s<sup>-1</sup>·kg<sup>-1</sup>. Darker curves correspond to higher repulsions. The salmon curve is the case in the absence of repulsion. Inset: zoom at low inter-drops distance. The blue vertical dashed line corresponds to the drop size in *λ*<sub>*F*</sub> units.</figcaption></figure>

[30] M. Labousse, S. Perrard, Y. Couder, and E. Fort, New Journal of Physics 16, 113027 (2014).

[31] P. J. Sáenz, T. Cristea-Platon, and J. W. M. Bush, Science Advances 6 (2020), 10.1126/sciadv.aay9234.

[32] M. Hubert, S. Perrard, M. Labousse, N. Vandewalle, and Y. Couder, Physical Review E 100, 032201 (2019).

[33] V. Bacot, S. Perrard, M. Labousse, Y. Couder, and E. Fort, Physical Review Letters 122, 104303 (2019).

[34] M. Durey, P. A. Milewski, and J. W. M. Bush, Chaos:672
An Interdisciplinary Journal of Nonlinear Science 28, 673
096108 (2018).

[35] M. Durey, Chaos: An Interdisciplinary Journal of Non-675
linear Science 30, 103115 (2020).

[36] M. Durey and J. W. M. Bush, Chaos: An Interdisci-681
plinary Journal of Nonlinear Science 31, 033136 (2021).

[37] M. Durey, S. Turton, and J. W. M. Bush, Proceedings685
of the Royal Society A 476, 2239 (2020).

[38] O. Devauchelle, E. Lajeunesse, F. James, C. Josserand,689
and P. Lagrée, Comptes Rendus. Mécanique 438, 591682
(2020).

[39] R. N. Valani, A. C. Slim, D. M. Paganin, T. P. Simula,684
and T. Vo, Physical Review E 104, 015106 (2021).

[40] M. Hubert, S. Perrard, N. Vandewalle, and M. Labousse,686
Nature Communications 13, 4357 (2022).

[41] S. Turton, M. Couchman, and J. Bush, Chaos 28, 096111 (2018).

[42] C. Borghesi, J. Moukhtar, M. Labousse, A. Eddi, E. Fort,690
and Y. Couder, Physical Review E - Statistical, Non-691
linear, and Soft Matter Physics 90, 063017 (2014),692
arXiv:1412.7701.

[43] R. N. Valani and A. C. Slim, Chaos 28, 096114 (2018). 694

[44] J. Arbelais, A. U. Oza, and J. W. M. Bush, Physical Review Fluids 3, 013604 (2018).

[45] M. M. P. Couchman, S. E. Turton, and J. W. M. Bush, Journal of Fluid Mechanics 871, 212 (2019).

[46] A. Eddi, J. Moukhtar, S. Perrard, E. Fort, and Y. Couder, Phys. Rev. Lett. 108 (2012).

[47] A. U. Oza, E. Siéfert, D. M. Harris, J. Moláček, and J. W. M. Bush, Physical Review Fluids 2, 053601 (2006).

[48] K. Papatryfonos, M. Ruelle, C. Bourdiol, A. Nachbin, B. J. W. M., and M. Labousse, Communication Physics 5, 142 (2022).

[49] V. Frumkin, J. M. Bush, and K. Papatryfonos, Phys. Rev. Lett. (2023).

[50] K. Papatryfonos, L. Vervoort, A. Nachbin, M. Labousse, and J. W. M. Bush, ArXiv preprint (2022).

[51] A. Eddi, A. Decelle, E. Fort, and Y. Couder, Europhys. Lett. 87 (2009).

[52] A. Eddi, A. Boudaoud, and Y. Couder, Euro. Phys. Lett. 94 (2011).

[53] M. M. P. Couchman, D. J. Evans, and J. W. M. Bush, Symmetry 14, 1524 (2022).

[54] B. Filoux, M. Hubert, and N. Vandewalle, Physical Review E 92, 041004 (2015).

[55] B. Filoux, M. Hubert, P. Schlagheck, and N. Vandewalle, Phys. Rev. Fluids 2, 013601 (2017).

[56] N. Vandewalle, B. Filoux, and M. Hubert, arXiv preprint (2019).

[57] M. M. P. Couchman and J. W. M. Bush, Journal of Fluid Mechanics 903, A49 (2020).

[58] S. J. Thomson, M. M. P. Couchman, and J. W. M. Bush, Physical Review Fluids 5, 083601 (2020).

[59] S. J. Thomson, M. Durey, and R. R. Rosales, Pro-
ceedings of the Royal Society of London Series A **476**,
20200155 (2020).

[60] S. J. Thomson, M. Durey, and R. R. Rosales, Phys. Rev. E 103, 062215 (2021).

[61] L. Barnes, G. Pucci, and A. U. Oza, Comptes Rendus Mécanique 348, 573 (2020).

[62] P. J. Sáenz, G. Pucci, S. E. Turton, A. Goujon, R. R. Rosales, J. Dunkel, and J. W. M. Bush, Nature 596, 58 (2021).

[63] C. A. Galeano-Rios, P. A. Milewski, and J.-M. Vandenbergheek, Journal of Fluid Mechanics 873, 856 (2019).

[64] L. Rayleigh, Phil. Mag. 16 (1883).

[65] S. Douady, J. Fluid Mech. 221, 383 (1990).

[66] K. Kumar and L. S. Tuckerman, Journal of Fluid Mechanics 279, 49 (1994).

[67] K. Kumar, Proceedings of the Royal Society A 452, 1113 (1996).

[68] S. Protière, A. Boudaoud, and Y. Couder, Journal of Fluid Mechanics 554, 85 (2006).

[69] L. Tadrist, J.-B. Shim, T. Gilet, and P. Schlagheck, Jour-
nal of Fluid Mechanics 848, 906 (2018).

[70] M. Labousse, *Etude d'une dynamique mémoire de chemin: une expérimentation théorique*, Ph.D. thesis, Université Pierre et Marie Curie-Paris VI (2014).

[71] L. Tambasco and J. Bush, Chaos (Focus Issue: Hydrodynamic Quantum Analogs) 28, 096115 (2018).

[72] F. Olver, D. Lozier, B. R.F., and C. Clark, NIST Handbook of Mathematical Functions (Cambridge University Press, 2010) Chap. 10, p. 247.