Energy spreading in strongly nonlinear disordered lattices

We study the scaling properties of energy spreading in disordered strongly nonlinear Hamiltonian lattices. Such lattices consist of nonlinearly coupled local linear or nonlinear oscillators, and demonstrate a rather slow, subdiffusive spreading of initially localized wave packets. We use a fractional nonlinear diffusion equation as a heuristic model of this process, and confirm that the scaling predictions resulting from a self-similar solution of this equation are indeed applicable to all studied cases. We show that the spreading in nonlinearly coupled linear oscillators slows down compared to a pure power law, while for nonlinear local oscillators a power law is valid in the whole studied range of parameters.


Introduction
The general understanding of the relation between chaos in classical systems, and ergodicity and thermalization is still far from complete nowadays. Intuitively, one expects from high-dimensional, non-integrable complex systems to demonstrate strong chaos and thus it seems reasonable to expect thermalization. This is essentially the fundamental assumption of classical thermodynamics [1]. The conditions under which this assumption can be safely made, however, is still an open question. It is not known what level of "chaoticity" or "complexity" is required to ensure thermalizing behavior. Chaos can be often seen as a consequence of nonlinear perturbations of an integrable system. The solutions of the unperturbed, integrable part of such a system are called modes. Starting with only a few initially excited modes, one can view thermalization as spreading in the mode-space, i.e. the excitation of new modes, due to the nonlinear chaotic interactions.
However, already the first attempts to follow such a thermalization of modes initiated by Fermi, Pasta, and Ulam revealed many extremely nontrivial effects, still not completely understood (see Refs. [2,3] for recent progress of the FPU problem). A very important case is when the integrable modes are spatially localized. Then the thermalization process is a spatial diffusion where more and more modes get excited. This allows to connect the rather abstract concept of thermalization in the mode-space with the very intuitive phenomenon of spatial diffusion. A prominent example where this has been studied very extensively in the past, is the interplay between nonlinearity and disorder. In this case, due to Anderson localization, linear eigenmodes are exponentially localized and the spectrum is purely discrete [4]. Recent numerical experiments with nonlinear disordered lattices have demonstrated that the initially localized wave packets spread in a very weak, subdiffusive manner [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. A complete theoretical understanding of the subdiffusive behavior has not been presented, but it is mostly agreed on that the spreading in these models is induced by weak chaos. However, the true asymptotic behavior is still discussed in some of these models with recent claims that spreading might stop due to an extinction of chaos [21,17,22].
In this paper we follow the scaling approach to this problem, first formulated in [18] and recently extended to two-dimensional systems [23]. We will try to establish and to check numerically the scaling relations for the properties of spreading, in dependence on the total energy of the initial wave packet. In these works, the nonlinear diffusion equation (NDE) was proposed to describe the spreading process and this assumption was verified by several numerical simulations. Here, we will generalize this model by introducing the fractional nonlinear diffusion equation (FNDE), and we will present new numerical results that will show that in some cases indeed only the FNDE gives a correct scaling description of the spreading process. We formulate the scaling relations based on this equation, and check their validity for nonlinear lattices.
We will start with formulating the object of our study, strongly nonlinear Hamiltonian lattices. We present a phenomenology of energy spreading and define the statistical quantities characterizing it in section 2. Next, we introduce a phenomenological model that we use to describe the properties of the spreading process, namely the fractional nonlinear diffusion equation (section 3). From its scaling properties, we derive spreading predictions for the strongly nonlinear Hamiltonian lattices. In section 4 we present extensive numerical calculations for different classes of the nonlinearity. These results are compared with the predictions from the FNDE and we identify different degrees of confirmation for different nonlinear classes. We end with concluding remarks where the found "universality classes" are summarized.

Strongly nonlinear lattices
The main goal of this paper is to study properties of energy spreading in strongly nonlinear lattices. By strongly nonlinear we understand lattices where the coupling is described by nonlinear functions that disappear in the linear limit. So there are no linear waves (phonons) in such lattices and energy transport can solely be induced by the nonlinear coupling. Such lattices can be introduced in the framework of equations for the complex amplitudes (and then one obtains a strongly nonlinear generalization of the nonlinear Schrödinger equation) or as a generalization of a Hamiltonian Klein-Gordon lattice. In this work we follow the latter way. Moreover, we restrict ourselves to pure power-law nonlinearities.

Hamiltonian
In one dimension we formulate a strongly nonlinear lattice in terms of a Hamilton function for positions q k and momenta p k of oscillators labelled by site index k: Here κ ≥ 2 and λ > 2 denote the powers of the on-site potential and the coupling term, respectively. For κ = 2 we have a chain of nonlinearly coupled linear on-site oscillators; for κ > 2 the on-site oscillators are nonlinear as well. Below we study situations with and without disorder, the latter is introduced via the variations of the parameters of the local potential ω k (these are linear frequences of the oscillators if κ = 2 and parameters of the nonlinear on-site potential if κ > 2). Note, that the integrable part of this system are uncoupled oscillators (α → 0), which means that the modes of the integrable system are extremely localized on one site.

Rescaling
Hamiltonian (1) contains two parameters W and α that determine the time scale and the ratio of local to coupling potentials. For different local and coupling nonlinearities, i.e. for κ = λ, we can get rid of these two parameters by rescaling the canonical variables and time as follows: with b = 1/(λ−κ). W and α disappear from the equations and we are left with the total energy E as the only relevant parameter depending on the initial state. Additionally, the distribution of local "frequencies" ω k is relevant, while the width of this distribution is rescaled together with W . We will consider three cases: (i) no disorder, ω k = 1; (ii) "soft" local oscillators, in this case the "frequencies" ω k are chosen iid. from [0, 1]; (iii) "hard" local disorder, here the "frequencies" ω k are chosen iid. from [0.5, 1.5]. In the rescaled coordinates the Hamiltonian reads One special and highly interesting case occurs when the on-site and coupling terms have the same nonlinearity κ = λ. As now all terms in q have the same power, one cannot set both parameters W and α to one by rescaling as before. Instead, one can use the remaining freedom to set the total energy E to unity: Particularly, this means that the energy is not a free parameter of the system but can rather be scaled to, say, E = 1, what also involves an appropriate change of the time scale. We note that the only remaining parameter is the ratio of strengths of on-site and coupling terms β = α/W . The rescaled Hamiltonian now reads

Phenomenology of energy spreading
For the Hamiltonian systems (1) we state the following question: How does an initially localized field spread over the lattice? We focus on very large systems, where boundary effects are not so important (we will discuss their relevance in some cases below). The distribution of energy is characterized with its density We start typically with non-zero values of w k in a small interval (in most runs 10 sites), by chosing initial momenta from a Gaussian distribution, and follow the distribution w k (t) in time.
In the case of a lattice without disorder (ω k = 1), regular waves can propagate along the lattice. Such localized solitary waves -compactons -have been thoroughly studied   in [24] for a lattice with W = 0 (i.e. without local potential). An initially localized perturbation emits compactons that dominate the process of energy spreading. For W = 0 it is not known if exact compactons exist in such lattices. In numerics, we quite often observe "quasi-compactons" that propagate ballistically over large distances but lose energy and therefore eventually stop. We illustrate this in Fig. 1a. In panel 1b we show the same initial conditions in a disordered lattice, here the propagation of "quasi-compactons" is blocked by disorder and one observes a slow spreading of the wavepacket, which will later be quantified as subdiffusive.
Additionally to the observation of "quasi-compactons", we see that in disordered strongly nonlinear lattices at any finite time the distribution of energies is strongly localized, and has sharp edges (we expect that generally the field at the edges decays superexponentially fast, as for breather solutions in such lattices [25]). This sharpeness is illustrated in Fig. 2a.

Measures of spreading
In a statistical context, the spread of a distribution w k can be quantified via entropies, most suitable are the Rényi entropies: that allow one to characterize also the spikeness/flattness of the distribution (in the context of energy spreading in disordered lattices this approach was introduced in [16]). We restrict here to the the entropies I 1 and I 2 , which are nothing else than the usual Boltzmann entropy and the logarithm of the participation number P :  The participation number is a characteristic of the width of the wave packet rather popular in the context of Anderson localization studies [16,10,26]. Both entropies define the effective width of the wave packet as L 1,2 = exp(I 1,2 ) (in particular, L 2 = P ). As individual dependencies L(t) demonstrate enormous fluctuations, we perform an averaging of the entropies I 1,2 (t) over many realizations of disorder, thus obtaining smoothly growing widths L(t).
For the strongly nonlinear lattices another approach [18] is even superior to the calculation of the entropies. Here we determine the width L of the wave packet as the distance between its sharp edges as seen in Fig. 2a (independently on the distribution of the energy between these edges). After determining the spatial extend L, we measure the time δT required to excite one new lattice site. So suppose we have L lattice sites being excited, then δT (L) is the time required to pass from L to L+1 lattice sites (so this quantity is in fact a first passage time). We define a lattice site as excited when its local energy exceeds some border E B = 10 −50 . The actual value of E B was chosen arbitrarily, but any other value, e.g. E B = 10 −100 would produce similar results. The quantity δT can be interpreted as a propagation time for L → L + 1, it can be determined for each particular realization of disorder and initial condition. After having the ensemble of these propagation times at given L, we calculate the average propagation time ∆T as the geometric average of δT (equivalently, we average the logarithms ∆T = exp[ log δT ]).
This second approach is superior to the measurement of the effective spatial extent L for two main reasons: (i) First, ∆T has no explicit time dependence, hence any prehistory does not appear in later measurements of ∆T . It is therefore easier to compare different realizations and simulations for different parameter values and different initial conditions using ∆T . (ii) The second advantage relates to the procedure of averaging over many realizations of the spreading trajectories. By averaging ∆T (L) for fixed L, we average over situations with the same energy density w = E/L. If, in contrast, different realizations of L(t) are averaged for a fixed time t, situations with different energy densities w are averaged together, which is not reasonable if the density w is the crucial parameter on which the properties of the propagation should depend. This is schematically sketched in Fig. 2b.

Fractional nonlinear diffusion equation
We study spreading which is induced by nonlinear chaotic interactions between oscillators, it disappears in the integrable (linear) limit. It is known for Hamiltonian systems that chaos might lead to diffusive behavior which can be understood as the result of "intrinsic stochasticity" induced by the chaotic motion [27,28]. In former works, the nonlinear diffusion equation (NDE) was introduced and remarkable similarities between its scaling properties and the spreading behavior and numerical results for strongly nonlinear lattices were found [18,23,29]. The nonlinear diffusion equation describes the spatio-temporal evolution of a density ρ(x, t) with a density dependent diffusion "constant" D(ρ) ∼ ρ a : The main idea for introducing such a macroscopic description is the hypothesis, that the average spreading of the energy in nonlinear Hamiltonian systems of type (1) follows this NDE. Thus, one identifies ρ(x, t) = w k (t) where k is understood as a discretized spatial coordinate and the averaging · is typically taken over ensembles of trajectories and time intervals. The essential prediction from the NDE is the oneparameter scaling of spreading with the nonlinear exponent a as the only parameter [18], which has been successfully tested by numerical studies in several cases of strongly nonlinear Hamiltonian systems [18,23,29]. The motivation for assuming a density dependent diffusion constant D(ρ) in the NDE above was that the strength of chaos in the Hamiltonian lattices decreases with the energy density. This also leads to a reduced stochasticity and thus also the diffusion constant should decrease when the energy density gets smaller. From the purely power-law nonlinearities in the Hamiltonian system it is natural to assume a power-law dependence for the diffusion constant D(ρ) ∼ ρ a . Diffusive behavior in the phase space, induced by chaos, has been studied also in low-dimensional Hamiltonian systems. There, anomalous transport might occur due to the mixed phase space structure with regular island in a chaotic sea. Chaotic trajectories might feel remainders of the destroyed integrability close to such regular islands which leads to so-called "accelerator modes" [30,31]. By analyzing the self-similarity of the structure of regular islands it was found that the diffusion process should be more precisely described by the fractional diffusion equation (FDE) [32]: where ∂ γ /∂t γ denotes the fractional derivative of order γ > 0 in the Caputo sense, defined later. This fractional time derivative introduces a memory effect and thus accounts for the sticking of trajectories to surviving integrable tori in the mixed phase space.
There is no general reason why such an effect should not be seen in the strongly Hamiltonian lattices discussed here. The phase space of coupled harmonic or nonlinear oscillators might also exhibit islands with integrable trajectories and thus possibly give rise to phenomenon describable by a fractional diffusion equation. To account for both effects, the reduction of chaoticity due to a decreasing density and the possibly mixed phase space, we introduce here the fractional nonlinear diffusion equation (FNDE) as a phenomenological model to describe the spreading process in nonlinear Hamiltonian systems (1): As above, ∂ γ t denotes the Caputo fractional derivative, defined as: with n = ⌈γ⌉ ∈ N being the smallest integer with n > γ.

Scaling properties of the FNDE
In the following, we will analyze the FNDE to deduce its scaling predictions for spreading states. Our analysis will closely follow previous considerations of the normal NDE (8), where the source-type solution can be found explicitly from a self-similar ansatz [33]. First, we look at the scaling properties related to a change of the conserved quantity E. Therefore, we assume that ρ(x, t) is a solution of (10). We rescale this solution to find a new solutionρ(x,t) using a scaling parameter b: with a scaling exponent α such thatρ is again a solution of (10). A straight forward substitution ofρ into the FNDE, defining △ x := ∂ 2 x , gives the terms: and This is a scaled version of (10) if α = − a γ . If we now set b = 1/E we can find the scaling relation of the time that is implied when reducing a solution with arbitrary energy E to the normalized caseẼ = 1, namely:t Note, that this result is compatible with previous findings for the usual NDE (γ = 1), where one indeed findst = E a t. It means that for any solution ρ(x, t) with arbitrary energy E, time and energy always have to appear in the combination above (15). For both, the nonlinear diffusion equation (γ = 1) and the linear fractional diffusion equation (a = 0), one finds source-type solutions by using a self-similar ansatz. It is therefore natural to expect that this ansatz would also be successful for the FNDE (10) considered here. Thus, we use the self-similar ansatz: to identify some scaling properties of the nonlinear fractional diffusion equation. We start with demanding the conservation of energy: Hence, we conclude µ = ν in the self-similar ansatz, because the r.h.s. has to be independent of time. Considering the FNDE directly, one finds for the r.h.s. of (10): The fractional derivative can be evaluated as: where we use x = yt µ and introduce the integral F (y): Using these expressions, the FNDE for this self-similar ansatz gives: Thus one is left with a closed integro-differential equation for f (y), if the scaling exponent is set to: Here, we will not look further at solutions f (y), but rather suppose that such a solution exists. We note that this result is consistent with self-similar solutions for the linear fractional NDE (a = 0) [34], where the scaling was found to be µ = γ/2. The scaling properties of such a solution then imply predictions on the spreading, namely L ∼ t µ , where L is some length scale of the spreading state, e.g. the width. Using the result from above (15) one can also deduce the correct energy scaling of this spreading law L ∼ (E a/γ t) µ , which gives the following scaling prediction for spreading: Solving for t and taking the derivative with respect to L, one also finds a scaling prediction for the excitation times: Both results resemble the relations for the NDE with γ = 1 reported earlier [18,29] and summarized in the next section.

Self-similar solution of the NDE
For the FNDE it is at the moment unclear if the profile of the self-similar solution f (y) can be found analytically by solving (19). For the usual NDE where γ = 1, the ordinary differential equation for the scaling function f (y) is much simpler [29]: with µ = 1/(a + 2) as above. This ODE can indeed be solved explicitly which leads, going back to the original variables ρ(x, t), to the following self-similar solution of the NDE [35]: with c being a constant of integration which follows from the energy conservation: . X(t) is the sharp front of the field and has the following time dependence:  The solution has sharp edges (see Fig. 3) and its spatial extension is given by X. The size of the wave packet grows in time as a power law, which can be represented in a scaling form as In order to get rid of the undetermined time offset t 0 , we calculate, following [18], the local inverse velocity of the spreading as dt/dX and obtain for it the following scaling law: Identifying the excitation edge X with the spatial extent L before, one sees that this indeed corresponds with the scaling result for the FNDE above (21), (22) for γ = 1.

Implications for spreading in lattices
In our numerical simulations of strongly nonlinear lattices, one approach is to calculate the characteristic size of the wave packet by appropriate averaging of the entropies (7).
In particular, we can directly attribute the size of the field at a given time to the participation number, so that in (21) L ∼ L. Thus, if we assume that the NDE provides an adequate description of the spreading in strongly nonlinear lattices, the spreading data for different energies should fulfill scaling (21), where the constants γ and a depend generally on the powers κ, λ: Similarly, in the second method we calculate the average time ∆T needed for spreading, in dependence of the field spatial extend L. This quantity directly corresponds to the inverse velocity in (22): dt/dL ∼ ∆T . Thus, we expect that the times ∆T behaves as: ∆T We note that in both predictions, (28) and (29), the two influences from the fractional derivative γ < 1 and the nonlinearity a > 0 can be nicely separated. The energy scaling in the l.h.s. of (29) is solely dependent on γ, so first by identifying the energy scaling in the numerical results one can determine γ. The power law of the subdiffusive spreading than determines the nonlinearity parameter a. We already note here that in some cases we numerically find a density dependent nonlinear exponent a(w) where w is the energy density w = E/L [18].

Results
In the following sections we report on extensive numerical simulations of strongly nonlinear lattices, trying to check the predictions of the NDE framework. For the numerical time evolution we used a 4th-order symplectic Runge-Kutta scheme [36,37], mostly with step-size ∆t = 0.1.

Homogeneous nonlinearity
Scaling induced spreading prediction. We start with the case of homogeneous nonlinearity κ = λ in (1), where the local and coupling potential have the same nonlinear power: Here, β is the parameter determining the relative strength between local and coupling potential and the total energy can be rescaled to unity as described in section 2.2 and is thus not a free parameter in the equations. We can find a relation between the order of the fractional derivative γ, the nonlinearity of the FNDE a and the parameter κ for this homogeneous case. Indeed, from the scaling invariance of the Hamiltonian in (4), we obtain that the time scales with energy as: On the other hand, the FNDE obeys the scaling relation: t−t 0 ∼ E −a/γ (15). Motivated from previous results, we assume that the FNDE gives a correct macroscopic description of the spreading process. If this assumption is true, then the spreading states have to fulfill both scaling relations above, which gives the nonlinearity parameter a as a function of γ and κ: To get an exact spreading prediction, one still has to obtain the parameter γ of the fractional derivative that is introduced to account for the mixed phase space of the system. However, in the following we will consider situations of large perturbations where it is reasonable to assume that the phase space is fully chaotic [38]. Large perturbation means that the coupling parameter is of the order β ≈ 1. In this case we expect that the perturbation is strong enough to destroy all remainders of the integrability for β = 0, and thus we make the assumption that γ = 1. This then gives the following spreading predictions: These exact relations will serve as a test for our assumption that the NDE adequately describes the spreading in nonlinear Hamiltonian lattices (1).
Comparison with numerical results. To test the theoretical predictions (33) we follow the evolution in a one-dimensional lattice with ω k ∈ [0, 1], started from a single site (or several sites for κ = 6) excitation in the middle. For several values for the nonlinear strength β = 0.25, 0.5, 1, 2 we integrated the system up to T = 10 6 and analyzed the spreading in terms of P (t) and ∆T (L). This was repeated for M = 100 realizations of disorder. Fig. 4 shows the averaged results of these runs for the excitation time ∆T (L) for κ = 4 (left panels) and κ = 6 (right panels). In both cases we find a quite nice agreement of the numerical results with the analytic predictions of the NDE (33). We also performed simulations choosing the disorder to be ω k ∈ [0.5, 1.5] and got similar results (not presented here). Additionally, results for the direct spreading measure P (t) were obtained which show the same agreement with prediction (33) and are also omitted here [18]. To our opinion, the agreement between numerics and prediction is a rather   The dashed lines shows the NDE prediction P ∼ t 4/9 and P ∼ t 3/7 respectively. In the inset we plot the numerical spreading exponent ν obtained from finite differences of the method above, the dashed line there also corresponds to the expectation from the NDE.
convincing evidence that the NDE provides the proper framework to model the average energy spreading in nonlinear Hamiltonian chains.
Remarkably, the NDE scaling (33) holds also for one-dimensional lattices without disorder. Fig. 5 shows the participation number evolution for lattice (5) with ω k = 1 and κ = 4, 6. In this case the excitation times ∆T are not the proper measures as they are dominated by quasi-compactons (cf. Fig. 1(a)), but the participation number calculations are insensitive to such modes. For them we expect from (33) the scaling P ∼ (t − t 0 ) 2κ 5κ−2 which is confirmed in Fig. 5. This is also supported by the direct comparison of spreading states from the numerical simulation to the self similar solution, as can be seen in right panel of Fig. 3. Note that in the end of the simulation for β = 1, 2 in Fig. 5a, the spreading state has hit the lattice boundaries leading to a saturation of the participation number and a decrease of the spreading exponent.
Summarizing these results, we have found that from the assumption of the validity of the NDE we derived an exact spreading predictions for a fully chaotic phase space in Hamiltonian lattices with homogeneous nonlinearity κ = λ. These predictions were to a high accuracy verified as the asymptotic behavior in numerous numerical simulations. We note that this spreading process can also be observed in the case of a regular onsite potential were disorder is completely absent, Fig. 5. This shows that disorder is not required for the spreading phenomenon, an observation already made for 2D lattices in [23]. Hence the subdiffusive spreading is not a result of the interplay between nonlinearity and disorder, but rather a more general phenomenon lately called "Chaotic Diffusion" [23,29]. To further verify the assumption of γ = 1 above due to the fully chaotic phase space it would be very interesting to study the behavior for smaller β. If our argument is correct one would expect some dependence γ(β) where γ also decreases     for smaller values of β. This will be the subject of future studies.

Nonlinear Oscillator, Nonlinear Coupling
Numerical Results. After having found that the NDE provides a good description of the spreading for the case of homogeneous nonlinearity, we turn now to a general situation with κ = λ. In this section we study the case of fully nonlinear oscillators, hence we choose κ = 4 and λ = 6, 8. In this case, the disorder parameter ω k in (3) does not have the meaning of an oscillator frequency, but is the coefficient determining the nonlinear strength. The real frequency of oscillations depends on the local energy at the site. In Fig. 1 we show an exemplary time evolution of an initially localized state in such a lattice. For this non-homogeneous case, the energy E is the crucial parameter in the Hamiltonian. That allows us to independently determine the parameter γ and a from numerical simulations by first identifying the energy-scaling of the spreading and then computing the slope of the subdiffusive process. That means we will compare the numerical results with the spreading predictions from (28) and (29). We start with the case κ = 4 and λ = 6 and investigate the excitation times ∆T (L) for different total energies. The results of our simulations for ω k ∈ [0.5, 1.5] are shown in Fig. 6. In the right panel the scaling as suggested by the FNDE (29) is applied and we found the best overlap of the individual curves for the parameter value γ = 1.08 that gives the scaling exponent 1 − 2/γ ≈ 0.85. That means we find only a slight deviation from the pure NDE case where γ = 1, hence the influence of the mixed phase space is rather small, but clearly identifiable as for γ = 1 the curves do not overlap as perfectly (cf. [18]). The numerical data also nicely follow a straight line as seen in Fig. 6b. This indicates subdiffusive behavior with a slope (a + 2 − β)/β ≈ 2.6 from a numerical fit   and we thus calculate the nonlinear exponent in the FNDE as a ≈ 1.8.
In a second simulation we studied the case κ = 4 and λ = 8. The results are shown in Fig. 7 where two scalings with γ = 1 (Fig. 7a) and γ = 1.18 (Fig. 7b) are compared. It is clear from these graphs that the normal NDE with γ = 1 does not predict the correct scaling as the curves for different energies do not overlap in Figure 7a. For γ = 1.18, the scaled variables according to the FNDE are ∆T /E 0.7 vs. L/E and Fig. 7b shows that for this choice indeed a convincing overlap of the individual curves is observed. In contrast to the case λ = 6, the scaled curve for λ = 8 does not follow a straight line. We explain this by the fact that we have not reached the asymptotic regime yet in this study. Indeed, the numerically accessible parameter range for the energy density E/L goes only down to E/L ≈ 10 −3 in Figure 7b, while for λ = 6 we were able to go almost two orders of magnitudes lower. To still quantify the slope in this case we performed a polynomial fit of the scaled data and plotted the derivative of this fitted curve in the inset in Figure 7b. The result indicates a convergence of this slope and hence an asymptotically constant value for a ≈ 3.5. The dashed line in this inset represents our theoretical prediction for this asymptotic value to be explained in the next subsection.
In summary, we have found here that for fully nonlinear oscillators, κ = 4 and λ = 6, 8, the spreading of initially localized excitations can be nicely described by the FNDE. With the scaling approach we were able to separate the two parameters γ and a of the FNDE and determine their values from the numerical results on the excitation times. The power of the fractional derivative was obtained as γ = 1.08 for λ = 6 and γ = 1.18 for λ = 6. For λ = 6 the scaled spreading was found to behave as a power law with some slope (a + 2 − γ)/γ ≈ 2.6 which gives the nonlinearity parameter of the FNDE as a ≈ 1.8. For λ = 8 we could not reach the asymptotic behavior and hence found a density dependent slope, but a numerical estimation of this slope indicates for a convergence against the value (a + 2 − γ)/γ ≈ 3.7 which means a ≈ 3.5. We conclude that the FNDE is a good model to describe spreading in fully nonlinear Hamiltonian systems.
Microscopic spreading dynamics. In [23] a microscopic model of spreading was developed that lead to an exact prediction of the spreading exponent for a regular twodimensional lattice (ω k = 1) of harmonic oscillators with nonlinear coupling. Here, we will follow this idea and try to find a reduced system that describes the dynamics at the excitation edge. The idea is to understand the mechanism of how a new oscillator is excited from the chaotic forcing induced by its already excited neighbor.
For harmonic oscillators with ω k = 1, the situation was particularly easy because all oscillators were in resonance due to the absence of disorder. Therefore, in [23] it was enough to consider only two coupled oscillators at the edge, one excited and one at rest, to obtain a correct spreading prediction. For the nonlinear oscillators with κ = 4 studied here it is immediately clear that considering only two oscillators will not be sufficient. The Hamiltonian for two coupled oscillators is: where q 1 , p 1 denote the already excited oscillator with a local energy density w ≈ p 2 1 /2 + q 4 1 /4 while the second oscillator is at rest: q r = p r = 0 with a zero energy density w r = 0. Because the second oscillator is subject to a non-resonant forcing it will, for small energy densities w, only become excited up to an energy density according to standard perturbation expansion which means w r ∼ w λ/4 ≪ w. Hence, for small densities there is almost no energy transport from the excited to the resting oscillator which would imply that spreading should stop because no new oscillators get excited. This prediction is clearly wrong as is seen from numerical spreading results presented above. The reason is that the two oscillator model is too simple to describe the spreading process. Thus, we consider more complex situations with N oscillators, were the first N − 1 oscillators are excited with some energy density w, while the last oscillator is at rest. From examining the geometric properties of the resonances of such coupled nonlinear oscillators it can be argued that only for N ≥ 5 energy transport to the last oscillator through a global chaotic layer is expected [29].
Here, we will verify this conjecture by a numerical simulation. Consider a situation with N − 1 excited oscillators with an energy density w as described above. One way to quantify the energy transport to the last, resting oscillator is by measuring the time T that is required for this oscillator to become excited to some critical energy density above the perturbative description. This time is very similar to the excitation times introduced earlier to quantify spreading. Here, we will fix the number of excited oscillators and just measure the time as a function of the energy density T (w). If T diverges then no energy transfer beyond the perturbative excitation is taking place. In Figure 8  we show the results from a Monte-Carlo study on T (w) for an ensemble of M = 100 random initial conditions and different numbers of oscillators N = 2 . . . 7. The bold lines correspond to the maximum times max T from this ensemble of initial conditions for each N and density w. In both cases, λ = 6 ( Figure 8a) and λ = 8 (Figure 8b), one definitely observes a divergence of T for N < 5. For N ≥ 5, however, we found an asymptotic power-law dependence T (w) ∼ w χ with χ ≈ −1.7 for λ = 6 and χ ≈ −3.0 for λ = 8. So firstly we note that the microscopic model with N ≥ 5 predicts spreading with an asymptotic power-law behavior. To connect these numerical results from the microscopic dynamics to the macroscopic spreading one can identify the microscopic and the macroscopic excitation times ∆T ∼ T . Noting that the number of oscillators in the microscopic remains constant and only the energy density changes one finds the following prediction for the macroscopic excitation time ∆T ∼ E χ . Translating this into the scaled variables used earlier one finds that (a + 2 − γ)/γ = 2/γ − 1 − χ and thus a/γ = −χ. For λ = 6 the nonlinear exponent was calculated from the numerical spreading as a/γ ≈ 1.7, which is in a very good agreement with the microscopic result −χ ≈ 1.7 shown in Figure 6b. For λ = 8 the asymptotic behavior of the macroscopic spreading is also in very good agreement with these microscopic results as −χ ≈ 3 appears to be very close to the asymptotic sprading behavior where a/γ ≈ 3 in Figure 7b. Thus we conclude that a microscopic model of N = 5 oscillators is enough to understand the macroscopic spreading properties in long, macroscopic chains of such oscillators. However, at this point the exponent χ was only obtained from numerical results and analytical treatments remain a challenge for future work.

Harmonic Oscillators, Nonlinear Coupling
Finally, we turn to the most complicated situation of harmonic oscillators with random frequencies and nonlinear coupling. Therefore, we assume the on-site potential to be quadratic κ = 2, for the coupling we chose λ = 4 and λ = 6. This case corresponds to a rather general situation of nonlinear disordered lattices, where in the representation of linear eigenmodes one can also interpret the system as an ensemble of nonlinearly coupled linear modes. The most prominent example of such a situation is the Discrete Anderson Nonlinear Schrödinger Equation (DANSE-model) [7] which, if treated in the eigenmode basis, consists of localized harmonic modes with nonlinear coupling. The main difference between this setup and the strongly nonlinear lattices considered here is that in the DANSE-model the coupling between the modes has random coefficients and is exponentially decaying in space due to the overlap integrals between the localized modes. In contrast, the strongly nonlinear lattices studied here have only a nearest neighbor coupling without a random coupling coefficient. Similar to the studies before, we analyze the excitation times ∆T (L) as function of excitation length L for different energies E to check the predictions of the FNDE scaling (22). At first, we report the results for κ = 2 and λ = 4. Note that some of these results have already been presented in [18]. All results are again averaged over different realizations of disorder and we studied two ways of choosing the random frequencies. Figure 9 shows the results for ω k ∈ [0, 1] and Figure 10 for ω k ∈ [0.5, 1.5]. Both cases are qualitatively very similar. At first, we identify the energy scaling to seemingly follow    the prediction of the FNDE with γ = 1 as seen from the good overlaps in Figures 9b  and 10b. The resulting curves, however, are not straight lines but rather exhibit a clear curvature bending upwards. This has already been reported earlier [18] and is not yet fully understood. Phenomenologically this behavior can be quantified by introducing a density dependent nonlinearity index a: From (22) one finds that for γ = 1 the slope of the rescaled curves is simply given by a(w)+1. Thus we evaluate this slope by means of a polynomial fit and plot the resulting numerical value for a(w) in the insets in Figures 9b and 10b. Qualitatively, there is no difference between the two choices of disorder in Figures 9 and 10, but quantitatively the increase of the nonlinearity index a(w) is faster for ω k ∈ [0.5, 1.5].
In Figures 11 and 12 we show the results of a similar study with the coupling nonlinear exponent λ = 6. The results are again qualitatively the same as above in that we find scaling with γ = 1 and a density dependent nonlinearity index a(w) shown in the insets of Figures 11b and 12b. Hence, this seems to be a universal picture for spreading in lattices of harmonic oscillators with random frequencies and nonlinear nearest neighbor coupling. It should be noted that the density dependent spreading can not be described by introducing a density dependent parameter of the fractional derivative γ(w), because this would mean a density dependent energy scaling which is not observed here. We also note that by introducing a density dependent nonlinearity index a(w) into the FNDE (or NDE as we have γ = 1 here) destroys the self-similar solution and even the scaling prediction. However, the density dependence is found to be very weak a(w) ∼ log 10 w    and thus the rate of change of a is much slower then the spreading time scale. Thus, it is reasonable to treat the energy spreading in a first approximation using a = const and then analyze the slow deviations afterwards. The question of the asymptotic behavior remains, however, open: from the data presented here we cannot judge whether the spreading effectively stops, or continues with an increasing index a, or some transition to another law of spreading (e.g., a logarithmic one) occurs.

Conclusions
Motivated by previous observations of subdiffusive behavior in nonlinear disordered systems and anomalous diffusion in chaotic Hamiltonian systems, we introduced the fractional nonlinear diffusion equation as a phenomenologic model to describe the spreading process in disordered one-dimensional Hamiltonian lattices of nonlinearly coupled oscillators. We have found that with the FNDE it is possible to explain in a consistent way the subdiffusive spreading behavior and the energy scaling of spreading states. Analysis of self-similar solutions of the FNDE not only predicts a subdiffusive spreading, but also induces a scaling of time and energy of the spreading process according to relations (21,22), which depend on parameters γ and a, responsible for the index of the fractional time derivative and of the nonlinearity, respectively. We tested these scaling laws on a class of nonlinearly coupled oscillators with different values of the nonlinear indices κ (local nonlinearity) and λ (coupling nonlinearity). Our main result is that there are three qualitatively different "universality classes" in regard of relations between γ, a and κ, λ. Specifically, we have found the following three cases of nonlinearities that demonstrate different scaling of spreading: (i) For homogeneous nonlinear potentials, where κ = λ, we were able to deduce an exact spreading prediction from the scaling property of the Hamiltonian equations and the FNDE when assuming a fully chaotic phase space. We argued that here the nonlinear diffusion equation with γ = 1, i.e. with normal time derivative and the nonlinearity index a = κ−2 2κ should be applied. This analytic prediction has been confirmed numerically as the asymptotic spreading behavior. As an important result we again note that subdiffusive spreading was also found in the regular case without disorder. This further supports the claim that disorder is not required for the spreading and it indeed seems reasonable to call this process "chaotic diffusion" [23,29].
(ii) In the fully nonlinear case with local nonlinearity index of the oscillators κ = 4 and the nonlinearity indexes λ = 6, 8 in the coupling, we have found that the numerical spreading results follow the energy scaling as predicted from the FNDE with the fractional time derivative of order γ = 1.08 (for λ = 6) and γ = 1.18 (for λ = 8). This is compatible with previous findings on anomalous diffusion in low-dimensional Hamiltonian systems were the mixed phase space also leads to a fractional diffusion equation with γ > 1 [39]. Furthermore, for this case we were able to construct a microscopic model of the dynamics at the excitation edge that predicts the correct spreading behavior verified in direct numerical simulations.
(iii) In the case of nonlinearly coupled harmonic disordered oscillators, we have verified that the energy scaling follows nicely the prediction of the normal nonlinear diffusion equation (fractional order γ = 1). However, the spreading does not follow a pure power law as predicted by the NDE. Instead, we have identified a remarkable dependence of the effective index of nonlinearity of the FNDE on the energy density a(w). In all cases considered we have observed that a increases as w becomes smaller, although the particular profiles of a(w) depend on the nonlinearity in coupling and on the disorder. As the effective nonlinearity increases in the course of spreading, this means a slowing down of the spreading process compared with the perfect power law, as in this case L ∼ t 1 a(w)+2 . Unfortunately, we are not able to present a microscopic model of the edge dynamics at this point, mainly due to the highly complicated resonance structure that emerges when considering nonlinearly coupled harmonic oscillators with random frequencies. Consequently, it is also not possible to judge from the data what is the asymptotic behavior of the spreading for times beyond those available in our numerics.
Our findings rely to a large extent on the novel quantity characterizing the spreading, the averaged excitation time introduced in [18]. This quantity is defined for a particular size of the wave packet, and thus for a particular value of the density. It thus allows us to reveal the density dependence of the spreading characteristics, what is hardly possible with old approaches where, e.g., averaged participation numbers have been followed. Unfortunately, the calculation of averaged excitations is relying on the sharp edges of the field, so its application to linearly coupled lattices where eigenmodes are exponentially (but not sharp) localized, remains a challenge for future studies. We stress once more that in our study we consider the fractional nonlinear diffusion equation as a phenomenological model guiding the scaling relations of the problem. Its derivation from the microscopic model appears, at the present stage, as a complex, not yet resolved problem. In this respect we refer to paper [40], where an attempt to derive a nonlinear diffusion equation for the two-dimensional disordered nonlinear Schrödinger equation is presented; the resulting conclusion on the linear in time growth of the variance of the wave packet (like in normal diffusion) does not, however, correspond to numerical findings of subdiffusion in this model [41]. Further attempts are necessary to resolve this problem.