Anomalous Diffusion in Systems with Concentration-Dependent Diffusivity

We show analytically that there is anomalous diffusion when the diffusion constant depends on the concentration as a power law with a positive exponent or a negative exponent with absolute value less than one and the initial condition is a delta function in the concentration. On the other hand, when the initial concentration profile is a step, the profile spreads as the square root of time. We verify our results numerically using particles moving stochastically.

Here C 0 is a constant reference concentration and D 0 is the diffusivity at that reference value. We e.g., find such behavior when compressible gases flow through porous media [1,2] or in filtration processes [3]. Another example is heat diffusion at high temperature [4,5]. Population dynamics gives rise to this kind of behavior [6][7][8], as does water ingress in zeolites as studied by where δϕ 0 is a reference excess porosity level. Comparing with Eqs 1 and 3, we see that c 1 here. There are two classes of examples in the list above: systems where there is a physical quantity that is diffusing due to random motion, and systems that are only described by the diffusion equation. The diffusion of excess porosity in the packing of beads is an example of the first class. On the other hand, the motion of a fluid-fluid interface under gravity in a porous layer [13] is an example of a system where there is no random diffusion.
In either case, there being physical diffusion or the system may be mapped onto equations having the form of a diffusion equation, we may consider an equivalent system consisting of diffusing particles. The concentration-dependent diffusivity may then be interpreted as mean-field type interaction between the particles, where the behavior of each particle is determined by the number of neighbors it has.
A key quantity to characterize diffusive behavior is the root mean square (rms) displacement traveled by the particlesi.e., random walkersas a function of time. If a random walker in one dimension starts at position x 0 when time t 0, then the subsequent rms displacement is Refs. 17, 18 When the diffusion exponent 0 < τ < 1/2, we are dealing with subdiffusion and when 1/2 < τ ≤ 1, we are dealing with superdiffusion. When τ 1, the diffusion is ballistic. When τ > 1, the term hyper-ballistic is used [19,20]. These are all examples of anomalous diffusion. Normal diffusion occurs when τ 1/2. There are several possible ways the behavior of the diffusivity in the diffusion equation may lead to anomalous diffusion. It may do so by having a position dependent diffusivity, see e.g., Refs. 21-24. The diffusivity we are considering in this work also depends on the position, but only through the concentration as in Eq. 3.
There is a relation between the exponent γ defined in Eq. 3 and the random walk exponent τ in Eq. 5, This equation was first worked out by Pattle [25] for negative γ indicating sub diffusion. We will in this paper extend this equation to positive γ, using the self-similarity approach [26][27][28]. We find that Eq. 6 is valid for the range c < 2/3, and that indeed there is super-diffusion for positive γ. This result comes from a closed form analytic solution that exists for the full range c < 1. In the 2/3 < c < 1 range the integral defining the second moment 〈x 2 〉 diverges. We furthermore note that hyper-diffusion cannot occur in one dimension if the concentration-dependent diffusivity Eq. 3 is at play. In spite of the relation Eq. 6 between γ and the exponent that defines anomalous diffusion, τ, it is sometimes believed that the Boltzmann transformation [29] demonstrates that there is no anomalous diffusion, even when the diffusivity depends on the concentration. This is, e.g., demonstrated in the famous textbook by Crank [30]. However, the random walk exponent τ, depends on the initial condition Eq. 2.
We support our analytical work by numerical simulations through the introduction of a particle model where random walkers move in discrete steps that depend on the local concentration. We show that the average description of this model, that is given by the Chapman-Kolmogorov equation, reduces to Eq. 7. Hence, we are able to reproduce the analytical results by particle simulations. These particle simulations are easily generalized to higher dimensions and different kinds of non-linearities in the diffusion equation that describes them.

ANALYTIC SOLUTIONS
The analytic results are most conveniently obtained by rewriting Eq. 1 in the form [31] Frontiers in Physics | www.frontiersin.org December 2020 | Volume 8 | Article 519624 Hence, we see that we need c < 1 for the equation to be defined when C(x, t) 0.

Boltzmann Transformation and the Step
In order for the Boltzmann transformation to be applicable we take the initial conditions to be C(x, 0) C 0 Θ(−x), where Θ(x) is the Lorentz-Heaviside function. The Boltzmann transformation consists in introducing the variable When the t and x derivatives are transformed to y-derivatives the diffusion Eq. 7 becomes the ordinary differential equation, with the x and t dependence through y only. Now, the initial condition too may be written in terms of y alone: For t > 0: C(y) C 0 Θ(−y). For this reason the solution of the diffusion Eq. 7 takes the form C(x, t) q x t 1/2 (10) for some function q that satisfies Eq. 9. This immediately leads to the conclusion i.e., that the diffusion is normal with τ 1/2 in Eq. 6. In other words, the step function initial condition cannot lead to the anomalous diffusion behavior defined by Eq. 6 and Pattle's solution. In the following we shall se that this conclusion is qualitatively changed by the introduction of a localized, and thus normalizable, δ-function inital condition.

Exact Solutions From Delta-Function Initial Condition
Solving Eq. 7 with a point source initial condition turns out to change the diffusion exponent τ. Hence, we start with where N p is the particle number, and δ(x) is the Dirac δ-function, which obviously, yields the normalization +∞ There is no intrinsic length or time-scale in Eq. 7 since the diffusivity depends on C through the power law of Eq. 3. This means that as long as boundary-or initial conditions do not introduce such scales either, the solutions C(x, t) must be scalefree too. More precisely, if x → λx, then there must be some rescaling of time t → t/f −1 (1/λ) so that the probability of finding the particle remains unchanged, that is, This ensures that the normalization of Eq. 13 remains constant with time. We now choose λ so that t/f −1 (1/λ) 1. That is, we set Combined with Eq. 14, this gives , where we have set p(y) ≡ C(y, 1). Since C(x, t) has dimension 1/ length, we will take p(y) to be dimensionless and f (t) to have dimension of length. We introduce the reduced variable and so get zy zx and Equation 7 may then be transformed into This means that for some dimensionless separation constant c d 2 and The property expressed by Eq. 14 is that the right-hand side is independent of λ. For this reason replacing f → c 1/(2−c) f will leave the right-hand side of Eq. 16 invariant. So, the c factor will cancel out in the final expression for C(x, t), and we might as well set c 1. Then Eq. 22 for f (t) is easily integrated assuming f (0) 0 since we are assuming Eq. 12, i.e., point-like initial conditions. This yields Our solution thus has the scaling form Frontiers in Physics | www.frontiersin.org December 2020 | Volume 8 | Article 519624 for some function g and with τ given by Eq. 6. Note that this form immediately gives This is in contrast to the Boltzmann transformation, which assumes step-like initial conditions, thus leading to 〈x 2 (t)〉 ∝ t.
We now turn to finding an expression for the probability density p(y). Equation 21 yields which may be integrated to where K is an integration constant.
On physical grounds we require that Fick's law be valid throughout the domain, and so C must be continuously differentiable. By symmetry it must also be symmetric around x 0 and so p ′ (0) 0 and p(0) finite. This implies that the integration constant K 0 so that we get the equation which is integrable. Rewriting it as dp y it may be integrated to yield where k is yet an integration constant. For c < 0, the factor multiplying y 2 is negative, c/2(1 − c) < 0, so that the solution is restricted to |c| < 2k(c − 1)/c , while for 0 < c < 1, the solution is valid for all y values. The constant k is determined by the normalization which follows from Eq. 13. Since the two cases of positive and negative γ-values give either finite or infinite support for p , they give rise to different functions k(c) which we shall denote k ± for c > 0 (+) and c < 0 (−). They are where the integrals These integrals may be calculated numerically, or evaluated through the analytic expressions and where Γ(z) is the Euler gamma function.
We may now reconstruct the normalized concentration field C(x, t) using eq. 16. For c < 0 we get where and for c > 0 We now calculate x RMS as, where for c < 0, For 0 < c < 2 3 , Frontiers in Physics | www.frontiersin.org December 2020 | Volume 8 | Article 519624 For c > 2/3, the integral in Eq. 25 does not converge. This means that as γ approaches 1 there is no continous transition to the ballistic regime as may be indicated by the τ 1 value at c 1. For c < 2/3 the diffusion exponent is given by Eq. 6 which is the result found by Pattle for c < 0. The fact that there is a gap of γ-values where the integral for 〈x 2 (t)〉 is undefined does not mean that the solution of Eq. 30 is invalid; it only means that the C(x, t) distribution becomes too broad. In fact, the one-sided 〈x(t)〉 for x > 0 remains well-defined with 〈x(t)〉 ∼ t τ .

DERIVATION OF THE DIFFUSION EQUATION FROM PARTICLE DYNAMICS
We now turn to stochastic modeling of the process we so far have described using the non-linear diffusion Eq. 7. Following the discussion in the van Kampen book Stochastic Processes in Physics and Chemistry [32], we derive the non-linear diffusion equation, which is an example of a Fokker-Planck equation.

Particle Model and the Fokker-Planck Equation
A population of N p particles are propagated by a sequence of random steps of zero mean using a concentration dependent step length. For every time the particle positions, which take on continuous values, are updated, the concentration field C(x, t) is updated onto a discrete onedimensional lattice of unit lattice constant. The value of C(x, t) is simply set to the number of particles between x int − 1/2 and x int + 1/2 where x int is the closest integer to x. The particle positions x i are updated according to the following algorithm where is a Wiener process and η is a random variable with 〈η〉 0 and 〈η 2 〉 1, where g(C) may be the normalized concentration. Now, following the discussion in the van Kampen book, Chapter VIII.2 we derive the corresponding Fokker-Planck equation as follows. The Chapman-Kolmogorov or master equation describing the above stochastic process is where W(x, r) is the number of particles per time and length that jump a distance r starting from x. We note that the formalism remains valid also when W(x, t) depends on x via C(x, t) itself.
Using the standard assumption that W(x, r) varies slowly with x, but rapidly with r, we may Taylor expand around x. This yields the Fokker-Planck equation where a 2 (x) is the mean squared jump length per time, according to Eq. 43. Setting g(C) bC −c/2 gives zC zt b 2 2 and requiring equivalence with Eq. 7 thus implies that b 2 2/(1 − c).

Itô-Stratonovitch Dilemma
However, the presence of a C-dependence in the diffusivity D introduces an ambiguity in the implementation of Eq. 42, since now Δx also depends on C, which in turn depends on all the Δx's. So, the question is whether one should use C(x) or C(x + Δx) or perhaps something in between? Since Δx ∼ Δt √ these choices are not equivalent.
Stratonovitch read Eq. 42 as [32] x while Itô read it as It turns out that it is the choice opted for by Itô that gives Eq. 45, while the Stratonovitch choice gives zC zt see the van Kampen book [32], Chapter VI.4 for a derivation of this result. By setting g(C) bC −c/2 again we can write the above equation as and equivalence with Eq. 7 now implies that b 2 2/(1 − c/2). It is also possible to read read Eq. 42 as This is known as the Hänggi-Klimontovich interpretation [22-24, 33, 34] Setting g(C) bC −c/2 again leads to the equation and equivalence with Eq. 7 now implies that b 2 2.
Frontiers in Physics | www.frontiersin.org December 2020 | Volume 8 | Article 519624 5 This means that the only difference between the different implementations of Eq. 42 is the magnitude of the random step. In the Itô case the step length has to be a bit smaller than in the Stratonovitch case, which in turn is smaller than in the Hänggi-Klimontovich case, in order to correspond to the same macroscopic descriptions for c ≠ 0. When c 0 the C-dependence of D goes away, and the three interpretations give the same b, as one would expect. In the simulations it is convenient to use the Itô implementation and thus The other implementations are complicated by the fact that the step length of the particles depends on the local concentration after the steps are taken. Numerically, this problem is similar to that of implicit solvers for partial differential equations.

SIMULATIONS
The numerical simulations are carried out using a population of N p random walkers that take on continous x-values but are sampled to yield a concentration profile C(x, t) on a discrete set of x-values, as discussed above. First, the particles are initialized on the negative x-axis to produce a step profile, then they are initialized at the same location, so that the initial concentration is a δ-function. The time step is dt 5 10 − 5 , unless otherwise stated, and N p 2000 particles are used.

The Boltzmann Case
We start by considering the step initial conditions in the simulations by setting C(x, 0) Θ(−x). In the simulations, which must be carried out over a finite domain, we truncate the initial condtion at some x min < 0. This value, which introduces a second step in C(x, 0), is chosen to be sufficiently large, so as to keep it from affecting the behavior for positive x. In all the simulations D 0 1, C 0 1. Figure 1A shows the concentration profile C(x, t) for different times plotted against the reduced variable y x/ t √ using c ± 0.5. There is data collapse in accordance with Eq. 10, allthough, interestingly, the master curves do in fact depend on γ.
In Figure 1B we show the mean square displacement Np i 1 x 2 i /N p where the sum runs over all particles with positions x i > 0. This quantity is easily calculated as the motion of each particle is traced. The black line shows a slope of 1, thus showing that the step profile leads to normal diffusion.

The Initial Delta-Function
In this case all the particles start at the origin for t 0, thus fulfilling the initial condition of Eq. 12. Figure 2 shows simulations of the concentration and how it evolves after a given time for different γ-values. It is seen that the positive γ-values yield a smaller diffusive spread but with more pronounced tails at larger x. The main effect is that the negative γ-values yield higher diffusivity in the high C regions,  while the the positive γ-values suppress the spreading from these regions while enhancing the spread at small C-values. Figure 3 show simulations compared to the analytic prediction of Eq. 30. We plot p(y) C(x, t)/f (t) against y x/f (t). The figures show a good agreement with theory as well as the expected data collapse between curves sampled at different times. Figure 4A, which shows 〈x 2 (t)〉, show that, unlike the Boltzmann case, the slopes are different for different γ values. Using a range of γ-values Figure 4B compares the measured values of τ to the prediction in Eq. 6 over the range of γ-values [−1, 2/3]. Good agreement between simulations and theory is observed.

SUMMARY AND CONCLUSION
With a power law diffusivity and a delta-function initial condition, there is no intrinsic length scale in the problem. In this case normalizability leads to the scaling form of Eq. 24 leading to the exponent relation Eq. 6 which is the defining characteristic of anomalous diffusion. With a step function initial condition however, as first studied by Boltzmann [26] the solution extends to x −∞ and cannot be normalized. Hence, in this case, Eq. 24 is replaced by Eq. 10 which gives normal τ 1/2 diffusion. If the step function were modified to a normalizable profile, it would necessarily imply the introduction of a length scale, the width of the profile. In this  case one would expect a late-time crossover to the anomalous diffusion of the δ-profile. We have reviewed the Boltzmann result and demonstrated, as did Pattle 60 years ago [25], that when D(C) ∼ C −c where c < 0, the non-linear diffusion equation is analytically solvable and leads to anomalous diffusion. We have proceeded, however, to consider the case when 0 < c < 1 which is indeed analytically solvable and yields super-diffusion.
We have also constructed a stochastic particle dynamics that satisfies the non-linear diffusion equation at hand, and implemented it computationally. Using this approach, we were able to verify the central results that we were derived analytically, and shown agreement between simulations and theory. In future applications, the particle model may also serve as a stepping stone for generalizations of our present model, that are not analytically tractable.
Both Küntz and Lavallée [15] and Gosh et al. [16] report subdiffusion when the diffusivity decreases with increasing concentration. Figure 2, where we show the concentration profile C(x, t) at a given time for different γ seems consistent with this: the smaller (i.e., less positive) the value of γ, the broader the concentration profile. However, Eq. 6, which is verified numerically in Figure 4, concludes the opposite: For positive γ, the exponent τ > 1/2, indicating super-diffusion. The value of τ controls the shape of the concentration profile, not how fast it spreads. There is thus a lack of consistency between the results presented here and those of refs. [15] and [16]. A possible explanation is that the more general approach of Refs. [21][22][23][24] employing a position-dependent diffusivity that does not depend on the position through the concentration.
As mentioned in the introduction, anomalous diffusion originating from a concentration dependent diffusivity may have been seen in diffusion in granular media [11,12]. These observations are based on rotating a bi-disperse composition of smaller and large glass beads in a horizontal cylindrical mixer. The mixer is filled with the larger beads except for a small disk of smaller beads. As the cylinder turns, the smaller beads diffuse into the larger beads and the concentration of smaller beads as a function of time and position along the cylinder is recorded. This setup mimics closely the initial conditions that we have studied here, except for Section 2.1, where we assumed a step initially. The connection with the present work is the proposal that the diffusivity of the smaller beads is larger when they are surrounded by other smaller beads than when they are surrounded by the larger beads; the higher the concentration of smaller beads, the larger their diffusivity is. We propose here to prepare the packing in a different way initially. Fill (say) the left half of the cylinder with the smaller beads and the right half with the larger beads. The system is therefore initiated with a step function in the concentration. One would then expect normal diffusion where the front evolves as x 2 ∼ t, i.e., the parabolic law as predicted by Boltzmann.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.