Dispersive Fractalization in Linear and Nonlinear Fermi-Pasta-Ulam-Tsingou Lattices

We investigate, both analytically and numerically, dispersive fractalization and quantization of solutions to periodic linear and nonlinear Fermi-Pasta-Ulam-Tsingou systems.

as expected, but rather exhibited an unexpected recurrence, in which energy from the low frequency modes would initially spread out into some of the higher modes but, after a certain time period, the system would almost entirely return to its initial configuration, as eloquently described in the above quote from Ulam's autobiography. It can be argued that the FPUT calculation inaugurated the contemporary fields of computational physics, [15], and experimental mathematics, [40], meaning "computer-based investigations designed to give insight into complex mathematical and physical problems that are inaccessible, at least initially, using more traditional forms of analysis" -thereby "spark[ing] a revolution in modern science".
In an attempt to understand this intriguing and unexpected phenomenon, Zabusky and Kruskal, [61], proposed a continuum model 3 that, in its unidirectional manifestation, turned out to be the Korteweg-deVries (KdV) equation, originally derived by Boussinesq, [8], in his pioneering studies of surface water waves. Zabusky and Kruskal's numerical integration of the periodic initial-boundary value problem for the Korteweg-deVries equation, starting with a smooth initial profile, led to their discovery of the soliton 4 and the consequent creation of an entirely new branch of mathematics -integrable nonlinear partial differential equations -whose remarkable repercussions continue to this day, [16]. The overall impact of the FPUT numerical experiment on modern mathematics and physics cannot be understated.
Much of the subsequent extensive research into the FPUT problem has concentrated on understanding their original surprising observation of the non-ergodic (almost) recurrence of the initial state; see [22] for historical remarks up to 1992. One key issue is to determine whether thermalization or ergodicity occurs if one extends the time range to be sufficiently long. An initial conjecture was that, while on a relatively short time scale the system returned close to its original configuration, subsequent behavior at each such "period" would take it farther and farther away. This was disproved in [53], where Tuck and Menzel (Tsingou) observed "superperiods" after which the initially increasing deviations decrease to a much smaller value. On the other hand, for sufficiently large nonlinearities, the FPUT recurrence no longer exists, and the systems exhibits a "strong stochasticity threshold", [28,29]. See [56] and [22] for reviews.
Explanations of the observed behavior have also appealed to Kolmogorov-Arnold-Moser (KAM) theory, [41,56], in which one views FPUT as a perturbation of either the linearized system, or the integrable Toda lattice, [52], or an integrable Birkhoff normal form, [22,42], or even some nearby as yet unknown integrable system. Other research directions include the derivation of explicit solutions; in [23], Friesecke and Wattis proved the existence of solitary waves solutions for a broad range of nonlinearities. This was followed by a series of papers investigating traveling wave solutions of lattices and their approximation by Korteweg-deVries solitons; see [39] for details. See also [1] for results on periodic solutions. Another productive line of research has been to investigate what happens when not all the masses and springs are identical, particularly those with alternating masses and/or spring stiffnesses; see, for example, [9,24,39].
In an unrelated but also surprising development, in the 1990s, Konstantin Oskolkov, [38], and, independently, Michael Berry and collaborators, [2,3,4], discovered the Talbot effect, which the latter named after a 1835 optical experiment of the Victorian scientist, inventor, and photography pioneer William Henry Fox Talbot, [51], who, in particular, invented the photographic negative. The Talbot effect arises in quantum mechanics through the behavior of rough solutions to the free space linear Schrödinger equation on a circular domain, i.e., subject to periodic boundary conditions. The evolution of a piecewise smooth but discontinuous initial profile, e.g., a step function, produces a fractal profile at irrational times (relative to the circumference of the circle) but "quantizes" into piecewise smooth but discontinuous profiles at rational times. Moreover, the fundamental solution, induced by an initial delta function, exhibits "revivals" at rational times, localizing into a finite linear combination of delta functions. This has the astonishing consequence that, at rational times, the solution to any periodic initial value problem is a finite linear combination of translates of the initial data and hence its value at any point on the circle depends only upon finitely many of the initial values! The effect underlies the experimentally observed phenomenon of quantum revival, [4,58,55], in which an electron, say, that is initially concentrated at a single location of its orbital shell is, at rational times, re-concentrated at a finite number of orbital locations.
The subsequent rediscovery of this remarkable phenomenon by the first author, in the context of the periodic linearized Korteweg-deVries equation, [34,35], showed that fractalization and quantization phenomena appear in a wide range of linear dispersive (integro-)differential equations, including models arising in fluid mechanics, plasma dynamics, elasticity, DNA dynamics, and elsewhere. Such linear systems exhibit a fascinating range of as yet poorly understood dynamical behaviors, whose qualitative features are tied to the large wave number asymptotics of the underlying dispersion relation. These studies were then extended, through careful numerical simulations, [11], to show that fractalization and quantization also appear in a variety of nonlinear dispersive equations, including integrable models, such as the nonlinear Schrödinger, Korteweg-deVries, and modified Korteweg-deVries equations, as well as their non-integrable generalizations with higher degree nonlinearities, [6]. (It is interesting to speculate as to how the history of the subject might have changed were Zabusky and Kruskal to have conducted their original investigations with discontinuous initial data!) Some of these numerical observations were subsequently rigorously confirmed in papers of Erdogan and collaborators, [12,18,19,20]; see also earlier analytical work of Oskolkov, [38], and Rodnianski, [43].
Given that the Korteweg-deVries equation and its generalizations arise as continuum models for FPUT chains, the question naturally arises as to whether dispersive fractalization and quantization effects appear in the discrete FPUT system. Resolving this question is the aim of this paper and its sequel(s). We initially focus our attention on the much simpler linear system, which can be analytically integrated. We find that, on an appropriately long time scale, the solutions to the periodic linear FPUT chain subject to a step function initial displacement do exhibit a suitably interpreted discrete version of fractalization. Although FPUT does not exhibit quantization in the sense that KdV does, we do observe coarse-scale similarity between the FPUT and KdV profiles at quantized times; this might be interpreted as approximate quantization, where the FPUT profile is the sum of a quantized profile and fractal "noise". On the other hand, we were unable to detect any trace of revival in the discrete linear system, which remains not entirely understood.
In the final section, we describe some numerical investigations in an initial attempt to extend our analysis to fully nonlinear FPUT systems. We briefly survey the use of geometric integrators that have been introduced to numerically integrate FPUT systems, including symplectic integrators such as the Störmer/Verlet scheme and trigonometric integrators.
However, we find that these are insufficiently accurate at high wave numbers, essential to our investigations, and so turn to a higher-order Hamiltonian splitting scheme to effect the computations. Currently we have some numerical data, for a relatively small number of masses, that indicate these phenomena also appear in nonlinear mass-spring chains. However, owing to the time scales involved and consequent large amount of computation required, we are currently unable to treat chains with a sufficient number of masses that would allow us to definitively address the basic question for nonlinear FPUT systems and generalizations thereof, and we defer the further development of analytical and numerical tools to study this intriguing problem to subsequent investigations.

Fermi-Pasta-Ulam-Tsingou Chains and their Continuum Models.
The Fermi-Pasta-Ulam-Tsingou (FPUT) system consists of a one-dimensional chain of masses that are connected by springs with nonlinear restoring forces, [21,22,39]. We will only consider the case when all the masses and springs are identical. The dynamics of the mass-spring chain follow immediately from Newton's Laws, taking the form of a system of second order ordinary differential equations for the mass displacements u n (t) for n ∈ Z at time t: where µ is the resonant frequency of the linear spring. The forcing function has the form where y indicates the elongation of an individual spring. The nonlinear intermass forcing term is prescribed by N (y) = W (y). As in [59,61], we will focus attention on the quadratic case when N (y) = α y 2 , (2.3) although higher degree polynomials in y, particularly cubic, are also of great interest. Another important system is the integrable Toda lattice, [52], where V (y) = α e β y . (

2.4)
Other examples include the Calogero-Moser integrable system and its trigonometric, hyperbolic, and elliptic generalizations, [10,33,50], where with P the Weierstrass elliptic function, as well as the Lennard-Jones potential, [27], which is used to model interactions between atoms and molecules.
In this note, we will concentrate on the periodic problem, viewing the system as a circular chain consisting of M masses, which are labelled so that u n 1 (t) = u n 2 (t) whenever n 1 ≡ n 2 mod M . Alternatively, one can consider an infinite chain, with the displacements at large distances subject to suitable decay conditions. Since the continuum models have smooth evolutionary behavior on the line -dispersive quantization being intimately tied to the periodic boundary conditions -we do not anticipate any unexpected effects in an infinite FPUT chain and so do not pursue it here. The Dirichlet problem, in which the first and last masses are pinned down, so that u 0 (t) = u M (t) = 0, with n = 0, . . . , M , is also of interest, [59], but its analysis will be deferred to subsequent investigations.
Following [61,59,44,45], we endeavor to better understand the discrete FPUT dynamics by passing to a continuum model. To this end, we assume the masses lie on a circle of fixed radius, say the unit circle of circumference 2 π. As the number of masses M → ∞, the equilibrium intermass spacing h = 2 π/M → 0. To maintain consistency, the time must be correspondingly rescaled, t → h t, and so we consider the system where c = µ h will be the wave speed of the limiting scalar wave equation.
We can view the individual displacements as the sample values of an interpolating function u(t, x) that is 2 π periodic in x, so that are the nodes or reference positions of the masses. To produce a continuum model, we apply Taylor's theorem to expand where the right hand side is evaluated at (t, x n ). Substituting into (2.1) and replacing x n → x, we arrive at the dispersive partial differential equation Thus, assuming the linear wave speed c and nonlinear scale parameter α are both O(1), we obtain, to second order in h, the bidirectional continuum model a potential form of the integrable (nonlinear) Boussinesq equation, [16], which can be obtained by differentiating (2.11) with respect to x and replacing u x → v. Note that, to leading order, the continuum model (2.11) coincides with the standard linear wave equation u tt = c 2 u xx with wave speed c > 0. The Korteweg-deVries equation is obtained through a standard "unidirectional factorization" of the preceding bidirectional system, [57], i.e., assuming the waves are only propagating in one direction, say in the direction of increasing x, producing Remark: For the cubically forced FPUT system, the unidirectional model is the integrable modified Korteweg-deVries equation, [45,59], in which the nonlinear term is a multiple of u 2 u x . Higher degree forcing polynomials produce generalized Korteweg-deVries equations with higher degree nonlinearities, which are no longer integrable, and, in fact, can produce blow up of solutions, [6].
To initiate our investigations, let us ignore the nonlinear contributions and concentrate on the linear FPUT system and its continuum models. The rescaled linear system becomes what is known as the discrete wave equation, [39]: (2.14) which, to the same order in h, has bidirectional continuum model known as the linearized "bad Boussinesq equation", [44], owing to the fact that it is an ill-posed partial differential equation. Indeed, its dispersion relation is found by the usual method, [57], of substituting the exponential ansatz producing the algebraic equation relating the temporal frequency ω to the wave number (spatial frequency) k. Because p 4 (k) < 0 for k 0, the bad Boussinesq model (2.15) is not purely dispersive since the high wave number modes induce complex conjugate purely imaginary values of ω and hence exponentially growing modes, [13], that underly the ill-posedness of the initial value problem. Interestingly, the corresponding linearized unidirectional Korteweg-deVries model does not suffer from this instability, since its dispersion relation is everywhere real, and coincides with the Taylor expansion, at k = 0, of one of the two branches of the bidirectional dispersion relation (2.17). One way of regularizing the bad Boussinesq model is to replace it by the linearized bidirectional Korteweg-deVries equation which agrees to order h 2 . It has solutions that are an exact linear combination of right-and left-moving linear KdV solutions. Positivity of the right hand side of the corresponding dispersion relation implies well-posedness of this sixth order model.
Remark: It remains somewhat mysterious to the authors how an ill-posed bidirectional wave model can have right-and left-moving unidirectional constituents that are both wellposed. This is clearly a consequence of the use of low wave number Taylor expansions of the dispersion relation near k = 0 in the approximation procedure, but we would argue that this seeming paradox warrants further study. Another means of regularizing the linear, and hence the nonlinear, continuum model is to retain the order h 4 terms in the preceding Taylor expansion. This produces the sixth order linear partial differential equation with dispersion relation Since p 6 (k) > 0 for all k, the regularized model (2.22) is purely dispersive, and hence well-posed, in that all Fourier modes maintain their form under translation. An alternative regularization procedure that avoids increasing the order of the differential equation is to replace two of the x derivatives in the fourth order term in the bad Boussinesq equation (2.15) by t derivatives, using the fact that, to leading order, u xx = c −2 u tt + O(h 2 ), thereby producing the continuum model known as the linear Boussinesq equation, [57, pp. 9, 462]. It arises as the linearization of Boussinesq's bidirectional model for shallow water waves, and has also been proposed as a model for DNA dynamics, [46]. In this case, the dispersion relation is (2.25) and hence the equation is purely dispersive and well-posed. Let us next derive the analogous "dispersion relation" for the discrete linearized FPUT system (discrete wave equation) (2.14); see also [60]. Substituting the usual exponential ansatz (2.16), evaluated at the node x = x n = n h, into (2.14) produces (2.26) We thus deduce the discrete FPUT dispersion relation that determines the temporal frequencies ω in terms of the wave numbers k. Since ω(k) is real for all k = 0, . . . , M , the FPUT system can be regarded as dispersive, in that the different Fourier modes propagate unchanged at different wave speeds. This implies "well-posedness" or, more accurately since we are dealing with a system of ordinary differential equations, stability of the zero equilibrium solution. Moreover, observe that the continuum model dispersion relations (2.17), (2.23), (2.21), and (2.25) all have the same order h 4 Taylor expansion at k = 0 as (2.27), and hence approximate it well at low or even moderately large wave numbers. However, they exhibit different high wave number asymptotics, which, as noted in [11], is the key property that governs the dispersive fractalization of rough solutions.

The Riemann Problem.
As in [11,34], we are particularly interested in the Riemann problem -a particular initial value problem of fundamental importance in the study of hyperbolic wave equations and shock waves, [47]. Here, the initial displacement is a (periodically extended) step function: whose values at integer multiples of π are specified in accordance with the convergence properties of its Fourier series We will also, for simplicity, impose zero initial velocity, concentrating on the pure displacement problem: By an elementary trigonometric identity, we can split the standing wave summands in (3.4) into right-and left-moving unidirectionally propagating waves: where the factor of 1 2 is introduced for comparative purposes, so that all three solutions have the same initial displacement: The right-moving constituent is and its left-moving counterpart is obtained by replacing t by − t. A key feature of these series solutions that produces the dispersive fractalization and quantization effects is the slow decay of their Fourier coefficients, which implies that they are conditionally but not absolutely convergent. The corresponding step function initial data for the discrete FPUT problem is obtained by sampling (3.1) at the nodes. To make the connection, we will take the number of masses to be even, M = 2 m, and the nodes to be x = x n = n h = π n/m, n = − m, . . . , m, identifying x −m = x m by periodicity. Thus, the initial data for the Riemann problem for the FPUT system is given by In other words, we displace each mass lying on the "right semicircle" by 1 unit, while those on the left remain at their equilibrium position, except the two masses lying at the interface that are displaced by only half a unit. As in (3.3), the masses are assumed to be at rest initially: We use the Discrete Fourier Transform to write the solution as a Fourier sum over the fundamental periodic modes, [36,Section 5.6], where the symbol ∼ will mean that the left and right hand sides agree at the nodes, i.e., when x = x n . The Discrete Fourier Transform applied to the sampled step function produces the interpolating discrete Fourier sum meaning that the displacement of the n-th mass is given by sampling the right hand side at the nodes: u n (t) = u(t, x n ), x n = n h = π n m . (3.12) Again, the solution (3.11) is a linear combination of standing waves, and can be decomposed into left and right moving constituents, as in (3.5). The right-moving constituent has the explicit form As above, its left-moving counterpart is obtained by replacing t by − t.
Remark: Note that if we omit the constant term, the bidirectional solutions constructed in (3.4), (3.11) also satisfy Dirichlet boundary conditions, with a half-size signum function as initial condition: u(t, x) = 1 2 sgn x, − π < x < π. Thus, all our subsequent remarks on their behavior also apply to this Dirichlet initial-boundary value problem. On the other hand, the unidirectional constituents do not individually satisfy Dirichlet boundary conditions. As shown in [34,38], the canonical linearized Korteweg-deVries equation u τ + u ξξξ = 0 (3.14) with step function initial data and periodic boundary conditions on − π ≤ ξ ≤ π exhibits dispersive fractalization and quantization in the following sense. At irrational times τ > 0, meaning τ /π ∈ Q, the solution profile u(τ, ξ) is a continuous but non-differentiable fractal.
On the other hand, at rational times, τ /π ∈ Q, the solution is discontinuous, but piecewise constant! Indeed, if τ = 2 π p/q where p, q ∈ Z have no common factors, then the solution is constant on the intervals 2 πj/q < ξ < 2 π(j + 1)/q for j ∈ Z. Thus, the larger the denominator q, the shorter the intervals of constancy. (It is possible that the solution achieves the same constant value on one or more adjacent intervals, and so an interval of constancy may be larger than specified above. See [37] for a number-theoretic investigation into when this occurs.) A rigorous proof of the fractal nature of the solution at irrational times, including the estimate that its fractal dimension d is bounded by 3 2 ≤ d ≤ 7 4 , can be found in [18].
Indeed, the results of [18] imply that a linear evolutionary integro-differential equation with dispersion relation that is (in an appropriate sense) asymptotic to a power of k at large wave numbers, ω(k) ∼ k α as k → ∞, for 1 = α > 0, exhibits fractalization at almost all times, provided the initial data is of bounded variation, but not too smooth, meaning it does not lie in any Sobolev space H β for β > 1 2 , i.e., its Fourier coefficients c n decay sufficiently slowly so that the series (1 + n 2 ) β | c n | 2 diverges. Moreover, if the asymptotic exponent α is integral, 2 ≤ α ∈ Z, numerical experiments, [11], indicate that the solution profiles quantize at other times, in the sense that they take a different form from the "generic" fractal profiles: piecewise smooth with either jumps or cusps, possibly with some much smaller fractal modulation superimposed. However, being so far based on numerical calculations, it is not yet known if the observed small scale fractals on the quantized profiles are genuine or just a manifestation of numerical error. See the forthcoming paper [7] for further developments.
Warning: The fractal dimension of the graph of a function has the potential to be misleading. For example, the graph of the sinusoidal function f (x) = sin(1/x) has fractal dimension 2 even though it is perfectly smooth, even analytic, except at the singularity at x = 0. For this reason, the results in [18], while striking, are, on a deeper level, unsatisfying. It would be preferable to know, or at least have estimates on, the Hausdorff dimension of (sections of) such solution profiles; however this seems to be beyond current analytic capabilities.
Turning our attention to the linearized Korteweg-deVries model (2.18), the leading first order term c u x represents linear transport moving at speed − c, and only affects the solution by an overall translation. We can map (2.18) to the preceding canonical form (3.14) by a Galilean shift to a moving coordinate frame, coupled with a rescaling of the time variable: Due to the temporal scaling, the dispersive quantization pattern occurring in the solution to the canonical KdV equation (3.14) at a rational time τ = 2 π p/q will now appear (suitably translated) at a much later time, namely t = 48 π p/(c h 2 q). Since the normalized model (3.14) exhibits quantization at every rational time, the same is true (modulo the scaling factor used to distinguish rational from irrational) of the FPUT model version (2.18). However, in the latter model, at a rational time t = O(1), the denominator q will be very large, of order O(h −2 ), hence the intervals of constancy are extremely small, O(h 2 ), and thus undetectable at the physical level, which has spatial scale ∆x = O(h). At such scales, both physical and graphical, it will be practically impossible to distinguish such profiles from fractals. Now let us compare the solution to the periodic Riemann initial value problem for the discrete linear Fermi-Pasta-Ulam-Tsingou system (2.14) with those of the three well-posed linear model equations: the bidirectional Korteweg-deVries model (2.20), the sixth order model (2.22), and the regularized Boussinesq model (2.24). For our numerical comparisons, we sum the same modes in the discrete and continuous Fourier series, truncating at k = m, and plot the resulting profile; in other words, we are performing exact (well, modulo floating point round off) computations on the truncated (discrete) Fourier series and not a numerical approximation. The initial data is a (periodically extended) step function. In the continuum models, one can work with either the continuous Fourier series representation of the initial data (3.2), or the corresponding discrete Fourier sum (3.10). However, in the given situations, we observe no appreciable differences between the solution profiles, and hence will use the discrete version in all figures representing solutions to the Riemann initial value problem. We fix the wave speed c = 1, and, for most of our numerical investigations, work with m = 512, so there are M = 1024 masses and h = π/m ≈ .006136. Solutions for other numbers of masses have been calculated, and the overall conclusions are similar, although the fewer their number, the less pronounced some effects tend to be.
The first observation is that, on the time scale and resolution under consideration, there is almost no noticeable difference between the sixth order and regularized Boussinesq models, and so we choose only to display the results for the latter. We will plot both the bidirectional solution u(t, x), as given in (3.4), (3.11) and its unidirectional right-moving constituent (3.6), (3.13). We consider the effects at times that are selected from three regimes: what we will call short times, where t = O(1), medium times, where t = O(h −1 ), and long times, First, on short time scales, the solutions to all four models exhibit little appreciable difference. For example, consider the profiles at t = 1 5 π graphed in Figure 1 -the top row being the full bidirectional solution and the bottom row its right-moving unidirectional constituent. Since all profiles remain rather similar at short times, in Figure 2 we just graph the FPUT solution profiles. What we observe is that, on the short time scale, the solution is an oscillatory perturbation of the traveling wave solution to the corresponding limiting bi-  In particular, at t = 1 2 π, the right-and left-moving waves have cancelled each other out, leaving only a constant solution profile for the traveling wave solution, with a superimposed fractal residue in the FPUT system, as well as its continuum models, all three of which take on a comparable form.
At medium times, of order O(h −1 ), the fractal nature of the oscillations superimposed upon the traveling wave solution profile has become more pronounced. Again, both the FPUT system and its continuum models exhibit similar behavior; Figure 3 graphs the former at some representative medium times.
Once we transition to the long time scale, of order O(h −2 ), significant differences arise in the observed behaviors. First let us consider the solution profiles at the irrational (meaning that h 2 t/π ∈ Q) times t = 1/h 2 ≈ 26561 and t = 400000, plotted in Figures 4 and 5. All 6 Recall that we have set c = 1. The unidirectional constituents are more uniformly fractal, while the bidirectional solutions exhibit some semi-coherent regions, perhaps indicating some remnant of the intervals of constancy of a nearby rational profile. However, at long rational times, the solution profiles differ dramatically, as illustrated in Figures 6 and 7 for two representative such times. The linearized KdV solution has quantized into a piecewise constant profile, whereas the FPUT system and the Boussinesq models retain a common fractal form. On the other hand, the latter profiles exhibit an observable adherence to the underlying piecewise constant KdV solution, albeit with a superimposed fractal modulation. Recall that the initial data is the discrete Fourier representation (3.10) of the step function. Interestingly, if we use the continuous version (3.2) instead, which only differs in its higher frequency modes, the graphs do not appreciably change, and so are not displayed. The only noticeable difference is that piecewise constant KdV profile exhibits a more pronounced Gibbs phenomenon at the discontinuities. Now, one might argue that the differences between the quantized KdV profiles and the fractal FPUT ones is due to discrepancies in their dispersion relations at the high frequency modes. So let us try eliminating the higher frequency terms by truncating the Fourier sum in order to bring the solutions closer in spirit. It is surprising that one must eliminate a large majority of the high frequency modes before their respective truncated solution profiles begin to align at the rational quantized times; on the other hand, at the irrational fractalized times they are quite similar no matter how one truncates.
Keeping in mind that we are working with m = 512 total modes, the first plots in Figure 8 are at the same quantized time illustrated in Figure 6, and show the results of summing the first 1/8, 1/16, and 1/32 of the terms in the discrete Fourier summation (i.e., 64, 32, and 16 terms). The top row shows the resulting truncated profiles for the unidirectional FPUT solution (3.13), while the bottom row shows the corresponding truncated KdV profile (3.6). Somewhat surprisingly, even retaining 1/8 of the modes leads to significant differences; these differences persist (albeit more subtly) at the 1/16 scale, and only at the very coarse 1/32 scale do they look very close. On the other hand, in Figure 9, which illustrates the   same results at the irrational time that were shown in Figure 6, all three pairs of truncated profiles exhibit very similar features, while, as expected, the overall local fractal nature of the profile is curtailed as the number of terms decreases. Of course, the FPUT mass-spring chain is not a continuum, and so the values of the trigonometric solution (3.11) -or its unidirectional counterpart (3.13) -only have physical meaning at the mass nodes. For the above cases of M = 1024 masses, the differences are imperceptible. To better illustrate, in Figure 10 we plot solution profiles for M = 64 masses, at selected times, comparing the discrete mass displacement profiles, the corresponding continuum FPUT bidirectional solution (3.11), and the continuum bidirectional solution (3.4) with KdV dispersion (2.19). Observe that the effects are quite similar to what was presented above, but less pronounced owing to the relatively small number of masses. Finally, let us investigate whether the Talbot revival phenomenon mentioned in the introduction appears in the FPUT system. Somewhat surprisingly, given the noticeable effects of dispersive quantization at long rational times, numerical experiments have failed to reveal any observable trace of revival.
To model the delta function initial displacement, we displace the center mass by a unit 7 . (This is equivalent to equipartitioning the initial energy into all the Fourier modes.) Figure 11 plots the resulting solutions, in the case of M = 64 masses, whose initial data is the (Fourier series for) the delta function, at the indicated long rational times. These, as always, are obtained by explicitly summing over the first m = 32 modes. Keep in mind that the Fourier series of the delta function and resulting fundamental solution to the continuum model is highly oscillatory, and only converges weakly to the distributional revival profile, consisting of a finite linear combination of delta functions, at rational times. The first column plots the solutions to the bidirectional KdV model; the discrete oscillatory peaks indicate the appearance of a revival. The second column plots the corresponding FPUT solution; here, there is no appreciable sign of concentration of the solution profiles and hence no apparent revival. Similar behaviors have been observed at other (long) times, with varying number of masses. The KdV profiles are fractal at irrational times and concentrated in accordance with t = 24 π/(5 h 2 ) t = 24 π/h 2 a revival at rational times, whereas the FPUT profiles are more or less uniformly oscillatory at all times.

Numerical Investigation of Nonlinear FPUT Chains.
For the linear problems featured in the previous sections, Fourier series techniques make it possible to compute exact solutions (up to floating point error) for extremely long times at rather low computational expense-essentially the cost of a Fast Fourier Transform (FFT). However, once nonlinearity is introduced, the resulting systems of differential equations can usually only be solved approximately, by using a numerical integrator with sufficiently small time step size. Since the FPUT system is Hamiltonian, it is natural to consider the class of symplectic integrators, which have excellent long-time numerical properties for Hamiltonian systems, including near-conservation of energy and slow growth of global error, [25].
There has been substantial investigation of the application of symplectic integrators (and other geometric numerical integrators) to the variant of the FPUT problem appearing in [24], with alternating stiff linear and soft nonlinear springs. Much of this work has focused on methods, such as trigonometric and modified trigonometric integrators, which can take large time steps in order to simulate the slow-scale nonlinear dynamics without needing to resolve the fast-scale linear oscillations. See [25, Chapter XIII] for a survey and [48,31] for more recent work on modified trigonometric integrators, including the IMEX method.
However, the phenomenon of dispersive quantization is quite delicate and requires the accurate resolution of high wave number oscillations, [11]. To illustrate the challenge this poses, we begin by considering the application of two widely-used symplectic integrators, the explicit Störmer/Verlet method and implicit midpoint method, to the harmonic oscillator y = −ω 2 y, ω ≥ 0. Let y j ≈ y(j∆t) denote the approximation produced by the method at the jth time step, where ∆t is the time step size. The Störmer/Verlet method gives y j+1 − 2y j + y j−1 = −(ω∆t) 2 y j , Hence, the Störmer/Verlet method produces harmonic oscillations with modified frequency ω. Note that | 1 2 ω∆t| ≤ 1 is necessary for ω to be defined, and this is precisely the linear stability condition for Störmer/Verlet. Similarly, the midpoint method gives, and substituting the exponential ansatz yields tan 2 1 2 ω∆t = ( 1 2 ω∆t) 2 =⇒ ω = 2 ∆t arctan 1 2 ω∆t = ω 1 − 1 12 (ω∆t) 2 + O (ω∆t) 4 . (4.4) In contrast to Störmer/Verlet, this modified frequency is defined without restrictions on ω∆t, which reflects the unconditional linear stability of the midpoint method. Figure 12 illustrates the effect of replacing ω by the modified frequencies ω in the dispersion relation for the bidirectional KdV model with m = 512 and t = 24 π/(5 h 2 ), showing that quantization does not become visible unless ∆t is very small, much smaller than needed for numerical stability. At ∆t = 10 −4 , the solution profile is qualitatively indistinguishable from the fractal profiles of the FPUT and Boussinesq models in Figure 6. The first hints of quantization are visible at ∆t = 10 −5 , and only by ∆t = 10 −6 does the solution appear to have converged sufficiently to the true, quantized KdV profile. Since t ≈ 4 · 10 5 such a simulation would require on the order of 10 11 time steps to observe quantization, even before the effects of nonlinearity are taken into account. Figure 13 repeats this experiment for a shorter chain with m = 32 and t = 24 π/h 2 ≈ 8 · 10 3 (compare Figure 10), where quantization occurs earlier and the time step size restriction is less severe, requiring on the order of 10 7 steps.
To overcome the computational obstacle of small step size, we turn to higher-order Hamiltonian splitting methods, cf. [25,32], which converge more quickly as ∆t → 0 while still preserving symplectic structure. Write the FPUT system (2.7) in the first-order form where F (y) = y + N (y). There are two natural ways to split this into two Hamiltonian systems, each of which can be integrated exactly. The first is where each of these can be integrated exactly since v is constant in the first system and u is constant in the second. The second splitting is where the first system is simply linear FPUT, which can be integrated exactly using the Fourier series techniques applied previously. Splitting methods approximate the time-∆t flow of the full system by alternating between the flows of the two subsystems, which we denote by ϕ A a i ∆t and ϕ B b i ∆t , where i a i = i b i = 1. One of the simplest splitting methods is ϕ B ∆t/2 • ϕ A ∆t • ϕ B ∆t/2 , called Strang splitting, [49]. For the splitting (4.6), this results in a first-order formulation of the Störmer/Verlet method, sometimes called velocity Verlet. Alternatively, for the splitting (4.7), where the linear flow ϕ A ∆t is computed using Fourier series techniques, the Strang splitting gives the so-called split-step Fourier method. Whenever the composition of flows is symmetric, the resulting method has even order, and in particular, the Strang splitting is order-2.
For the linear FPUT system, the split-step Fourier method gives the exact solution in a single step. However, once nonlinearity is introduced, we find that the split-step Fourier method has only modest benefits over Störmer/Verlet, which are outweighed in practice by the additional cost of performing an FFT and inverse FFT at every step. Table 1 illustrates this for the quadratic nonlinearity N = αy 2 , showing that the error in Störmer/Verlet is O (∆t) 2 uniformly as α → 0, while the split-step Fourier method is O α(∆t) 2 . Therefore, we focus our attention on methods based on the splitting (4.6), whose steps are less expensive to compute since they do not require transforming to Fourier space and back. Runge-Kutta-Nyström (RKN) methods are designed specifically for splittings of the form (4.6), i.e., for second-order Newtonian systems written in first-order form using a velocity variable. Of these, we chose the optimal 14-stage order-6 RKN method of Blanes and Moan, [5], which has the symmetric form  Blanes and Moan calculated these coefficients to minimize the constant in the order-6 error estimate. Although a step of this method is 14 times as expensive as a step of Störmer/Verlet, it requires vastly fewer steps owing to its faster convergence. For the linear FPUT problem with m = 32 and t = 24 π/h 2 , we observe that the Störmer/Verlet method gives an error on the order of 10 −3 for 10 8 steps and 10 −5 for 10 9 steps; by contrast, the RKN method gives an error on the order of 10 −3 for only 10 5 steps and 10 −9 for 10 6 steps. Figure 14 shows discrete solution profiles at t = 7500 and t = 24 π/h 2 ≈ 7822 for the FPUT model with m = 32 and quadratic nonlinearity N (y) = αy 2 , computed using the Blanes-Moan RKN method with 10 6 time steps. For α = 0.005, the solutions are nearly identical to the linear discrete FPUT profiles observed in Figure 10. As the strength of the nonlinearity increases, we observe noticeably different profiles at α = 0.05 and 0.5. However, unlike with the KdV model, we still do not observe any dispersive quantization, and there does not appear to be any qualitative difference between the profiles at t = 7500 and t = 24 π/h 2 , just as with the linear FPUT model.

Discussion.
In conclusion, we have shown that the solution to the periodic linear Fermi-Pasta-Ulam-Tsingou chain, with a step function as initial displacement and zero initial velocity, exhibits , where M is the number of masses, and h their spacing around the unit circle. Of course, being purely discrete, the solution cannot be genuinely fractal, even when extended into a continuous trigonometric interpolating function, because it only involves a sum over a finite number of Fourier modes. Moreover, it does not become fractal in the final h → 0 limit since the limiting equation is merely the very basic linear second order wave equation (3.16), whose solution is a combination of traveling waves, and hence piecewise constant at all times. Indeed, as h → 0, all of the observed behavior on medium and long time scales moves off to infinity, and the solution converges (weakly) to the corresponding solution to the simple limiting wave equation, with all times now being classified as "short". On the other hand, all of the regularized bidirectional continuum models have genuinely fractal solutions at a dense set of times, which closely follow the FPUT solution at the given resolution. In contrast, the bi-and uni-directional Korteweg-deVries models mimic the FPUT and Boussinesq solutions at irrational times, but exhibit a very different dispersive quantization profile at rational times. Be that as it may, the latter solutions retain an observable trace of the overall quantized character within their fractal profiles. Finally, the lack of any noticeable form of revival in the FPUT system is, in light of the previous results, not well understood. The next stage of this project will be to investigate which of these properties, if any, carry over to the other nonlinear FPUT systems and other nonlinear lattices of interest. Keeping in mind the numerical observations of dispersive quantization in the Korteweg-deVries equation and its generalizations, [11], we expect that this will indeed be the case. Numerical schemes that retain accuracy over long times, while allowing for efficient computation of longer chains, will be essential to this endeavor. This may require the construction of novel numerical methods, e.g., exponential-type Fourier integrators in the spirit of [26], in addition to the splitting methods considered above.