## Abstract

The dynamic evolution of linearly dispersive waves on periodic domains with discontinuous initial profiles is shown to depend remarkedly upon the asymptotics of the dispersion relation at large wavenumbers. Asymptotically linear or sublinear dispersion relations produce slowly changing waves, while those with polynomial growth exhibit dispersive quantization, a.k.a. the Talbot effect, being (approximately) quantized at rational times, but a non-differentiable fractal at irrational times. Numerical experiments suggest that such effects persist into the nonlinear regime, for both integrable and non-integrable systems. Implications for the successful modelling of wave phenomena on bounded domains and numerical challenges are discussed.

## 1. Introduction

The most common device used to model physical wave phenomena is to approximate the fully nonlinear fluid mechanics problem in a targeted physical regime by first suitably rescaling the original variables and then expanding the rescaled system in powers of an appropriate small parameter. Truncation at low (typically first) order produces a mathematically simpler system, whose solutions, one hopes, closely track the originals, at least over an initial time span.

At the linear level (at least in the homogeneous, isotropic regime), conservative wave dynamics is entirely determined by the dispersion relation, which relates the temporal frequency of an oscillatory complex exponential solution to its wavenumber (spatial frequency). The standard modelling paradigm can thus, at the linear level, be viewed as replacing the full dispersion relation by a suitable approximation. However, in most perturbative models, while the approximate dispersion relation is designed to fit its physical counterpart well at small wavenumbers, e.g. they have the same Taylor polynomial at the origin, they often have little in common at the high-frequency range. For example, the free boundary problem governing water waves has a dispersion relation that is asymptotically proportional to a square root of the wavenumber. On the other hand, in the shallow water regime, the unidirectional Korteweg–de Vries model's dispersion relation is cubic, and of the wrong sign, and the competing regularized long wave (RLW) or Benjamin–Bona–Mahony (BBM) model's is asymptotically zero, while that of some standard bidirectional Boussinesq systems is asymptotically constant. Provided one restricts attention to smooth solutions, or to unbounded domains, the discrepancies in the dispersion of high-frequency modes do not seem to play a very noticeable role. However, on bounded domains, rougher solutions, which have more substantial high-frequency components, are affected in a significant manner owing to subtle and as yet poorly understood resonances among the varyingly propagating Fourier modes.

In Olver [1], the second author observed that, for integral polynomial dispersion relations, i.e. those whose coefficients are integer multiples of a common real number *λ*, the solution in question has the remarkable behaviour that, at irrational times (relative to the ratio *β*=ℓ/*λ* between the length ℓ of the interval and the dispersion constant *λ*), it is a continuous, but fractal, non-differentiable function, whereas at rational times, the solution is piecewise constant! While previously unnoticed^{1} by experts in dispersive wave theory, the phenomenon of dispersive quantization had, in fact, been first discovered in the early 1990s, in the contexts of optics and quantum mechanics, by Berry and various co-workers [2–6], who named it the *Talbot effect* in honour of a striking 1836 optical experiment by the photographic pioneer Talbot [7]. Experimental confirmations of the Talbot effect in optics and in quantum revival are described in Berry *et al.* [5]. Rigorous analytical results and estimates justifying the Talbot effect can be found in the work of Kapitanski & Rodnianski [8], Rodnianski [9], Oskolkov [10,11] and Taylor [12].

The initial aim of the present paper was to extend these studies to basic linear differential and integro-differential equations, depending on a single spatial variable, which have non-polynomial dispersion relations. Our main conclusion is that the large wavenumber asymptotics of the dispersion relation plays the dominant role governing the qualitative features exhibited by their rough solutions on periodic domains. In analogy with the basic Riemann problem of gas dynamics [13], we focus on the specific initial data provided by a step function. The resulting Fourier series solution is easily constructed, and of independent interest. For example, when the dispersion relation is polynomial, the solution at rational times has the form of a Gauss or Weyl exponential sum, of fundamental importance in analytic number theory since the work of Hardy & Littlewood [14] and Vinogradov [15]. Indeed, as explained in Ford [16], because the behaviour of such exponential sums has implications for the size of the zero-free region of the Riemann zeta function in the critical strip, further refinement of known estimates could provide a route to proving the Riemann hypothesis!

We will also initiate the investigation of such effects in the nonlinear regime. The revolutionary discovery of the *soliton* was sparked by the original movie of Zabusky & Kruskal [17], which displayed a numerical simulation of the solution to a particular periodic initial-boundary value problem for the Korteweg–de Vries equation with small nonlinearity. (By contrast, the celebrated studies of Lax & Levermore [18], Lax *et al.* [19], Venakides [20,21] are concerned with the small dispersion regime and convergence to shock wave solutions of the limiting nonlinear transport equation.) Because Zabusky and Kruskal's selected initial data were a smoothly varying cosine, the Talbot fractalization effect was (perhaps fortunately!) not observed. (And, technically, the elastically interacting waves that emerge from the initial cosine profile are not true solitons, in that these exist only for the full line problem, but, rather, cnoidal waves embedded in a hyperelliptic finite gap solution [22,23].) Our numerical experiments with step function initial data strongly indicate that the dispersive quantization/fractalization effect is present in the periodic boundary value problems for both the integrable nonlinear Korteweg–de Vries and Schrödinger equations, as well as non-integrable models of a similar nature, and extends well beyond the weakly nonlinear regime.

One consequence of these studies is that, contrary to the conventional wisdom, when dealing with nonlinear wave models on bounded intervals, the principal source of analytic difficulty may be, counterintuitively, not the nonlinear terms, but rather the poorly understood behaviour of linearly dispersive partial differential equations. Our investigations imply that the qualitative behaviour of the solution to the periodic problem depends crucially on the asymptotic behaviour of the dispersion relation at large wavenumbers, reinforcing Benjamin *et al.*'s critique [24] that the Korteweg–de Vries equation, say, is an unsatisfactory model for surface waves, because its cubic dispersion relation effectively transmits the high-frequency modes in the wrong direction, with unboundedly negative phase velocity and group velocity, inciting unphysical interactions with other solution components. And indeed, in the periodic problem, this shortcoming is observed as the apparently number-theoretic resonant interaction of high-frequency modes spawned by the initial data serve to produce a-physical fractalization and quantization effects in the engendered solution.

## 2. Dispersion relations and wave models

Consider a linear, scalar, constant coefficient partial differential equation for a function *u*(*t*,*x*) of time *t* and a single spatial variable *x*. Recall [13] that the *dispersion relation* *ω*(*k*) relates temporal frequency *ω* to wavenumber (spatial frequency) *k* of an oscillatory complex exponential solution of the form
2.1The differential equation will be called *purely dispersive* if the resulting function *ω*(*k*) is real. (Complex dispersion relations indicate the presence of dissipative or viscous effects that damp out solutions, or, alternatively, ill-posedness, depending on the sign of their imaginary part.) Given a real dispersion relation, the *phase velocity* *c*_{p}=*ω*(*k*)/*k* prescribes the speed of an individual oscillatory wave, while, as a consequence of the method of stationary phase, wave packets and energy move with the *group velocity* *c*_{g}=d*ω*/d*k* [13]. These differ when the dispersion relation is nonlinear, causing initially localized disturbances to spread out in a dispersive fashion—a physical effect that can be seen, for example, by throwing a rock into a pond and observing that the individual waves move at a different speed, namely the phase velocity, than the overall disturbance, which moves at the group velocity. The dispersion relation of nonlinear equations is that of their linearization.

For example, the basic *transport equation*
2.2has linear dispersion relation *ω*(*k*)=*ck*, whereas the one-dimensional *wave equation*
2.3has bidirectional dispersion relation *ω*(*k*)=±*ck*. Linearity of these dispersion relations has the consequence that all Fourier modes travel with a common speed, thereby producing travelling waves of unchanging form. The simplest quadratic dispersion relation, *ω*(*k*)=*k*^{2}, appears in the free space *Schrödinger equation*
2.4which is fundamental to quantum mechanics [25]. The linear *beam equation*
2.5which models small vibrations of a thin elastic beam, and has bidirectional quadratic dispersion: *ω*(*k*)=±*k*^{2}. The third-order linear wave model
2.6which is sometimes known as the *Airy partial differential equation*, because it admits solutions in terms of Airy functions [26] or, alternatively, the *linear Korteweg–de Vries equation*, has cubic dispersion relation *ω*(*k*)=*k*^{3}.

A fertile source of interesting wave models is the free boundary problem for water waves, meaning the dynamics of an incompressible, irrotational fluid under gravity. Let us review, following [13,27,28], the derivation of models in the shallow water regime. For simplicity, we restrict our attention to two-dimensional flow, with *x* indicating the horizontal coordinate and *y* the vertical. (See Yuen & Lake [29] for the deep water regime, and [27,30] for extensions to three-dimensional models.)

Suppose that the water sits on a flat, impermeable horizontal bottom, which we take as *y*=0, and has undisturbed depth *h*. We assume further that the surface elevation at time *t* is the graph of a single-valued function *y*=*h*+*η*(*t*,*x*), i.e. there is no wave overturning. The periodic problem of interest here takes *η* to be a periodic function of *x* with mean 0, while, in contrast, solitary waves require *η*(*t*,*x*)→0 as . The motion of the fluid is prescribed by the velocity potential *ϕ*(*t*,*x*,*y*), which is defined within the domain occupied by the fluid *D*_{η}={0<*y*<*h*+*η*(*t*,*x*)} (figure 1). In the absence of surface tension, the free boundary problem governing the dynamical evolution of the water is then given by
2.7in which *g* is the gravitational constant. As shown in Whitham [13], §13.4, the dispersion relation for the linearized water wave system, obtained by omitting the nonlinear terms from the free surface conditions, is
2.8

In the shallow water regime, one introduces the rescaling
2.9in which ℓ represents a typical wavelength, *h* the undisturbed depth, the leading order wave speed and *a* the wave amplitude. Substituting (2.9) into the free boundary value problem (2.7), one formally solves the resulting rescaled Laplace equation for the potential *ϕ*(*t*,*x*,*y*), as a power series in the vertical coordinate *y*, in terms of the horizontal velocity *u*(*t*,*x*)=*ϕ*_{x}(*t*,*x*,*θ*) at a fraction 0≤*θ*≤1 of the undisturbed depth *h*. Substituting the result into the free boundary conditions results in a system of equations of Boussinesq type, which is then truncated to, say, first order in the small parameters
2.10representing the ratio of undisturbed depth to height of the wave, and the square of the ratio of height to wavelength. We are interested in the regime of long waves in shallow water, in which both *α* and *β* are small, but of comparable magnitudes. As in [27,28,31], the resulting *Boussinesq systems* are all of the general form
2.11where the parameters *a*,*b*,*c*,*d*, depend on the depth, and are subject to the physical constraints
2.12The particular system
2.13valid for *θ*=1, i.e. on the free surface, was derived by Whitham [32] and Broer [33], and shown to be completely integrable (in fact, tri-Hamiltonian) by Kaup [34] and Kupershmidt [35]. Bona *et al.* [36] subsequently found that this system is, in fact, linearly ill-posed—although the question of nonlinear ill-posedness is, apparently, still open. The corresponding second-order models can found in earlier studies [27,28].

The second stage is to restrict the bidirectional Boussinesq system (2.12) to a submanifold of approximately unidirectional solutions, building on the usual factorization of the second-order wave equation (2.3) into right- and left-moving waves. Substituting the right-moving ansatz
2.14and truncating to first order reduces the Boussinesq system to the *Korteweg–de Vries equation*
2.15for unidirectional propagation of long waves in shallow water. In the long wave regime, *β*≪1, the nonlinear effect is negligible, and the equation reduces to the linear equation
2.16which, in turn, can be mapped to the third-order Airy equation (2.6) by passing to a moving coordinate frame and then rescaling. Alternatively, because to leading order *u*_{t}≈−*u*_{x}, one can replace the third derivative term, *u*_{xxx}≈−*u*_{xxt}, leading to the *RLW/BBM model* [24,37],
2.17Observe that both first-order models (2.15) and (2.17) are independent of the depth parameter *θ*, which appears only in the second- and higher-order terms in *α*,*β*. Inverting the unidirectional constraint (2.14) leads to the same Korteweg–de Vries equation and RLW/BBM equations for the surface elevation *η*. Again, this is only true up to first order, and the higher-order terms are slightly different [28]. Because the preceding unidirectional models do not admit the exact water wave dispersion relation (2.8), Whitham [32] proposed the alternative model
2.18in which the Fourier integral operator
exactly reproduces the rescaled water wave dispersion relation given by (2.8) when *g*=*h*=1. However, Benjamin *et al.* [24] argue that the linear dispersion in Whitham's model is insufficiently strong to counteract the tendency to shock formation.

Let us also list a few important integrable nonlinear evolution equations, from a variety of physical sources. The *nonlinear Schrödinger equation*
2.19arises in nonlinear optics as well as in the modulation theory for water waves [29,38]. Its linearization has quadratic dispersion relation *ω*(*k*)=*k*^{2}. Integrability was established by Zakharov & Shabat [39], and the fact that it admits localized soliton solutions that preserve their shape under collision plays an important role in transmission of signals through optical fibres [40]. The integrable *Benjamin–Ono equation*
2.20in which is the *Hilbert transform*, given by the Cauchy principal value integral
2.21arises as a model for internal waves [41]. Its linearization has quadratic dispersion relation *ω*(*k*)=*k*^{2} sign *k*, which is the odd counterpart of the quadratic Schrödinger dispersion.

Another integrable model is the *Boussinesq equation*
2.22also derived by Boussinesq from the water wave problem. While the Boussinesq equation admits waves propagating in both directions, it is not, in fact, equivalent to any of the bidirectional Boussinesq systems (2.11), and, indeed, its derivation from the water wave problem relies on the unidirectional hypothesis. The plus sign is the ‘bad Boussinesq’, which is unstable, whereas the minus sign is the stable ‘good Boussinesq’. However, both versions admit solutions that blow up in finite time; see McKean [42] for a detailed analysis of solutions to the periodic problem, and also [43,44] for further results on stability and blow-up of solutions. An alternative, regularized version
2.23was subsequently investigated by Whitham [13]. The integrable Boussinesq equation (2.23) and its regularization (2.2) were proposed as model for DNA dynamics by Muto [45] and Scott [46]. The associated periodic boundary value problems can thus be viewed as a model for the dynamics of DNA loops, whereas individual DNA strands require suitably adapted boundary conditions at each end.

In table 1 below, we summarize dispersion relations arising from water waves (in units so that *g*=*h*=1) and some of the more important shallow water models, rescaled so that *α*=*β*=1. In the first three cases, taking the positive square root corresponds to right-moving waves. The final column lists their leading-order asymptotics (omitting constant multiples) at large wavenumber.

## 3. Linear dispersion in periodic domains

In general, let *L* be a scalar, constant coefficient (integro-)differential operator with purely imaginary Fourier transform (FT) . The associated scalar evolution equation
3.1has real dispersion relation .

We subject (3.1) to periodic boundary conditions on the interval −*π*≤*x*≤*π*, (or, equivalently, work on the unit circle: *x*∈*S*^{1}) with initial conditions *u*(0,*x*)=*f*(*x*). To construct the solution, we expand the initial condition in a Fourier series:
3.2The solution to the periodic initial-boundary value problem then has time-dependent Fourier series
3.3

In particular, the *fundamental solution* *u*=*F*(*t*,*x*) takes a delta function as its initial data, *u*(0,*x*)=*δ*(*x*), and hence has the following Fourier expansion:
3.4The solution (3.3) can then be rewritten as a convolution integral with the fundamental solution:
3.5

The following result underlies the dispersive quantization effect for equations with ‘integral polynomial’ dispersion relations.

### Definition 3.1

A polynomial *P*(*k*)=*c*_{m}*k*^{m}+⋯+*c*_{1}*k*+*c*_{0} will be called *integral* if all its coefficients are integers: , *j*=0,…,*m*.

### Theorem 3.2

*Suppose that the dispersion relation of the evolution equation (*3.1*) is a multiple of an integral polynomial: ω(k)=λP(k), for some* *. Let β=2π/λ. Then, at every rational time t=βp/q, with p and* *the fundamental solution (*3.4*) is a linear combination of q periodically extended delta functions concentrated at the rational nodes x*_{j}*=2πj/q for* *.*

As in Olver [1], this result is an immediate consequence of the following well-known lemma:

### Lemma 3.3

*Let* . *The coefficients of a complex Fourier series* (3.2) *are* *q*-*periodic in their indices; so* *c*_{k+q}=*c*_{k} *for all* *k*, *if and only if the series represents a linear combination of the* *q* periodically *extended delta functions concentrated at the rational nodes* *x*_{j}=2*πj*/*q* *for*
*for certain* .

Applying theorem 3.2 to the superposition formula (3.5) produces the following intriguing corollary.

### Corollary 3.4

*Under the assumptions of theorem 3.2, at a rational time* *t*=*βp*/*q*, *the solution profile to the periodic initial-boundary value problem is a linear combination of* ≤*q* *translates of its initial data* *u*(0,*x*)=*f*(*x*), *i.e*.

In the case of the ‘Riemann problem’, the initial data are the unit step function:
3.6Under the assumption that the dispersion relation is odd, *ω*(−*k*)=−*ω*(*k*), the Fourier expansion of the resulting solution is
3.7whose spatial derivative is the difference of two translates of the fundamental solution, namely ∂*u*^{☆}/∂*x*=*F*(*t*,*x*)−*F*(*t*,*x*−*π*). Non-odd dispersion relations produce complex solutions, whose real part is given by the more complicated expression
3.8

In the following figures, we display graphs of the solution (3.7) at some representative times, for various *odd* dispersion relations, most of which are associated with water waves models. The figures are ordered by their asymptotic order *α* at large wavenumber: *ω*(*k*)=*O*(|*k*|^{α}) as . We did not use numerical integration to obtain these plots; rather, they were produced in Mathematica by summing the first 150 terms in the Fourier series (3.7), which was adequate to capture the essential detail. (In all the cases we looked at, summing the first 1000 terms does not appreciably change the solution graphs.) In some plots, a tiny residual Gibbs effect can be seen at the jump discontinuities. By comparing the various profiles, we conclude that the qualitative behaviour of the solution (3.7) depends crucially on the asymptotic behaviour of the dispersion relation *ω*(*k*) for large wavenumber. However, a full understanding of the various qualitative behaviours cannot be gleaned from these few static solution profiles, and so the reader is strongly encouraged to watch the corresponding movies posted at http://www.math.umn.edu/∼olver/mvdql.html. Each movie shows the dynamical behaviour of the solutions (3.7) for various dispersion relations over a range of times; the associated figures at the end of the paper are selected frames extracted from the corresponding movie. In some cases, the solution evolves at a glacially slow pace; so a representative collection of sub-movies are posted. Again, we emphasize that all dispersion relations used in the figures and movies are odd functions of *k*. More general dispersion relations, which produce complex oscillatory solutions whose real part is given by the more complicated formula (3.8), exhibit additional qualitative features that we will investigate more thoroughly elsewhere.

Here is an attempt to convey in words the behaviours that can observed in the movies, but are sometimes less evident in the plots. In all cases considered, the large wavenumber asymptotics is fixed by a certain power of the wavenumber
3.9In our experiments, the overall qualitative behaviour seems to be entirely determined by the asymptotic exponent *α*, although the intricate details of the solution profile do depend upon the particularities of the dispersion relation. As *α* varies, there is a continual change from one type of qualitative behaviour to the next. Over the range , there appear to be five main regions, which are roughly demarcated as follows.

—

*α*≤0.*Large-scale oscillations*: the two initial discontinuities, at 0 and*π*, remain stationary and begin to produce oscillating waves moving to the left. After this initial phase, the solution settles into what is essentially a stationary up- and down-motion, with very gradually accumulating waviness superimposed. Two examples are the dispersion relations associated with the RLW/BBM equation (2.17), with*α*=−1, shown in figure 2, and the regularized Boussinesq equation (2.23) as well as some versions of the Boussinesq system (2.11), for which*α*=0, shown in figure 3. Interestingly, the former becomes wavier sooner, whereas in the latter, even after time*t*=1000, the smaller scale oscillations are still rather coarse, and numerical limitations prevented us from following it far enough to see very high-frequency oscillations appear.— 0<

*α*<1.*Dispersive oscillations*: each of the two initial discontinuities changes into a slowly moving ‘seed’ that continually spawns a right-moving oscillatory wave train. After a while, the wave train emanating from one seed interacts with the other, and the seeds follow along these increasingly rapid oscillations. Much later, the solution has assumed a slightly fractal wave form superimposed over a slowly oscillating ocean, similar to ripples on a swelling sea that moves up and down while gradually changing form. Examples include the full water wave dispersion relation, plotted in figure 4, and the square root dispersion graphed in figure 5. Observe that both have the same large wavenumber asymptotics, and the corresponding solutions have very similar qualitative behaviour, while differing in their fine details.—

*α*=1.*Slowly varying travelling waves*: as*α*approaches 1, the oscillatory waves acquire an increasingly noticeable motion to the right. Indeed, if the dispersion relation is exactly linear,*ω*(*k*)=*k*, the solution is a pure travelling wave*u*(*t*,*x*)=*f*(*x*−*t*), which moves unchanged with unit speed. Having asymptotically, but not exactly linear, dispersion relation superimposes slowly varying oscillations on top of the travelling wave. An example appears in figure 6, based on the asymptotically linear dispersion relation associated with, for example, the equation^{2}3.10— 1<

*α*<2.*Oscillatory becoming fractal*: in this range, each of the initial discontinuities produces an oscillatory wave train propagating to the right. After the wave train encounters the other discontinuity, the oscillations become increasingly rapid over the entire interval, and, after a while, the entire solution acquires a fractal profile, which displays an overall motion to the right while its small scale features vary rapidly and seemingly chaotically. An example is provided in figure 7, based on the dispersion relation*ω*(*k*)=|*k*|^{3/2}sign*k*corresponding to the equation 3.11—

*α*≥2.*Fractal/quantized*: once*α*reaches 2, the solution exhibits the fractalization/quantization phenomenon, which provably happens when the dispersion relation is an integral polynomial, e.g.*ω*(*k*)=*k*^{m}, where . In such cases, theorem 3.2 implies that the solution quantizes at rational times, meaning that, when*t*=*βp*/*q*, it is piecewise constant on intervals of length 2*π*/*q*, and, as proved by Oskolkov [11], continuous but nowhere differentiable at irrational times. For non-polynomial dispersion with integral asymptotic exponent , most frames in the movie show a fractal profile, but occasionally the profile will abruptly quantize, with jump discontinuities punctuating noticeably less fractal subgraphs. We presume that, as in the polynomial case, the times at which the solution (approximately) quantizes are densely embedded in the times at which it has a continuous, fractal profile. However, a rigorous proof of this observation appears to be quite difficult, requiring delicate Fourier estimates. Examples include the linearized Korteweg–de Vries equation (2.6), and linearized Benjamin–Ono equation (2.20), both of which precisely quantize, as seen in figures 8 and 11, and the linearization of the integrable Boussinesq equation (2.2), which appears to do so approximately, as in figure 9. On the other hand, if is not an integer (or at least not very close to an integer), only fractal solution profiles are observed; see figure 10 for the case*ω*(*k*)=|*k*|^{5/2}sign*k*corresponding to the equation 3.12Observe that, as an immediate consequence of corollary 3.4, for integral polynomial dispersion relations, and, presumably, more general relations in this region, the quantization effect persists under the addition of noise to the initial data, although the small jumps at rational numbers with large denominators would be overwhelmed by the noise, whereas at irrational steps one ends up with a ‘noisy fractal’.

Let us outline a justification of our observation that the quantization and fractalization phenomena hold for dispersion relations that are asymptotically polynomial at large wavenumber. Assume that
3.13where , and *P*(*k*) is an integral polynomial. Let us rewrite the corresponding complex Fourier series solution in the form
We assume that the coefficients *b*_{k}, which are prescribed by the initial data, are uniformly bounded, |*b*_{k}|≤*M*, as is the case with the step function. With this assumption, *r*_{k}=*O*(*k*^{−2}), and hence the final series represents a uniformly and absolutely convergent Fourier series, whose sum is a continuous function. Thus, the discontinuities of *u*(*t*,*x*) will be determined by the initial series, which, by theorem 3.2, exhibits jump discontinuities and quantization at rational times. We conclude that solutions to linear equations with such asymptotically integral polynomial dispersion relations will exhibit quantization at the rational times *t*=*βp*/*q*, where *β*=2*π*/*λ*. A rigorous proof of fractalization at irrational times will require a more detailed analysis.

A similar argument explains why the discontinuities in the solution remain (almost) stationary, when the asymptotic exponent *α*<0. Formally, suppose *ω*(*k*)=|*k*|^{α}*μ*(*k*), where *α*<0 and *μ*(*k*)=*O*(1). Then, the Fourier series solution has the form
Again, assuming that the coefficients *b*_{k} are uniformly bounded, the remainder terms are of order *k*^{−1+α} with *α*<0, and hence represent an uniformly convergent series whose sum is continuous. We conclude that the discontinuities of *u*(*t*,*x*) are completely determined by those of the leading term, which are stationary.

Finally, in figure 12, we display graphs of the temporal evolution of the solution at the origin, *u*(*t*,0), for a representative subset of the dispersion relations we've considered. In the first two figures, for the RLW/BBM and water wave dispersion, we graph on the long time interval 0≤*t*≤100, while in the other plots, the time interval is 0≤*t*≤2*π*. When the asymptotic power (3.9) governing the high wavenumber dispersion satisfies *α*≤1, the time plot appears smooth. Following the non-rigorous analysis presented in Berry [2], §2, for *α*≥1, the fractal dimension of the time plot should be 2−1/(2*α*). This seems correct as stated when *α*≥2, as these plots become increasingly fractal with increasing asymptotic dispersion power *α*. However, for *α* near 1, the plot appears to be initially smooth for a short time, which is followed by the appearance of increasingly rapid oscillations. After some further time has elapsed, the graph eventually assumes a fractal form.

## 4. Fractalization and quantization in nonlinear systems

In this final section, we turn our attention to the nonlinear regime. Basic numerical techniques will be used to approximate the solutions to the Riemann initial value problem on a periodic domain. While we do not have any rigorous theoretical results to justify the resulting observations, our numerical studies strongly indicate that the variety of qualitative solution features observed in the underlying linear problem persist, at the very least, into the weakly nonlinear regime. In particular, the fractalization and quantization of solutions to the Airy evolution equation are observed in the numerical solutions to both the integrable Korteweg–de Vries equation and its non-integrable quartic counterpart. Similar phenomena have been found in nonlinear Schrödinger equations, but these will be reported on elsewhere.

Our numerical algorithms are based on a standard operator splitting method [47]. Consider the initial value problem for an evolution equation
4.1where *K* is a differential operator (linear or nonlinear) in the spatial variable. Here, the system is always supplemented by periodic boundary conditions, but the numerical scheme is of complete generality. Assuming the basic existence and uniqueness of solutions to the initial value problem, at least over a certain time interval 0<*t*<*t*^{☆}, we denote the resulting solution to (4.1), i.e. the induced flow, by *u*(*t*, ⋅ )=*Φ*_{K}(*t*)*u*_{0}.

Operator splitting relies on writing the spatial operator as a sum
4.2of two simpler operators, which can each be more readily solved numerically. In all cases treated here, *L*[*u*] represents the linear part of the evolution equation (4.1), whereas the nonlinear terms appear in *N*[*u*]. For example, writing the Korteweg–de Vries equation in the reduced form^{3}
4.3we split the right-hand side into the sum of the linear part *L*[*u*]=−*u*_{xxx}, with cubic dispersion, and the nonlinear part given by the *inviscid Burgers' operator* *N*[*u*]=*uu*_{x}.

For numerical purposes, we are interested in approximating the solution *u*_{n}(*x*)=*u*(*t*_{n},*x*) at the mesh points 0=*t*_{0}<*t*_{1}<*t*_{2}<⋯ which, for simplicity, we assume to be equally spaced, with *Δt*=*t*_{n}−*t*_{n−1} fixed. (Of course, we will also discretize the spatial variable.) We use to denote the numerical approximation to the solution to the original evolution equation at times *t*_{n}=*nΔt*, which is obtained by successively applying the flows *Φ*_{L} and *Φ*_{N} corresponding, respectively, to the linear and nonlinear parts in the splitting (4.2). The simplest splitting algorithm is the *Godunov scheme*:
4.4In favourable situations, for a suitable choice of norm and appropriate conditions on the initial data, it can be proved that the Godunov scheme is first-order convergent:
The convergence can be improved to second order by use of the *Strang splitting scheme*:
4.5In this case,
again under appropriate assumptions. In our numerical experiments, the results of the Godunov and Strang splitting are very similar, and so we display only the Godunov versions in the figures.

In particular, to solve the initial value problem for the Korteweg–de Vries equation (4.3), we apply the Godunov splitting scheme (4.4) to its linear and nonlinear parts 4.6each of which is subject to periodic boundary conditions. Periodicity and discretization of the spatial variable enables us to apply the fast Fourier transform (FFT) to exactly solve the discretized linear equation. On the other hand, the nonlinear inviscid Burgers' equation is first written in conservative form as above. As in the linear case, we use FFT to write the numerical approximation in the form of a Fourier series, apply convolution to square the numerical solution and then differentiate by multiplication by the frequency variable. The resulting system of ordinary differential equations is mildly stiff, and is integrated by applying the backward Euler method, based on the implicit midpoint rule, using fixed point iteration to approximate the solution to the resulting nonlinear algebraic system to within a prescribed tolerance. Because the time step is small, difficulties with the formation of shock waves and other discontinuities do not complicate the procedure.

Our results, at rational and irrational times, are plotted in figures 13 and 14. These clearly suggest the presence of both quantization and fractalization phenomena in the nonlinear system. It is striking that these appear, not just, as one might expect, in the weakly nonlinear regime, but in a fully nonlinear framework in which the coefficients of both the linear and nonlinear terms in (4.3) have comparable magnitude. Indeed, recent results of Erdoğan *et al.* [48], imply that, for periodic boundary conditions, the temporal evolution of high-frequency data for the Korteweg–de Vries equation and other nonlinear water wave models is nearly linear over long time intervals. So, while it may happen that these effects may eventually be overcome by the nonlinearity, any serious investigation will require very accurate long time numerical integration schemes. The corresponding movies have been posted at http://www.math.umn.edu/∼olver/mvdqn.html. Our numerical results are also in agreement with a theorem of Erdoğan & Tzirakis [49], corollary 1.5, that, when the initial data *u*(0, ⋅ ) lies in the Sobolev space *H*^{s} with , then the difference between solutions of the Korteweg–de Vries equation and its linearization can be controlled, meaning that the difference lies in *H*^{r} for any , and has at most polynomially growing *H*^{r} norm.

To verify that these phenomena are not restricted to integrable models such as the Korteweg–de Vries equation, in figures 15 and 16, we plot the solution to the non-integrable *quartic Korteweg–de Vries equation*
4.7The numerical solution algorithm is again based on Godunov splitting, using the same Fourier-based algorithms to integrate the individual linear and nonlinear parts.

Keep in mind that the coefficients in both the Korteweg–de Vries and quartic Korteweg–de Vries equations can be rescaled to any convenient non-zero values, and thus quantization and fractalization of the solution can be expected to appear in any version, perhaps with the overall magnitude of the effect governed by the relative size of the initial discontinuities. Thus, the Talbot effect appears to be rather robust, and hence of importance to a wide range of linear and nonlinear dispersive equations on periodic domains.

## 5. Further directions

We conclude by listing a few of the possibly fruitful directions for further research.

— Our numerical integration schemes are fairly crude, and it would be worth implementing more sophisticated algorithms. See Zuazua [50] for a recent survey on the numerics of dispersive partial differential equations. In the linear case, the complicated behaviours exhibited by the (partial) Fourier series sums would serve as a good testing ground for the application of numerical methods for dispersive waves to rough data.

— While the explicit formulas for the solutions are rather simple Fourier series, the amazing variety of observed behaviours, solution profiles, and fascinating spatial and temporal patterns indicates that formulation of theorems and rigorous proofs will be very challenging. Progress may well rely on delicate estimates of the type used in the analysis of number-theoretic exponential sums of Gauss and Weyl type [15]. One is also reminded of Bourgain's celebrated analysis of the periodic problems for the nonlinear Schrödinger and Korteweg–de Vries equations [51,52], which involves similarly subtle number-theoretic analysis.

— Zabusky and Kruskal's numerical experimentations that led to the discovery of the soliton were motivated by the fact that the Korteweg–de Vries and Boussinesq equations are continuum limits of the nonlinear mass-spring chains, whose surprising lack of thermalization came to light in the seminal work of Fermi

*et al.*[53]. It would be interesting to see how the foregoing dispersive effects might be manifested in the original Fermi–Pasta–Ulam chains when the initial displacement includes a sharp transition.— An interesting and apparently difficult question is to determine the fractal dimension of the graphs of solutions at fixed times that are given by such slowly converging Fourier series. The wide variety of individual space plots in the accompanying figures indicates that this is considerably more subtle than the Fourier series that were rigorously analysed by Chamizo & Córdoba [54], or those of Weierstrass–Mandelbrot form that were investigated by Berry & Lewis [55]. Indeed, for fixed

*t*, all of our series have the same slow*O*(1/*k*) rate of decay of their Fourier coefficients, and hence the fractal nature of their graphs depends essentially on the detailed dispersion relation asymptotics. Further, when the large wavenumber asymptotic power is in the range 1<*α*<2, the graphs appear to be (mostly) smooth over some initial time interval, then increasingly oscillatory, only later apparently achieving a fractal nature. Using Besov space methods, Rodnianski [9] rigorously proves that the graphs of real and imaginary parts of rough solutions to the linear Schrödinger equation have fractal dimension at a dense subset of irrational times, reconfirming some of Berry's observations [2]. We do not know how difficult it would be to extend Rodnianski's techniques to more general dispersion relations.— We have concentrated on the periodic boundary value problem for linearly dispersive wave equations. The behaviour under other boundary conditions, e.g.

*u*(*t*,0)=*u*_{x}(*t*,0)=*u*(*t*,2*π*)=0, is not so clear because, unlike the Schrödinger equation, these boundary value problems are not naturally embedded in the periodic version. Fokas & Bona [56,57,58] have developed a new solution technique for initial-boundary value problems for linear and integrable nonlinear dispersive partial differential equations based on novel complex integral representations. It would be instructive to investigate how the effects of the dispersion relation in the periodic and other problems are manifested in this approach.— Another important direction would be to extend our analysis to linearly dispersive equations in higher space dimensions. Berry [2] gives a non-rigorous argument that, at least for the linearized Schrödinger equation, a fractal Talbot effect persists for higher dimensional boundary value problems. Other important examples worth investigating include the integrable Kadomtsev–Petviashvili and Davey–Stewartson equations [41], and non-integrable three-dimensional surface wave models found in [27,30].

## Acknowledgements

We thank Michael Berry, Jerry Bona, Adrian Diaconu, Dennis Hejhal, Svitlana Mayboroda and Konstantin Oskolkov for supplying references and helpful advice. We also thank the anonymous referees for helpful and thought-provoking comments. The work of second author is supported in part by NSF grant no. DMS–1108894.

## Footnotes

↵1 The work of Oskolkov [11], which analyses both the linear Schrödinger equation and the linearized Korteweg–de Vries equation on periodic domains, is a notable exception.

↵2 The figure can be viewed as graphing the real part of the right-moving component of the solution to the periodic Riemann initial value problem.

↵3 The model version (2.15) can be changed into the reduced form by applying a combination of scaling and Galilean boost.

- Received July 9, 2012.
- Accepted September 17, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.