The eigenspectrum of the Koopman operator enables the decomposition of nonlinear dynamics into a sum of nonlinear functions of the state space with purely exponential and sinusoidal time dependence. For a limited number of dynamical systems, it is possible to find these Koopman eigenfunctions exactly and analytically. Here, this is done for the Korteweg–de Vries equation on a periodic interval using the periodic inverse scattering transform and some concepts of algebraic geometry. To the authors’ knowledge, this is the first complete Koopman analysis of a partial differential equation, which does not have a trivial global attractor. The results are shown to match the frequencies computed by the data-driven method of dynamic mode decomposition (DMD). We demonstrate that in general, DMD gives a large number of eigenvalues near the imaginary axis and show how these should be interpreted in this setting.

## I. INTRODUCTION

The Koopman operator was introduced by Koopman^{1} to describe the nonlinear behavior of a dynamical system as the linear evolution of nonlinear observables of that system. It is well known that a finite dimensional nonlinear system can be converted into an infinite dimensional linear system; Koopman analysis extends this to infinite dimensional nonlinear systems. Recent interest in the Koopman operator was initiated by Mezić,^{2,3} and it has become popular for its close connection to the computational method of dynamic mode decomposition.

DMD was originally invented by Schmid^{4} as a method for simultaneously extracting the important spatial and temporal features of a timeseries. It has been successfully applied to a wide array of numerical, observational, and experimental data, most notably in fluid dynamics.^{5–7} Under certain conditions,^{8} the results of DMD can be seen as a numerical approximation to the Koopman modes and eigenvalues of the underlying dynamical system, and therefore, understanding the Koopman operator aids interpretation of the results of DMD. In particular, to understand the spatial patterns called DMD modes, it would be helpful to have analytic results on nonlinear partial differential equations (PDEs).

Several authors^{9–11} have successfully performed Koopman analysis on the Burgers equation, a dissipative nonlinear PDE, which can be transformed into a linear heat equation with a suitable change of variables. Nakao and Mezić^{12} considered the Burgers equation as well as a non-trivial transformation of it to the phase-diffusion equation. Though insightful, these PDEs admit only a single steady global attractor and are strongly dissipative, which precludes much interesting nonlinear behavior relevant to applying DMD to situations involving sustained waves, a common use-case.

Parker and Page^{13} considered the Korteweg–de Vries (KdV) equation, an integrable partial differential equation of one variable. This Hamiltonian dynamical system behaves very differently from the dissipative systems mentioned above. In that work, some Koopman eigenfunctions were found for soliton solutions of the KdV equation on the real line, which excluded a large class of solutions; in particular, this excluded the spatially periodic solutions, which are the natural choice for computer simulations of solitons. It was argued that the simplest periodic solutions, cnoidal waves, give purely imaginary Koopman eigenvalues, in sharp contrast to the isolated solitons, which have purely real Koopman eigenvalues, despite being a naturally limiting case of the former.

A periodic domain is the natural setting for numerical solutions of one-dimensional PDEs, such as the KdV equation. Indeed, early pioneering numerical simulations of the KdV equation on a periodic domain^{14} were used to shed light on the famous Fermi–Pasta–Ulam–Tsingou problem,^{15} in which the initial state of a system recurs to arbitrary precision after complex nonlinear dynamics. It was later proven^{16} that this is because the solutions lie on quasiperiodic invariant tori. We will show that it is possible to use DMD to determine the underlying frequencies. In fact, when the dynamics are confined to an invariant torus, DMD is equivalent to a Fourier transform in time, though it generalizes when other modes grow or decay.

In the present paper, we study a particular but very general class of solutions of the KdV equation on a periodic interval, for which we are able to define Koopman eigenfunctions. These eigenfunctions require the evaluation of contour integrals on Riemann surfaces, which must be performed numerically in all but the simplest of cases. This method is not, therefore, recommended as a general approach for nonlinear PDEs, where DMD could easily be applied. However, the semi-formal mathematical treatment presented here gives instructive results: we will see that in this case, Koopman eigenvalues necessary for the decomposition of the state of the system densely fill the imaginary axis, and therefore, the results of DMD are subtle to interpret. This paper proceeds as follows: in Sec. II, the KdV equation is introduced, and its history and significance are briefly described; in Sec. III, the Koopman operator and its spectrum are defined; in Sec. IV, we then derive Koopman eigenfunctions for the KdV equation; Sec. V presents the results of applying these to an example solution, which is compared to the numerical results of DMD in Sec. VI. Concluding remarks are given in Sec. VII.

## II. THE KdV EQUATION

^{17}The equation was derived by Korteweg and de Vries

^{18}to describe the weakly nonlinear evolution of shallow water waves propagating in one direction. The pioneering computational results of Zabusky and Kruskal

^{14}demonstrated the existence of soliton solutions when this equation is solved on a finite periodic domain. For an infinite domain, the celebrated inverse scattering method of Gardner

*et al.*

^{19}gives a straightforward procedure to solve the equation and provides intuitive interpretations for the existence of solitons as conserved quantities. Analytical results on a periodic interval have proven much more complicated, despite the early computation successes and the well known cnoidal wave solution.

In addition to shallow water waves, the KdV equation naturally arises from weakly nonlinear theory in many physically relevant flows.^{20–22} Our interest, however, derives from the fact that the equation admits complex, nonlinear but non-chaotic solutions amenable to analytic treatment. Unlike other PDEs for which Koopman spectra have been derived, it is a Hamiltonian system with an infinite number of conserved quantities, rather than having one unique attractor.

The solution to (1) is well-posed on $ T= R/2\pi Z$, the periodic domain of length $2\pi $, for initial conditions in Sobolev spaces $ H s( T, R)$ with $s\u2265\u22121$,^{23} with a well-defined evolution operator $ S t: H s( T, R)\u2192 H s( T, R)$ for each time $t$. A single periodic wave of sufficient amplitude breaks down into a spatially periodic solution, with quasiperiodic behavior in time. That is to say, solutions lie on an invariant torus, as is usual for integrable Hamiltonian systems, and these invariant tori foliate phase space. In fact, any solution of (1) can be approximated to arbitrary precision as an invariant $M$-torus,^{16}^{,} $M\u2208 N$. These are the so-called finite gap solutions, which will be the focus of the present study. Practically, the evolution of an arbitrary initial condition can be approximated by truncating the scattering data at a judiciously chosen $M$.^{24}

Let us, therefore, define our solution space $ \Omega M$ to be the subset of $ L 2( T, R)$ for which there are $M$-gap solutions (this terminology should become clearer in Sec. IV). This is a well-posed invariant subspace. For convergent Koopman decompositions, it will be necessary to further restrict this space, in a manner analogous to that of Balabane *et al.*,^{11} in Sec. V.

^{25}

^{,}

^{9–11}though in this case, the transformation seems at first glance to have made the equation more complicated. The utility comes from the fact that $\u03d1$ can be expressed as a Riemann theta function, as explained in Sec. IV.

## III. THE KOOPMAN OPERATOR

*Koopman eigenvalues*$\lambda \u2208 C$ and

*Koopman eigenfunctions*$\phi \u2208 O M$ satisfy

*Koopman modes*, first defined by Mezić,

^{2}which encode spatial information for each eigenvalue and are independent of the particular choice of initial condition $ u 0$, whose contribution is included in the value of $ \phi \nu ( u 0)$. If this is possible, it means that the dynamics of (1) can be decomposed as a sum over nonlinear functions whose temporal behavior is purely oscillatory in time.

## IV. KOOPMAN EIGENFUNCTIONS OF THE KdV EQUATION

It would take a whole textbook to fully explain the periodic inverse scattering transform. We refer readers to Novikov *et al.*,^{26} Belokolos *et al.*,^{17} and Osborne^{27} for accessible introductions, including the necessary background of Riemann surfaces and theta functions, though note the differing notations and conventions between these (we follow the notation of Belokolos *et al.*^{17}). Here, we give a summary of the relevant results for the KdV equation, which are implemented in the Mathematica notebook given in the supplementary material.

^{28}a pair of linear operators $L(u),A(u): L 2( R)\u2192 L 2( R)$ such that

^{29}The admissible eigenvalues $\lambda $ for a bounded eigenfunction $\psi $ reside in intervals $[ E 1, E 2]$, $[ E 3, E 4]$, …, where $\u2212\u221e< E 1< E 2\u2264 E 3< E 4\u2264 E 5<\cdots $. Outside these regions, only unbounded solutions are possible, and these are termed forbidden gaps. The $ E k$ are the values of $\lambda $ for which $F(\lambda )$, defined as half the trace of the monodromy matrix of (11), is $\xb11$

^{29}(see Fig. 1). Though the monodromy matrix is not invariant under the dynamics (1), its trace, and, therefore, also the $ E k$ are invariant. We explicitly consider only the case when there is a finite number $g$ of non-degenerate forbidden gaps $( E 2, E 3)$, $( E 4, E 5)$, …, so that $ E 2 k= E 2 k + 1$ for all $k>g$.

*genus*of this surface is simply the number of gaps $g$. It is then possible to define a basis of contours $ a j$ and $ b j$ ( $1\u2264j\u2264g$) for the Riemann surface such that any contour can be expressed, up to continuous deformations, as a sum of $ a j$ and $ b j$. Such a choice of basis is not unique and will have implications for the final results—see discussions of the wave basis and the soliton basis in Osborne.

^{27}Our convention is shown in Fig. 3. We also define a basis of holomorphic differentials on this surface,

*period matrix*of the Riemann surface

^{17}

As the frequencies $ W j$ are incommensurate, for genus $g=2$ and greater, these Koopman eigenvalues densely fill the imaginary axis. This is a significant complication over previously studied PDEs. Note that $ \phi m 1(u) \phi m 2(u)= \phi m 1 + m 2(u)$, and $ \phi 0(u)=1$. Indeed, rather generally, integer linear combinations of Koopman eigenvalues are eigenvalues and products of Koopman eigenfunctions are Koopman eigenfunctions,^{30} though these eigenfunctions are not necessarily useful for decompositions in all cases.

## V. KOOPMAN DECOMPOSITIONS

## VI. DYNAMIC MODE DECOMPOSITION

In the original and most basic form of the algorithm, the number of DMD modes that come from the spectral decomposition will be equal to the spatial dimension. We can increase both the robustness and number of discovered of DMD modes found with delay embedding, so-called Hankel-DMD,^{31} a higher-order extension in which temporal resolution is substituted for spatial resolution.^{7} We found much better results when employing this method; we used $20$ delays.

As DMD is designed to detect the important temporal frequencies of the dynamics, it should be possible to use it to reconstruct an approximation $ W ~$ for $W$ from a time series. However, as discussed in Sec. IV, the Koopman eigenvalues densely fill the imaginary axis, and therefore, the results of DMD are obscured. For example, if we believe that the solution will be well-represented by an invariant 2-torus—and the DMD eigenfrequencies $ \omega k$ are sufficiently well-resolved—we expect to see frequencies $i\omega = m 1 W 1+ m 2 W 2$ for many different $ m 1, m 2\u2208 Z$.

However, due to the exponential term in Eq. (24), we expect the amplitudes to be very low for large values of $ m k$ (in absolute terms), and thus, we concentrate on the largest amplitude DMD modes, associating these with small values of $ m 1$ and $ m 2$. We can guess the two fundamental frequencies $ W ~ 1$ and $ W ~ 2$ to be high amplitude modes, and furthermore, the smaller of these should be equal to the gap between the modes. See Fig. 6 for a visual representation of this. We can extend this argument in an obvious way for $g>2$.

## VII. DISCUSSION

We have performed a Koopman decomposition of the periodic KdV equation. This is almost immediate once the convoluted but well-defined process of periodic inverse scattering is performed. Additionally, we have shown how this result relates to the output of DMD for such a system. DMD gives a very large number of near-imaginary eigenvalues associated with the different harmonics of the nonlinearly interacting waves, which correspond to the Koopman eigenvalues found analytically, which densely fill the imaginary axis.

Furthermore, by exploiting the $\theta $-function representation of the solution, we are able to use DMD to approximately recover the necessary parameters.

We note in passing that since we expect purely imaginary eigenvalues in our system, it is a potential use-case for the physics-informed DMD method^{32} of finding a unitary matrix $A$ to fit the data. However, we found that this obfuscates the results, as it prevents the use of the real part of the eigenvalue as a measure for how well-resolved a given mode is.

The analytic results of this paper could be extended to other integrable PDEs, which admit Lax pairs, such as the nonlinear Schrödinger equation, the sine-Gordon equation, or the Kadomtsev–Petviashvili equation. The latter could be particularly insightful since it describes two-dimensional wave fields, a significant increase in complexity over the one-dimensional PDEs, which have been studied heretofore.

## SUPPLEMENTARY MATERIAL

See the supplementary material that consists of Mathematica notebooks reproducing the results of Secs. IV and V.

## ACKNOWLEDGMENTS

This work started life at the Geophysical Fluid Dynamics (GFD) summer school at the Woods Hole Oceanographic Institution, and Sec. VI represents a small part of the fellow’s project of C.V. The authors wish to thank Peter Schmid for his help with this project and everyone at GFD for many fruitful discussions. J.P.P. would like to thank Al Osborne for some useful pointers. J.P.P. was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant No. 865677). C.V. was supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1839302.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jeremy P. Parker:** Conceptualization (lead); Investigation (equal); Methodology (lead); Software (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal). **Claire Valva:** Conceptualization (supporting); Investigation (equal); Methodology (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the supplementary material.

### APPENDIX: RIEMMAN THETA FUNCTION DETERMINATION USING DMD

Here, we give a brief overview of the procedure to recover the parameters for the Riemann theta function, i.e., $U$, $W$, and $ B$, from a time series of $\u03d1$.

We apply our DMD analysis to KdV data with initial condition $ u 0(x)=sin\u2061x$, where we analyze the value of $\u03d1$, rather than $u$. To three decimal places, we recover the same values for the frequencies (29). Figure 7 shows the eigenvalues, along with the three modes corresponding to frequencies $ W ~ 1$ , $ W ~ 2$, and $ W ~ 1+ W ~ 2$, which were used to determine $ B ~$. The DMD spectrum for $\u03d1$ is much cleaner than for $u$, which shows that the Hirota transform has in some sense simplified the dynamics.

## REFERENCES

*Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems*

*Algebro-Geometric Approach to Nonlinear Integrable Equations*

*Non-Linear Waves in Dispersive Media*

*The Direct Method in Soliton Theory*

*Theory of Solitons: The Inverse Scattering Method*

*Nonlinear Ocean Waves and the Inverse Scattering Transform*