State space geometry of the chaotic pilot-wave hydrodynamics

TL;DR

This paper develops a symmetry-reduced dynamical systems framework to analyze the chaotic behavior of a bouncing droplet in a central potential, revealing bifurcations, chaos onset, and angular momentum reversals.

Contribution

It introduces a rotation symmetry reduction method applied to pilot-wave hydrodynamics, enabling detailed analysis of bifurcations and chaos in the system.

Findings

Identification of local bifurcations leading to chaos

Description of chaotic region emergence and merging bifurcations

Observation of spontaneous angular momentum sign changes

Abstract

We consider the motion of a droplet bouncing on a vibrating bath of the same fluid in the presence of a central potential. We formulate a rotation symmetry-reduced description of this system, which allows for the straightforward application of dynamical systems theory tools. As an illustration of the utility of the symmetry reduction, we apply it to a model of the pilot-wave system with a central harmonic force. We begin our analysis by identifying local bifurcations and the onset of chaos. We then describe the emergence of chaotic regions and their merging bifurcations, which lead to the formation of a global attractor. In this final regime, the droplet's angular momentum spontaneously changes its sign as observed in the experiments of Perrard et al. (Phys. Rev. Lett., 113(10):104101, 2014).

Tables

Table 1. Table 1: Summary of invariant pilot-wave solutions that are studied in this paper. Each solution has a reflection-copy with the exception of PO2l PO2l {\rm PO{2l}} , which itself is symmetric under the reflection ( 53 ). The last column shows the Me value at which the trajectories are plotted.

Name	Type	Shape	Trajectory	Reduced trajectory	at Me
$\mathrm{REQ0c}$	Relative equilibrium	Circle			$10.0$
$\mathrm{RPO1o}$	Relative periodic orbit	Oval			$15.0$
$\mathrm{RPO4os}$	Relative periodic orbit	Oval			$18.0$
$\mathrm{RPO4ou}$	Relative periodic orbit	Oval			$18.0$
$\mathrm{PO2l}$	Periodic orbit	Lemniscate			$14.0$
$\mathrm{RPO2l}$	Relative periodic orbit	Lemniscate			$16.0$
$\mathrm{RPO1t}$	Relative periodic orbit	Trefoil			$14.0$

Equations

$$a(\tau )=\left(\mathbf{r}(\tau ),\mathbf{v}(\tau ),h(r,\theta ;\tau )\right)$$

$$a(\tau )=f^{\tau }(a(0)).$$

$$\dot{a}=v(a)=\underset{h\to 0}{\lim }\left(f^h(a)-a\right)/h$$

$$a^{\prime }$$

$$R(\varphi )=\left(\begin{array}{cc}\cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{array}\right).$$

$$f^{\tau }(g(\varphi )a)=g(\varphi )f^{\tau }(a).$$

$$a_{\mathrm{REQ}}=g(-c\tau )f^{\tau }(a_{\mathrm{REQ}}).$$

$$a_{\mathrm{RPO}}=g(-{\varphi }_{\mathrm{RPO}})f^{T_{\mathrm{RPO}}}(a_{\mathrm{RPO}}),$$

$$\hat{a}(\tau )=g(-\hat{\varphi }(\tau ))a(\tau )$$

$$\hat{\varphi }=\text{mbox}Arg(v_x(\tau )+i{v}_y(\tau )).$$

$$\hat{a}=\left(\hat{x},\hat{y},{\hat{v}}_x,0,\hat{h}(r,\theta ;t)\right).$$

$$t(a)$$

$$T=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)$$

$${\hat{a}}^{\prime }=(\hat{x}=0,\hat{y}=0,{\hat{v}}_x=1,{\hat{v}}_y=0,\hat{h}(r,\theta ;\tau )=0),$$

$$⟨\hat{a}-{\hat{a}}^{\prime },t^{\prime }⟩=0\text{mbox}and⟨t(\hat{a}),t^{\prime }⟩>0,$$

$$⟨a_m,a_n⟩$$

$$\delta \hat{a}=H(\hat{a})\delta a,$$

$$H(\hat{a})=\mathrm{I}-\frac{t(\hat{a})\otimes t^{\prime }}{⟨t(\hat{a}),t^{\prime }⟩},$$

$$\hat{v}(\hat{a})$$

$$\dot{\hat{\varphi }}(\hat{a})$$

$$\hat{v}({\hat{a}}_{\mathrm{REQ}})=0.$$

$$\hat{A}({\hat{a}}_{\mathrm{REQ}}){\hat{e}}_i={\lambda }_i{\hat{e}}_i,$$

$$\hat{A}(\hat{a})=A(\hat{a})-\frac{t(\hat{a})\otimes A(\hat{a}{)}^T{t}^{\prime }}{⟨t(\hat{a}),t^{\prime }⟩}+\dot{\hat{\varphi }}\frac{t(\hat{a})\otimes {\mathcal{T}}^T{t}^{\prime }}{⟨t(\hat{a}),t^{\prime }⟩}-\dot{\hat{\varphi }}\mathcal{T},$$

$${\hat{a}}_{\mathrm{RPO}}={\hat{f}}^{T_{\mathrm{RPO}}}({\hat{a}}_{\mathrm{RPO}}),$$

$${\hat{J}}^{T_{\mathrm{RPO}}}({\hat{a}}_{\mathrm{RPO}}){\hat{V}}_i={\Lambda }_i\widehat{V_i}.$$

$$m\ddot{\mathbf{r}}+D\dot{\mathbf{r}}=-mg\nabla h(\mathbf{r},\tau )-k\mathbf{r},$$

$$h(\mathbf{r},\tau )=A\sum _{n=-\infty }^{\lfloor \tau /T_F\rfloor }{J}_0(k_F|\mathbf{r}-\mathbf{r}(n{T}_F)|)e^{-(\tau -n{T}_F)/T_F\text{Me}},$$

$$h(\mathbf{r},\tau )=\frac{A}{T_F}{\int }_{-\infty }^{\tau }{J}_0(k_F|\mathbf{r}-\mathbf{r}({\tau }^{\prime })|)e^{-(\tau -{\tau }^{\prime })/T_F\text{Me}}d{\tau }^{\prime }.$$

$$h(\mathbf{r},\tau )=A\sum _{n=0}^{\infty }$$

$$C_n(\tau )$$

Full text

hydrodynamic quantum analogs, symmetry reduction

During the quantum physics’ infancy, Louis de Broglie deBroglie1 imagined it as the outcome of the dynamics of a point-like particle that is interacting with a continuous background field. This viewpoint was to a large extent forgotten during the second half of the twentieth century due to the great success of the Copenhagen interpretation. In recent years, a resurgent interest in de Broigle’s “wave-particle duality” has developed as a result of its discovery in a completely different field: fluid mechanics. In a series of experiments pioneered by Couder et al. CFGB2005 ; CPFB2005 , several phenomena that were once thought to be exclusive to quantum physics were demonstrated in the macroscopic setting of a droplet of silicon oil bouncing on vertically vibrating bath of the same fluid. We examine a class of these systems with rotational symmetry and formulate a novel method for their analysis. We demonstrate the utility of our approach in a numerical study where we describe the chaotic dynamics of a hydrodynamic pilot-wave model through a geometrical approach, in which we identify solutions with qualitative differences and intermittent transitions in between.

I Introduction

In 2005, Couder et al. CFGB2005 showed that a droplet of silicon oil can indefinitely bounce on the bath of the same fluid when the bath is vibrating vertically with an acceleration close to but below the Faraday instability point, at which surface waves appear spontaneously. In the same year, Couder et al. CPFB2005 also reported that a transition from “bouncing” to “walking” takes place when the bath’s oscillation amplitude is large enough. In this latter state, the droplet becomes a “walker”, and the system exhibits a macroscopic “wave-particle duality”: The force that the droplet experiences at each bounce is determined by the shape of the bath’s surface, which itself is determined by the waves generated at the droplet’s previous bounces. A year later, Couder and Fort reported CouderFort2006 that when these walkers were realized in single- or double-slit geometries, their final positions obeyed a statistical distribution akin to the single-particle diffraction experiments of quantum particles111 The single- and double-slit experiments with walkers were revisited by Andersen et al. AMRRLT2015 and Pucci et al. PHFB2018 and the outcomes of these experiments appear to be very sensitive to the experimental conditions. .

Early experiments with bouncing droplets were followed by others that demonstrated various analogies to the quantum mechanics Bush2015 . Fort et al. FEBMJC2010 and Harris and Bush HarrisBush2014 studied walkers in rotating frames and reported quantization of droplet trajectories. Similarly quantized orbits were also found PLMFC2014 when the two-dimensional pilot-wave system was realized with a central harmonic potential, which was generated by applying magnetic field on a droplet filled with a ferromagnetic fluid. Furthermore, the dynamics of bouncing droplets under central harmonic force was shown to become chaotic PLFC2014 . These systems are of particular interest to us due to their rotation symmetry: Assuming the bath’s walls are of circular shape (or far away), governing laws of motion preserve their shape under the rotation of coordinates. As we shall illustrate in this paper, the rotation symmetry translates to a redundant degree of freedom in the system’s state space.

Experimental developments in bouncing droplets were accompanied by modeling efforts. Couder et al. CPFB2005 already predicted the transition from bouncing to walking as a pitchfork bifurcation by postulating a sinusoidal force to be exerted on the droplet by the bath during the bounces. Through a detailed analysis, Oza et al.  ORB2013 reached an integro-differential equation of motion for the bouncing droplets and coined the term “pilot-wave hydrodynamics”. Numerical investigations of this model with confining central-potential terms (Coriolis, harmonic, or Coulomb type) demonstrated different routes to chaos THORB2016 and rich subsequent dynamics. Similar to their experimental counterparts, these systems also exhibit rotation symmetry and are at the focus of the current paper.

The references that we cited so far are the key developments in the pilot-wave hydrodynamics literature which are relevant to our work. For a thorough review, we refer the reader to the recent article by Bush Bush2015 . We are now going to turn our attention to a different branch of literature, where chaotic dynamics take place in the presence of continuous symmetries.

In recent years, numerical studies of nonlinear partial differential equations such as the complex Ginzburg–Landau equation Luce95 , the Kuramoto–Sivashinsky equation SCD07 ; BudCvi15 , and the Navier–Stokes equations ACHKW11 ; WFSBC15 , demonstrated that degeneracies due to the continuous symmetries in these systems tend to obscure their dynamics. Take, for example, the fluid flow through a cylindrical pipe. This system is symmetric under axial translations and azimuthal rotations. Each generic (non-symmetric) solution of pipe flow has infinitely many copies that can be obtained by continuous symmetry transformations, all of which correspond to the same physical solution. Budanur et al. BudCvi14 showed that a simple polar-coordinate transformation in the Fourier space representation of these systems can be used to eliminate these redundant degrees of freedom and such a transformation can be interpreted as a “slice”, that is a codimension-1 hyperplane in the state space of the system. This reformulation allowed for the straightforward use of the established geometrical methods of the slicing literature, which were mostly developed within the equivariant bifurcation theory ChossLaut00 . Budanur et al. BudCvi14 named this method “first Fourier mode slice” and applied it to the Kuramoto–Sivashinsky equation. Later on, this method was successfully adapted for the simulations of two-dimensional Kolmogorov flow Faraz15 and three-dimensional pipe flow BudHof17 ; BudHof18 . In this paper, we formulate this method for pilot-wave systems and demonstrate its utility in a numerical study.

The rest of this paper is organized as follows. In the next section, we formulate a symmetry-reduced representation of the pilot-wave hydrodynamics with rotation symmetry. At this stage, our treatment will be general and applicable to both experimental and numerical data. In section III, we describe the trajectory equation of Oza et al. ORB2013 with an additional central harmonic potential and apply our symmetry-reduction scheme to this system. We analyze the system’s local bifurcations, the transition to chaos via period doubling cascades, and global bifurcations of the system’s chaotic sets. We discuss our results in section IV.

II Pilot-wave hydrodynamics as an infinite-dimensional dynamical system

We are going to investigate the dynamics of a point-like droplet as it bounces on the vertically vibrating bath of the same fluid in the presence of a central force. Let us denote the two-dimensional position and velocity of the droplet at time $\tau$ with ${\bf r(\tau)}=(x(\tau),y(\tau))$ and ${\bf v(\tau)}=(v_{x}(\tau),v_{y}(\tau))$ respectively, and the fluid height on the bath surface as a function of polar coordinates $r$ and $\theta$ at time $\tau$ with $h(r,\theta;\tau)$. We assume that the future time-evolution of the system is completely determined by the current position and velocity of the droplet and the shape of the bath’s surface. Under this assumption, we can define a state space vector

[TABLE]

and a finite-time flow mapping ${f^{\tau}(a)}$, which maps an initial state $a(0)$ to a final state $a(\tau)$ as

[TABLE]

Definitions (1) and (2) make up an infinite-dimensional dynamical system since $h(r,\theta;\tau)$ takes values in a function space. In numerical applications, $h(r,\theta;\tau)$ would be expressed as a finite-dimensional approximation such as a discretization on a grid, or a spectral expansion. In order to illustrate what we mean by these definitions, we show two snapshots of a simulated (details will be presented in section III) circular trajectory of the pilot-wave system with harmonic potential in figure 1. Snapshots in this figure are separated by $20$ nondimensional units of time; and each panel contains the instantaneous information necessary to determine the future evolution of the system. According to our definitions, panel (a) ($a(0)$) is mapped to panel (b) ($a(20)$) by the flow as $a(20)={f^{20}(a(0))}$.

If we further assume that the flow (2) is smooth, then we can also define the state space velocity

[TABLE]

and represent the system as an infinite set of ordinary differential (ODE) equations. In general, especially for experimental data, it may not be possible to obtain an explicit ODE form (3). Nevertheless, if the dynamics is sufficiently smooth, we can assume that the state space velocity $v(a)$ (3) exists and we are going to use its definition (3) in what follows.

II.1 The rotation symmetry and the relative invariant solutions

While we have not specified the exact physics that governs the motion of a bouncing droplet yet, the conditions that we described in the previous section along with the additional assumption that the boundaries are far away (or circular) imply rotation symmetry. Since the shape of the bath’s surface is determined by the past trajectory and the forcing is only radially dependent; the physics of the problem stays the same under the rotation of coordinates. This, however, does not imply that the solutions would be symmetric under rotations. In fact, it is already visible on the example of figure 1 that the individual snapshots of the circular trajectory do not have rotation symmetry. However, for each generic solution of the system, there exist infinitely many physically equivalent copies that can be obtained by continuous rotations. In the dynamical systems theory, this is called “equivariance”, which we will now formalize for the pilot-wave system.

Let $g(\phi)$ be the rotation operator which acts on a state vector (1) as

[TABLE]

where $R(\phi)$ are the standard $2\times 2$ rotation matrices

[TABLE]

If for every trajectory $a(t)$ with $t\in[0,t_{f}]$ there exists another one $g(\phi)a(t)$ obtained by the symmetry action, then the finite-time flow (or equivalently the state space velocity, when it exists) commutes with the symmetry action, i.e.

[TABLE]

We are going to refer to (6) as the “equivariance condition”. Note that the equivariance condition does not imply that the flow is “invariant” under the action of the symmetry, but rather, it transforms as does the state space coordinates themselves. In our treatment, we are going to assume that the pilot-wave systems that we consider satisfy the equivariance condition (6) with the rotation operator defined as in (4).

One way of elucidating the state space geometry of chaotic systems is identifying the time-invariant solutions (equilibria, periodic orbits, invariant tori …) that are embedded in the system’s chaotic set DasBuch . In systems with continuous symmetry, roles of these solutions are taken over by their higher-dimensional equivalents, which we will refer to as “relative invariant solutions”. Specifically, a “relative equilibrium” is a solution that evolves only in the symmetry direction at a constant phase speed $c$. Let $a_{{\rm REQ{}}}$ be a state space point on a relative equilibrium, then its time-evolution satisfies

[TABLE]

Note that a time-invariant set is formed by all solutions on the relative equilibrium, i.e. the one-dimensional curve $\{g(c\tau)a_{{\rm REQ{}}}:\tau\in[0,2\pi/c)\}$. The second type of relative invariant solutions we are going to consider here are “relative periodic orbits”, which have trajectories that recur at a symmetry-transformed location after a finite period $T_{{\rm RPO{}}}$, i.e.

[TABLE]

where $a_{{\rm RPO{}}}$ is a point on a relative periodic orbit, $T_{{\rm RPO{}}}$ is its period, and $\phi_{{\rm RPO{}}}$ is its phase-shift. Similar to the relative equilibrium, all points on a relative periodic orbits, i.e. the 2-torus $\{g(\phi){f^{\tau}(a_{{\rm RPO{}}})}:\phi\in[0,2\pi)\,,\tau\in[0,T_{{\rm RPO{}}})\}$, form a time-invariant set. In other words, the relative invariant solutions that appear as a consequence of a single-parameter continuous symmetry correspond to manifolds with dimensions 1 higher than those of their standard counterparts. Notice that the special case $c=0$ in (7) corresponds to an equilibrium solution, and similarly $\phi_{{\rm RPO{}}}=0$, in (8) defines a periodic orbit.

It is possible to define the linear stability of relative invariant solutions by defining Jacobians of (7) and (8) ChossLaut00 ; DasBuch . However, we are going to skip these definitions here as they will be redundant once we introduce the symmetry reduction in the next section.

II.2 Continuous symmetry reduction

Symmetry reduction is a state space coordinate transformation $a\rightarrow\hat{a}$ such that each symmetry-related (physically equivalent) set of states $\{g(\phi)a^{*}\,|\,\phi\in[0,2\pi)\}$ is represented by a single solution $\hat{a}^{*}$ in the symmetry-reduced state space. An obvious choice of such coordinates for the droplet system with the rotation symmetry is the time-dependent polar coordinate transformation that fixes the polar-angle in $(x,y)$ plane to a certain arbitrary value. This angle, however, is not defined when the droplet is at the origin, and would experience very rapid phase fluctuations if the droplet is close to the origin. Observing that the droplet’s horizontal velocity is always non-zero CPFB2005 , we define the following operation

[TABLE]

where

[TABLE]

This transformation fixes the polar angle in the $(v_{x},v_{y})$ plane to [math], hence maps all symmetry-equivalent solutions to the one with $v_{y}=0$ as illustrated in figure 2(a), as long as $||{\bf v}||\neq 0$. Thus, the symmetry-reduced state space coordinates $\hat{a}$ have $\hat{v}_{y}$ coordinate identically equal to [math], i.e.

[TABLE]

Budanur et al. BudCvi14 showed that a polar coordinate transformation similar to (9), when defined for a partial differential equation with translation symmetry, can be interpreted as a projection onto a codimension-1 hyperplane, named a “slice”, in the state space. We are now going to take analogous steps and redefine the transformation (9) in this terminology in order to exploit the well-established tools of this method.

Let $\mathcal{T}$ be the generator of infinitesimal rotations satisfying $g(\phi)=e^{\phi\mathcal{T}}$. Its action on the state space coordinates can be expressed as

[TABLE]

where $T$ is the generator of the $2\times 2$ rotations (5)

[TABLE]

and $t(a)$ is called the “group tangent” evaluated at $a$, since it gives the direction of an infinitesimal rotation at $a$. Let us now define the “slice template” as the state vector

[TABLE]

whose only nonzero component is $\hat{v}_{x}=1$. It is straightforward to verify that the symmetry-reduced state space coordinates (9) satisfy the half-hyperplane condition

[TABLE]

where $t^{\prime}=t(\hat{a}^{\prime})$ and $\left\langle.,\,.\right\rangle$ is the inner product defined as

[TABLE]

where w(r) is an appropriately chosen weight function and the integral is computed over the whole domain $\mathcal{D}$. Projection of the full state space dynamics onto this codimension-1 half-hyperplane is illustrated in figure 2 (b). In the slicing method, one looks for the slice phases $\hat{\phi}(\tau)$, such that $\hat{a}(\tau)=g(-\hat{\phi}(\tau))a(\tau)$ satisfies (15). The main advantage of slicing framework is that it allows for a straightforward symmetry-reduction of the tangent space. Let $\delta a$ be a small perturbation (in the full state space) to $\hat{a}$; it is projected onto the slice as

[TABLE]

where

[TABLE]

${\rm I}$ is the identity matrix and $\otimes$ denotes the outer product. For a derivation, see the appendix of Budanur & Hof (2017) BudHof17 . In particular, we can obtain the symmetry reduced state space velocity $\dot{\hat{a}}=\hat{v}(\hat{a})$ as

[TABLE]

By taking the time derivative of $\hat{a}=g(-\hat{\phi}(\tau))$ and requiring the slice hyperplane condition (15) to be satisfied, one can show that the multiplier of the group tangent in (19) is in fact the slice phase velocity. In other words, we can write

[TABLE]

and integrate the slice phase $\hat{\phi}$ along with the symmetry-reduced evolution. Rowley and Marsden rowley_reconstruction_2000 referred to (20) and (21) as “reconstruction equations” since by integrating them simultaneously, one can reconstruct the full state space evolution. In other words, symmetry reduction does not lose any information, it rather separates the time-evolution into the “physical” (in-slice) and “redundant” (phase) parts. Another observation to make in (21) is that it diverges when $\left\langle t(\hat{a}),\,t^{\prime}\right\rangle=0$. This, in general, implies that a hyperplane slice that is constructed around an arbitrary template $\hat{a}^{\prime}$ is “local”, and applicable in a finite neighborhood where the denominator of (21) is nonzero. Notice, however, that for our particular choice of the slice template (14), this condition becomes $\left\langle t(\hat{a}),\,t^{\prime}\right\rangle=\hat{v}_{x}=\sqrt{v_{x}^{2}+v_{y}^{2}}$; thus, our slice is applicable as long as the droplet’s velocity does not vanish.

In the symmetry-reduced state space, a relative equilibrium (7) becomes an equilibrium, which satisfies

[TABLE]

Linear stability of this solution is determined by the eigenvalue equation

[TABLE]

where $\hat{A}_{ij}(\hat{a}^{*})=\partial\hat{v}_{i}(\hat{a})/\partial\hat{a}_{j}|_{\hat{a}=\hat{a}^{*}}$. By straightforward algebra, this matrix is obtained as

[TABLE]

where $A_{ij}(a^{*})=\partial v_{i}(a)/\partial a_{j}|_{a=a^{*}}$ . We are going to refer to $A(a)$ and $\hat{A}(\hat{a})$ respectively as the stability matrix and the reduced stability matrix. $\lambda_{i}$ are the linear stability eigenvalues and $\hat{e}_{i}$ are the corresponding stability eigenvectors within the slice. ${\rm Re}\,\lambda_{i}<0$ and ${\rm Re}\,\lambda_{i}>0$ respectively correspond to the stable and unstable directions associated with the relative equilibrium and a bifurcation takes place when ${\rm Re}\,\lambda_{i}=0$.

Similar to the relative equilibria, relative periodic orbits (8) satisfy the periodicity condition

[TABLE]

where $\hat{a}(\tau)={\hat{f}^{\tau}(\hat{a}(0))}$ is the symmetry-reduced finite-time flow induced by the dynamics. Consequently, the linear stability of a relative periodic orbit is determined by the eigenvalue equation

[TABLE]

where $\hat{J}^{\tau}_{ij}(\hat{a}^{*})=\partial\hat{f}_{i}^{\tau}(\hat{a})/\partial\hat{a}_{j}|_{\hat{a}=\hat{a}^{*}}$. Numerically, the reduced Jacobian $\hat{J}^{\tau}_{ij}(\hat{a}^{*})$ can be obtained as the time-ordered integral of the reduced stability matrix $\hat{A}(\hat{a})$ along a symmetry-reduced state space trajectory. $\Lambda_{i}$ are known as the Floquet multipliers of the relative periodic orbit and $\hat{V}_{i}$ are the corresponding Floquet vectors in the symmetry-reduced state space. $|\Lambda_{i}|<1$ and $|\Lambda_{i}|>1$ correspond to the stable and unstable directions respectively and $|\Lambda_{i}|=1$ corresponds to a marginal direction. Note that the symmetry-reduced state space velocity $\hat{v}(\hat{a}_{{\rm RPO{}}})$ is an eigenvector of $\hat{J}^{T_{{\rm RPO{}}}}(\hat{a}_{{\rm RPO{}}})$ with the eigenvalue $1$.

This concludes our general introduction of the symmetry-reduction, which can – in principle – be applied to experimental, as well as numerical data, independent of the specifics of system in consideration, as long as the rotation symmetry is present. In the next section, we are going to demonstrate the utility of our transformation by applying it to the simulations of a pilot-wave model with central harmonic force.

III A pilot-wave model with central harmonic force

Following the efforts of Moláček and Bush MB2013a ; MB2013a , Oza et al. ORB2013 derived a trajectory equation for the bouncing droplets. With the addition of a central harmonic force $-k\mathbf{r}$, this equation reads

[TABLE]

where $\mathbf{r}=(x,y)$ is the position of the droplet, $m$ is its mass, $D$ is the viscous damping coefficient, $g$ is the acceleration due to gravity, $k$ is the effective spring constant of the harmonic force, and $h$ is the wave field on the bath surface determined by the previous bounces of the droplet as

[TABLE]

where $A$ is the initial amplitude of the waves generated, $J_{0}$ is the Bessel function of the first kind, $T_{F}$ is the time between the bounces, known as the ‘Faraday time’, $k_{F}=2\pi/\lambda_{F}$ is the wave number corresponding to the least stable Faraday wave, and $\lambda_{F}$ is the corresponding wavelength. The tunable parameter Me is called “memory” and it determines the damping speed of the waves generated in previous bounces. One of the approximations that were made by Oza et al. ORB2013 is already apparent from (28), where each bounce instantly generates a Bessel wave, with no time needed for propagation.

With the wave field given as in (28), the trajectory equation (27) is a delay-differential system, which is very hard to study both analytically and numerically. Therefore, Oza et al. ORB2013 and several subsequent studies LOPB2016 ; THORB2016 ; KOB2017 approximated the sum in (28) by the integral

[TABLE]

This approximation is justified since the time scale of the horizontal motion is much greater than the time scale of bouncing ORB2013 . At this stage, (27) is a four-dimensional integro-differential system. Labousse labousse2014thesis ; PerrardLabousse2018 showed that this system can further be manipulated into an infinite-dimensional ODEs and this is the form we are going to use in the following. It follows from Graf’s addition theorem watson1944 that (29) can be written as

[TABLE]

where

[TABLE]

We can obtain the time derivatives of the mode amplitudes $C_{n}$ and $S_{n}$ as

[TABLE]

Rearranging (27), the acceleration of the walker in $x$ and $y$ directions can be written explicitly as

[TABLE]

where the partial derivatives are

[TABLE]

with $K_{n},L_{n}$ and $M_{n}$ defined as

[TABLE]

Finally, we nondimensionalize trajectory equations (35) and (36) by scaling time as $\tau/T_{F}\rightarrow\tau$ and lengths as $(k_{F}x,k_{F}y)\rightarrow(x,y)$ and rewrite everything compactly as

[TABLE]

where,

[TABLE]

With the addition of trivial equations $\dot{x}=\dot{x}$ and $\dot{y}=\dot{y}$, equations (44–47) yields a set of first order ODEs (3) for the state space vector

[TABLE]

Formally, this state space is infinite dimensional. In the numerical work to follow, the Fourier-Bessel expansion (30) is truncated at a finite number of $N+1$ modes (starting counting from the [math]th mode), yielding a $2N+5$ dimensional representation. We define the $L_{2}$ inner product in this representation as

[TABLE]

For the state space coordinates ordered as in (49), we can write the rotation-operator (4) explicitly as

[TABLE]

with the subsequent $2\times 2$ blocks along its diagonal populated with the rotation matrices (5) $R(2\phi),R(3\phi),\ldots$ The corresponding generator of infinitesimal rotations then becomes

[TABLE]

In addition to the rotational symmetry, the pilot wave model with a central force is also equivariant under the reflection $y\rightarrow-y$, or equivalently $\theta\rightarrow-\theta$. The action of this transformation on the state space coordinates is represented by the diagonal matrix

[TABLE]

and the finite-time flow mapping implied by the time-evolution under (44–47) satisfies

[TABLE]

The model parameters (48) are determined by the experimental conditions as explained by Oza et al. ORB2013 and Perrard et al. PLMFC2014 . In the following numerical work, we are going to adopt the numbers reported by Tambasco et al. THORB2016 and set $m=0.25\times 10^{-6}kg$, $f=80Hz$, $D=2.0\times 10^{-6}kg/s$, $A=3.5\times 10^{-6}m$, and $k=3.2\times 10^{-6}N/m$. These choices yield the nondimensional system parameters

[TABLE]

In the next section, we are going to explore the bifurcations of this model and the onset of chaos numerically for $\textit{Me}\in[10,19]$. We observed for this parameter range that a truncation of the Fourier-Bessel expansion (30) at $N=25$ ensures at least $8$ order of magnitude drop in the mode amplitudes $||C_{n}+S_{n}||$, hence the numerical results that we are going to report adopts this resolution. Note that the resulting state space representation is $55$-dimensional. In the numerical results to follow, we simulate this $55$-dimensional ODE system (44–47) using odeint function of scipy.integrate module scipy of Python programming language. Note that odeint adapts the time-steps in order to keep estimated numerical integration errors below $10^{-8}$. In the results of the following sections, the trajectories are sampled with a time-step $\delta\tau=0.01$. For simulations in the symmetry-reduced state space, we integrate (19), which we obtain explicitly by projecting (44–47).

III.1 Local bifurcations and the onset of chaos

As we shall now demonstrate, the explicit ODE form of the dynamical equations (44–47) along with the symmetry reduced dynamics (20–21) allow us to use standard numerical analysis methods for a bifurcation study of the pilot-wave system. In this section and the rest of the paper, we are going to use the following naming convention for the solutions we describe. Initial capital letters will refer to the type of the solution: EQ: equilibrium, REQ: relative equilibrium, PO: periodic orbit, RPO: relative periodic orbit, C: chaotic. The following integer (except for chaotic solutions) will refer to the number of times the solution intersects the Poincaré section (60), which we are going to define and use for our illustrations later. The following lowercase letter will refer to the shape of the solution on the $(x,y)$-plane: c: circle, o: oval, l: lemniscate, t: trefoil. The solutions that are not reflection-symmetric will appear in pairs of positive and negative average angular momenta. We will refer to the ones with negative average angular momentum with an additional $\sigma$ in front. For example, the circular solution is a relative equilibrium that does not intersect the Poincaré section (60); thus, the circular solution with positive angular momentum will be ${\rm REQ{0c}}$ and its reflection copy with negative angular momentum will be $\sigma{\rm REQ{0c}}$. When needed, additional letters will be used to further specify the solutions and these will be explained in place. For reference, table 1 shows a summary of invariant solutions in the numerical bifurcation study to follow. In table 1 and the rest of the trajectory plots of this paper, the lengths are scaled by $2\pi$ so that one unit of distance is equal to one Faraday wavelength.

III.1.1 Bifurcations from the circular solution

As shown by Labousse et al. LOPB2016 and Tambasco et al. THORB2016 , when $\textit{Me}=10$, the dynamical system described by (44–47) has a stable circular orbit with radius $r\approx 0.8\lambda_{F}$, which we have visualized in figure 1. This solution is a relative equilibrium (7), which we are going to refer to as ${\rm REQ{0c}}$. ${\rm REQ{0c}}$ intersects the slice hyperplane (15) at a single point $\hat{a}_{{\rm REQ{0c}}}$ corresponding to the snapshot on figure 1(a). This solution satisfies the equilibrium condition (22) in the symmetry-reduced state space. Starting from the stable $\hat{a}_{{\rm REQ{0c}}}$ at $\textit{Me}=10.0$, we increased the Me in steps of $0.1$, find the root of (22) numerically using fsolve function from scipy.optimize scipy and compute the spectrum of (24) at this point for each Me value. In these computations, we discarded the $4^{th}$ column of $\hat{v}(\hat{a})$, which corresponds to the $\hat{v}_{y}$ direction that is identically set to [math] in the slice. Similarly, when computing the eigenvalues of (24) we discard the $4^{th}$ column and row.

The leading (with the largest ${\rm Re}\,\lambda$) two pairs of linear stability eigenvalues of the circular trajectory computed this way are plotted in figure 3(a). As shown, the leading pair of complex conjugate eigenvalues $\lambda_{1,2}$ cross the imaginary axis upon increasing Me. This corresponds to the interval $\textit{Me}\in(13.7,13.8)$ and at $\textit{Me}=13.8$, the circular trajectory is unstable. For verification, we translated this Me value to the to the forcing acceleration $\gamma$ which is related to Me by $\textit{Me}\ T_{F}=T_{d}/(1-\gamma/\gamma_{F})$, where $T_{d}$ is the decay time of waves in the absence of forcing and $\gamma_{F}$ is the Faraday instability threshold. With $T_{F}=2/f=0.025s$ and $T_{d}\approx 0.0182$ LOPB2016 , $\textit{Me}=13.8$ corresponds to $\gamma/\gamma_{F}\approx 0.9472$. This value is in agreement with Tambasco et al. THORB2016 who reported the destabilization of the circular orbit at $\gamma/\gamma_{F}=0.948$.

The destabilization of the circular orbit as Me is increased corresponds to a supercritical Hopf bifurcation as we confirm by the appearance of a relative periodic orbit in its vicinity. This is illustrated in figure 3(b), where we show an orbit at $\textit{Me}=15$ starting approximately on the unstable manifold of ${\rm REQ{0c}}$ with the initial condition

[TABLE]

In order to find this relative periodic orbit to the numerical precision, we implemented Newton’s method (Appendix A) for periodic orbits to find the roots of (25). Once again, discarding the $4$-th element of the equation, which is set to [math] by the symmetry reduction. We visualized the relative periodic orbit’s trajectory on $(x,y)$ plane along with the corresponding wave field at $\textit{Me}=15$ as three snapshots in figure 4. Perrard et al.  PLMFC2014 referred to the orbits with similar trajectories as “ovals”, which we will also adopt for the family of solutions to follow. At $\textit{Me}=15$, ${\rm RPO{1o}}$ is stable and its period is $T_{{\rm RPO{1o}}}\approx 40.04$, thus, snapshots in figure 4(a) and (c) are approximately one period apart. Notice that the wave fields in figure 4 (a) and (c) are almost the same up to a rotation of the coordinates.

As we did with the circular solution, we varied $\textit{Me}\in[14,19]$ with steps of $0.1$ and numerically solved (25) using Newton’s method and obtained the Floquet multipliers of ${\rm RPO{1o}}$ by computing the eigenvalues of (26) at each step. We plotted the leading (with largest $|\Lambda_{i}|$) three Floquet multipliers of ${\rm RPO{1o}}$ for $\textit{Me}\in[14,19]$ on figure 5(a). For each Me value, there is a marginal Floquet multiplier with $\Lambda=1$ corresponding to the perturbations along the periodic orbit, as expected. Leading non-marginal pair of complex conjugate Floquet multipliers cross the unit circle at $\textit{Me}\in(16.8,16.9)$ (annotated in figure 5(a)) and after this point, ${\rm RPO{1o}}$ is unstable.

In order to study the dynamics on the $({\rm Re}\,V_{1},{\rm Im}\,V_{1})$-plane in the vicinity of $\hat{a}_{{\rm RPO{1o}}}$ after it becomes unstable, we define a Poincaré section as the half-hyperplane

[TABLE]

where $\hat{a}_{{\rm RPO{1o}}}$ is an arbitrary point on ${\rm RPO{1o}}$ and $\hat{a}_{\cal P}$ denotes the state space points on the Poincaré section. We started a simulation with the initial condition

[TABLE]

and tracked its intersections with the Poincaré section (57). In order to visualize the dynamics on this Poincaré section, we projected it onto the orthonormal bases $(\hat{\xi}_{1},\hat{\xi}_{2})$ that span the $({\rm Re}\,V_{1},{\rm Im}\,V_{1})$-plane in figure 5(b) at $\textit{Me}=17$. As suggested by the appearance of an invariant 2-torus in the vicinity of ${\rm RPO{1o}}$, its leading complex conjugate Floquet multipliers’ crossing of the unit circle corresponds to a supercritical Neimark-Sacker bifurcation Kuznetsov2010 .

The invariant 2-torus surrounding ${\rm RPO{1o}}$ disappears at $\textit{Me}\in(17.1,17.2)$ and leaves its place to a stable and unstable pair of period-4 (on the Poincaré section) orbits. In order to elucidate this process, we analyzed the system before ($\textit{Me}=17.1$) and after ($\textit{Me}=17.2$) the bifurcation. At $\textit{Me}=17.1$, we parameterized the invariant torus surrounding ${\rm RPO{1o}}$ with $\theta_{{\cal P}}\in[0,2\pi)$ and generated a return map $\theta_{{\cal P}}[n+1]=f_{{\cal P}}(\theta_{{\cal P}}[n])$ as shown in figure 6(a). Upon our investigation of the fourth iterate of this map, we have found that at $4$ locations it comes close to being tangent to the identity map as plotted in figure 6(b). This is a typical situation as the system proceeds towards a saddle-node bifurcation that would give rise to a stable-unstable pair of period-$4$ orbits. In order to elucidate the dynamics after this bifurcation, we approximated ${\rm RPO{1o}}$’s unstable manifold by forward integrating the initial conditions

[TABLE]

where we used $4$-equidistant values for $\delta\in[0,1)$, $36$-equidistant values for $\phi\in[0,2\pi)$ and set $\epsilon=10^{-4}$. Time-forward dynamics of the initial conditions (59) approximately cover the linearized unstable manifold of ${\rm RPO{1o}}$, thus we expect their further evolution to take the shape of the nonlinear unstable manifold BudCvi15 . We visualized these orbits at their intersections with the Poincaré section (57) as projections onto $({\rm Re}\,V_{1},{\rm Im}\,V_{1})$-plane on figure 7. We found that all of these trajectories converged to a new stable periodic orbit (marked magenta on figure 7) that intersects the Poincaré section 4 times. We will refer to this orbit as ${\rm RPO{4os}}$ with the additional “s” signifying that this solution appears as the stable one of a saddle-node pair. We have also observed that this manifold was separated in two parts at 4 different locations, where nearby trajectories went in opposite directions when they arrived at the edge of the manifold as illustrated by red and black orbits on figure 7. By utilizing a Newton search, we confirmed that these separation points corresponded to an unstable period-4 orbit ${\rm RPO{4ou}}$ (“u” standing for “unstable”), which we marked yellow on figure 7. These observations together confirm that a saddle-node bifurcation on the invariant torus surrounding ${\rm RPO{1o}}$ took place at $\textit{Me}\in(17.1,17.2)$.

We then continued both ${\rm RPO{4ou}}$ and ${\rm RPO{4os}}$ up to $\textit{Me}=19.0$. For this computation, we needed to implemented pseudoarclength continuation DKK1991 since ${\rm RPO{4ou}}$ underwent several fold bifurcations and became more and more unstable. ${\rm RPO{4os}}$, on the other hand, underwent a period doubling bifurcation at $\textit{Me}\in(17.8,17.9)$, marking the beginning of the cascade which leads to chaos.

Until the period doubling point, the bifurcation sequence we described here qualitatively222The exact bifurcation points are slightly different, which could be due to differences in the numerical methods. agrees with the observations of Tambasco et al. THORB2016 , who studied the same system with a different numerical method. Notably, they also observed that the relative two-torus is replaced by a period-4 relative periodic orbit, which they refer to as “a frequency-locked wobbling orbit”. Upon further increase of Me, Tambasco et al. THORB2016 observes a new incommensurate frequency in the Fourier spectrum of the wobbling orbit, which does not agree with the period doubling bifurcation we find here. Consequently Tambasco et al.’s description to transition to chaos disagrees with ours, which we will present in section III.2. In section IV, we discuss the differences between two transition scenarios in more detail.

We summarized the bifurcations that we analyzed up to this point in figure 8. The vertical axis of figure 8 shows the $\hat{y}$ coordinate of $\hat{a}_{{\rm REQ{0c}}}$ and for the rest of the solutions, the $\hat{y}$ coordinate of their intersections with the Poincaré section defined as the half-hyperplane

[TABLE]

where $\hat{t}_{{\cal P}_{y}}=(0,1,0,\ldots)$ and $\hat{n}_{{\cal P}_{y}}=(1,0,0,\ldots)$. The Poincaré section (60) is equivalent to marking the reduced state space trajectories’ intersections with the $\hat{y}$-axis in $+\hat{x}$ direction.

At the period doubling bifurcation of ${\rm RPO{4os}}$, annotated with PD in figure 8, a period-8 cycle appears in its vicinity, which itself undergoes another period doubling upon further increase in Me. We postpone the analysis of the chaotic attractor, which forms following this cascade of period doublings to section III.2, and turn our focus to another set of solutions that coexist with the ones we presented so far.

III.1.2 The lemniscate solution and the symmetry breaking

In the previous section, we followed primary and subsequent bifurcations from the counterclockwise rotating circular solution. All the solutions we presented so far had positive angular momentum $L=xv_{y}-yv_{x}>0$ at all times. Since the system is equivariant under the reflection symmetry (53), all of these solutions have reflection-symmetric counterparts with $L<0$. A different type of stable periodic solution, which we will refer to as ${\rm PO{2l}}$ coexist in the state space at $\textit{Me}=14.0$. This solution can be found by initiating a simulation with a random initial condition at $\textit{Me}=14.0$: We observed at this regime that an initial condition $a(0)$ that is populated by standard normally distributed numbers eventually lands on ${\rm PO{2l}}$, or ${\rm RPO{4}}$ (or its reflection copy $\sigma{\rm RPO{4}}$). Both of these orbits are stable at $\textit{Me}=14.0$.

We visualized three snapshots from one full period of ${\rm PO{2l}}$ in figure 9. Perrard et al. PLMFC2014 called similar solutions “lemniscates”, referring to their figure-8 shape. From figure 9, it can be seen that the droplet’s angular momentum is reversed as it traverses through one full period. In fact, the angular momentum along ${\rm PO{2l}}$ is exactly zero as the second half of the orbit is the reflection symmetry of the first half. Similar solutions of the Kuramoto-Sivashinsky system are sometimes referred to as “pre-periodic” SCD07 ; BudCvi15 since they can be treated as relative periodic orbits satisfying

[TABLE]

In other words, after half period $T_{{\rm PO{}}}/2$, the periodic orbit comes to the reflection of its initial position. Such solutions generically do not undergo period doubling bifurcations  SwiWie84 and when the Jacobian of (61) has a negative real eigenvalue outside the unit circle, this corresponds to the symmetry-breaking.

A symmetry-breaking bifurcation of ${\rm PO{2l}}$ takes place at $\textit{Me}\in(15.6,15.7)$ and there appears drifting lemniscate solutions in its vicinity. We visualized these dynamics in figure 10 by showing the trajectories of small perturbations $\pm 10^{-2}V_{1}$ to ${\rm PO{2l}}$, where $V_{1}$ is the leading (unstable) Floquet vector. Figure 10 (a) shows the symmetry-broken dynamics in full state space along with the unstable ${\rm PO{2l}}$. In figure 10 (b), when the symmetry is reduced, it becomes apparent that two new relative periodic orbits appear in ${\rm PO{2l}}$’s vicinity. We are going to refer to these drifting lemniscates as ${\rm RPO{2l}}$ and $\sigma{\rm RPO{2l}}$.

At $\textit{Me}\in(16.5,16.6)$, ${\rm RPO{2l}}$ undergoes a period doubling bifurcation, which leads to a cascade and a chaotic attractor formation. Before examining these dynamics in section III.2, we must present another type of solution, which coexist with the ones we studied so far.

III.1.3 The trefoil

In the previous section, we have stated that at $\textit{Me}=14.0$, the pilot-wave system has two stable solutions, ${\rm RPO{1o}}$ and ${\rm PO{2l}}$, and a random initial condition eventually lands on one or the other. A natural question to ask at this stage is the following. What separates the basins of attraction of these two solutions? Experience from previous studies GOY1982 ; GOY1983 ; ASY1997ch10 of low dimensional systems suggests that such disconnected regions might be separated by the stable manifold of another invariant solution. One of the most prominent examples of this kind of state space structure is manifested in the subcritical transition to turbulence in shear flows, where laminar and turbulent flows are coexisting attractors. In numerical simulations of channel flow, Itano and Toh IT01 were first to find a traveling wave solution (a relative equilibrium) by a shooting method which bisects between trajectories whose energy fluctuations eventually reach turbulent levels and those that vanishes. This suggested that the turbulent and laminar regions of the state space were separated by the stable manifold of this solution. Following numerical studies TI03 ; SchEckYor07 ; SGLDE08 ; MMSE09 ; ScMaEc10 ; ZamEck14a ; KKSDEH2016 of shear flows found similar solutions333The exact type of the invariant solution depends on the system under study and the parameters. with shooting methods akin to the one used by Itano and Toh  IT01 . Such separating solutions in subcritical transition literature are now generally referred to as “edge states”, following the terminology introduced by Schneider et al. SchEckYor07 and the bisection methods for locating them is usually called “edge tracking”. We are going to adopt this terminology in what follows.

In order to investigate the edge between the initial conditions that eventually land on ${\rm RPO{1o}}$ and those that go to ${\rm PO{2l}}$, we formulated the following edge tracking method. For tracking the laminar-turbulent boundary in shear flows, the kinetic energy of turbulent fluctuations is usually used as a proxy to determine whether the solution eventually becomes turbulent. In our case, we found the time-averaged angular momentum of the final state to be an appropriate quantity and formualted the edge tracking algorithm in appendix B.

Let $\hat{a}_{{\rm PO{2l}}}$ and $\hat{a}_{{\rm RPO{1o}}}$ be random points on ${\rm RPO{1o}}$ and ${\rm PO{2l}}$. In order to search for the basin boundary between ${\rm PO{2l}}$ and ${\rm RPO{1o}}$, we run edge tracking algorithm (appendix B) with $\epsilon=10^{-9}\,,\,\tau_{f}=10000\,,\,\tau_{w}=1000,L_{th,1}=1.8\,,\,L_{th,2}=0.1$. We plotted the time-series of the window-averaged angular momentum in figure 11. As expected, the average angular momentum eventually settles around $\langle L\rangle\approx 0$ (${\rm PO{2l}}$) or $\langle L\rangle\approx 2.1$ (${\rm RPO{1o}}$). Notice, however, that before landing on these final states, the average angular momentum of the bisection trajectories oscillates around $\langle L\rangle\approx 1.65$. We observed approximately periodic symmetry-reduced dynamics in these episodes; and upon a Newton search, found a new relative periodic orbit, which we will refer to as ${\rm RPO{1t}}$.

In Figure 12, we visualize three snapshots from ${\rm RPO{1t}}$’s trajectory simulated in the full state space. Experimentally observed solutions with similar shapes were named “trefoils” by Perrard et al. PLMFC2014 . By numerical continuation, we confirmed that when $\textit{Me}\in[12.91,19.0]$, ${\rm RPO{1t}}$ has one positive Floquet multiplier outside the unit circle, which renders ${\rm RPO{1t}}$’s stable manifold codimension-1. This is in agreement with our initial intuition since a codimension-1 manifold divides the state space into two distinct regions. Thus, we conclude that at $\textit{Me}=14.0$, the edge between the basins of attraction of ${\rm RPO{1o}}$ and ${\rm PO{2l}}$ is set by the stable manifold of the trefoil solution ${\rm RPO{1t}}$.

${\rm RPO{1t}}$ appears in a saddle-node bifurcation at $\textit{Me}\approx 12.91$ along with a stable branch, trajectory of which is also trefoil-shaped. We observed this stable solution for $\textit{Me}\in[12.91,12.95]$, however, our efforts in continuing this solution to higher Me have failed. Trefoil-shaped trajectories have been observed in both laboratory PLMFC2014 and numerical KOB2017 experiments, which could only be possible if the associated solution is stable. Our observation of the stable trefoil orbit following the saddle-node bifurcation provides one explanation for experimental findings by fine-tuning the system parameters and initial conditions. It is conceivable that other trefoil-shaped solutions might appear at saddle-node bifurcations at different parameter values.

III.2 Chaotic attractors and their global bifurcations

In the previous sections, we showed that the droplet system has three different families of solutions: (1) the circular solution and the subsequent ovals, (2) the lemniscate solutions, and (3) the trefoil. We have already mentioned that two of these solutions, namely ${\rm RPO{4os}}$ and ${\rm RPO{2l}}$, undergo period-doubling bifurcations. In this section, we are going to illustrate that both of these period-doublings yield chaotic dynamics, which subsequently merge and generate full complexity of the chaotic pilot-wave dynamics.

The first set of period doublings take place in the neighborhoods of ${\rm RPO{2l}}$ and its reflection symmetry $\sigma{\rm RPO{2l}}$. In order to illustrate these bifurcations, we generated the orbit diagram in figure 13, where we plotted the $\hat{y}$ coordinate of the very long trajectories’ intersections with the Poincaré section (60). Starting at $\textit{Me}=16.0$ with initial conditions on ${\rm RPO{2l}}$ and $\sigma{\rm RPO{2l}}$, we increased Me in steps of $0.001$, used the final state of the previous step as the initial condition of the next step, and generated a long data set. Figure 13 shows the intersections of these trajectories with the Poincaré section (60) after initial transients are discarded. Green markers correspond to orbits that started from ${\rm RPO{2l}}$ and the blue ones correspond to the ones that started from $\sigma{\rm RPO{2l}}$. Additionally, the symmetric lemniscate ${\rm PO{2l}}$ is shown as black dashed lines, and the trefoil solutions ${\rm RPO{1t}}$ and $\sigma{\rm RPO{1t}}$ are drawn as yellow dashed lines.

The first period doubling of ${\rm RPO{2l}}$ and $\sigma{\rm RPO{2l}}$ takes place at $\textit{Me}\approx 16.437$. This bifurcation is followed by a cascade of period doublings, which yields chaotic sets ${\rm C{l}}$ and $\sigma{\rm C{l}}$ in the vicinity of ${\rm RPO{2l}}$ and $\sigma{\rm RPO{2l}}$, respectively. As illustrated in figure 10, these two reflection-related regions are separated by ${\rm PO{2l}}$ (drawn dashed black in figure 13) and we observed that the two chaotic regions merge at $\textit{Me}\approx 16.580$, as they cross the border set by ${\rm PO{2l}}$. At this point, the chaotic attractors, which are symmetry-copies of one another merge through a symmetry-increasing bifurcation ChoGol1988 and result in a reflection-symmetric attractor ${\rm C{l}}=\sigma{\rm C{l}}$.

Upon further increase of Me, we observed that when $\textit{Me}>16.616$, the long trajectories eventually leave this chaotic neighborhood and settle on ${\rm RPO{1o}}$ or $\sigma{\rm RPO{1o}}$, which are stable until $\textit{Me}\approx 16.9$, see figure 5. This suggests that the chaotic attractor crosses the stable manifold of ${\rm RPO{1t}}$, which sets the boundary between the neighborhoods of ${\rm RPO{2l}}$ and ${\rm RPO{1o}}$. At the first glance, approach of the chaotic attractor towards ${\rm RPO{1t}}$ is not obvious, however, when we zoomed into the parameter region $\textit{Me}\in[16.6,16.616]$ (the inset of figure 13), we observed that the chaotic dynamics came close to ${\rm RPO{1t}}$. These visits are rather rare: We observed only a few instances in more than $3500$ intersections with the Poincaré section (60) at $\textit{Me}=16.616$. The sudden transition from a chaotic attractor to a chaotic transient is called a “boundary crisis” GOY1982 ; ASY1997ch10 . It is not clear from our data whether or not the expansion of the attractor towards ${\rm RPO{1t}}$ is continuous. Nevertheless, we believe that around $\textit{Me}\in(16.616,16.617)$, a boundary crisis causes the destruction of this chaotic attractor. For further evidence, we carried out a bisection search similar to that of figure 11 at $\textit{Me}=16.6$ and found, once again, that ${\rm RPO{1t}}$ to be the edge state in between the neighborhoods of ${\rm RPO{2l}}$ and ${\rm RPO{1o}}$.

In section III.1.1, we have shown that upon the increase of Me, ${\rm RPO{1o}}$ first undergoes a Neimark-Sacker bifurcation and as we further increased the Me, the invariant torus in the neighborhood of ${\rm RPO{1o}}$ is replaced the stable and unstable manifolds of a pair of period-4 orbits ${\rm RPO{4os}}$ and ${\rm RPO{4ou}}$. We ended section III.1.1 by mentioning that ${\rm RPO{4os}}$, which is initially stable, undergoes a period doubling bifurcation at $\textit{Me}\in(17.8,17.9)$. We illustrated the following cascade of period-doublings and the emergence of chaotic dynamics on the orbit diagrams of figure 14

Similar to figure 13, we produced figure 14 by recording the intersections of long trajectories with the Poincaré section (60) after discarding the initial transients for $\textit{Me}\in[17.8,19.0]$. At $\textit{Me}=17.8$, we started with an initial condition on ${\rm RPO{4os}}$ and for the following simulations, we used the final point of the last step as the initial condition. Figure 14 illustrates the results of these simulations, where we show the $\hat{y}$ values where trajectories intersect (60) at different Me. For clarity, we did not show the bifurcations of $\sigma{\rm RPO{4os}}$ on figure 14, but one should keep in mind that the same cascade also takes place in the vicinity of $\sigma{\rm RPO{4os}}$.

Three consecutive period doublings are visible in figure 14(a) and more can be found by zooming into smaller parameter regions. This result shows that the transition to chaos in this system follows a typical period-doubling route as opposed to Ruelle-Takens-Newhouse scenario nhouse78 , which was previously suggested by Tambasco et al. THORB2016 . We are going to refer to the resulting chaotic attractors as ${\rm C{o}}$ and $\sigma{\rm C{o}}$.

As visualized in figure 14(b), dynamics of ${\rm C{o}}$ for $\textit{Me}>18.2$ is chaotic with intermittent periodic windows, which is typical for attractors that follow a period doubling cascades strogb . Notice the sudden expansion of the chaotic attractor at $\textit{Me}\approx 18.455$. This is another symmetry-increasing bifurcation where ${\rm C{o}}$ connects to its reflection symmetry $\sigma{\rm C{o}}$ yielding the final chaotic attractor of the pilot-wave system. As we shall now illustrate, this global bifurcation is indeed the result of the merging of ${\rm C{o}}$ and $\sigma{\rm C{o}}$ with ${\rm C{l}}$.

As we explained earlier, at $\textit{Me}\approx 16.617$, ${\rm C{l}}$ becomes unstable by merging with ${\rm RPO{1t}}$. After this boundary crisis, ${\rm C{l}}$ includes ${\rm RPO{1t}}$ and becomes the chaotic edge state between the attractors ${\rm RPO{1o}}$ and $\sigma{\rm RPO{1o}}$. ${\rm C{l}}$ carries this role up until $\textit{Me}\approx 18.455$, when the distinct attractors are chaotic. We illustrated this in figure 15 at $\textit{Me}=16.8$ and $\textit{Me}=18.45$ where the red trajectories were obtained by edge tracking (appendix B) between solutions with positive and negative angular momenta. In both cases, red trajectories in figure 15 are chaotic transients that would fall onto one of the attractors if integrated long enough. Figure 15 (b) illustrates the state space picture just before the merging of two chaotic regions. At this point, blue and red regions are still dynamically disconnected and their apparent overlap is due to the lower-dimensional projection.

Finally, we are going to present the chaotic attractor following the symmetry-increasing global bifurcation at $\textit{Me}\approx 18.455$. After this point, two chaotic regions ${\rm C{o}}$ and $\sigma{\rm C{o}}$ are connected through ${\rm C{l}}$ and the trajectories intermittently switch between these neighborhoods. We illustrate this dynamics in figure 16 where we plotted a long chaotic trajectory along with ${\rm REQ{0c}}$, ${\rm RPO{1o}}$, ${\rm RPO{1t}}$, ${\rm PO{2l}}$, ${\rm RPO{2l}}$ in the full state space (figure 16 (a)) and in the symmetry-reduced state space (figure 16 (b)) at $\textit{Me}=18.5$.

The chaotic attractor in figure 16 (b), which is revealed after the symmetry-reduction, is qualitatively similar to the well-known Lorenz attractor  lorenz63 . Akin to the counterclockwise-rotating circular solution ${\rm REQ{0c}}$ and its reflection-symmetry $\sigma{\rm REQ{0c}}$, the Lorenz system at standard parameter values has two non-trivial equilibria that are related by a rotation by $\pi$ and the dynamics chaotically switches between their neighborhoods. In the droplet system, the switchings between the neighborhoods of the circular solutions correspond to the reversal of the droplet’s angular momentum. We illustrated this in the time-series of figure 16(c), where we show the chaotically moving droplet’s angular momentum and its distance from ${\rm REQ{0c}}$ as functions of time. Since ${\rm REQ{0c}}$’s angular momentum is positive (counterclockwise rotation), when the trajectory is close to it, the droplet’s angular momentum is positive and vice versa. This, in fact, is also visible from the projection of figure 16(b). Since the symmetry reduction (9) sets $\hat{v}_{y}=0$ and $\hat{v}_{x}>0$, the angular momentum in the symmetry-reduced state space is $L=-\hat{y}\hat{v}_{x}$ and its sign is the opposite of $\hat{y}$’s. Hence, the $\hat{y}=0$ hyperplane separates solutions with positive angular momenta from those with negative angular momenta in the symmetry-reduced state space.

Besides the chaotic attractor, we also visualized the invariant solutions of the pilot-wave system in figure 16 (a,b). It is clear from figure 16(b) that the oval (${\rm RPO{1o}}$) and the trefoil (${\rm RPO{1t}}$) solutions’ angular momenta does not change sign along one period whereas the angular momenta of the lemniscates reverse twice in one full period. The chaotic trajectories’ angular momentum reversals correspond to the episodes during which they are shadowed by the lemniscates.

IV Conclusion and discussion

In this paper, we introduced a continuous symmetry-reduction method for hydrodynamic pilot-wave systems with rotational symmetry. In essence, our method is fixing the polar-angle in the velocity plane $(v_{x},v_{y})$ (as opposed to $(x,y)$-plane) in order to obtain a transformation that is nonsingular for generic dynamics. We formulated this transformation in such a way that it can be used to simplify both experimental and numerical data by eliminating the symmetry degeneracy of the solutions. We then proceeded to reformulate this transformation in the framework of the method of slices. This formulation brought us a set of geometrical tools, which are useful mainly for theoretical and computational undertakings. In order to demonstrate the utility of the introduced symmetry-reduction technique, we applied it to the numerical study of Oza et al.’s  ORB2013 trajectory equation with central harmonic forcing. With the help of the symmetry-reduction, we analyzed the bifurcations and the subsequent chaotic dynamics in this system.

The main message we would like to deliver is that the continuous symmetry reduction substantially simplifies the pilot-wave dynamics by eliminating infinitely many copies of each generic solution. This is perhaps best illustrated in figure 16, where the symmetry-reduction reveals a relatively simple chaotic attractor with qualitatively different regions.

Moreover, the symmetry reduction allowed us to numerically identify the unstable invariant solutions and their roles in organizing the state space of the system.

In the previous experimental PLMFC2014 ; PLFC2014 and numerical KOB2017 studies of hydrodynamic pilot-wave systems with central harmonic forcing, invariant solutions such as circles, ovals, lemniscates, and trefoils, were observed by tuning system parameters to the values at which these solutions were stable. These papers also demonstrated that the chaotic solutions transiently resembled the invariant solutions, hence they conjectured that in the chaotic regime, the dynamics can be decomposed into episodes approximated by different invariant solutions. Our results essentially confirm this conjecture: The chaotic pilot-wave dynamics is a union of two chaotic state space regions that are formed in the vicinity of the oval and lemniscate solutions and these solutions persist in the chaotic regime, albeit they are unstable. This observation can potentially be brought to a quantitative level: Long-term averages of the observables associated with a chaotic system can be predicted from its periodic orbits if the system and periodic orbits satisfy certain properties such as ergodicity and hyperbolicity. The collection of methods to carry out such computations is known as the periodic orbit theory DasBuch ; Cvitanovic1992 and our results suggest that such a calculation can be done in this context. Hence, the periodic orbit analysis of chaotic pilot-wave systems is one of our future research directions.

Symmetries of the radially-confined pilot-wave systems were recognized in previous literature at varying levels. For example, Perrard et al. PLMFC2014 noted that the symmetries of the droplet trajectories were reflected in the Fourier-Bessel decomposition of the bath’s surface height. Durey & Milewski DurMil2017 observed that the perfect lemniscate solution had [math] mean angular momentum due to its reflection symmetry. Durey & Milewski DurMil2017 also used successive maxima of radius as the averaging window in their data analysis. Similarly, Perrard et al. PLFC2014 used successive maxima of the radius in order to construct a return map. The choice of radius as a diagnostic variable is inherently symmetry-related since the radius is an invariant of rotation and reflection. The main difference of our work against the previous techniques is that we do not only use symmetry-invariant variables for extracting information but also express the dynamics in terms of them. Thus, our symmetry-reduced representation retains all dynamical information.

Transition to chaos in the model studied here were also examined by Tambasco et al. THORB2016 who conjectured that the chaotic dynamics appear following a Ruelle-Takens-Newhouse transition scenario nhouse78 . Their conjecture was based on the observation of a three-frequency orbit (fig. 6d of ref. THORB2016 ) prior to chaos. However, as noted by the authors themselves, two frequencies $f_{1}$ and $f_{2}$ that are visible in the spectrum of this orbit are in fact commensurate, i.e. $f_{1}=4f_{2}$. Thus, even if the third frequency is incommensurate with the first two, the associated symmetry-invariant dynamics takes place on a two-torus which does not meet the necessary criterion for a Ruelle-Takens-Newhouse transition scenario nhouse78 . Here, we found that the transition from the stable period-4 orbit ${\rm RPO{4os}}$ to chaos is a result of a period-doubling cascade, through which a new commensurate frequency appears at each bifurcation. We would like to note that the numerical method OWHRB2014 adopted in ref. THORB2016 is very different from ours: They integrate the integro-differential trajectory equation (27) with (29) using an explicit fourth-order Adams-Bashfort method with a fixed time-step. We, on the other hand, simulate a truncated ODE representation. In order to check our results’ robustness against truncation at this regime, we quadrupled number of Fourier-Bessel modes and confirmed that we still observe period-doubled orbits prior to the chaos.

While we focused on pilot-wave systems in this paper, our symmetry-related methods can be extended to a larger class of systems in which a finite-size (or point-like) object interacts with a continuous environment. Similarly, it could be possible to generalize our method to three-dimensional space by introducing Euler angles for polar coordinate transformations.

In summary, we have shown that the continuous symmetry reduction drastically simplifies the data generated by a hydrodynamic pilot-wave system with chaotic dynamics. To this end, we formulated a symmetry-reduction method and applied it to the numerical simulations of a rotation-equivariant pilot-wave model with a central harmonic force. We identified local bifurcations, the onset of chaos, and the global attractor merging bifurcations in the system; as well as a qualitative description of its state space geometry.

Acknowledgement

We acknowledge stimulating discussions with Predrag Cvitanović and Matthieu Labousse. We are grateful to George Choueiri for his critical reading of an early version of this manuscript.

Appendix A Newton’s method

Let $(\hat{a}_{{\rm RPO{}}}^{(i)},T_{{\rm RPO{}}}^{(i)}$ be the guesses for an initial condition on a (relative) periodic orbit and its period, which approximately satisfy (25). In order to bring them to a desired numerical precision, we utilize Newton’s method as follows. Our goal is to find the small corrections $(\Delta\hat{a},\Delta T)$ which yield

[TABLE]

Expanding (62) to the linear order in $\Delta\hat{a}$ and $\Delta T$, we obtain

[TABLE]

Defining $\hat{J}^{(i)}\equiv\hat{J}^{T_{{\rm RPO{}}}^{(i)}}(\hat{a}_{{\rm RPO{}}}^{(i)}$), $\hat{v}^{(i)}\equiv\hat{v}({\hat{f}^{T_{{\rm RPO{}}}^{(i)}}(\hat{a}_{{\rm RPO{}}}^{(i)})})$, and $E^{(i)}={\hat{f}^{T_{{\rm RPO{}}}^{(i)}}(\hat{a}_{{\rm RPO{}}}^{(i)})}-\hat{a}_{{\rm RPO{}}}^{(i)}$ , and rearranging the terms in (63), we obtain the linear system of equations

[TABLE]

for the approximate corrections $(\Delta\hat{a}^{(i)},\Delta T^{(i)}))$. (64) provides $d$ equations for $d+1$ unknowns, where $d$ is the dimension of the symmetry-reduced state space. We supplement (64) with the constraint

[TABLE]

which disallows for updating our guesses in the flow direction. We can now solve (64) and (65) for $(\Delta\hat{a}^{i},\Delta T^{i})$ and update our guesses as

[TABLE]

where we introduced the step size $\alpha^{k}\in(0,1]$. We set $\alpha=0.5$ and choose the smallest $k=0,1,2\ldots$ such that $||E_{i+1}||<||E_{i}||$. We iterate this procedure until the error $||E_{i+1}||$ is less than a tolerance $\epsilon$.

At each Newton step, we sample the relative periodic orbit with $\lfloor T_{{\rm RPO{}}}^{i}/\delta\tau\rfloor$ steps in time, hence the final converged time-step $\delta\tau_{{\rm RPO{}}}$ of the periodic orbit is close to the rest of the simulations.

Appendix B Edge-tracking algorithm

Let $\langle L\rangle(\tau)$ denote the moving average of the angular momentum with a time-window $\tau_{w}$ and $\hat{a}_{1}$ and $\hat{a}_{2}$ be two state vectors that belong to attractors $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, respectively. We assume that $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ have mean angular momenta $L_{1}$ and $L_{2}$ respectively such that $L_{1}>L_{th,1}$ and $L_{2}<L_{th,2}$ for some thresholds $L_{th,1}$ and $L_{th,2}$. We iterate the following algorithm in order to locate the edge that separates $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$’s basins of attraction.

while $||\hat{a}_{1}-\hat{a}_{2}||>\epsilon$ do

    $\hat{a}_{3}\leftarrow(\hat{a}_{0}+\hat{a}_{1})/2$

    Set initial condition to $\hat{a}_{3}$

    Compute $\hat{a}(\tau)$ and $\langle L\rangle(\tau)$ for $\tau\in[0,\tau_{f}]$

    if $\langle L\rangle(\tau_{f})>L_{1}$ then

         $\hat{a}_{1}\leftarrow\hat{a}_{3}$

    else if $\langle L\rangle(\tau_{f})<L_{2}$ then

         $\hat{a}_{2}\leftarrow\hat{a}_{3}$

    else

         Terminate

    end if

end while