We describe a simple mechanical system, a ball rolling along a specially-designed landscape, which mimics the well-known two-bounce resonance in solitary wave collisions, a phenomenon that has been seen in countless numerical simulations but never in the laboratory. We provide a brief history of the solitary wave problem, stressing the fundamental role collective-coordinate models played in understanding this phenomenon. We derive the equations governing the motion of a point particle confined to such a surface and then design a surface on which to roll the ball, such that its motion will evolve under the same equations that approximately govern solitary wave collisions. We report on physical experiments, carried out in an undergraduate applied mathematics course, that seem to exhibit the two-bounce resonance.

Solitary waves are solutions to partial differential equations that maintain their spatial profile while moving at constant speed. A wide variety of systems display a behavior called chaotic scattering when two such waves collide. The two waves may bounce off each other one or more times before escaping to infinity, or they may capture each other and never escape. The number of collisions and the final speed of those that separate depends in a complex way on their initial speeds. To our knowledge, this process has never been observed in a laboratory experiment involving solitary waves. The problem is described well by a small finite-dimensional system of ordinary differential equations. We describe an experiment, a ball rolling on a specially designed surface, that obeys the same system of ordinary differential equation**. We report on laboratory experiments demonstrating that the experimental system has similar dynamics to the solitary wave collisions.**

## I. INTRODUCTION

*Solitary waves*, localized structures that translate at constant velocity while maintaining their spatial profile, are among the most important concept in nonlinear physics. *Solitons* are a type of solitary waves which possess an additional underlying mathematical structure, one which renders the equations solvable. This enforces strikingly simple behavior upon colliding: solitons survive a collision with their shape and velocity unchanged, but with their positions shifted by a finite, computable amount.

As long as scientists have known about solitary waves and had computers capable of simulating them, we have been colliding solitary waves together numerically and observing what comes out. An important example is the collision of so-called *kink* and *antikink* solutions to the $\phi 4$ equation,

The first studies in the 1970s gave a fleeting hint of rich structure and the first step toward understanding it. This was explored more deeply and given a name, the *two-bounce resonance*, in the 1980s. In the 1990s, more detailed numerical experiments revealed chaotic scattering. Finally, in the early 2000s, the mechanism behind this chaotic scattering was explained fully using techniques from dynamical systems. Figure 1 shows the speed at which a kink and antikink separate as a function their speed of approach, demonstrating a rich structure. Numerical simulations of partial differential equations (PDE) arising in diverse areas of physics have revealed the same phenomenon. In all that time, however, *this behavior has never been reported in a physical experiment in a real nonlinear wave system*.

This paper describes a simple experiment—a ball set rolling on a manufactured landscape—which mimics the dynamics of solitary wave collisions and reproduces some features of Figure 1. The main tool for understanding solitary wave interaction has been the derivation and analysis of *collective coordinate* (CC) models. The landscape is constructed such that the ball's equations of motion strongly resemble the ordinary differential equation system (ODE) governing the evolution of solution parameters described by a CC reduction describing solitary wave collisions in the $\phi 4$ equation.

Chaotic scattering is an appealing phenomenon to study because images such as Figure 1 cry out for explanation. The mathematical structure that has been developed is equally appealing, as is well-described by Ott.^{1,2} A scattering process is a physical phenomenon in which an object's trajectory begins and ends in free motion (i.e., with zero acceleration) but spends a finite time in a region where it is subject to forces. Such a process is called chaotic scattering if the final state depends in a complex, i.e., fractal or multi-fractal, manner on its initial state.

The paper is organized as follows. Section II discusses where the problem comes from and what approaches have been most useful in understanding it. It begins with a brief recap of the history of solitons, followed by an explanation of the difference between solitons and garden-variety solitary waves. It then describes the long history of chaotic scattering of solitary waves, the progress that has been made, and the methods that have been used to explain the observations. In Sec. III, we describe the laboratory experiment and the mathematical model that describes it, while in Sec. IV, we present the experimental results, which we interpret as the two-bounce resonance, the aspect of chaotic scattering that is most robust in the presence of dissipation. We end in Sec. V by reviewing the physical systems that have motivated many of earlier studies of the two-bounce/chaotic scattering phenomenon and expressing our hope that someone may perform this experiment to display this phenomenon in an extended physical system.

## II. HISTORICAL MOTIVATION

### A. Solitary waves and solitons

*Solitary waves* are solutions to a nonlinear evolutionary partial differential equation (PDE) that translate at a constant speed while maintaining their spatial profile. *Solitons* are solitary waves which exhibit unexpectedly simple dynamics upon collision. This is due to a hidden mathematical structure that was discovered beginning in the 1960's, first for the Korteweg-de Vries equation,^{3,4}

and later for many other equations, such as the sine-Gordon equation,

Such soliton equations can be solved exactly using a method called the inverse scattering transform (IST).^{5} Good brief histories of solitary waves and solitons can be found elsewhere, e.g., in encyclopedia articles by Scott,^{6} and by Zabusky and Porter.^{7}

The IST can be used to explicitly show that after colliding, solitons continue propagating with the same shape and speed, but with a phase shift and a time delay. Most nonlinear wave equations do not possess the IST, however, and solitary wave collisions in such systems display more complicated dynamics. Beginning in the 1970s, numerical simulations hinted this picture, but it took the development of more powerful computer hardware and numerical methods for the true complexity of this dynamics to fully emerge. The simultaneous development of *collective coordinate methods* would prove useful for untangling this complexity.

### B. Collective coordinates methods

The term “collective coordinates” refers to a variety of methods that are used to obtain simple ODE models that approximate the behavior of a PDE, at least as long as the solution stays in some small region of solution space. A very useful approach is the “variational method,” which applies to partial differential equations that arise as the Euler-Lagrange equations for a system with Lagrangian density

For example, solutions of the $\phi 4$ equation (1) minimize the action due to

over all *C*^{1} functions satisfying appropriate boundary conditions at infinity.

The variational method, due originally to Bondeson,^{8} works by constructing a solution ansatz whose spatial profile depends on a small number of time-dependent parameters and minimizing the action (3) with respect to these parameters. The resulting equation is a finite-dimensional Lagrangian system of ODE governing the parameters' evolution. Many examples and some generalizations are given in the review paper of Malomed.^{9}

### C. Numerical and analytical studies of solitary wave collisions

Solitary waves exist in the $\phi 4$ equation (1), which is non-integrable, despite its clear similarity to the sine-Gordon equation (2). These “kink” waves are given by

The “antikink” solution is just $\phi K\xaf(x,t,v)=\u2212\phi K(x,t;v)$.

Many papers have been written about kink-antikink ($KK\xaf$) collisions in this system. In 1975, Kudryavtsev^{10} found numerically that when *v* = 0.1, the $KK\xaf$ pair coalesce into a single bound state at the origin which oscillates irregularly and is slowly damped to zero. He explained this with a formal argument based on potential energy which can be considered a first step toward an explanatory collective coordinates model.

In 1979, Sugiyama^{11} performed additional experiments for a few values of *v* between 0.1 and 0.6. He found that those with speeds below a critical value *v*_{c} ≈ 0.25, the $KK\xaf$ pair merge into a localized bound state, while for *v* > *v*_{c}, the pair collides inelastically. Careful examination of his simulations showed that, upon collision, an oscillatory mode (sometimes called a “shape mode,” since its effect is to alter the shape of the kink) is excited, which removes energy from the translation component. If the initial kinetic energy is less than the amount lost to the oscillatory mode, the $KK\xaf$ pair cannot escape to infinity. If the *v* > *v*_{c}, they will escape to infinity but with reduced speed and with oscillations superimposed. He also derived the first CC model for this phenomenon, and used it in a formal calculation to determine *v*_{c} in agreement with his numerical observations.

The shape mode is an eigenfunction for the linearization of Eq. (1) about the kink solution (4), and is given by

and oscillates with frequency $\omega 1=32$. Sugiyama based his CC ODE for this system on the ansatz

It leads to a two degree-of-freedom Hamiltonian system for the evolution of *X*(*t*) and *A*(*t*). The ODE system that arises is somewhat complicated,^{11} but its fundamental features are captured by the simpler model system

The term *U*(*X*) is a potential energy describing the interaction of the two kinks, and was essentially derived by Kudryavtsev.^{10} For large *X*, both $F\u2032(X)$ and $U\u2032(X)$ vanish, so that the kink and antikink centers ±*X*(*t*) move at constant speed. The most important feature, for our purposes, is the phase diagram for *X* when *c* = 0, shown in Figure 2(b). The trajectories are level sets of the energy $E=12X\u03072+U(X)$, with *E* < 0 on bounded orbits, a separatrix orbit with *E* = 0, and unbounded orbits for *E* > 0. We denote the separatrix orbit by *X*_{S} (*t*), and note that it is an even function of time.

Figure 1 shows that this model qualitatively captures much of the dynamics of solitary wave collisions in (1). Observe, however, a few fundamental differences. While both systems conserve energy, in extended system (1), radiation (phonons) can carry energy away from the immediate vicinity of the kinks, effectively adding dissipation. This can be seen in the figure, where only solutions that escape after four or fewer collisions are plotted. In ODE (6), all solutions escape except for a set with measure zero, while for the PDE, a significant fraction of solutions are trapped, due to energy loss to radiation. The radiation also steals some energy from solutions that do escape, so that *v*_{out} < *v*_{in} even at the maxima of the resonance windows.

Also in 1979, Ablowitz, Kruskal, and Ladik reported on a series of numerical experiments for a few different nonlinear Klein-Gordon equations of the form $\phi tt\u2212\phi xx+f(\phi )=0$, including the $\phi 4$ equation and the sine-Gordon equation.^{12} Their numerical results were similar to Sugiyama's, but the paper ends with an observation which was to prove important: On certain velocity intervals below *v*_{c}, the kink and antikink eventually separate, instead of forming a bound state as seen by Kudryavtsev (from their References, it seems that they had not seen Sugiyama's work at this point). When the initial velocity is 0.3—greater than *v*_{c}, the final velocity is 0.135. However, when the initial velocity is 0.2—less than *v*_{c}, the final velocity is 0.155. They conclude: “The reason for the apparent “resonance” between these interacting aperiodic waves and the radiation is not yet fully understood.”

In 1983, Campbell *et al.*^{13} performed a more systematic numerical sweep of initial velocities and found significantly more detailed structure. Their calculation revealed the black curve and the leftmost nine blue curves in Figure 1. The black curve shows the final velocities of all the solutions that escape after exactly one collision, i.e., those with initial speed above *v*_{c}, which they estimate numerically to be 0.26. The blue curve shows the final velocities of those with exactly two collisions. They called this phenomenon the *two-bounce resonance*.

The first question is: “What is the difference between solutions in one two-bounce window and those in the next?” A good way to understand the solutions is to fit the numerical solution to the ansatz (5), and plot the approximate values of *X*(*t*) and *A*(*t*). This is done in Figure 3 for five increasing initial velocities. Subfigure (a) shows a collision leading to capture; subfigures (b) and (d) show solutions from the first two two-bounce windows; (c) shows a solution from a three-bounce window; and (e) shows a one bounce solution with *v*_{in} > *v*_{c}. The number of oscillations of *A*(*t*) is the same for all initial velocities in a given window and increases by one from one window to the next. Let *T*_{shape} be the period of the shape oscillation and *T _{n}* be the time interval between the two collisions for the

*v*

_{in}in the

*n*th window at which

*v*

_{out}is maximized, with

*T*

_{2}and

*T*

_{3}indicated in subfigures (b) and (d). We remark that the nonlinear projection (5) is singular when

*X*(

*t*) = 0, because the coefficient of

*A*(

*t*) then vanishes,

^{14}in which case the fitting algorithm fails.

Treating the numerical results as laboratory data, Campbell *et al.* reasoned there should be a relation of the form

which they confirm by least-squares fitting the data to a line. They conjectured that at the first collision, the kink loses energy to the shape-mode oscillation and that, if the timing is right on the second collision, the shape mode returns enough energy to the propagating mode to allow escape. What remained was to explain the mechanism.

A few years later, Anninos *et al.*^{15} showed, via more detailed experiments, the existence of higher bounce windows arranged in a fractal structure. That is, they found chaotic scattering, in the solitary wave collisions, including now the three and four bounce windows of Figure 1.

Over the next two decades, these numerical experiments were reproduced for many other nonlinear wave systems, and the reasoning of Campbell *et al.* was applied to these systems as well. Campbell, as well as Peyrard, looked at a variety of $\phi 4$-like equations. They demonstrated that, at least for the systems they considered, the resonance phenomenon was only present when the linearization around the kink possesses an internal mode.^{16–18} Of particular interest is a study by Fei *et al.*^{19} which considers the sine-Gordon equation perturbed by a small localized defect,

They find that kink solutions impinging on the defect show a behavior very similar to that seen in Figure 1. They derive a CC system similar to system (6) where *A*(*t*) now measures the amplitude of a mode localized in neighborhood of the defect. They use conservation of energy to derive an implicit formula for *v*_{c}(*ϵ*) and apply the reasoning of Campbell *et al.* to fit the time between collisions to the period of the secondary oscillator.

### D. Reduction of collective coordinates to an iterated map

Beginning in 2004, Goodman and Haberman published a series of papers explaining the chaotic scattering in detail by reducing the CC ODE system to a discrete-time iterated map, called a separatrix map or scattering map.^{20,21} This was first done^{22} for the system studied by Fei,^{19} because the the CC ODE system for equation (8) has a small parameter *ϵ* that can be used in a perturbation analysis. Explicit formulas were found that approximate the critical velocities and the two- and three-bounce resonance, eliminating the need for data fitting as in equation (7). Over subsequent papers,^{23–25} the derivation of the map was streamlined and applied to other solitary waves systems and the model system (6). It was subsequently put in a more explicit form.^{26} We outline this last approach here.

Figure 4 depicts the results of one numerical simulation of the initial value problem (6) with the solitons initially far apart, propagating toward each other, and with the shape mode unexcited, i.e., *X* ≫ 1, $X\u0307<0$, and $A=A\u0307=0$. Additionally, $X\u0307$ is small enough for capture to occur. In this simulation, *E*(*t*) > 0 before the first collision, this computed solution has energy *E* > 0 so that the solution begins outside the separatrix in Figure 2(b). It crosses to the inside of the separatrix at the first collision time *t*_{1}. At each subsequent collision, *E*(*t*) jumps, reaching a plateau *E _{j}* between collisions at

*t*

_{j}_{−1}and

*t*, and escaping to infinity when

_{j}*E*(

*t*) > 0 once again. Upon each collision, the amplitude and phase of

*A*(

*t*) jump. We represent the solution before the collision at time

*t*by the energy level

_{j}*E*and by assuming

_{j}for *t* between *t _{j}*

_{−1}and

*t*, where setting the constant

_{j}allows us to scale the variables to be *O*(1) in the final form of the system. We seek a map (*E _{j}*

_{+1},

*c*

_{j}_{+1},

*s*

_{j}_{+1}) =

*M*(

*E*,

_{j}*c*,

_{j}*s*).

_{j}The map is constructed by building a matched asymptotic expansion which alternates between “outer expansions” where *X*(*t*) is approximated by *X*_{S} (*t* − *t _{j}*), Figure 2(b), and “inner expansions.” On the outer solution, valid when

*X*(

*t*) ≫ 1, the modes exchange energy. We find, by variation of parameters, that outer solutions with “before-collision” condition (9) as

*t*−

*t*→ −

_{j}*∞*satisfy the “after-collision” condition

and that

This change of energy is computed using a *Melnikov integral*.^{20} Condition (10) is written in terms of (*t* − *t _{j}*), whereas the form of the map requires (

*t*−

*t*

_{j}_{+1}). We can compute the time between collisions by matching between two the outer expansions connected by an inner expansion, which gives $tj+1\u2212tj=\alpha Ej+1\u22121/2+o(1)$ for some

*α*determined by matched asymptotics. This gives, approximately,

where *θ _{j}*

_{+1}=

*ω*(

*t*

_{j}_{+1}−

*t*).

_{j}The discrete map possesses a conserved quantity $H=Ej+12\omega 2C2(cj2+sj2)$ that allows us to eliminate *E _{j}* from the system. Defining a complex variable

*z*=

_{j}*c*+

_{j}*is*, the new map may be written

_{j}Figure 1 shows that the map reproduces the qualitative and many quantitative features of the ODE system and enables the determination of many of the features of the dynamics of system (6). For example, in the experiments depicted in Figure 1, the initial condition is *z _{j}* = 0. Therefore, we may estimate the critical velocity: the solution is captured if

*E*

_{1}< 0, i.e., if $v<vc=\omega C$, making use of Eq. (11). The analysis also allows us to calculate approximate formulas for the centers of the two- and three bounce resonance windows, allowing the computation of the constants in Eq. (7) without data fitting.

^{26}

## III. A PHYSICAL MODEL OF THE COLLECTIVE COORDINATE EQUATIONS

We wish to design a surface such that a ball confined to roll along that surface satisfies equations of motion similar to system (6). We derive the evolution equations corresponding for motion along a general surface and then design a surface that gives the desired dynamics.

### A. The mathematical model

Consider the behavior of a point particle of mass *m* confined to a surface *z* = *h*(*x*, *y*) and moving under the influence of constant gravity in the *z*-direction, and ignoring any friction or other dissipative mechanisms. We do not model in detail the rotation of the ball, instead noting that the rotational kinetic energy of a sphere rolling with speed *v* is $25mv2$. Adding this to the translational component gives a total kinetic energy

The gravitational potential energy is just *U* = *mgz* = *mgh*(*x*, *y*). The evolution of a Lagrangian system with coordinates *q _{j}* is governed by the Euler-Lagrange equations

For this system, this gives, after some algebra,

where $g\u0303=57g$. Under the assumption that the height *h*(*x*, *y*) varies slowly (e.g., *h*(*x*, *y*) = *h*(*δx*, *δy*), *δ* ≪ 1), this is approximately

Thus, when

Eq. (14) is identical to Eq. (6). A contour plot of such a surface is shown in Figure 5(a). Numerical simulations of both Eqs. (13) and (14) show qualitatively the same chaotic scattering behavior (not shown here).

We denote by *x* and *y*, respectively, the longitudinal and transverse, directions. To interpret the ODE, consider the case *ϵ* = 0, in which case the subspace $y=y\u0307=0$ is invariant. A ball that begins rolling down the center of the surface with *y* = 0 stays along that line. When *ϵ* > 0, this symmetry is broken, deflecting the longitudinal trajectory in the transverse direction. We will say that a “bounce” occurs when the ball reaches a minimum in the longitudinal direction.

### B. Effects of dissipation

To account for experimental energy loss, we include a viscous damping force $F\u2192damp=\u2212\mu (x\u0307,y\u0307)$, which modifies the equations of motion to:

In Goodman *et al.*,^{27} a dissipative correction to system (6) is derived to account for the effect of radiative damping in nonlinear wave collisions. It adds a term like $F(X)2A2A\u0307$ to the left-hand side of Eq. (6b). This damping term is strongly localized in *X*, only applying when the kink and antikink are close together, and nonlinear in *A* so that small radiation damps very slowly. Both types of dissipation alter the fractal structure displayed by the ODE system, as seen in Figure 6, but their effect is rather different. The localized damping preserves much more of the fine structure of the chaotic scattering, whereas the viscous damping destroys almost everything but the two-bounce resonance windows. Their behavior for *v* near *v*_{c} is also quite different, with viscous damping reducing the final speed of solutions in the two-bounce solutions much more than the localized damping. *We therefore expect to see mostly the two-bounce resonance phenomenon in our physical experiments, rather than the full chaotic scattering picture.*

Numerical experiments were used to both verify that system (13) could reproduce the dynamics demonstrated by solitary wave collisions. These produced the expected results, so we proceeded with laboratory experiments.

## IV. THE LABORATORY EXPERIMENT

To test experimentally whether the physical system modeled by (13) displays chaotic scattering, we needed to fabricate a surface satisfying (15). The surface was milled out of high-density urethane foam using a three-axis mill in the Fabrication Laboratory of the NJIT Department of Architecture. Dimensions of the surface are given in the Appendix. The rough milled surface was sanded and painted, as shown in Figure 5(b). A ramp was placed at one end, and a rubber coated steel ball (from a computer mouse) was rolled down this ramp and allowed to move along the surface until either (1) it became clear that it would not escape, or (2) the ball returned close enough to the starting point that we deemed it to have escaped.

We computed an effective friction constant by the following procedure. The ball was placed at the edge of the channel near the left edge of Figure 5(b) and allowed to move under the influence of gravity. Coupling to the *x*-direction is weak here, so the motion remains largely confined to the *y*-direction over the time scale of observation. From the video, we find the sequence of successive maxima $| y(tk) |$ from which we can find the frequency and decay rate. Nondimensionalizing using these scales gives *y* motion satisfying the nondimensional equation

This nondimensional system has a time unit of 0.15 seconds. The video trials shown below last on the order of 5–10 s, which is 30–60 nondimensional time units. Over this time, dissipation may be significant.

In the initial experiments, we used a large marble, which lost energy rapidly. The rubber-coated steel ball, being heavier and quieter, kept rolling for much longer. Dissipation was not negligible, however, so any trial in which the ball did not exit after three “bounces” was considered to be trapped.

The motion was recorded at 125 frames/second with a Photron FASTCAM 1024 PCI high-speed camera, positioned above the surface and pointing downward. Its high frame rate was not necessary, but its short exposures prevented the ball's image from being smeared into a snake-like shape. The videos were analyzed using the MATLAB Image Processing Toolbox, giving a time series for its *x* and *y* coordinates. A screenshot of the video and the running program are shown in Figure 7 (Multimedia view).

The initial velocity was varied by varying the height of the ball's release. Due to the size of our foam block, the surface could only be milled to a fairly shallow depth. This meant that a ball with sufficiently large initial velocity would simply fly off the milled portion of the surface and, unfortunately, left us unable to determine the critical velocity *v*_{c} experimentally. Figure 8 shows several interesting trajectories: First, five two-bounce trajectories with the number of transverse oscillations increasing from two in subfigure (a) to six in subfigure (e). Subfigure (f) shows a fit of the time between the two bounces vs. the number of transverse oscillations, showing the same approximate linear fit as given in Eq. (7), in this case, *T _{n}* ≈ 0.961

*n*+ 0.755. Unfortunately, the initial velocities that led to these solutions do not increase monotonically as predicted by the theory, due, perhaps, to our inability to control with sufficient precision the transverse component of the initial location and velocity vectors. The figure also shows (g) a trapped trajectory and (h) a three-bounce trajectory. These represent, as far as we knox, the first (non-numerical) experimental evidence of the two-bounce resonance phenomenon, and a trace of chaotic scattering.

## V. CONCLUSIONS AND PLEA

We have designed and built an experimental system that is approximately governed by the same reduced system that gives rise to the two-bounce resonance phenomenon in solitary wave collisions. Experiments run with this apparatus qualitatively reproduce the dynamics seen with solitary waves.

Nonetheless, we remain hopeful that someone will find a way to demonstrate this phenomenon in a laboratory setting with actual solitary waves. The basic requirement is a solitary wave that supports an additional degree of freedom to which it can transfer energy, and which can transfer energy back. Usually, this takes the form of a mode that is localized near a stationary potential^{19,22} or an internal mode of oscillation that moves with the wave,^{13,17} although this is not always necessary.^{28–30}

Combing through some of the literature on this phenomenon, we find that many papers discuss the physical systems described by their equations and mention that—perhaps—this phenomenon may be found in these systems. Experimental verification for some of these systems is clearly impossible: Anninos *et al.*^{15} describe the $\phi 4$ equation as describing the interactions of large-scale domain walls in the universe. They admit, “Because the possibility of head-on collisions is small in the real Universe, our findings are not expected to have a profound effect cosmologically.” At the other extreme of scales, Kudryavtsev describes the $\phi 4$ system as a model for the Higgs field.^{10,31} Other unpromising applications for the theory are information transport in brain microtubules,^{32} and various quantum field theories.^{32,33}

Other applications are perhaps closer to being experimentally realizable. In their original paper on the subject, Campbell *et al.* suggest excitations in polymeric chains and phase transitions in uniaxial ferroelectrics.^{13} Tan and Yang study the phenomenon in vector solitons in optical fiber.^{23,28,29} Numerical experiments of Fei *et al.* potentially describe the scattering of solitons off a defect in a long Josephson junction.^{19,22} Finally, Forinash *et al.* mention the denaturation of DNA.^{34} This list is far from exhaustive, and we encourage readers to think if an experiment is possible in their favorite system.

## ACKNOWLEDGMENTS

The experiment was the focus of the Spring 2010 NJIT “capstone” course in applied mathematics, a class in which students perform laboratory experiments, learn the math needed to model them, and apply analytical and numerical methods to that model to explain the results. The first author was the instructor, and the others were undergraduates enrolled in the course, who continued working on the project the following summer. The NSF Capstone Laboratory in the NJIT Department of Mathematical Sciences was supported by NSF DMS-0511514. Thanks to Richard Haberman, Shane Ross, and Lawrie Virgin for useful discussions and to Richard Garber and Gene Dassing in the FabLab at the New Jersey School of Architecture, NJIT, for fabricating the surface and donating the material. Special thanks to David Campbell for early discussions about the form of this article.

### APPENDIX: SPECIFICATION OF THE EXPERIMENTAL APPARATUS

To design the experimental surface, we performed numerical experiments of system (14), varying not just the parameters in the potential (6) but also the formulas for the potentials *U*(*x*) and *F*(*x*), trying to ensure we could observe interesting dynamics using the materials available to us. The surface chosen was

where the units of distance in all coordinate directions is centimeters. The parameters chosen were

## References

^{4}model

^{4}theory

^{4}models: Problems in the variational approach

^{4}equation: The n-bounce resonance and the separatrix map

^{6}model

^{4}field theory