Hyperballistic superdiffusion and explosive solutions to the non-linear diffusion equation

By means of a particle model that includes interactions only via the local particle concentration, we show that hyperballistic diffusion may result. This is done by findng the exact solution of the corresponding non-linear diffusion equation, as well as by particle simulations. The connection between these levels of description is provided by the Fokker-Planck equation describing the particle dynamics.

By means of a particle model that includes interactions only via the local particle concentration, we show that hyperballistic diffusion may result. This is done by findng the exact solution of the corresponding non-linear diffusion equation, as well as by particle simulations. The connection between these levels of description is provided by the Fokker-Planck equation describing the particle dynamics.

PACS numbers:
Superdiffusion is characterized by the fact that the root mean square displacement of some kind of particles, increases with time t as r rms ∼ t τ with the exponent τ > 1/2, the normal diffusion value being τ = 1/2. This behavior may arise in physical, biological or geological systems; examples include Levy flights [1,2], particle motion in random potentials or the seemingly random paths of objects moving in turbulent flows [3,4].
Biological examples may be found in the foraging movement of spider monkeys [5] and the flight paths of albatrosses [6,7]; in both cases τ ≈ 0.85. These movements are Levy walks, which are random walks of uncorrelated steps of length δx, that take their value from a distribution p(δx) ∼ 1/δx µ+1 . They result in superdiffusive behavior with r rms ∼ t 2/µ when 0 < µ < 2 [1].
However, the mere observation that the step length distribution has a fat tail, does not by itself provide any physical model to explain the superdiffusive behavior. The simplest physical example of superdiffusion is perhaps provided by the undamped Langevin equation which describes a random walk in momentum space and a corresponding real space displacement with τ = 3/2 [8]. This kind of behavior is termed hyperballistic as τ > 1. Quantum-or classical particles in random potentials behave much like those described by the undamped Langevin equation, and yield hyperballistic diffusion with τ = 3/2 [8] too, though Golubovic et al. [9] studied a case where τ = 9/8. In optical experiments [10,11] where the spatial coordinate in the direction of the light plays the role of the time coordinate, hyperballistic spreading has been observed as well. This effect is linked to Anderson localization [12], and comes from a transition where the light modifies its mean free path as it passes through the medium.
Hyperballistic diffusion seems almost a contradiction in terms, for how could a random walker move faster * Electronic address: flekkoy@fys.uio.no than a directed walker that never changes direction? The explanation lies in the fact that the velocity, and thus the step length, keeps increasing with time without limits. This behavior is of course unphysical in the context of the Langevin equation as there will always be dissipative forces that match the fluctuations, but has a physical basis in random potentials. Without diverging velocities or step lengths, long range time-correlations are required for superdiffusion, an example being the elephant random walk, so named because both the walkers and elephants have long memories, which in the model give rise to (subballistic) superdiffusion [13].
Generally, superdiffusion has been modelled by independent agents interacting with an environment, or possessing a long term memory [14]. The main question of the present letter is if superdiffusion, including the hyperballistic case, could result directly from a Markovian description of particle interactions. Such interactive systems could include crowds of people, bacteria swimmers competing for food [15,16] or the evolution of the porosity in a granular packing. For the purpose of addressing this question we investigate the potentially simplest description of particle interactions, namely, that where a conserved concentration C of particles is governed by Ficks law j = −D(C)∇C. Here the C-dependence in D reflects interactions between the particles; in many cases of interest these interactions are well captured by this type of mean field description. Already in 1959 Pattle [17] solved the diffusion equation where C = C(r, t) is the concentration and D is given by the power law D = D 0 (C(r, t)/C 0 ) −γ where C 0 is a constant reference concentration, D 0 is the diffusivity at that reference value, and the exponent γ < 0. Pattle found the root mean square displacement r rms (t) ∼ t τ with arXiv:2012.03748v1 [cond-mat.stat-mech] 3 Dec 2020 where d is the dimension. For negative γ this will always lead to sub-diffusion. We have recently shown that in d = 1 there are exact solutions with positive γ as well [18], which still satisfy Eq. (2), thus yielding superdiffusion with 1/2 < τ < 1 as γ < 1 always. In the present article we take this result further by deriving the solution for C(r, t)and r rms (t) for γ > 0 in any dimension. When d ≥ 2 the corresponding exponent τ will then take on any value, including those of the hyperballistic regime, implying that hyperballistic diffusion is a higherdimensional effect. We coin the term 'explosive' for the corresponding time dependence of C(r, t) because the decay of an intially localized C-profile is qualitatively faster than normal diffusive, or even superdiffusive, decay.
To validate the mean field description and provide it with a physical basis, we introduce a particle model that is described by Eq. (1). The step lengths in this model ∼ C −γ/2 , and therefore correspond to velocities that diverge as C → 0. This would correspond to an unlimited access to thermal energy. However, unlike the Langevin equation where τ = 3/2 [8], this model can produce any τ -value.
Following the same lines as in [18] we rewrite Eq. (1) as Hence, we see that we need γ < 1 for the equation to be defined when C(r, t) = 0. The initial condition at t = 0 is a point source pulse containing N p particles, C(r, 0) = N p δ(r). This means that there is no intrinsic length-or time scale in the problem, and the particle number N (r, t) inside a radius r should satisfy the scale invariance condition N (r, t) = N (λr, h(λ)t) for some h(λ). Differentiating this equation with respect to r, using the fact that dN (r, t) ∝ C(r, t)r d−1 dr leads to the scaling relation We are free to chose λ such that h(λ)t = 1, which yields λ as a function of t, say λ = 1/f (t), and where we have introduced p(y) ≡ C(y, 1) and the reduced variable for some dimensionless constant c, which can be absorbed in the definition of f (t). The point-like initial condition, implies f (0) = 0, and the left hand side of Eq. (6) can be easily integrated to give Note that this form immediately gives with τ given by Eq. (2). From equation (??), we also have an expression for p(y), which can be integrated to give, For Fick's law to be valid throughout the domain, C(r, t), and therefore, p(y), must be differentiable everywhere.
To avoid a spike at the origin we must have p (0) = 0 and also a finite p(0), which implies that K = 0. So, Eq. (10) may be integrated to yield where k is an integration constant. This expression is independent of the dimension d. The value of the constant k can be determined through the normalization, dV C(r, t) = N p , which gives and yields the concentration field by means of Eqs. (5), (7) and (11). The mean square displacement is given by which is limited to the range of γ-values where the integrals in Eq. (8) converge. Since r d+1 C(r, t) ∼ r d+1−2/γ for large r this range is 0 < γ < 2/(d + 2). However, in any particle simulation there will always be a largest particle position r max that will act as a cut-off. This means that the t 2τ factor in r 2 rms survives, but that its prefactor will fluctuate with the r max -value. The behavior with different d and γ is summarized in table I.
We will employ two simulation models, both in d = 3 with N p random walkers, labeled i, that have positions r i → r i + δr i . The particles interact only via the value of C, which is the local population density. The steps are chosen isotropically at each time step; to find their length we need to derive the appropriate Fokker-Planck equation and match it to Eq. (3). For every time step ∆t the walkers move τ > 1 rrms prediction converges C is normalizable d=1 never γ < 2/3 All γ < 1 d=2 γ > 1/2 γ < 1/2 All γ < 1 d=3 γ > 1/3 γ < 2/5 γ < 2/3   where α is a Cartesian index and the function g(C) is to be determined. This defines a Wiener process with η as a random variable with η = 0 and η 2 = 1. Now, following the same steps as in [18,19] we use the standard Chapman-Kolmogorov, or master equation, to derive the following Fokker-Planck equation for the particle concentration C(r, t) ∂C(r, t) ∂t = 1 2 ∇ 2 (a 2 (r)C(r, t)) .
Here a 2 (r) is the mean squared jump length per time, where W (r, x) is the probability per unit time that a walker jumps a distance x from r. Setting g(C) = bC −γ/2 gives and requiring equivalence with Eq. (3) thus implies that γ)). This leads to the step where the random variable η is given above. This defines the particle model that is described by Eq. (3). In the finite interaction range model C is calculated by assuming a maximum interaction range ∆x between particles. This is done by calculating C onto a lattice with lattice constant ∆x: The local value C(r n , t) at the discrete site r n is simply 1/∆x d times the number of particles at positions x i that satisfy |x iα − x nα | < ∆x/2. The step length for a particle that is located at x depends on the C-value at the nearest lattice site. The finite interaction range of this model has a discretization effect: Once C is so small that there is only one-or zero particles in each ∆x-cell, the step length will always be the same, and as a result, there will be a cross-over to normal r rms ∼ t 1/2 diffusion, an effect that is observed in the γ = 0.35 curve of Fig. 3 (a). The other, infinite interaction range model employs no lattice at all, but evaluates C at any particle position x as C(x, t) = N r /V r (x) where N r ∼10 is a fixed particle number and V r (x) is the volume of the sphere that contains N r nearest neighbors. There is no upper limit to the size of V r (x), and it is in this sense that the model has a potentially infinite interaction range. When γ = 0 this model will never cross over to normal diffusive behavior.
In Fig. 1 the analytic solution of Eq. (11) is plotted for different γ-values. The term 'explosive' seems an appropriate label for the behavior of the concentration for two reasons: First, as γ → 1/2 close to the critical value of 2/3, the initial concentration C(0, 0) drops by more than 10 orders of magnitude in the same time that the negative γ solutions (taken from Pattle [17]), drop by less than 2 orders. Second, the divergence of the integral in Eq. (8) defining r rms (t) signals a cross-over to a regime where the break-away particles dominate the r rms (t)-behavior at ever increasing step lengths.
In Fig. 2 the data collapse anticipated in Eq. (5) is seen to be satisfied. Figs. 3 (a) and (b) demonstrate that the particle displacement is in fact characterized by Eq. (13), the difference between Fig. 3 (a) and (b), being that the first figure compares simulations and the full analytic prediction of Eq. (13), while the hyperballistic transport shown in Fig. 3 (b), only confirms the prediction of the τ exponent, Eq. (2). Note that in Fig. 3 (a) the convergence to the prediction of Eq. (13), happens over a time that increases with γ, signalling the end of the regime where r rms (t) has an exact analytical expression. Figure 4 summarizes this comparison for the full range of relevant γ-values, using the finite-range model for the smaller-and the infinite range model for the larger γ-values.
In conclusion, we have shown that particle interactions described entirely in terms of their local concentration will yield superdiffusion, and even hyperballistic diffusion when d ≥ 2. In Pattles classical 1959 paper [17] the γ < 0 solution of Eq. (1) is not actually derived, but only written down. So, for completeness we derive it here along the same lines as those leading up to Eq. (11). In the solutions thus derived C(r, t) has a final support outside which it is strictly zero. For any γ the normalization where we have used the isotropic nature of the problem to perform the angular integration and thus introduced the geometric factor C d = 1, 2π, 4π when d = 1, 2, 3.
We see from equation (11) that, for γ < 0, the domain of the probability density, p(y), is limited to y < y c = 2k(γ − 1)/γ, so that the normalization condition is yc 0 dy y d−1 p(y) = N p C d .
yielding the normalization constant