Time-step restrictions for numerical approximations of the Poisson-Nernst-Planck (PNP) equations

The Poisson–Nernst–Planck (PNP) system is an accurate model of electrodiffusion of ionic species. It is commonly used in situations where nanoscale resolution is required, for instance close to ion channels in the membranes of biological cells. The inherent stiffness of the equations has made them challenging to solve and has limited the applicability of the system. In particular, the time step required for stable solutions has typically needed to be very short (nanoseconds), which makes simulations on the time scale of an action potential (milliseconds) difficult. Recently, it has been observed that avoiding operator splitting and instead solving the concentration equations and the electrostatic equation in a coupled manner relaxes the time-step limitation considerably. However, no theoretical explanation of this observation has been provided. Here, we aim to explain why the coupled scheme allows much larger time steps. We illustrate the mechanism by considering special cases that define necessary, but not sufficient, conditions for stability. We also show that these conditions remain relevant for the fully coupled PNP model in 3D.