PHYSICAL REVIEW E 108, 034213 (2023)

Pilot-wave dynamics: Using dynamic mode decomposition to characterize bifurcations, routes tochaos, and emergent statistics
J. Nathan Kutz ,1 André Nachbin ,2 Peter J. Baddoo ,3 and John W. M. Bush
1

3

Department of Applied Mathematics and Electrical and Computer Engineering, University of Washington, Seattle, Washington 98195, USA
2
Department of Mathematical Sciences, Worcester Polytechnic Institute, Worcester, Massachusetts 01742, USA
3
Department of Mathematics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA(Received 14 October 2022; revised 7 August 2023; accepted 29 August 2023; published 21 September 2023)
We develop a data-driven characterization of the pilot-wave hydrodynamic system in which a bouncing droplet
self-propels along the surface of a vibrating bath. We consider drop motion in a confined one-dimensional
geometry and apply the dynamic mode decomposition (DMD) in order to characterize the evolution of the
wave field as the bath’s vibrational acceleration is increased progressively. Dynamic mode decomposition
provides a regression framework for adaptively learning a best-fit linear dynamics model over snapshots ofspatiotemporal data. Thus, DMD reduces the complex nonlinear interactions between pilot waves and droplet toa low-dimensional linear superposition of DMD modes characterizing the wave field. In particular, it provides
a low-dimensional characterization of the bifurcation structure of the pilot-wave physics, wherein the excitation
and recruitment of additional modes in the linear superposition models the bifurcation sequence. This DMD
characterization yields a fresh perspective on the bouncing-droplet problem that forges valuable new links withthe mathematical machinery of quantum mechanics. Specifically, the analysis shows that as the vibrational
acceleration is increased, the pilot-wave field undergoes a series of Hopf bifurcations that ultimately lead to achaotic wave field. The established relation between the mean pilot-wave field and the droplet statistics allows us
to characterize the evolution of the emergent statistics with increased vibrational forcing from the evolution of thepilot-wave field. We thus develop a numerical framework with the same basic structure as quantum mechanics,
specifically a wave theory that predicts particle statistics.
DOI: 10.1103/PhysRevE.108.034213

I. INTRODUCTION

Pilot-wave theories have been developed to describe
particle-wave interactions on both the microscopic and macroscopic scales. In quantum mechanics, the double-solution
pilot-wave theory of de Broglie [1,2] was proposed in the
1920s on the premise that microscopic particles have an associated internal vibration at the Compton frequency thatgenerates waves. The resonant wave-particle interaction was
posited to result in the particle being propelled by its guiding
or pilot wave with the de Broglie wavelength, giving rise tostatistical behavior consistent with the standard formalism.
Hydrodynamic pilot-wave theory [3,4] has been developed todescribe the motion of millimetric droplets self-propelling on
the surface of a vibrating liquid bath, a system discovered
by Couder and Fort in 2005 [5,6]. By virtue of a resonantinteraction between the bouncing droplets and the underlying
wave field, the droplet is piloted by a quasimonochromaticwave field with the Faraday wavelength. The periodic bouncing of the droplet thus plays the role of the particle vibrationin de Broglie’s mechanics and the Faraday wavelength that ofthe de Broglie wavelength [3]. This walking-droplet system
is remarkable in that it represents a macroscopic realizationof the type of pilot-wave dynamics envisaged by de Broglie
and all the more remarkable in that it has yielded a growing
list of hydrodynamic quantum analogs [3,4,7]. These include
analogs of single-particle diffraction and interference [6,8,9],
2470-0045/2023/108(3)/034213(11)orbital quantization [10–14], tunneling [15–17], quantum
corrals [18,19], Friedel oscillations [20], spin lattices [21]and long-range correlations in bipartite systems [22–24]. In
several of these systems, wavelike statistical forms emergethat are strikingly similar to those arising in their quantum
counterparts [18–20]. Hydrodynamic quantum analogs thus
suggest that a pilot-wave dynamics of the form engendered in
the hydrodynamic system might plausibly underlie quantum
statistics [3,4].
The walking-droplet system is complex and strongly nonlinear, being governed by a wave equation for the fluid bathand a trajectory equation for the droplet. The waves are governed by a system of partial differential equations with anonlocal Fourier operator and a nonlinear forcing induced by
the walking droplets. The droplet bounces vertically, and sointeracts intermittently with the bath surface, and is propelled
by a force proportional to the local slope of the interface
during impact. Despite this dynamical complexity and nonlinearity, we proceed by showing that by leveraging the dynamic
mode decomposition (DMD), the mean pilot-wave field anddroplet statistics can be well characterized by a low-ranklinear dynamics.
Specifically, the DMD model shows that linear superposition allows for a reconstruction of the mean wave field, fromwhich we recover the probability density function of the particle. In quantum mechanics, the particle statistics is governedby a linear Schrödinger equation. Here DMD extracts from a034213-1©2023 American Physical Society

KUTZ, NACHBIN, BADDOO, AND BUSH

PHYSICAL REVIEW E 108, 034213 (2023)

nonlinear system of coupled dynamical equations the linear
features of the particle statistics. Dynamic mode decomposition also allows for a characterization of the corresponding
bifurcation diagram, specifically the period-doubling route tochaos arising as the vibrational forcing is increased progressively. Our analysis of the walking-droplet system focuseson the dynamics of the waves rather than the particle andprovides a simple spatiotemporal modal analysis of theunderlying low-dimensional wave interactions that characterize the overall pilot-wave physics and emergent particle
statistics.
The feature of the walking-droplet system responsible
for the emergent quantum features is path memory [25],which renders the drop dynamics non-Markovian and results from the persistence of the pilot wave on the bathsurface. The drop is propelled by its wave field, whose form
is prescribed by the droplet’s history and environment. The
critical control parameter is the bath’s vibrational acceleration , which prescribes the bath’s proximity to the Faraday
threshold F , above which waves would form even in theabsence of the droplet. While  is always less than F in
the laboratory, the closer  is to F , the more persistent
the waves generated by the droplet, and the longer the path
memory. The manner in which the droplet dynamics changes
as the vibrational acceleration (or memory) is increased
progressively towards the Faraday threshold has been characterized and rationalized for the free droplet [26–28], droplet
pairs [29–31], confined rings [32,33], and free rings [34].
Transitions from steady orbital motion to chaotic dynamics
have been reported and characterized in a number of settings,
including motion in a rotating frame [10–12,35], motion confined by a central force [13,14,36], and motion confined byboundaries [18,19,37]. Rich dynamical properties have been
revealed, including period-doubling cascades to chaotic trajectories [11,35,38,39]. The dynamical system aspects of thewalking-droplet system have been highlighted in the recent
review by Rahman and Blackmore [40].
In characterizing the bifurcation structure of the walkingdroplet system, prior work has focused primarily on thedroplet dynamics, specifically how the particle trajectorychanges with increasing memory. However, for confined
walker motion, connections between the mean pilot-wave
field and the emergent statistical forms have been both reported [19] and rationalized [37,41]. As the resulting mean
pilot wave is expressible in terms of the droplet’s statistical behavior, Bush and Oza [4] proposed that it may playa role similar to that of the quantum potential in Bohmianmechanics. While the mean walker dynamics is thus effectively nonlocal, as is Bohmian mechanics, the time-resolved
walker dynamics is entirely local. The present work focuseson the bifurcation structure of the pilot wave as the systemmemory is increased progressively. By leveraging the DMD,
the underlying pilot-wave field can be characterized, revealing its own period-doubling cascade to chaos. The analysis
provides a spatiotemporal modal analysis of the underlying
low-dimensional wave interactions that characterize the overall pilot-wave physics.
Dynamic mode decomposition originated as a modal analysis method in the fluid dynamics community. Introduced
as an algorithm by Schmid and Sesterhenn [42,43], it has

rapidly become a commonly used data-driven analysis tooland the standard algorithm to approximate the Koopman
operator from data [44]. Specifically, DMD was used to identify dominant spatiotemporal coherent fluid structures fromhigh-dimensional time-series data. The DMD analysis offeredan alternative to standard dimensionality reduction methods
such as the proper orthogonal decomposition (POD), whichhighlighted low-rank features in fluid flows using the computationally efficient singular value decomposition (SVD) [45].
The advantage of using DMD over SVD is that the DMD
modes are linear combinations of the SVD modes that have acommon linear (exponential) behavior in time, given by oscillations at a fixed frequency with growth or decay. Specifically,
optimized DMD [46,47] is a regression to solutions of theformx(t ) =rb j φ j eω j t ,

(1)

j=1

where x(t ) is a rank-r approximation to a collection of state
space measurements xk = x(tk ) (k = 1, 2, . . . , n). The algorithm finds the DMD eigenvalues ω j , DMD modes φ j , andtheir loadings b j . The ω j determines the temporal behaviorof the system associated with a modal structure φ j , thus
giving a highly interpretable representation of the dynamics. Such a regression can also be learned from time-series
data [48]. Dynamic mode decomposition may be thought ofas a combination of the SVD and POD in space with theFourier transform in time, combining the strengths of eachapproach [49,50].
Dynamic mode decomposition also provides a natural
mathematical connection with pilot-wave hydrodynamics.
Specifically, the emergent steady-state statistics of the walking
droplets have been shown in certain instances to be comparable to solutions of the time-independent Schrödinger equation of quantum mechanics [4,18,19]. The time-dependent
Schrödinger equation takes the general formiψt = −ψxx + V (x)ψ,(2)

where ψ (x, t ) is the wave function, V (x) is a potential, andthe coefficients have been normalized to unity. The solutionof the Schrödinger equation is given by a linear superposition of modes, each of which has a characteristic oscillation
frequency. In quantum mechanics, these frequencies are theenergy levels of the quantum system. Although the interpretation of these solutions is quite different from that of thefluidic system under consideration, the representation of thesystems in terms of a linear superposition of modes oscillating
at characteristic frequencies is the same. Thus, despite thecomplexity of the fluidic system, its approximate DMD solution generates the canonical Sturm-Liouville solution form.
Specifically, the state-space variable x(t ) (1) is the vectorizedversion of the Schrödinger wave field ψ (x, t ). The generalsolution to the discretized Schrödinger equation is given by
the DMD approximation (1) where the modes and frequencies
are determined by regression in order to model the complex
nonlinear and nonlocal interactions between wave field anddroplet. In this case, the DMD modes and eigenvalues correspond to the eigenstates and their corresponding energy levels,
respectively, of the quantum system. Thus, DMD provides not

034213-2

PILOT-WAVE DYNAMICS: USING DYNAMIC MODE …

PHYSICAL REVIEW E 108, 034213 (2023)

only a valuable approximation method but a regression to thesolution form that is standard in quantum mechanics.
II. PILOT-WAVE PHYSICS

The walking-droplet system is modeled theoreticallythrough a trajectory equation for the droplet and a partial differential equation that describes the spatiotemporal dynamics
of the accompanying waves [26–28]. The wave model may
be obtained from an asymptotic simplification of the linear
free-surface Navier-Stokes equations [51]. The Helmholtz decomposition of the velocity field yields, in the weakly viscous
regime, an irrotational velocity potential perturbed by a weak
shearing component, arising through the stream function. The
pilot wave along the undisturbed free surface z = 0 is expressed through a weakly diffusive Bernoulli law as well as aweakly diffusive kinematic condition, as detailed below. The
velocity potential is denoted by φ(x, z, t ), where we note thatin the formulation that follows, the dynamics are constrained
to the x-z plane. The velocity field components in the bulkof the fluid are given by (u, v) = ∇φ. The wave elevation isdenoted by η(x, t ). The fluid parameters are the fluid density
ρ, the surface tension σ , and the kinematic viscosity ν. In thereference frame of the fluid bath, which oscillates at frequency
ω0 , gravity takes the form g(t ) = g[1 +  sin(ω0t )]. The freesurface wave equations are given by [16,51]∂φ
σ1(x, 0, t ) = −g(t )η + ηxx + 2ν φxx − Pd [x − X (t )],
∂t
ρ
ρ
(3)
∂η
(x, t ) = ODtN [φ] + 2νηxx .
(4)∂t
The diffusive terms, in both the Bernoulli law and the kinematic condition, are the leading-order terms from the vorticalcomponent of the Helmholtz decomposition. The presenceof the droplet is felt through the pressure term Pd , centeredat the droplet’s position X (t ), which acts as a wave maker.
This pressure term is compactly supported in space, overthe droplet’s diameter, and is discontinuous in time, being
activated periodically at each bounce. The velocity potential satisfies Laplace’s equation, which enables one to define
the Dirichlet-to-Neumann operator that maps the Dirichlet data φ(x, 0, t ) onto the free surface’s normal speed, attime t:ODtN [φ] = φz (x, 0, t ).

(5)

The operator ODtN is defined as a Fourier integral operator andis computed in a straightforward manner using a conformal
mapping and the fast Fourier transform [16].
To complete the wave-particle model, the above wave
system is coupled to the droplet’s horizontal trajectory equation [27]
∂η
dX
d 2X
= −F (t ) (X (t ), t ).
(6)
+ cF (t )
dt 2
dt
∂x
The drop is propelled by the wave force, which is proportional
to the local gradient of the wave field, and resisted by a linear
drag. The magnitude of the propulsive wave force transmitted during the contact time, here prescribed as Tc = TF /4,
where TF is the Faraday period, is denoted by F (t ). Thismtime-dependent coefficient also appears in the drag term, sincedrag is also imparted during impact. Further modeling details
can be found elsewhere [16,51].
The evolution of the PDE model is shown in Fig. 1. As
the forcing amplitude is increased, the walking-droplet system
undergoes a bifurcation sequence. For low-amplitude forcing
 = 4.8, the droplet bounces in place in the middle of the well[Fig. 1(a)]. As the forcing amplitude is increased, the droplet
begins to oscillate periodically in the well [Fig. 1(b)]. A further increase of  generates larger oscillations until eventuallyperiod doubling occurs, which is detailed in Sec. V. The dynamics of the droplet can also be plotted in phase space using
the droplet position, denoted by X , and the droplet speed,
denoted by V = Ẋ . Figure 2 shows the underlying attractors
associated with the dynamics illustrated in Fig. 1. Increasing
 leads to a chaotic motion of the droplet, as further detailed inFig. 3. The center of oscillation varies in a chaotic fashion. As
will be detailed in Sec. V, the bifurcation sequence exhibitsa period-doubling bifurcation sequence to chaos, which is
canonical for damped-driven systems [52–55].
The corresponding pilot-wave dynamics of Fig. 2 is highlighted in Fig. 4. As  takes on the values 4.8, 5.0, 5.3, and 5.4,
the pilot wave evolves from a steady field to a periodic wave
to a doubly periodic wave to spatiotemporal chaos. In what
follows we will illustrate that the spatiotemporal dynamics
can be completely characterized by the DMD. Indeed, DMD
exploits the low-rank structure of the wave field as shownin Fig. 5. Truncation of the SVD modes is typically donein a heuristic fashion. Often, a prescribed variance, such as99%, is sought for the truncation criterion. Here the truncatedmodes selected are prescribed by cross validation of the reconstruction error. Thus, modes aside from the dominant zerobackground mode are selected as complex conjugate pairs
until the reconstruction error is only marginally improved by
the inclusion of additional modes. Ultimately, the smallest
number of modes capable of accurately reconstructing thesolution is sought. Thus, in Fig. 5 the numbers of modes
used as  takes on the values 4.8, 5.0, 5.3, and 5.4 are one,
three, five, and nine, respectively. The addition of more modes
in any of these dynamical regimes has little impact on thereconstruction error.
III. DYNAMIC MODE DECOMPOSITION

The DMD algorithm can be best understood from theso-called exact DMD [56], which is simply a least-square
fitting procedure. Specifically, the exact DMD algorithm seeks
a best-fit linear operator A that approximately advances thestate of a system, x ∈ Rn , forward in time according to thelinear dynamical system
xk+1 = Axk ,(7)

where xk = x(k t ), with t denoting a fixed time step that issmall enough to resolve the highest frequencies in the dynamics. Thus, the operator A is an approximation of the Koopman
operator K restricted to a measurement subspace spanned bydirect measurements of the state x [44,57]. One can work
directly with the measured variables x or construct alternative observables via a number of different methods [57].
Time-delay embeddings, for instance, often offer an improved

034213-3

KUTZ, NACHBIN, BADDOO, AND BUSH

PHYSICAL REVIEW E 108, 034213 (2023)

FIG. 1. Wave-droplet dynamics over 150 Faraday periods in a cavity of width 1 cm. The vibrational forcing takes the values (a)  = 4.8,(b)  = 5.0, (c)  = 5.3, and (d)  = 5.4. For this geometry the Faraday threshold is approximately F = 6.15. All other parameters are kept
unchanged, including the droplet’s initial position at X = −0.35. The wave field η(x, t ) is denoted by blue while the horizontal droplet position
X (t ) is indicated by red.

characterization of nonlinear systems [58–61]. However, thefluidic wave-particle system considered here is well characterized by a linear DMD model from measurements of the

(a)

(b)

wave field alone. Thus there is no need to seek improved
observables for the system.
Bagheri [62] first highlighted that DMD is particularly
sensitive to the effects of noisy data, with systematic biases
introduced to the eigenvalue distribution [63–66]. The optimized DMD algorithm of Askham and Kutz [46], which uses
a variable projection method for nonlinear least squares toη(x, t)

(c)

X(t)

(d)

FIG. 2. Phase-space dynamics. The droplet position is denoted
by X and the droplet speed by V = Ẋ . Panels correspond to extensions of the sequences illustrated in Fig. 1, for which (a)  = 4.8,(b)  = 5.0, (c)  = 5.3, and (d)  = 5.4, and the initial particle
position X (0) = −0.35. The total time for the simulations is (a) and(b) t = 4000TF and (c) and (d) t = 8000TF .

FIG. 3. Simulation of particle-wave dynamics with  = 5.4, displayed over a longer time interval than in Fig. 1(d). The wave field
η(x, t ) is denoted by blue while the droplet position X (t ) is indicatedby red.

034213-4

PILOT-WAVE DYNAMICS: USING DYNAMIC MODE …

PHYSICAL REVIEW E 108, 034213 (2023)

Thus, the data matrix X may be reconstructed asX ≈  diag(b)T(ω)
⎡
⎤⎡
|
| b1
⎢
...
= ⎣φ1 · · · φr⎦⎣||⎤⎡eω1 t1⎥⎢ ..
⎦⎣ .
br eωr t1

···
...
···

⎤
eω1 tm.. ⎥.
. ⎦eωr tm(9)

The variable projection method for this approximation [46]determines the parameters of the matrix components. Thus the
optimization is formulated asarg minω,b X − b T(ω)F ,FIG. 4. Pilot-wave dynamics corresponding to Fig. 1, where takes on the values (a) 4.8, (b) 5.0, (c) 5.3, and (d) 5.4. The pilotwave field is shown for t ∈ [2000, 2080] after the droplet behavior
has settled into its long time dynamics, as illustrated in Figs. 1 and 2.

compute the DMD for unevenly timed samples, provides thebest performance of any algorithm currently available. This isbecause it directly solves the exponential fitting problem of
DMD. Consider the data matrix⎡⎤⎢
X = ⎣x(t1 )

x(t2 )···⎥x(tm )⎦.

(8)

where b =  diag(b). Thus optimized DMD provides a direct approximation to the solution (1), which is a regression
to the form of the modal solution expansion to the vectorized
Schrödinger equation. This has been shown to provide a superior decomposition due to its ability to optimally suppress
bias and handle snapshots collected at arbitrary times. Dynamic mode decomposition approximations do not constrain
the DMD eigenvalue to the imaginary axis. However, theoptimized DMD regression can be modified to fit eigenvalues on the imaginary axis or to constrain eigenvalues to theleft half plane for stability [46]. These constraints are notenforced in what follows. Moreover, with statistical bagging,the BOP-DMD algorithm (bagging, optimized DMD) also
provides robust models with uncertainty quantification in thepresence of noise, corruption, and outliers [47]. Physicallymotivated constraints can also be embedded into the DMD
architecture [67], further enhancing the algorithm. The disadvantage of optimized DMD is that one must solve a nonlinear
optimization problem, which can fail to converge.

The optimized DMD algorithm directly solves the exponential
fitting problem in order to produce the matrix decomposition.
100

(a)

50

100

(b)

50

0

0
0

10

100

20
(c)

50

0

10

100

20
(d)

50

0

0
0

10

20

0

10

20

FIG. 5. Low-rank structure of the pilot-wave dynamics corresponding to Fig. 4, where  takes on the values (a) 4.8, (b) 5.0,(c) 5.3, and (d) 5.4. The singular value decomposition shows thepercentage of variance in each SVD mode, i.e., the jth mode varianceis equal to 100σ j / σk . The variance of only the first 20 modes
is shown, with the magenta modes depicting the number of modes
required for accurate DMD reconstruction.

(10)

IV. INFERRING PARTICLE STATISTICS FROMWAVE DYNAMICS

Durey et al. [37] deduced a simple relation between themean pilot-wave field and the droplet statistics for walker
motion in a statistically steady state, either periodic or ergodic.
For a particle moving in a bounded region (for example,
when constrained by a central force), the mean wave fieldη may be expressed as the convolution of the wave field ofa stationary bouncer ηB and the droplet’s histogram μ(x),normalized to have total mass one. Note that in this instance,
where boundary effects are negligible, ηB is independent of
the droplet position but does depend on the system memory.
For a domain with variable bottom topography relevant for thecavity considered here, the result was generalized by Dureyet al. [41] through incorporation of the spatial dependenceof the bouncer wave field resulting from the influence ofboundaries. The convolution is thus generalized to an integral
operator, where the kernel ηG (x, y) is given by the wave fieldof a bouncer located at position y. In this case, the statisticsof the particle and its pilot wave are related by (see Eq. 3.2
in [41])
η(x) =∞−∞ηG (x, y)μ(y)dy,(11)

where η(x) = limN→∞ N1 Nk=1 η(x, tk ) is the time-averaged
wave field and tk = kTF . We note that numerical verification

034213-5

KUTZ, NACHBIN, BADDOO, AND BUSH

PHYSICAL REVIEW E 108, 034213 (2023)

FIG. 6. Pilot-wave dynamics for the lowest forcing value considered,  = 4.8. The simulation was run for 4000 Faraday periods(or droplet bounces). Figure 5 shows that a rank-1 decompositionis sufficient to model the steady-state wave dynamics. (a) Full PDE
evolution. (b) One-mode DMD approximation. (c) DMD eigenvalue.
(d) DMD eigenfunction used for reconstruction of the spatiotemporal
dynamics. Note that the DMD eigenvalue is at the origin, which isconsistent with the steady-state behavior in this parameter regime.
(e) Time-averaged pilot-wave field calculated directly (blue solidline) and deduced from the measured particle probability distribution
via (11) (red dashed line). (f) Particle probability distribution as
computed via particle tracking (blue solid line) and inferred fromthe mean pilot wave by inverting Eq. (11) (red dashed line).

FIG. 7. Pilot-wave dynamics for  = 5.0. The simulation was
run for 4000 Faraday periods. Figure 5 shows that a rank-3 decomposition is adequate to model the spatiotemporal periodic behavior
observed. (a) Full PDE evolution. (b) Three-mode DMD approximation. (c) DMD eigenvalues. (d) DMD eigenfunctions used forreconstruction of the spatiotemporal dynamics. (The black dottedline shows the steady-state background mode of Fig. 6.) Note thatthe DMD eigenvalues are manifest in complex conjugate pairings
typical of an underlying Hopf bifurcation. (e) Time-averaged pilotwave field (blue solid line) and that deduced from the measured
particle probability distribution via (11) (red dashed line). (f) Particle
probability distribution as computed via particle tracking (blue solidline) and inferred from the mean pilot wave by inverting Eq. (11) (reddashed line).

of (11) was not presented in [41]. In Figs. 6(e), 7(e), 8(e),and 9(e) we compare the time-averaged wave field, computed
directly from the numerical simulations, with that predictedby the expression (11). The agreement is excellent, even forthe cases where the phase-space orbits indicate complex cycles. To be explicit, the time-averaged wave field is the DMD
mode associated with the zero eigenvalue. Moreover, DMD
shows that the inferred particle statistics are expressible interms of a very few eigenmodes of the respective spectrum.
We proceed by exploring the relation (11) in a novel fashion. Inverting the integral operator allows us to infer particle
statistics from the mean wave field (as may be viewed as aneffective potential) computed over the time interval of interest.
The particle statistics is thus obtained with no specific knowledge of the particle dynamics, a step that evokes quantum
mechanics. Numerical aspects of the operator inversion arepresented in the Appendix. Figures 6(f), 7(f), 8(f), and 9(f)compare the particle’s probability density function, as deduced directly from particle tracking, with that produced byinverting equation (11). The agreement between the two issatisfactory in all cases. Thus, without having recorded theparticle dynamics, we obtain a good approximation to theparticle statistics from the time-averaged wave field η.
V. DYNAMICS AND BIFURCATIONS

The optimized DMD algorithm is used on the wave data
shown in Fig. 4 using a low-rank truncation suggested by
Fig. 5. Specifically, as the forcing parameter  takes on thevalues 4.8, 5.0, 5.3, and 5.4. the low-rank structure is well
captured by one, three, five, and nine modes, respectively. As
will be detailed below, the number of modes used was dictatedby the minimum number of modes required to reproduce anaccurate representation of the dynamics.
Figures 6–9 illustrate the dynamics as the dynamics transitions from the lowest value of  = 4.8 to the highest
value of  = 5.4. Each of the four figures have six panels
which illustrate the characteristics dynamics as the forcing parameter  takes on the values 4.8, 5.0, 5.3, and 5.4.
Figures 6(a), 7(a), 8(a), and 9(a) show the wave dynamics

034213-6

PILOT-WAVE DYNAMICS: USING DYNAMIC MODE …

PHYSICAL REVIEW E 108, 034213 (2023)

FIG. 8. Pilot-wave dynamics for  = 5.3. The simulation was
run of 8000 Faraday periods. Figure 5 shows that a rank-5 decomposition is adequate to model the spatiotemporal periodic behavior
observed. (a) Full PDE evolution. (b) Five-mode DMD approximation. (c) DMD eigenvalues. (d) DMD eigenfunctions used forreconstruction of the spatiotemporal wave dynamics. (The black
dotted line shows the steady-state background mode.) Note thatthe DMD eigenvalues appear in complex conjugate pairs whose
frequencies are approximately harmonic, leading to the secondary
period-doubling dynamics manifest in the pilot-wave field. (e) Timeaveraged pilot-wave field as computed directly (blue solid line)
and deduced from the particle probability distribution via (11) (reddashed line). (f) Particle probability distribution as computed via
particle tracking (blue solid line) and inferred from the mean pilot
wave by inverting Eq. (11) (red dashed line).

after the decay of all transients in the system (Fig. 2).
Recall that the wave field has been strobed to remove oscillations at the Faraday frequency. Figures 6(a), 7(a), 8(a),and 9(a) show the evolution of the full PDE simulations, while
Figs. 6(b), 7(b), 8(b), and 9(b) show the low-rank DMD reconstruction. The accompanying DMD eigenvalue and DMD
eigenfunction are shown in Figs. 6(c), 7(c), 8(c), and 9(c) and
Figs. 6(d), 7(d), 8(d), and 9(d), respectively. Note that the
eigenvalue at the origin captures a stationary behavior, while
the complex conjugate pairs of eigenfunctions capture theperiodic dynamics of the wave field. Figures 6(e), 7(e), 8(e),and 9(e) show the mean pilot-wave field (blue solid line) andthat deduced from Durey’s convolution theorem [37,41] (reddashed line). Figures 6(f), 7(f), 8(f), and 9(f) show the particle
probability distribution as computed directly from particle
tracking (blue solid line) and inferred from the mean pilot
wave by inverting Eq. (11) (red dashed line).

FIG. 9. Pilot-wave dynamics for the highest forcing value considered,  = 5.4. The simulation was run for 8000 Faraday periods.
Figure 5 shows that a rank-9 decomposition is adequate to model
the spatiotemporal chaotic behavior observed. (a) Full PDE evolution. (b) Nine-mode DMD approximation. (c) DMD eigenvalues.
(d) DMD eigenfunctions used for reconstruction of the spatiotemporal wave dynamics. (The black dotted line shows the steady-statebackground mode as computed over the limited time interval shown.)
Note that the DMD eigenvalues appear in complex conjugate pairings whose frequencies are incommensurate, which is responsible
for the relatively complex dynamics. (e) Time-averaged pilot-wave
field computed directly (blue solid line) and that deduced from themeasured particle probability distribution via (11) (red dashed line).
(f) Particle probability distribution as computed via particle tracking(blue solid line) and inferred from the mean pilot wave by inverting
Eq. (11) (red dashed line).

Note that as the bifurcation sequence progresses, a perioddoubling bifurcation is observed with the generation ofharmonics in the DMD spectra (see Figs. 7 and 8). Specifically, we note that the simulations reported were executed justbeyond the bifurcation point at which the period doubling is
expected to hold exactly; thus, the frequencies in normalized
units are 0.17 and 0.32. The period doubling of the pilot wave
corroborates the well-characterized period-doubling behaviorof the droplet itself as reported in a number of settings involving orbital dynamics [11,35,38,39]. Once the forcing issufficiently large, the complex conjugate pairs are no longer
harmonic, thus producing the observed spatiotemporal chaos.
The overall period-doubling route to chaos is thus evident in
the simulations. Moreover, as the forcing becomes sufficiently
large, the DMD spectrum generated from the optimized DMD
regression produces eigenvalues with small but nonzero real

034213-7

KUTZ, NACHBIN, BADDOO, AND BUSH

PHYSICAL REVIEW E 108, 034213 (2023)

TABLE I. Comparison of the quantum description of a particle in a one-dimensional well and the pilot-wave system considered here.

driving parameter
waveform
statistical inferenceaQuantum mechanics

Pilot-wave hydrodynamics

particle energy
quantum wave function
Born’s rule

memory parameter γ
Faraday pilot wave h(x, t )integral operator inversiona

From [41].

parts in the spectrum of (9). This is expected since the regression is not forced to find eigenvalues constrained to theimaginary axis.
The evolution from Fig. 6 to Fig. 9 shows the underlying
onset of instability in this pilot-wave hydrodynamic system.
This period doubling is canonical for damped-driven systems,
from optics to combustion [52–55]. The optimized DMD
algorithm provides a data-driven algorithm that provides aregression to the exponential solution form (1) commonlyused to model physical systems, from quantum mechanics toelectrodynamics. Thus, DMD provides interpretable models
and a clear quantification of the evolution of the pilot-wave
field accompanying the droplet dynamics.

VI. CONCLUSION

Particle-wave interactions arise throughout the physicalsciences. The specific example considered here of the pilotwave hydrodynamic system has generated significant interest
due to its connection to quantum mechanics. Indeed, pilotwave hydrodynamics has provided a compelling example ofwell-resolved classical particle-wave interactions producing
quantumlike statistics. To date, the particle dynamics hasreceived the bulk of the mathematical attention, and the emergent statistics have been seen as a compelling feature of thesystem that is not always simply rationalized. Here we haveinstead focused on the evolution of the pilot wave and the inference of system statistics. By leveraging the dynamic mode
decomposition, the pilot-wave dynamics has been shown
to execute a period-doubling cascade typically observed in
damped-driven systems ranging from detonation waves tomode-locked lasers. Dynamic mode decomposition provides
a low-rank approximation of the pilot-wave dynamics into
a set of spatial modes with associated temporal frequencies.
The regression framework of DMD provides a best-fit linear
dynamics model over snapshots of spatiotemporal data. The
DMD characterization of the wave field yields a different
perspective on the walking-droplet problem that forges valuable links with quantum mechanics. In particular, it naturally
decomposes the wave field into modes of the form prevalent
in standard quantum theory.
Our analysis has shown that as the vibrational acceleration
is increased progressively, the pilot wave undergoes a series ofHopf bifurcations in which new modes at approximately harmonic frequencies emerge, culminating in a period-doublingcascade to spatiotemporal chaos. Such a period-doublingroute to chaos is a canonical feature of damped-driven systems, which can be related to the logistic map in which thesame canonical bifurcation structure is evident. The simplicityof the logistic map belies the rich and complicated behavior
that it captures [68]. Given the diversity of models capableof producing this same bifurcation structure [52–55], it ishighly suggestive that energy balance considerations alone
can dictate the overall physics in such damped-driven systems.
Specifically, the gain and loss dynamics in the hydrodynamic
pilot-wave system produce a mapping between the driving
energy input and the damping losses that generate a perioddoubling cascade [39].
Our study has demonstrated how classical pilot-wave dynamics, like quantum mechanics, may yield predictions forthe statistical behavior of particles on the basis of a wave theory. The map between the quantum and classical pilot-wave
descriptions of particle motion in a one-dimensional well ispresented in Table I. In quantum mechanics, the relevant wave
form is the wave function , while in our system it is the pilot
wave. In quantum mechanics, the number of modes exciteddepends on the particle energy. In our pilot-wave system, thesystem memory plays an analogous role. Note that as either ofthese control parameters is increased progressively, a discrete
set of new wave modes is introduced. In quantum mechanics, the waves correspond to the complex wave function ;in our system, the waves correspond to the pilot wave. In
quantum mechanics, the statistical behavior of the particles isprescribed by Born’s rule: The density of states is prescribedby the square of the wave function | |2 . In our pilot-wave
system, the density of states is inferred from the mean pilot
wave via the operator inversion suggested by the theorem of
Durey et al. [41].
In summary, the DMD algorithm is a regression to exponential solutions of the form (1). In quantum mechanics,
the solution (1) is typically constructed using analytical (orsemianalytical) techniques where appropriate boundary andmatching conditions are imposed for a given potential. Since
Schrödinger’s equation is linear, such an eigendecomposition
provides a set of linear basis modes which can be superimposed to express any solution. Dynamic mode decomposition
is a regression directly to the solution form (1), using variable
projection to compute amplitudes b j , modes φ j , and frequencies ω j of the exponential solutions directly from data. In thequantum scenario, one must know the potential in order toconstruct the solutions. In contrast, DMD can be used in acompletely data-driven manner without knowledge of the underlying potential. Thus DMD is advantageous for scenarios
where an unknown effective potential determines the underlying dynamics. Moreover, it provides the best approximation of
the form (1) even when the underlying dynamics is nonlinear,such as is the case with the hydrodynamic analog system.

034213-8

PILOT-WAVE DYNAMICS: USING DYNAMIC MODE …

PHYSICAL REVIEW E 108, 034213 (2023)

FIG. 10. Effects of terminating the NN-FCGLS algorithm at different numbers of iterations. If the algorithm is terminated too early thenthe least-squares residual in (A2) is large; if the algorithm is terminated too late then the ensuing PDF is too jagged. The plotted PDFs in the
top row are computed with 20, 80, and 200 iterations, respectively.
ACKNOWLEDGMENTS

P.J.B. and J.N.K. acknowledge support from the National
Science Foundation AI Institute in Dynamic Systems (Grant
No. 2112085). J.W.M.B. gratefully acknowledges support
from the National Science Foundation through Grant No.
CMMI-2154151. A.N. acknowledges support from CNPq
under (PQ1D) Grant No. 307078/2021-3 and FAPERJ Cientistas do Nosso Estado Project No. E-26/201.156/2021.
APPENDIX: CALCULATING THE PARTICLE
PROBABILITY DENSITY FUNCTIONFROM THE AVERAGE WAVE FIELD(A1)

where {μ}i = μ(xi ), {η}i = η(xi ), and {G}i, j = ηG (xi , x j ).
Our goal is now to find a μ that solves (A1) for a given
Green’s function matrix G and an observed mean wave field
η. The wave field is observed over a finite time horizon, so theresulting time average invariably contains noise. As such, theobserved average wave differs from the true average wave, so
we writeη = ηtrue + ,μ = arg minμ Gμ − η22 .

(A2)

Additionally, we are only interested in solutions that represent
probability density functions. In other words, μ must be nonnegative and have unit mass, specifically,μ  0,In this Appendix we detail the practical steps involvedin computing the particle probability density function (PDF)
from the average wave field. We use the full time window of
the simulated system (3)–(6) to calculate the average wave
field and the true PDF. Discretizing (11) on a grid of equispaced points {x1 , . . . , x128 } produces the system
η = Gμ,where  is a vector of noise corresponding to the temporal
truncation error. Therefore, we consider a least-squares optimization problem to find the particle histogram: We seek the
optimal vector μ∗ defined by xμ 1 = 1,

(A3)

where x = x2 − x1 is the equispaced quadrature weight.
The above constraints effectively regularize the problem (A2)and remove unphysical solutions. However, the matrix G is
ill-conditioned [with condition number O(1010 )], so taking
an inverse or repeatedly applying G can be unreliable. Furthermore, G has a large approximate null space, so many
candidate PDFs produce small residuals in (A2).
We deployed several classical methods [69] to solve theinverse problem in (A2). Ultimately, we found that the
non-negative flexible conjugate-gradient least-squares (NNFCGLS) method of Gazzola and Wiaux [70] provided themost robust results. The NN-FCGLS method is a Krylov
subspace method and is implemented in the IR TOOLS package [71]. Although the NN-FCGLS method does not enforce
the unit mass constraint, the computed solutions have mass
very close to one.
It is important to control the number of iterations executedby the NN-FCGLS method, as illustrated in Fig. 10. Therein,

034213-9

KUTZ, NACHBIN, BADDOO, AND BUSH

PHYSICAL REVIEW E 108, 034213 (2023)

we plot the least-squares residual (Gμ − η2 ) and the error
between the calculated and the true PDFs (μ − μtrue 1 ) whenμ is calculated after the given number of iterations. In practice, the true PDF is unavailable, so we include the PDF error
here only for the sake of illustration. We select the iteration
number that produces a PDF satisfying two criteria: It must
attain a small residual error in (A2) and it must be sufficiently smooth. Too few iterations can result in large residualsin (A2) and too many iterations result in a jagged (nonsmooth)solution. Since the data are generated from a continuous dynamical system, we can exclude solutions that are not smooth.
Chapter 2 of [69] suggests alternative criteria for selecting thenumber of iterations of conjugate gradient methods.
Figure 10 shows three candidate PDFs that we evaluate
according to these two criteria. The first PDF is produced after

20 NN-FCGLS iterations. It is smooth, but has a low leastsquares residual, so we discard it. The third PDF is produced
after 200 iterations. It has a small residual, but is jagged, so
we discard it. The second PDF is produced after 80 iterations.
It has a small residual and is smooth. Accordingly, we selectthis PDF as our solution μ . The green solid curve in Fig. 10
confirms that this is a good choice of PDF. Indeed, selecting
anywhere between 70 and 100 NN-FCGLS iterations would
have produced a satisfactory PDF.
Finally, an advantage of using Krylov methods is that they
can be applied when one has access only to a linear operator
that evaluates matrix-vector products Gv, as opposed to G
itself. Accordingly, the NN-FCGLS method scales well andwill be suitable for higher-dimensional problems, such as aparticle in a circular corral [18,41].

[1] L. de Broglie, Ondes et quanta, Comptes Rendus 177, 507(1923).
[2] L. de Broglie, Interpretation of quantum mechanics by the double solution theory, Ann. Fond. L. de Broglie 12, 1 (1987).
[3] J. W. M. Bush, Pilot-wave hydrodynamics, Annu. Rev. Fluid
Mech. 47, 269 (2015).
[4] J. W. M Bush and A. U. Oza, Hydrodynamic quantum analogs,
Rep. Prog. Phys. 84, 017001 (2020).
[5] Y. Couder, S. Protiere, E. Fort, and A. Boudaoud, Walking andorbiting droplets, Nature (London) 437, 208 (2005).
[6] Y. Couder and E. Fort, Single-Particle Diffraction and Interference at a Macroscopic Scale, Phys. Rev. Lett. 97, 154101(2006).
[7] J. W. M. Bush, Y. Couder, T. Gilet, P. A. Milewski, and A.
Nachbin, Introduction to focus issue on hydrodynamic quantum
analogues, Chaos 28, 096001 (2018).
[8] G. Pucci, P. J. Sáenz, L. M. Faria, and J. W. M. Bush, Nonspecular reflection of walking droplets, J. Fluid Mech. 804, R3(2016).
[9] C. Ellegaard and M. T. Levinsen, Interaction of wave-driven
particles with slit structures, Phys. Rev. E 102, 023115(2020).
[10] E. Fort, A. Eddi, J. Moukhtar, A. Boudaoud, and Y. Couder,Path-memory induced quantization of classical orbits, Proc.
Natl. Acad. Sci. USA 107, 17515 (2010).
[11] D. M. Harris and J. W. M. Bush, Drops walking in a rotating
frame: From quantized orbits to multimodal statistics, J. Fluid
Mech. 739, 444 (2014).
[12] A. U. Oza, D. M. Harris, R. R. Rosales, and J. W. M. Bush,
Pilot-wave dynamics in a rotating frame: On the emergence oforbital quantization, J. Fluid Mech. 744, 404 (2014).
[13] S. Perrard, M. Labousse, M. Miskin, E. Fort, and Y. Couder,Self-organization into quantized eigenstates of a classical wavedriven particle, Nat. Commun. 5, 3219 (2014).
[14] M. Labousse, S. Perrard, Y. Couder, and E. Fort, Build-up of
macroscopic eigenstates in a memory-based constrained system, New J. Phys. 16, 113027 (2014).
[15] A. Eddi, E. Fort, F. Moisy, and Y. Couder, Unpredictable Tunneling of a Classical Wave-Particle Association, Phys. Rev.
Lett. 102, 240401 (2009).
[16] A. Nachbin, P. A. Milewski, and J. W. M. Bush, Tunneling witha hydrodynamic pilot-wave model, Phys. Rev. Fluids 2, 034801(2017).

[17] L. Tadrist, T. Gilet, P. Schlagheck, and J. W. M. Bush, Predictability in a hydrodynamic pilot-wave system: Resolution ofwalker tunneling, Phys. Rev. E 102, 013104 (2020).
[18] D. M. Harris, J. Moukhtar, E. Fort, Y. Couder, and J. W. M.
Bush, Wavelike statistics from pilot-wave dynamics in a circularcorral, Phys. Rev. E 88, 011001(R) (2013).
[19] P. J. Sáenz, T. Cristea-Platon, and J. W. M. Bush, Statisticalprojection effects in a hydrodynamic pilot-wave system, Nat.
Phys. 14, 315 (2018).
[20] P. J. Sáenz, T. Cristea-Platon, and J. W. M. Bush, A hydrodynamic analog of Friedel oscillations, Sci. Adv. 6, 20 (2020).
[21] P. J. Saenz, G. Pucci, S. E. Turton, A. Goujon, R. R. Rosales,
J. Dunkel, and J. W. M. Bush, Emergent order in hydrodynamic
spin lattices, Nature (London) 596, 58 (2021).
[22] A. Nachbin, Walking droplets correlated at a distance, Chaos
28, 096110 (2018).
[23] K. Papatryfonos, M. Ruelle, C. Bourdiol, A. Nachbin, J. W. M.
Bush, and M. Labousse, Hydrodynamic superradiance via
wave-mediated cooperative tunneling, Commun. Phys. 5, 142(2022).
[24] A. Nachbin, Effect of isolation on two-particle correlations inpilot-wave hydrodynamics, Phys. Rev. Fluids 7, 093604 (2022).
[25] A. Eddi, E. Sultan, J. Moukhtar, E. Fort, M. Rossi, and Y.
Couder, Information stored in Faraday waves: The origin of apath memory, J. Fluid Mech. 674, 433 (2011).
[26] J. Moláček and J. W. M. Bush, Droplets bouncing on a vibrating
fluid bath, J. Fluid Mech. 727, 582 (2013).
[27] J. Moláček and J. W. M. Bush, Drops walking on a vibrating bath: Towards a hydrodynamic pilot-wave theory, J. Fluid
Mech. 727, 612 (2013).
[28] A. U. Oza, R. R. Rosales, and J. W. M. Bush, A trajectory equation for walking droplets: Hydrodynamic pilot-wave theory,
J. Fluid Mech. 737, 552 (2013).
[29] J. Arbelaiz, A. U. Oza, and J. W. M. Bush, Promenading pairsof walking droplets: Dynamics and stability, Phys. Rev. Fluids
3, 013604 (2018).
[30] A. U. Oza, E. Siéfert, D. M. Harris, J. Moláček, and J. W. M.
Bush, Orbiting pairs of walking droplets: Dynamics and stability, Phys. Rev. Fluids 2, 053601 (2017).
[31] M. M. P. Couchman, S. E. Turton, and J. W. M. Bush,
Bouncing phase variations in pilot-wave hydrodynamics andthe stability of droplet pairs, J. Fluid Mech. 871, 212(2019).

034213-10

PILOT-WAVE DYNAMICS: USING DYNAMIC MODE …

PHYSICAL REVIEW E 108, 034213 (2023)[32] S. J. Thomson, M. M. P. Couchman, and J. W. M. Bush, Collective vibrations of confined levitating droplets, Phys. Rev. Fluids
5, 083601 (2020).
[33] S. J. Thomson, M. Durey, and R. R. Rosales, Collective vibrations of a hydrodynamic active lattice, Proc. R. Soc. A 476,
20200155 (2020).
[34] M. M. P. Couchman and J. W. M. Bush, Free rings of bouncing droplets: Stability and dynamics, J. Fluid Mech. 903, A49(2020).
[35] A. U. Oza, Ø. Wind-Willassen, D. M. Harris, R. R. Rosales, andJ. W. M. Bush, Pilot-wave hydrodynamics in a rotating frame:
Exotic orbits, Phys. Fluids 26, 082101 (2014).
[36] M. Durey, Bifurcations and chaos in a Lorenz-like pilot-wave
system, Chaos 30, 103115 (2020).
[37] M. Durey, P. A. Milewski, and J. W. M. Bush, Dynamics, emergent statistics and the mean-pilot-wave potential of walking
droplets, Chaos 28, 096108 (2018).
[38] L. D. Tambasco, D. M. Harris, A. U. Oza, R. R. Rosales,
and J. W. M. Bush, The onset of chaos in orbital pilot-wave
dynamics, Chaos 26, 103107 (2016).
[39] A. Rahman and J. N. Kutz, Walking droplets as adamped-driven system, SIAM J. Appl. Dyn. Sys. 22, 1219
(2023).
[40] A. Rahman and D. Blackmore, Walking droplets through thelens of dynamical systems, Mod. Phys. Lett. B 34, 2030009(2020).
[41] M. Durey, P. A. Milewski, and Z. Wang, Faraday pilotwave dynamics in a circular corral, J. Fluid Mech. 891, A3(2020).
[42] P. J. Schmid and J. Sesterhenn, in Proceedings of the 61st
Annual Meeting of the APS Division of Fluid Dynamics, San
Antonio, 2008 (American Physical Society, Ridge, 2008).
[43] P. J. Schmid, Dynamic mode decomposition of numerical andexperimental data, J. Fluid Mech. 656, 5 (2010).
[44] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S.
Henningson, Spectral analysis of nonlinear flows, J. Fluid
Mech. 641, 115 (2009).
[45] J. N. Kutz, Data-Driven Modeling & Scientific Computation:
Methods for Complex Systems & Big Data (Oxford University
Press, Oxford, 2013).
[46] T. Askham and J. N. Kutz, Variable projection methods for an
optimized dynamic mode decomposition, SIAM J. Appl. Dyn.
Syst. 17, 380 (2018).
[47] D. Sashidhar and J. N. Kutz, Bagging, optimized dynamic mode
decomposition (BOP-DMD) for robust, stable forecasting withspatial and temporal uncertainty-quantification, Phil. Trans. R.
Soc. A. 380, 20210199 (2022).
[48] H. Lange, S. L. Brunton, and J. N. Kutz, From Fourier to Koopman: Spectral methods for long-term time series prediction, J.
Mach. Learn. Res. 22, 1881 (2021).
[49] K. K. Chen, J. H. Tu, and C. W. Rowley, Variants ofdynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses, J. Nonlinear Sci. 22, 887(2012).
[50] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex
Systems (SIAM, Philadelphia, 2016).
[51] P. A. Milewski, C. A. Galeano-Rios, A. Nachbin, and J. W. M.
Bush, Faraday pilot-wave dynamics: Modelling and computation, J. Fluid Mech. 778, 361 (2015).

[52] F. Li, P. Kong A. Wai, and J. N. Kutz, Geometrical descriptionof the onset of multi-pulsing in mode-locked laser cavities,
J. Opt. Soc. Am. B 27, 2068 (2010).
[53] K. M. Spaulding, D. H. Yong, A. D. Kim, and J. N. Kutz,
Nonlinear dynamics of mode-locking optical fiber ring lasers,
J. Opt. Soc. Am. B 19, 1045 (2002).
[54] E. Ding and J. N. Kutz, Operating regimes, split-step modeling,
and the Haus master mode-locking model, J. Opt. Soc. Am. B
26, 2290 (2009).
[55] J. Koch, M. Kurosaka, C. Knowlen, and J. N. Kutz, Multiscale
physics of rotating detonation waves: Autosolitons and modulational instabilities, Phys. Rev. E 104, 024210 (2021).
[56] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton,
and J. N. Kutz, On dynamic mode decomposition: Theory andapplications, J. Comput. Dyn. 1, 391 (2014).
[57] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, Modern
Koopman theory for dynamical systems, SIAM Rev. 64, 229(2022).
[58] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, andJ. N. Kutz, Chaos as an intermittently forced linear system, Nat.
Commun. 8, 19 (2017).
[59] H. Arbabi and I. Mezić, Ergodic theory dynamic mode decomposition, and computation of spectral properties of theKoopman operator, SIAM J. Appl. Dyn. Syst. 16, 2096 (2017).
[60] M. Kamb, E. Kaiser, S. L. Brunton, and J. N. Kutz, Time-delay
observables for Koopman: Theory and applications, SIAM J.
Appl. Dyn. Syst. 19, 886 (2020).
[61] S. M. Hirsh, S. M. Ichinaga, S. L. Brunton, J. N. Kutz, and B. W.
Brunton, Structured time-delay models for dynamical systems
with connections to Frenet-Serret frame, Serret frame Proc. R.
Soc. A. 477, 20210097 (2021).
[62] S. Bagheri, Effects of weak noise on oscillating flows: Linking
quality factor, Floquet modes, and Koopman spectrum, Phys.
Fluids 26, 094104 (2014).
[63] D. Duke, J. Soria, and D. Honnery, An error analysis of the
dynamic mode decomposition, Exp. Fluids 52, 529 (2012).
[64] S. Bagheri, Koopman-mode decomposition of the cylinder
wake, J. Fluid Mech. 726, 596 (2013).
[65] S. T. M. Dawson, M. S. Hemati, M. O. Williams, and C. W.
Rowley, Characterizing and correcting for the effect of sensor
noise in the dynamic mode decomposition, Exp. Fluids 57, 42(2016).
[66] M. S. Hemati, C. W. Rowley, E. A. Deem, and L. N. Cattafesta,
De-biasing the dynamic mode decomposition for applied Koopman spectral analysis, Theor. Comput. Fluid Dyn. 31, 349(2017).
[67] P. J. Baddoo, B. Herrmann, B. J. McKeon, J. N. Kutz, andS. L. Brunton, Physics-informed dynamic mode decomposition,
Proc. R. Soc. A 479, 20220576 (2023).
[68] R. M. May, Simple mathematical models with very complicateddynamics, Nature 261, 459 (1976).
[69] J. P. Kaipio and E. Somersalo, Statistical and Computational
Inverse Problems, Applied Mathematical Sciences Vol. 160(Springer, New York, 2005).
[70] S. Gazzola and Y. Wiaux, Fast nonnegative least squares
through flexible Krylov subspaces, SIAM J. Sci. Comput. 39,
A655 (2017).
[71] S. Gazzola, P. C. Hansen, and J. G. Nagy, IR Tools: A MATLAB
package of iterative regularization methods and large-scale test
problems, Numer. Algor. 81, 773 (2019).

034213-11