Numerical solution for a non-Fickian diffusion in a periodic potential

Numerical solutions of a non-Fickian diffusion equation belonging to a hyperbolic type are presented in one space dimension. The Brownian particle modelled by this diffusion equation is subjected to a symmetric periodic potential whose spatial shape can be varied by a single parameter. We consider a numerical method which consists of applying Laplace transform in time; we then obtain an elliptic diffusion equation which is discretized using a finite difference method. We analyze some aspects of the convergence of the method. Numerical results for particle density, flux and mean-square-displacement (covering both inertial and diffusive regimes) are presented.

💡 Research Summary

The paper addresses the numerical solution of a non‑Fickian (hyperbolic) diffusion equation, often called the telegraph equation, for a Brownian particle moving in a one‑dimensional symmetric periodic potential. The authors start by presenting the governing equations: the classical Fickian diffusion equation with a potential term and its non‑Fickian counterpart that includes an inertial second‑time‑derivative term. Both can be derived from the Kramers phase‑space equation. The non‑Fickian model reduces to the ordinary diffusion equation for times much larger than the friction time 1/γ, but retains short‑time inertial effects that allow a meaningful definition of particle velocity.

A flexible periodic potential V(x;α)=J₀(iα)⁻¹e^{α cos x}−1 is introduced, where the single parameter α controls the shape from shallow sinusoidal (α≈1) to sharply peaked barriers (α≈16). After nondimensionalisation, the system consists of a second‑order time equation for the particle density n(x,t) and a first‑order equation for the flux j(x,t). Initial conditions are Gaussian for n and the corresponding flux, while the domain is taken effectively infinite with vanishing boundary values.

The numerical strategy is three‑fold. First, a Laplace transform in time eliminates the temporal derivatives, yielding for each complex Laplace variable s an ordinary differential equation in space: d²ĥn/dx²−λ_s ĥn−d/dx(P ĥn)=−(1+s)n(x,0), where λ_s=s²+s and P(x)=−(1/mγD) dV/dx. Second, the spatial domain is discretised on a uniform grid using central finite differences, giving a tridiagonal (three‑band) linear system K(s)·η(s)=b(s) for the transformed density values η_i(s). The matrix K(s) depends on s but remains sparse, allowing efficient solution for many s values, which is advantageous for parallel computation. Third, the inverse Laplace transform is performed numerically. The Bromwich integral is rewritten as a Fourier‑type integral and evaluated with the trapezoidal rule. Because the resulting series converges slowly, the authors accelerate it by representing the series as a continued fraction and truncating when successive approximations differ by less than a prescribed tolerance (TOL). The parameter β, which shifts the integration contour to the right of all singularities, controls the exponential decay of the trapezoidal error E_T; choosing β sufficiently larger than the growth rate σ of the solution makes E_T arbitrarily small.

Error analysis is carried out in detail. Spatial discretisation error E_S is shown to be O(Δx²) by a Taylor expansion of the exact solution in the finite‑difference operator. The Laplace inversion error consists of the trapezoidal error E_T (controlled by β) and the continued‑fraction truncation error E_MF (controlled by TOL). To validate the error estimates, the authors consider a test problem with constant drift P and pure Fickian diffusion, for which an analytical solution is known. Numerical experiments demonstrate that with β chosen as −ln(10⁻⁶)/(2T) and TOL≈1/N², the global error can be reduced below 10⁻⁶, and the number of continued‑fraction iterations M_i concentrates in regions with steep gradients, confirming the adaptivity of the algorithm.

Physical results are then presented. Simulations for α=1 and α=16 illustrate how the periodic potential influences the evolution of the density, flux, and mean‑square displacement ⟨x²(t)⟩. At very short times the inertial term dominates, leading to ballistic growth ⟨x²⟩∝t². As time progresses, the system transitions to a diffusive regime with ⟨x²⟩∝t, as predicted by the long‑time limit of the telegraph equation. Larger α (deeper barriers) suppress particle spreading, causing the density to accumulate near potential minima and the flux to develop non‑Fickian contributions (the term γP(x)n) that are absent in ordinary diffusion. These findings are consistent with recent experimental observations of instantaneous Brownian velocities and with models of ion transport in insulating media that also employ hyperbolic diffusion equations.

In the concluding section the authors emphasise that the Laplace‑transform‑based method is particularly well suited for problems requiring long‑time integration, because it avoids the step‑size restrictions of explicit time‑marching schemes and provides a clear framework for error control. They suggest that the approach can be extended to higher spatial dimensions, to non‑linear potentials, and to more complex boundary conditions. Future work may involve higher‑order spatial discretisations, adaptive mesh refinement, and coupling with stochastic particle methods to further enhance accuracy and efficiency. Overall, the paper delivers a robust, theoretically grounded numerical technique for hyperbolic diffusion problems and demonstrates its applicability to physically relevant scenarios involving periodic potentials.