An asymptotic preserving and accurate scheme for multiscale Poisson-Nernst-Planck (MPNP) system

In this paper, we propose and validate a two-species multiscale model for a Poisson-Nernst-Planck (PNP) system, focusing on the correlated motion of positive and negative ions under the influence of a trap. Specifically, we {\color{red}aim} to model surface traps whose attraction range, of length $\delta$, is much smaller than the \textcolor{blue}{length} scale of the problem. The physical setup considered here is an anchored gas drop (bubble) surrounded by a flow of charged surfactants {(composed by positive and negative ions) that diffuses in water. When the diffusing surfactants reach the surface of the trap, the negative ions are adsorbed because of their hydrophobic tail that is attracted by the air bubble}. As in our previous works \cite{astuto2023multiscale,ASTUTO2023111880,astuto2025time,astuto2024high}, the effect of the attractive potential is replaced by a suitable boundary condition derived by mass conservation and asymptotic analysis. The main novelty of this work is the extension of the model proposed in \cite{astuto2023multiscale}, now incorporating the simultaneous influence of both carriers -- positive and negative ions -- which is often neglected in traditional approaches that treat ion species independently. {The two carriers interact through the Coulomb potential, which is computed from a Poisson equation. This leads to a multiscale model with two additional equations compared to the initial problem formulated in \cite{astuto2023multiscale}.} In the second part of the paper, we address the treatment of the Coulomb interaction. When the Debye length $\lambda_D$ ( {related} to a small parameter $\varepsilon$) is very small, one can adopt the so-called quasi-neutral limit, which significantly simplifies the system, reducing it to a {single} diffusion equation for {the sum of the two carriers} with effective diffusion coefficient \cite{jungel,CiCP-31-707}. {While this approach significantly simplifies the mathematical model, by reducing the system from three equations to a single one and eliminating the stiffness induced by a vanishingly small Debye length, it fails to capture the effects arising for non-negligible values of $\varepsilon$. In the regime where the Debye length is small but not negligible, capturing small deviations from the quasi-neutral limit may become computationally expensive when using standard methods in the literature.} One of the objectives of this work is to develop a second-order \textit{Asymptotic Preserving} (AP) numerical scheme that works for all Debye lengths and reduces to a consistent discretization of the quasi-neutral limit as $\varepsilon \to 0$, with no stability restriction on the time step. Furthermore, the numerical scheme we propose is also \textit{Asymptotic Accurate} (AA), which means that it preserves second-order accuracy in the quasi-neutral limit.