# Local Relaxation Fast Poisson Methods on Hierarchical Meshes
**Authors**: Zhenli Xu, Qian Yin, Hongyu Zhou
> School of Mathematical Sciences, MOE-LSC, and CMA-Shanghai, Shanghai Jiao Tong University, Shanghai 200240, China.
> Department of Applied Mathematics, the Hong Kong Polytechnic University, Hong Kong, China.
> Corresponding author. School of Mathematical Sciences, MOE-LSC, and CMA-Shanghai, Shanghai Jiao Tong University, Shanghai 200240, China.
## Abstract
The local relaxation algorithm is promising for fast solution of Poisson’s equations, which computes the electric field distribution in a stepwise manner via local curl-free updates while strictly enforcing Gauss’s law. We propose a novel hierarchical local relaxation (HLR) method for speeding up the convergence of curl-free iterations. The local algorithm reformulates the Poisson’s equation into the electric-field form and sweeps each cell to minimize the associate electric energy, avoiding the solution of linear systems. The updates with hierarchical meshes significantly accelerate the slow convergence of low-frequency components of the residual in the local curl-free update process. Convergence analysis is performed to obtain the convergence of the hierarchical relaxation approaches. Numerical results show that the HLR methods have the nice properties in accuracy and efficiency and the hierarchical construction leads to an overall computational complexity of $\mathcal{O}(N\log N)$ with respect to the number of grid points. Particularly, the applications in solving the Poisson–Boltzmann and Poisson–Nernst–Planck equations demonstrate the attractive performance for problems which frequently solve the Poisson’s equations.
Key words and phrases: Poisson’s equation; Inhomogeneous dielectric permittivity; Hierarchical structure; Local curl-free relaxation.
## 1 Introduction
The Poisson’s equation is a fundamental model describing the relation between source distributions and potential fields and arises ubiquitously in many fields [3, 15, 27, 43]. In molecular dynamics (MD) simulations, accurate evaluation of long-range electrostatic forces relies on the repeated solution of the Poisson’s equation or its equivalent formulations [16, 38, 44]. Besides, it is coupled with the Nernst–Planck or Vlasov equations to describe ionic transport and the self-consistent evolution of charged particles under electric fields [19, 46]. In all cases, efficient and scalable solvers are essential for efficient simulations. Repeatedly solving the Poisson’s equation is a major computational for simulating large or spatially heterogeneous systems. The finite element methods (FEMs) [7, 40, 48] and boundary element methods (BEMs) [2, 6, 13] are two classes of traditional methods for the Poisson’s equation, but they are considered less efficient for applications of time-dependent problems. The FEM requires repeated matrix assembly when the dielectric permittivity evolves in time. The BEM relies on analytical fundamental solutions, which are generally unavailable for problems with spatially varying dielectric permittivity. Local algorithms provide an efficient alternative by treating the electric field as the fundamental variable and enforcing Gauss’s law as a constraint, rather than directly solving for the electrostatic potential [26, 37]. Their local formulation accommodates time-dependent dielectric permittivity without extra cost compared to the constant permittivity case.
The local algorithm was first introduced for systems of discrete charges, for which Maggs and his collaborators [26, 25, 37] proposed to use curl-free relaxation for performing local Monte Carlo simulations of Coulomb systems. This idea was extended to MD methods [10, 33] with converting globally coupled Coulomb interactions into locally updated electromagnetic fields, enforcing Gauss’s law and quasi-static constraints while preserving thermodynamic properties. Such local approaches enable linear scaling and are naturally parallelizable, while maintaining physically consistent particle–field dynamics. Subsequently, the same principle was applied to the solution of the Poisson–Boltzmann (PB) equation, by reformulating the PB model as a variational problem that minimizes a free energy functional and solving it using a Maggs-type algorithm [1, 47]. Recently, Qiao et al. [34, 35] further generalized local algorithms to the Poisson–Nernst–Planck (PNP) model, developing the Maxwell–Ampère Nernst–Planck (MANP) framework to describe charge transport. Besides, local algorithms provide an efficient and physically consistent strategy for particle-in-cell (PIC) simulations [36], allowing electric fields and particle positions to be updated locally while strictly preserving Gauss’s law and supporting large-scale, parallel computations.
In existing literature, electric field is updated locally within each discretized cell. However, such relaxation update is less efficient in eliminating low-frequency components of the curl residual. Motivated by multigrid methods for addressing the slow convergence of Gauss–Seidel iterations for low-frequency errors in linear systems [14, 28, 30], we propose an efficient hierarchical local relaxation (HLR) method for solving Poisson’s equation. The HLR method decomposes the computational domain into a hierarchy of grids with different resolutions and performs curl-free relaxation at each level. Relaxation on fine grids efficiently damps high-frequency errors, while coarse-grid corrections effectively reduce low-frequency errors, leading to accelerated overall convergence. We show that this idea can be naturally extended to other algorithms based on the Maggs-type local formulation, including solvers for the PB equation and the MANP system. Moreover, by employing the method of image charges [29], we generalize the local algorithms to handle non-periodic boundary conditions.
Compared with traditional fast Poisson solvers such as FFT-based [8, 11, 32] and AMG-based [39, 41] methods, the HLR method benefits from its local formulation, which naturally accommodates spatially and temporally varying dielectric permittivity. Unlike AMG-based approaches that require discretization into large sparse linear systems and repeated stiffness matrix assembly for time-dependent media, the HLR method avoids global algebraic formulation and matrix assembly altogether. Furthermore, by incorporating a hierarchical structure into the local curl-free relaxation, the HLR method efficiently removes low-frequency curl residuals, significantly accelerating convergence and improving overall computational efficiency. Due to these advantages, the HLR is promising for broad applications which require a fast Poisson solver.
The rest of this paper is organized as follows. In Section 2, we reformulate Poisson’s equation and show that solving the Poisson’s equation is equivalent to computing the electric field. Section 3 introduces a class of local algorithms for the Poisson’s equation and presents the HLR method. In Section 4, we provide a rigorous theoretical analysis of the proposed HLR method. Section 5 presents numerical experiments that demonstrate the accuracy and efficiency of the proposed algorithm. Finally, conclusions are drawn in Section 6.
## 2 Electric-field formulation of the Poisson’s equation
Consider the Poisson’s equation of the form
$$
-\nabla\cdot(\varepsilon(\bm{x})\nabla\phi)=\rho(\bm{x}),
$$
where $\phi$ is a scalar field, $\varepsilon$ is a spatially varying coefficient, and $\rho$ is a given source term. In charged systems, the Poisson’s equation describes how the electrostatic potential is generated by a given charge distribution and mediated by the dielectric response of the medium with $\phi$ corresponding to the electric potential, $\varepsilon$ being the dielectric permittivity, and $\rho$ describing the charge density. Accordingly, the electric field is given by $\bm{E}=-\nabla\phi$ . With $\bm{E}$ , the Poisson’s equation can be written as,
$$
\displaystyle\nabla\cdot\left(\varepsilon\bm{E}\right)=\rho, \displaystyle\nabla\times\bm{E}=0.
$$
For a bounded simply connected domain, the equivalence between (2.1) and (2.2) can be constructed by Helmholtz decomposition [20, 31].
The electrostatic energy of a charged system can be expressed as a functional of the electric field,
$$
\mathcal{F}_{\text{pot}}[\bm{E}]:=\int_{\Omega}\frac{\varepsilon|\bm{E}|^{2}}{2}d\bm{x}.
$$
In electrostatics, the dielectric permittivity $\varepsilon$ is required to be positive to ensure a physically admissible and thermodynamically stable continuous medium, for which the electrostatic energy density remains positive definite. Proposition 1 states the convex property of electrostatic energy, which enables the electrostatic field to be computed by solving a constrained energy functional minimization problem.
**Proposition 1**
*The electrostatic energy functional $\mathcal{F}_{\text{pot}}:L^{2}(\Omega,\mathbb{R}^{n})\rightarrow\mathbb{R}$ is convex. Specifically, For any $\bm{E}_{1},\bm{E}_{2}\in L^{2}(\Omega,\mathbb{R}^{n})$ , and $\forall\theta\in(0,1)$ we have
$$
\mathcal{F}_{\text{pot}}(\theta\bm{E}_{1}+(1-\theta)\bm{E}_{2})\leq\theta\mathcal{F}_{\text{pot}}(\bm{E}_{1})+(1-\theta)\mathcal{F}_{\text{pot}}(\bm{E}_{2}).
$$
The equality holds if and only if $\bm{E}_{1}=\bm{E}_{2}$ .*
* Proof*
Since $\varepsilon>0$ , (2.4) is equivalent to the following inequality,
$$
|\theta\bm{E}_{1}+(1-\theta)\bm{E}_{2}|^{2}\leq\theta|\bm{E}_{1}|^{2}+(1-\theta)|\bm{E}_{2}|^{2}.
$$
It remains to show that (2.5) holds,
| | $\displaystyle|\theta\bm{E}_{1}+(1-\theta)\bm{E}_{2}|^{2}$ | $\displaystyle\leq(|\theta\bm{E}_{1}|+|(1-\theta)\bm{E}_{2}|)^{2}$ | |
| --- | --- | --- | --- |
Equality holds if and only if $\bm{E}_{1}=\bm{E}_{2}$ , which is evident. ∎
Consider the minimization problem of the above convex electrostatic energy with the Gauss’s law constraint (2.2a). Since the electrostatic energy functional (2.3) is convex and the Gauss’s law constraint is linear, the associated constrained minimization problem is a convex optimization problem. Consequently, the Karush–Kuhn–Tucker (KKT) conditions [17] are necessary and sufficient condition for minimizer. The Lagrangian functional of the optimization problem reads
$$
\mathcal{L}[\bm{E},\phi]:=\int_{\Omega}\frac{\varepsilon|\bm{E}|^{2}}{2}d\bm{x}-\int_{\Omega}\phi(\nabla\cdot(\varepsilon\bm{E})-\rho)\,d\bm{x},
$$
where $\phi$ is the Lagrange multiplier. By the the stationarity condition and the primal feasibility condition, one has
$$
\bm{E}=-\nabla\phi\quad\text{and}\quad\nabla\cdot(\varepsilon\bm{E})=\rho.
$$
This indicates that the Lagrange multiplier $\phi$ corresponds precisely to the electrostatic potential, and the electric field at the minimizer satisfies the curl-free condition $\nabla\times\bm{E}=0$ .
In fact, it can be shown that the electric field $\bm{E}$ resulting from the minimization of the above energy functional coincides with the negative gradient of the weak solution of the Poisson’s equation in the appropriate Sobolev space. For simplicity, we assume that $\Omega=(0,L_{x})\times(0,L_{y})$ is a domain with periodic boundary conditions. Let $L_{per}^{p}(\Omega)$ and $H_{per}^{k}(\Omega)$ denote the space of all $\overline{\Omega}$ -periodic functions on $\mathbb{R}^{n}$ whose restrictions onto $\Omega$ are in $L^{p}(\Omega)$ and $H^{k}(\Omega)$ respectively. Define
| | $\displaystyle\mathring{L}_{per}^{p}(\Omega)=\{f(\bm{x})\in L_{per}^{p}(\Omega):\mathscr{A}_{\Omega}(f)=0\},$ | |
| --- | --- | --- |
where $\mathscr{A}_{\Omega}(f)$ denote the average
$$
\mathscr{A}_{\Omega}(f):=\frac{1}{|\Omega|}\int_{\Omega}f(\bm{x})\ d\bm{x}.
$$
The space of all $\overline{\Omega}$ -periodic $C^{k}$ -functions are denote by $\ C_{per}^{k}(\overline{\Omega})$ for any $k\in\mathbb{Z}$ , further define
| | $\displaystyle H(\operatorname{div},\Omega)=\{\bm{E}\in L^{2}(\Omega):\nabla\cdot\bm{E}\in L^{2}(\Omega)\},$ | |
| --- | --- | --- |
Given $\rho\in L_{per}^{2}(\Omega)$ and $\varepsilon\in L^{\infty}_{per}(\Omega)$ , denote
$$
S_{\rho}=\{\bm{E}\in H_{per}(\operatorname{div},\Omega):\nabla\cdot\varepsilon\bm{E}=\rho\}.
$$
Assume that $\rho\in\mathring{L}^{2}_{per}(\Omega)$ and $\varepsilon$ satisfies
$$
0<\varepsilon_{\text{min}}<\varepsilon(\bm{x})<\varepsilon_{\text{max}},\quad\forall\bm{x}\in\mathbb{R}^{n},
$$
where $\varepsilon_{\text{min}}\ \text{and}\ \varepsilon_{\text{max}}$ are the lower and upper bounds. Theorem 2.1 was provided in [22].
**Theorem 2.1**
*([22]) Let $\rho\in\mathring{L}^{2}_{per}(\Omega)$ and $\varepsilon\in L^{\infty}_{per}(\Omega)$ satisfy (2.7). Define
$$
\displaystyle I[\phi]=\int_{\Omega}\bigg(\frac{\varepsilon}{2}|\nabla\phi|^{2}-\rho\phi\bigg)\ d\bm{x}\quad\forall\phi\in H_{per}^{1}(\Omega), \displaystyle F[\bm{E}]=\int_{\Omega}\frac{\varepsilon}{2}|\bm{E}|^{2}\ d\bm{x}\quad\forall\bm{E}\in S_{\rho}.
$$
Then one has,
1. There exists a unique $\phi_{min}\in\mathring{H}^{1}_{per}(\Omega)$ such that $I[\phi_{min}]=\operatorname{min}_{\phi\in\mathring{H}^{1}_{per}(\Omega)}I[\phi]$ , and it is the unique weak solution in $\mathring{H}^{1}_{per}(\Omega)$ to Poisson’s equation (2.1).
1. There exists a unique $\bm{E}_{min}\in S_{\rho}$ such that $F[\bm{E}_{min}]=\operatorname{min}_{\bm{E}\in S_{\rho}}F[\bm{E}]$ . Moreover, one has $\bm{E}_{min}=-\nabla\phi_{min}.$*
This theorem demonstrates that solving the Poisson’s equation (2.1) is equivalent to solving the minimization problem in (2.9). And the condition enforcing restriction $\bm{E}\in S_{\rho}$ is precisely the Gauss’s law constraint (2.2a).
## 3 Local relaxation method
For simplicity of presentation, all theoretical derivations, algorithmic steps, and notation definitions in this section are explicitly limited to two-dimensional configurations. It is crucial to emphasize that the proposed numerical algorithms are applicable to three-dimensional problems without introducing any new challenges and extra structural adjustments, and corresponding 3D numerical experiments are included in Section 5.
### 3.1 Discretization
The computational domain $\Omega=(0,L_{x})\times(0,L_{y})$ is covered by a uniform grid with grid spacings $\Delta x$ and $\Delta y$ in two coordinate directions. Let $\Delta\Omega=\Delta x\Delta y$ denote the area of a cell. Let $N_{x}$ and $N_{y}$ be the number of grid points in each dimension. Denote $\Delta x\mathbb{Z}\times\Delta y\mathbb{Z}=\{(i\Delta x,j\Delta y):i,j\in\mathbb{Z}\}$ and $\Delta x(\mathbb{Z}+1/2)\times\Delta y(\mathbb{Z}+1/2)=\{((i+\frac{1}{2})\Delta x,(j+\frac{1}{2})\Delta y):i,j\in\mathbb{Z}\}$ . For any grid function $f:\Delta x\mathbb{Z}\times\Delta y\mathbb{Z}\rightarrow\mathbb{R}$ , one denotes $f_{ij}=f(i\Delta x,j\Delta y)$ . For $\overline{\Omega}$ -periodic grid functions, one defines spaces
$$
\displaystyle V_{h}=\{\overline{\Omega}\text{-periodic grid functions }f:\Delta x\mathbb{Z}\times\Delta y\mathbb{Z}\rightarrow\mathbb{R}\}, \displaystyle\mathring{V}_{h}=\{f\in V_{h}:\mathscr{A}_{h}(f)=0\},
$$
where $\mathscr{A}_{h}$ denotes the discrete average
$$
\mathscr{A}_{h}(f)=\frac{1}{N_{x}N_{y}}\sum_{i=0}^{N_{x}-1}\sum_{j=0}^{N_{y}-1}f_{ij}.
$$
Given $\bm{E}=(E_{x},E_{y})$ , one defines a discrete electric field as a vector-valued function
$$
\mathbb{E}_{i+\frac{1}{2},j+\frac{1}{2}}=(E_{i+\frac{1}{2},j},E_{i,j+\frac{1}{2}})\quad\forall i,j\in\mathbb{Z},
$$
where $E_{i+\frac{1}{2},j}$ represents the value of $E_{x}$ at the grid point $((i+\frac{1}{2})\Delta x,j\Delta y)$ , and $E_{i,j+\frac{1}{2}}$ represents the value of $E_{y}$ at the grid point $(i\Delta x,(j+\frac{1}{2})\Delta y)$ . Denote
$$
Y_{h}=\{\overline{\Omega}-\text{periodic functions }\mathbb{E}:\Delta x(\mathbb{Z}+1/2)\times\Delta y(\mathbb{Z}+1/2)\rightarrow\mathbb{R}^{2}\}.
$$
For any $\mathbb{E}\in Y_{h}$ , define the discrete average
$$
\mathscr{A}_{h}(\mathbb{E})=\frac{1}{N_{x}N_{y}}\sum_{i=0}^{N_{x}-1}\sum_{j=0}^{N_{y}-1}(E_{i+\frac{1}{2},j},E_{i,j+\frac{1}{2}}),
$$
and the norm
$$
||\mathbb{E}||_{h}=\sqrt{\Delta\Omega\sum_{i=0}^{N_{x}-1}\sum_{j=0}^{N_{y}-1}\left(E_{i,j+\frac{1}{2}}^{2}+E_{i+\frac{1}{2},j}^{2}\right)}.
$$
The discrete divergence and curl on $Y_{h}$ are defined by
| | $\displaystyle(\nabla_{h}\cdot\mathbb{E})_{i,j}$ | $\displaystyle=\frac{E_{i+\frac{1}{2},j}-E_{i-\frac{1}{2},j}}{\Delta x}+\frac{E_{i,j+\frac{1}{2}}-E_{i,j-\frac{1}{2}}}{\Delta y},$ | |
| --- | --- | --- | --- |
Given $\varepsilon\in C_{per}(\overline{\Omega})$ satisfying (2.7), define a new average on $\Delta x(\mathbb{Z}+1/2)\times\Delta y(\mathbb{Z}+1/2)$ by
$$
\varepsilon_{i+\frac{1}{2},j}=\frac{\varepsilon_{i,j}+\varepsilon_{i+1,j}}{2},\quad\varepsilon_{i,j+\frac{1}{2}}=\frac{\varepsilon_{i,j}+\varepsilon_{i,j+1}}{2}.
$$
Given $\rho^{h}\in\mathring{V}_{h}$ the discrete charge density, Gauss’s law (2.2a) can be discretized as
$$
\nabla_{h}\cdot\left(\varepsilon\mathbb{E}\right)=\rho^{h}.
$$
And the electrostatic energy (2.3) admits an approximation of the form
$$
\mathcal{F}^{h}[\mathbb{E}]=\frac{\Delta\Omega}{2}\sum_{i=0}^{N_{x}-1}\sum_{j=0}^{N_{y}-1}(\varepsilon_{i,j+\frac{1}{2}}E_{i,j+\frac{1}{2}}^{2}+\varepsilon_{i+\frac{1}{2},j}E_{i+\frac{1}{2},j}^{2}),
$$
which is to be minimized under the constraint of the discrete Gauss’s law (3.6).
### 3.2 Single mesh relaxation
The local algorithm begins with initializing an electric field $\mathbb{E}^{(0)}$ that satisfies the discrete form of Gauss’s law (3.6). In certain applications where the Poisson’s equation must be solved at every time step, the electric field that satisfies Gauss’s law can be obtained through time evolution from the electric field at the previous time step. In contrast, in a broader range of scenarios, this electric field is unknown a priori and thus requires explicit computation. Notably, the procedure for constructing an electric field $\mathbb{E}^{(0)}$ is not unique, one feasible approach is presented in [1]. Let $\overline{\rho}_{j}=\sum_{i=0}^{N_{x}-1}\rho^{h}_{ij}/N_{x}$ denote the average charge over the $j$ th $x$ -aligned line. For $i=0,\cdots,N_{x}-1$ , the $y$ -component of the electric field is then computed by
$$
\varepsilon_{i,j+\frac{1}{2}}E_{i,j+\frac{1}{2}}=\varepsilon_{i,j-\frac{1}{2}}E_{i,j-\frac{1}{2}}+\Delta y\overline{\rho}_{j},\quad\ j=1,\cdots,N_{y}-1,
$$
with $E_{i,\frac{1}{2}}=0$ for all $i=0,\cdots,N_{x}-1.$ Subsequently, the $x$ -directional component of the electric field is updated by manipulating the discrete Gauss’s law (3.6). For $j=0,\cdots,N_{y}-1$ , the computation follows
$$
\varepsilon_{i+\frac{1}{2},j}E_{i+\frac{1}{2},j}=\varepsilon_{i-\frac{1}{2},j}E_{i-\frac{1}{2},j}+\Delta x(\rho^{h}_{i,j}-\overline{\rho}_{j}),\quad i=1,\cdots,N_{x}-1,
$$
with $E_{\frac{1}{2},j}=0,~j=0,\cdots,N_{y}-1.$ Crucially, the discrete electric field $\mathbb{E}^{(0)}$ constructed through the sequential application of (3.8) and (3.9) is guaranteed to satisfy the discrete Gauss’s law (3.6), aligning with the initial requirement for the algorithm.
The next step of the local algorithm involves iteratively updating the electric field to minimize the discrete energy functional $\mathcal{F}^{h}$ (3.7) with constraint of Gauss’s law. This minimization ensures the electric field to be curl-free, thus coinciding with the exact solution to Poisson’s equation. It is emphasized that the local algorithm naturally and rigorously preserves Gauss’s law, which is achieved by designing a specialized numerical update scheme for electric field. The minimization of $\mathcal{F}^{h}$ is implemented via local updates [26, 25] to the discrete electric field within each computational cell, so we call it as single mesh relaxation method. Starting from the initial field ${\mathbb{E}^{(0)}}$ , the electric field is updated by introducing a rotational flux $\eta$ along the four edges of each cell, illustrated in Fig. 1. Notably, this same flux $\eta$ is chosen uniformly across the four edges of a single cell — a specific design constraint that preserves the discrete Gauss’s law strictly. For details, the process is as follows:
$$
\displaystyle E_{i,j+\frac{1}{2}} \displaystyle\leftarrow E_{i,j+\frac{1}{2}}-\frac{\eta}{\varepsilon_{i,j+\frac{1}{2}}\Delta x}, \displaystyle E_{i+1,j+\frac{1}{2}} \displaystyle\leftarrow E_{i+1,j+\frac{1}{2}}+\frac{\eta}{\varepsilon_{i+1,j+\frac{1}{2}}\Delta x}, \displaystyle E_{i+\frac{1}{2},j} \displaystyle\leftarrow E_{i+\frac{1}{2},j}+\frac{\eta}{\varepsilon_{i+\frac{1}{2},j}\Delta y}, \displaystyle E_{i+\frac{1}{2},j+1} \displaystyle\leftarrow E_{i+\frac{1}{2},j+1}-\frac{\eta}{\varepsilon_{i+\frac{1}{2},j+1}\Delta y}.
$$ *[Error downloading image: 2603.02763v1/x1.png]*
Figure 1: Schematic of the electric-field update on a single grid cell $(i,j)$ with flux $\eta$ .
It can be simply proved that discrete Gauss’s law (3.6) defined at the four nodes of this updating cell still strictly holds. With this update, the energy change is given by
$$
\displaystyle\delta\mathcal{F}^{h}[\eta]= \displaystyle\frac{\eta^{2}}{2}\left[\frac{\Delta x}{\Delta y}\left(\frac{1}{\varepsilon_{i+\frac{1}{2},j}}+\frac{1}{\varepsilon_{i+\frac{1}{2},j+1}}\right)+\frac{\Delta y}{\Delta x}\left(\frac{1}{\varepsilon_{i,j+\frac{1}{2}}}+\frac{1}{\varepsilon_{i+1,j+\frac{1}{2}}}\right)\right] \displaystyle+\eta\bigg[\Delta x\bigg(E_{i+\frac{1}{2},j}-E_{i+\frac{1}{2},j+1}\bigg)+\Delta y\left(E_{i+1,j+\frac{1}{2}}-E_{i,j+\frac{1}{2}}\right)\bigg].
$$
This is a quadratic form of $\eta$ , thus the minimum value can be exactly obtained when
$$
\eta=-\frac{\Delta y(\Delta x)^{2}\left(E_{i+\frac{1}{2},j}-E_{i+\frac{1}{2},j+1}\right)+\Delta x(\Delta y)^{2}\left(E_{i+1,j+\frac{1}{2}}-E_{i,j+\frac{1}{2}}\right)}{(\Delta x)^{2}\left(\frac{1}{\varepsilon_{i+\frac{1}{2},j}}+\frac{1}{\varepsilon_{i+\frac{1}{2},j+1}}\right)+(\Delta y)^{2}\left(\frac{1}{\varepsilon_{i,j+\frac{1}{2}}}+\frac{1}{\varepsilon_{i+1,j+\frac{1}{2}}}\right)}.
$$
This formula provides the optimal flux for the change of electric field. Such a local update traverses all discrete grid cells, and iterative updates continue until the stopping condition is met.
The single mesh update is not sufficient [21] to ensure ergodicity in the sampling for Monte Carlo simulation. It may not force the electrostatic energy (3.7) reach to its minimum [22]. To address this issue, a line shift update [18, 22] has been introduced, which uses the principle of ensuring that the flux entering a discrete node is equal to that exiting it. This design enables the strict preservation of the discrete Gauss’s law, while the optimal flux value is derived from energy minimization calculations. Specifically, for $x$ -direction, the update of electric field at $j$ th $x$ line with flux $\eta_{x}$ is
$$
E_{i+\frac{1}{2},j}\leftarrow E_{i+\frac{1}{2},j}+\frac{\eta_{x}}{\varepsilon_{i+\frac{1}{2},j}},\qquad i=0,1,...,N_{x}-1.
$$
The energy change after this update is
$$
\delta\mathcal{F}^{h}[\eta_{x}]=\Delta x\Delta y\bigg[\frac{\eta_{x}^{2}}{2}\bigg(\sum_{i=0}^{N_{x}-1}\frac{1}{\varepsilon_{i+\frac{1}{2},j}}\bigg)+\eta_{x}\bigg(\sum_{i=0}^{N_{x}-1}E_{i+\frac{1}{2},j}\bigg)\bigg].
$$
The energy reaches its minimum value, when
$$
\eta_{x}=-\sum_{i=0}^{N_{x}-1}E_{i+\frac{1}{2},j}\bigg/\sum_{i=0}^{N_{x}-1}\frac{1}{\varepsilon_{i+\frac{1}{2},j}}.
$$
Similarly, for $y$ -direction, the process of update at $i$ th $y$ line with $\eta_{y}$ is
$$
E_{i,j+\frac{1}{2}}\leftarrow E_{i,j+\frac{1}{2}}+\frac{\eta_{y}}{\varepsilon_{i,j+\frac{1}{2}}},\qquad j=0,1,...,N_{y}-1.
$$
The optimal flux $\eta_{y}$ is
$$
\eta_{y}=-\sum_{j=0}^{N_{y}-1}E_{i,j+\frac{1}{2}}\bigg/\sum_{j=0}^{N_{y}-1}\frac{1}{\varepsilon_{i,j+\frac{1}{2}}}.
$$
To summarize, the first step of the local algorithm involves the construction of an initial electric field that strictly adheres to the discrete Gauss’s law. The subsequent step includes the minimization of the electrostatic energy functional, which is tailored to achieve a curl-free electric field configuration while rigorously preserving the discrete Gauss’s law throughout the entire process. The local algorithm with single mesh relaxation is outlined in Algorithm 1.
Algorithm 1 Local Algorithm with Single Mesh Relaxation
1: Step 1. Initialize field $\mathbb{E}^{(0)}$ to satisfy the discrete Gauss’s law (3.6). Set $m=0$ .
2:
3: Step 2. Curl-free relaxation $\mathbb{E}:=\mathbb{E}^{(m)}$ .
4: Update $\mathbb{E}$ along each cell via (3.10).
5: Update $\mathbb{E}$ by the line shift via (3.13) and (3.16).
6:
7: Step 3. Termination check:
8: if $\delta\mathcal{F}^{h}<\varepsilon_{tol}$ then
9: stop
10: else
11: Set $\mathbb{E}^{(m+1)}=\mathbb{E}$ , $m\leftarrow m+1$ .
12: Go to Step 2
13: end if
**Remark 1**
*For non-periodic boundaries, we employ the image charges method (ICM) [23] to handle Neumann and Dirichlet conditions. By mapping the system to an equivalent periodic configuration, the electric field can then be computed using the local method developed here.*
### 3.3 Hierarchical mesh relaxation
The single mesh relaxation updates the electric field sequentially in each cell of the finest computational grid. A critical limitation of this approach, however, is that it yields only slow reduction of rotational errors with particularly in low-frequency error components, thereby restricting the inefficient rate at which the electrostatic energy can be minimized. Recall that the core principle of such local algorithms is to design specialized update fluxes that drive the most rapid possible reduction in electrostatic energy, while rigorously preserving the discrete Gauss’s law at all times. In fact, without considering the curl-free condition (2.2b), the general solution of Eq. (2.2a) can be represented as
$$
\varepsilon\bm{E}=-\varepsilon\nabla\phi+\nabla\times\mathbf{Q},
$$
where $\phi$ satisfies $-\nabla\cdot(\varepsilon\nabla\phi)=\rho$ , and $\mathbf{Q}$ is the rotational degree of freedom. Based on the orthogonality of the two components of $\bm{E}$ , the electrostatic energy (2.3) can be written as
$$
\mathcal{F}_{\text{pot}}[\bm{E}]:=\int_{\Omega}\frac{\varepsilon}{2}[(\nabla\phi)^{2}+(\nabla\times\mathbf{Q})^{2}]d\bm{x}.
$$
From this, local algorithms essentially keep the gradient terms, namely, Gauss’s law unchanged, and the minimization of energy is equivalent to the minimization of the energy associated with rotational term. Inspired by geometric multigrid methods [14, 28], a well-established technique that accelerates the iterative solution of linear systems by leveraging a hierarchy of computational grids with different spatial scales, we propose a HLR method by integrating a multi-level hierarchical structure into the framework. This novel design is expected to effectively accelerate the propagation of rotational degrees of freedom across the entire computational domain while keeping the principle of local algorithms.
In the HLR, the computational domain $\Omega$ is partitioned into a set of nested grid blocks across multiple refinement levels $L_{k},~k=1,\ldots,M$ , where $L_{1}$ represents the coarsest level and $L_{M}$ the finest level. In specific, the grid partitioning at level $L_{k}$ is defined as:
$$
\Omega=\bigcup_{m,n=1}^{2^{k}}\Omega_{m,n}^{(k)},
$$
where $\Omega_{m,n}^{(k)}$ denotes the $(m,n)$ -th grid block at refinement level $L_{k}$ , with each block at level $L_{k}$ being partitioned into $2\times 2$ sub-blocks to form the grid at level $L_{k+1}$ (i.e. $\Omega_{m,n}^{(k)}=\cup_{a,b=1}^{2}\Omega_{2m-a,2n-b}^{(k+1)}$ ). Notice that $2^{M}=N_{x}=N_{y}$ if the same resolutions are used in $x$ - and $y$ - coordinates.
Based on this hierarchical grid structure, we propose two HLR methods, forward HLR and zigzag HLR, within the local algorithm framework, as illustrated in Fig. 2 (a) for a single iteration. Unlike the single mesh relaxation method, where curl-free relaxation is restricted to the finest grid ( $L_{M}$ ), the HLR method shifts the curl-free relaxation from individual finest grid cells to grid blocks at different levels. Performing an update on any grid block $\Omega^{(k)}_{m,n}$ implies executing a local process on the grid nodes covered by that block at its corresponding resolution.
The forward HLR method follows a strictly sequential, layer-by-layer update path from coarse scales to fine scales grids. Starting at the coarsest grid level $L_{1}$ (e.g., the initial $2\times 2$ block grid), the algorithm completes the field update on the current coarse level, then proceeds forward and irreversibly to the next finer level $L_{2}$ . After finishing updates on $L_{2}$ , it moves directly to the finer level $L_{3}$ , and repeats this process until reaching the finest grid level $L_{\max}$ . This strategy follows a one-way progressive pathway:
$$
L_{1}\to L_{2}\to\cdots\to L_{M}.
$$
No backward revisit or correction to coarser levels is performed during the entire procedure.
In contrast, the zigzag HLR scheme adopts a reciprocal, backtracking update pathway that repeatedly travels from coarse to fine, then back to coarse, and forward to fine again, forming a zigzag trajectory across grid levels. After updating the coarsest level $L_{1}$ and refining to the finer levels $L_{2}$ and $L_{3}$ , the algorithm intentionally revisits the coarser level $L_{2}$ to apply secondary corrections based on the fine-scale information. It then refines again to $L_{3}$ and proceeds toward $L_{4}$ , with similar back-and-forth corrections at each intermediate level. Its typical pathway can be summarized as:
$$
L_{1}\to L_{2}\to L_{3}\to L_{2}\to L_{3}\to L_{4}\to\cdots\to L_{M}.
$$
This zigzag backtracking enables two-way information exchange between coarse and fine grids, rather than the one-way propagation used in the forward scheme.
*[Error downloading image: 2603.02763v1/figures/finalpm.png]*
(a) (a)
*[Error downloading image: 2603.02763v1/x2.png]*
(b) (b)
Figure 2: (a) Schematic illustration of relaxation on hierarchical meshes with forward (left), and zigzag (right) HLR. (b) Schematic of the electric-field update on a coarse-grid block covering multiple finest cells.
The procedures for local curl-free relaxation on coarse grids are analogous to those on finest grids described in the previous section. At the grid level $L_{k}$ , the curl-free relaxation is performed on grid blocks $\Omega^{(k)}_{m,n},~m,n=1,\cdots,2^{k}$ , where each block is uniquely specified by its lower-left coordinate $(i_{\min},j_{\min})=\left(2^{M-k}(m-1),2^{M-k}(n-1)\right)$ and upper-right corner coordinate $(i_{\max},j_{\max})=(2^{M-k}m,2^{M-k}n)$ in the global finest-grid index system. For any given block $\Omega^{(k)}_{m,n},~m,n=1,\cdots,2^{k}$ , the curl-free relaxation update is applied exclusively to all electric field components defined on the block boundary edges, using the same rotational flux–based strategy established earlier. As illustrated in Fig. 2 (b), the update procedure for a block with lower-left coordinate $(i_{\min},j_{\min})$ and upper-right coordinate $(i_{\max},j_{\max})$ proceeds as follows:
$$
\displaystyle E_{i_{\min},j+\frac{1}{2}} \displaystyle\leftarrow E_{i_{\min},j+\frac{1}{2}}-\frac{\eta}{\varepsilon_{i_{\min},j+\frac{1}{2}}\Delta x},\quad j_{\min}\leq j<j_{\max}, \displaystyle E_{i_{\max},j+\frac{1}{2}} \displaystyle\leftarrow E_{i_{\max},j+\frac{1}{2}}+\frac{\eta}{\varepsilon_{i_{\max},j+\frac{1}{2}}\Delta x},\quad j_{\min}\leq j<j_{\max}, \displaystyle E_{i+\frac{1}{2},j_{\min}} \displaystyle\leftarrow E_{i+\frac{1}{2},j_{\min}}+\frac{\eta}{\varepsilon_{i+\frac{1}{2},j_{\min}}\Delta y},\quad i_{\min}\leq i<i_{\max}, \displaystyle E_{i+\frac{1}{2},j_{\max}} \displaystyle\leftarrow E_{i+\frac{1}{2},j_{\max}}-\frac{\eta}{\varepsilon_{i+\frac{1}{2},j_{\max}}\Delta y},\quad i_{\min}\leq i<i_{\max}.
$$
This update preserves the discrete Gauss’s law defined on all grids of the boundary of this block. Similarly, the energy change after the electric displacement update is
$$
\displaystyle\delta\mathcal{F}^{h}[\eta]= \displaystyle\frac{\eta^{2}}{2}\left[\frac{\Delta x}{\Delta y}\sum_{i=i_{\min}}^{i_{\max}-1}\left(\frac{1}{\varepsilon_{i+\frac{1}{2},j_{\min}}}+\frac{1}{\varepsilon_{i+\frac{1}{2},j_{\max}}}\right)+\frac{\Delta y}{\Delta x}\sum_{j=j_{\min}}^{j_{\max}-1}\left(\frac{1}{\varepsilon_{i_{\min},j+\frac{1}{2}}}+\frac{1}{\varepsilon_{i_{\max},j+\frac{1}{2}}}\right)\right] \displaystyle+\eta\bigg[\Delta x\sum_{i=i_{\min}}^{i_{\max}-1}\bigg(E_{i+\frac{1}{2},j_{\min}}-E_{i+\frac{1}{2},j_{\max}}\bigg)+\Delta y\sum_{j=j_{\min}}^{j_{\max}-1}\left(E_{i_{\max},j+\frac{1}{2}}-E_{i_{\min},j+\frac{1}{2}}\right)\bigg].
$$
The maximization of the energy change gives the exact solution of $\eta$ :
$$
\eta=-\frac{\Delta y(\Delta x)^{2}\sum_{i=i_{\min}}^{i_{\max}-1}\bigg(E_{i+\frac{1}{2},j_{\min}}-E_{i+\frac{1}{2},j_{\max}}\bigg)+\Delta x(\Delta y)^{2}\sum_{j=j_{\min}}^{j_{\max}-1}\left(E_{i_{\max},j+\frac{1}{2}}-E_{i_{\min},j+\frac{1}{2}}\right)}{(\Delta x)^{2}\sum_{i=i_{\min}}^{i_{\max}-1}\left(\frac{1}{\varepsilon_{i+\frac{1}{2},j_{\min}}}+\frac{1}{\varepsilon_{i+\frac{1}{2},j_{\max}}}\right)+(\Delta y)^{2}\sum_{j=j_{\min}}^{j_{\max}-1}\left(\frac{1}{\varepsilon_{i_{\min},j+\frac{1}{2}}}+\frac{1}{\varepsilon_{i_{\max},j+\frac{1}{2}}}\right)}.
$$
The local algorithms with HLR methods are summarized in Algorithm 2.
Algorithm 2 Local Algorithm with HLR Methods.
1: Step 1. Initialize displacement $\mathbb{E}^{(0)}$ satisfy the discrete Gauss’s law (3.6). Set $m=0$ .
2:
3: Step 2. Curl-free Relaxation $\mathbb{E}:=\mathbb{E}^{(m)}$ .
4: Forward:
5: for $\text{Level}:~L_{k}=L_{1}$ to $L_{M}$ do
6: for Each block: $m,n=1$ to $2^{k}$ do
7: Update $\mathbb{E}$ to get $\mathbb{E}^{*}$ via (3.18) and $\mathbb{E}\leftarrow\mathbb{E}^{*}$ .
8: end for
9: end for
10: Zigzag:
11: for $\text{Level}:~L_{\ell}=L_{1}$ to $L_{M-2}$ do
12: for $\text{Level}:~L_{k}=L_{\ell}$ to $L_{\ell+2}$ do
13: for Each block: $m,n=1$ to $2^{k}$ do
14: Update $\mathbb{E}$ to get $\mathbb{E}^{*}$ via (3.18) and $\mathbb{E}\leftarrow\mathbb{E}^{*}.$
15: end for
16: end for
17: end for
18: Update $\mathbb{E}$ by the line shift via (3.13) and (3.16).
19:
20: Step 3. Termination check:
21: if $\delta\mathcal{F}^{h}<\varepsilon_{tol}$ then
22: stop
23: else
24: Set $\mathbb{E}^{(m+1)}=\mathbb{E}$ , $m\leftarrow m+1$ .
25: Go to Step 2
26: end if
## 4 Convergence Analysis
In this section, we will present the proof of convergence for the new proposed HLR relaxation method, following the idea of [22]. For convenience, one assumes that $L_{x}=L_{y}=L$ and $N_{x}=N_{y}=N=2^{M}$ , such that $\Delta x=\Delta y=h=L/{2^{M}}$ . Thus, one can denote $h\mathbb{Z}^{2}=\Delta x\mathbb{Z}\times\Delta y\mathbb{Z}$ and $h(\mathbb{Z}+1/2)^{2}=\Delta x(\mathbb{Z}+1/2)\times\Delta y(\mathbb{Z}+1/2)$ . Given $\varepsilon\in Y_{h}$ and $\rho^{h}\in\mathring{V}_{h}$ , define
$$
S_{\rho,h}=\{\mathbb{E}\in Y_{h}:\nabla_{h}\cdot\left(\varepsilon\mathbb{E}\right)=\rho^{h}\ \text{on}\ h\mathbb{Z}^{2}\}.
$$
Thus, the discrete electrostatic energy (3.7) can be treated as a functional $\mathcal{F}^{h}:S_{\rho,h}\rightarrow\mathbb{R}$ . Lemma 4.1 and 4.2 were obtained in [22], which provide some equivalent conditions on the minimizer of $\mathcal{F}^{h}$ .
**Lemma 4.1**
*Let $\varepsilon\in Y_{h}$ satisfy (2.7) and $\rho^{h}\in\mathring{V}_{h}$ . There exists a unique minimizer $\mathbb{E}_{min}^{h}$ of $\mathcal{F}^{h}:S_{\rho,h}\rightarrow\mathbb{R}$ . For $\mathbb{E}\in S_{\rho,h}$ , $\mathbb{E}=\mathbb{E}_{min}^{h}$ if and only if:
1. Local equilibrium: $\mathbb{E}$ is curl free, i.e, $\nabla_{h}\times\mathbb{E}=0$ on $h(\mathbb{Z}+1/2)^{2}$ ;
1. Zero total field : $\mathscr{A}_{h}(\mathbb{E})=0$ in $\mathbb{R}^{2}$ .*
**Lemma 4.2**
*Let $\varepsilon\in Y_{h}$ satisfy (2.7), $\rho^{h}\in\mathring{V}_{h}$ and $\mathbb{E}\in S_{\rho,h}$ . Given $i,j\in\{0,1,\ldots,N-1\}$ , $\eta_{ij}$ is given by (3.12). Then $\mathbb{E}$ is curl free, i.e., $\nabla_{h}\times\mathbb{E}=0$ on $h(\mathbb{Z}+1/2)^{2}$ , if and only if $\eta_{ij}=0$ for all $i$ and $j$ .*
Lemma 4.3 gives some properties of local relaxation methods, which will be used to prove the convergence of the algorithms proposed in Section 3.3.
**Lemma 4.3**
*Let $\varepsilon\in Y_{h}$ satisfy (2.7), $\rho^{h}\in\mathring{V}_{h}$ and $\mathbb{E}\in S_{\rho,h}$ .
1. Given $i,j\in\{0,1,\ldots,N-1\}$ . Let $\hat{\mathbb{E}}$ be updated from $\mathbb{E}$ by (3.10) with $\eta_{i,j}$ given by (3.12). Then $\hat{\mathbb{E}}\in S_{\rho,h}$ and $\mathscr{A}_{h}(\varepsilon\hat{\mathbb{E}})=\mathscr{A}_{h}(\varepsilon\mathbb{E})$ .
1. Given $i_{\min},j_{\min},i_{\max},j_{\max}\in\{0,1,\ldots,N-1\}$ with $i_{\min}<i_{\max}$ and $j_{\min}<j_{\max}$ . Let $\hat{\mathbb{E}}$ be updated from $\mathbb{E}$ by (3.18) with $\eta$ given by (3.20). Then $\hat{\mathbb{E}}\in S_{\rho,h}$ and $\mathscr{A}_{h}(\varepsilon\hat{\mathbb{E}})=\mathscr{A}_{h}(\varepsilon\mathbb{E})$ .
1. Let $\hat{\mathbb{E}}$ be updated from $\mathbb{E}$ by (3.13) with $\eta_{x}$ given by (3.15) for all $j\in\{0,1,\ldots,N-1\}$ and (3.16) with $\eta_{y}$ given by (3.17) for all $i\in\{0,1,\ldots,N-1\}$ . Then $\hat{\mathbb{E}}\in S_{\rho,h}$ and $\mathscr{A}_{h}(\hat{\mathbb{E}})=0$ .*
* Proof*
The first two items can be verified by direct substitution. The proof of the third item is given below. Since the updates (3.13) for different $j$ are independent, given $j\in\{0,1,\ldots,N-1\}$ fixed,
| | $\displaystyle\sum_{i=0}^{N-1}\hat{E}_{i+\frac{1}{2},j}$ | $\displaystyle=\sum_{i=0}^{N-1}E_{i+\frac{1}{2},j}+\frac{\eta_{x}}{\varepsilon_{i+\frac{1}{2},j}}$ | |
| --- | --- | --- | --- |
Therefore,
$$
\sum_{j=0}^{N-1}{\sum_{i=0}^{N-1}\hat{E}_{i+\frac{1}{2},j}}=0.
$$
Similarly, based on the update (3.16), it follows that
$$
\sum_{i=0}^{N-1}{\sum_{j=0}^{N-1}\hat{E}_{i,j+\frac{1}{2}}}=0.
$$
Thus, we get $\mathscr{A}_{h}(\hat{\mathbb{E}})=0$ directly. ∎
By Lemma 4.1, $\mathscr{A}_{h}(\mathbb{E}_{min}^{h})=0$ is necessary for $\mathbb{E}_{min}^{h}$ to be a minimizer of the functional $\mathcal{F}^{h}$ . Moreover, Lemma 4.3 shows that the line shift acts as the mechanism that enforces the total zero field constraint on the discrete electric field. Without this step, the iterative scheme cannot ensure that the limiting electric field satisfies the total zero field condition. Therefore, the line shift is essential. We now present the convergence analysis of the algorithm in Theorem 4.1.
**Theorem 4.1**
*Let $\varepsilon\in Y_{h}$ satisfy (2.7), $\rho^{h}\in\mathring{V}_{h}$ and $\mathbb{E}_{min}^{h}\in S_{\rho,h}$ be given by the unique minimizer of $F^{h}:S_{\rho,h}\rightarrow\mathbb{R}$ . Let $\mathbb{E}^{(0)}\in\ S_{\rho,h}$ and $\mathbb{E}^{(t)}$ be the sequence generated by Algorithm 2.
1. If the sequence is finite ending at $\mathbb{E}^{(m)}$ , then $\mathbb{E}^{(m)}=\mathbb{E}_{min}^{h}$ .
1. If the sequence is infinte, then $\mathbb{E}^{(t)}\rightarrow\mathbb{E}^{h}_{min}$ on $h(\mathbb{Z}+1/2)^{2}$ .*
* Proof*
1. If the sequence is finite ending at $\mathbb{E}^{(m)}$ , then $\eta_{i,j}=0$ for all $i,j\in\{0,1,\ldots,N-1\}$ . Thus , $\mathbb{E}^{(m)}$ is curl-free by Lemma 4.2 and $\mathscr{A}_{h}(\mathbb{E}^{(m)})=0$ by Lemma 4.3. Therefore, we can get $\mathbb{E}^{(m)}=\mathbb{E}_{min}^{h}$ from Lemma 4.1. 2. Given $i,j\in\{0,1,\ldots,N-1\}$ , denote $\eta^{(t)}_{i,j}$ by (3.12) for $\mathbb{E}=\mathbb{E}^{(t)}$ . We claim that
$$
\lim_{t\rightarrow\infty}\eta_{i,j}^{(t)}=0\quad\forall i,j\in\{0,1,\ldots N-1\}.
$$ Suppose that (4.2) is true. We prove that $\mathbb{E}^{(t)}\rightarrow\mathbb{E}_{min}^{h}$ . It suffices to prove that for any convergent subsequence $\mathbb{E}^{(t_{r})}$ of $\mathbb{E}^{(t)}$ with $\mathbb{E}^{(t_{r})}\rightarrow\mathbb{E}^{(\infty)}$ , then $\mathbb{E}^{(\infty)}=\mathbb{E}_{min}^{h}$ . Since function $\mathscr{A}_{h}:Y_{h}\rightarrow\mathbb{R}$ depends continuously on $\mathbb{E}$ , we have $\mathscr{A}_{h}(\mathbb{E}^{(\infty)})=0$ . And with the assumption (4.2), $\mathbb{E}^{(\infty)}\in S_{\rho,h}$ and $\eta_{i,j}^{(\infty)}=0$ . Hence $\mathbb{E}^{(\infty)}$ is curl-free by Lemma 4.2. Thus, $\mathbb{E}^{(\infty)}=\mathbb{E}_{min}^{h}$ by Lemma 4.1. We now prove the claim (4.2). Note that the iteration from $\mathbb{E}^{(t)}$ to $\mathbb{E}^{(t+1)}$ consists of many cycles of hierarchical mesh updates, single mesh updates and line shift updates. For convenience, we redefine the sequence of updates. Still start from $\mathbb{E}^{(0)}$ , $\mathbb{E}^{(t+1)}$ is obtain form $\mathbb{E}^{(t)}$ by a hierarchical mesh update (3.18), a single mesh update (3.10) and a line shift update (3.13) or (3.16). Since the original sequence is a subsequence of the new one, we only need to prove claim (4.2) is true for the new sequence. Since $\mathcal{F}^{h}[\mathbb{E}^{(t)}]>0$ decreases as $t$ increases, the limit $\mathcal{F}^{h}_{\infty}:=\lim_{t\rightarrow\infty}\mathcal{F}^{h}[\mathbb{E}^{(t)}]$ exists. And the change between two steps has the following property,
$$
\lim_{t\rightarrow\infty}\delta_{t}:=\lim_{t\rightarrow\infty}(\mathcal{F}^{h}[\mathbb{E}^{(t)}]-\mathcal{F}^{h}[\mathbb{E}^{(t+1)}])=0.
$$
We still denote $\eta^{(t)}_{i,j}$ by (3.12) for $\mathbb{E}=\mathbb{E}^{(t)}$ . Let $\tau<t$ be the time of the last single mesh update at $(i,j)$ and $\eta^{(\tau)}_{i,j}$ denotes the corresponding value at that update, thus we have $\eta_{i,j}^{(\tau+1)}=0$ . Notice that $\eta^{(t)}_{i,j}$ is a linear combination of $E^{(t)}_{i+\frac{1}{2},j},E^{(t)}_{i+\frac{1}{2},j+1},E^{(t)}_{i,j+\frac{1}{2}}$ and $E^{(t)}_{i+1,j+\frac{1}{2}}$ , whose values is obtained from some previous updates. The update that determines all these values has three cases: the first case is a single mesh update (3.10) for some $(m,n)$ in $\mathbb{E}^{(t^{\prime})}$ with $\tau<t^{\prime}<t$ , the second one is a hierarchical mesh update (3.18) for some $(i_{\min},j_{\min},i_{\max},j_{\max})$ in $\mathbb{E}^{(t^{\prime})}$ with $\tau<t^{\prime}<t$ , the last one is line shift update (3.13) or (3.16) in $\mathbb{E}^{(t^{\prime})}$ with $\tau<t^{\prime}<t$ . Consider the first case. Denote the perturbation is $\eta_{m,n}^{(t^{\prime})}$ , by assumption (2.7), it follows that
$$
\frac{4}{\varepsilon_{max}^{2}}[\eta_{m,n}^{(t^{\prime})}]^{2}\leq||\mathbb{E}^{(t^{\prime}+1)}-\mathbb{E}^{(t^{\prime})}||_{h}^{2}\leq\frac{4}{\varepsilon_{min}^{2}}[\eta_{m,n}^{(t^{\prime})}]^{2},
$$
and by (4.3) and we have the fact that $t^{\prime}\rightarrow\infty$ as $t\rightarrow\infty$ ,
$$
\lim_{t\rightarrow\infty}||\mathbb{E}^{(t^{\prime}+1)}-\mathbb{E}^{(t^{\prime})}||_{h}^{2}=\lim_{t\rightarrow\infty}[\eta_{m,n}^{(t^{\prime})}]^{2}=0.
$$
It follows from equation (3.10) that
$$
|E^{(t^{\prime}+1)}_{i+\frac{1}{2},j}-E^{(t^{\prime})}_{i+\frac{1}{2},j}|\leq\frac{\eta_{m,n}^{(t^{\prime})}}{\varepsilon_{\min}h},
$$
and the variation of $E^{(t^{\prime})}_{i+\frac{1}{2},j+1},E^{(t^{\prime})}_{i,j+\frac{1}{2}}$ and $E^{(t^{\prime})}_{i+1,j+\frac{1}{2}}$ admit the same upper bound as that in (4.5). Since the single part update is calculated by (3.12), we have
$$
|\eta_{i,j}^{(t^{\prime}+1)}-\eta_{i,j}^{(t^{\prime})}|\leq\frac{\varepsilon_{\max}}{\varepsilon_{\min}}\eta_{m,n}^{(t^{\prime})}.
$$
Combining with (4.4), we obtain $\eta_{i,j}^{(t^{\prime}+1)}-\eta_{i,j}^{(t^{\prime})}\rightarrow 0$ as $t\rightarrow\infty$ . Next consider the second case. Denote the perturbation is $\tilde{\eta}^{(t^{\prime})}$ . Similarly, we have
$$
\frac{2(i_{\max}+j_{\max}-i_{\min}-j_{\min})}{\varepsilon_{\max}^{2}}\leq\frac{||\mathbb{E}^{(t^{\prime}+1)}-\mathbb{E}^{(t^{\prime})}||_{h}^{2}}{[\tilde{\eta}^{(t^{\prime})}]^{2}}\leq\frac{2(i_{\max}+j_{\max}-i_{\min}-j_{\min})}{\varepsilon_{\min}^{2}},
$$
and thus
$$
\lim_{t\rightarrow\infty}||\mathbb{E}^{(t^{\prime}+1)}-\mathbb{E}^{(t^{\prime})}||_{h}^{2}=\lim_{t\rightarrow\infty}[\tilde{\eta}^{(t^{\prime})}]^{2}=0.
$$
Similar to the proof of Case 1, it holds that $\eta_{i,j}^{(t^{\prime}+1)}-\eta_{i,j}^{(t^{\prime})}\rightarrow 0$ as $t\rightarrow\infty$ . Last, consider Case 3. We still have $||\mathbb{E}^{(t^{\prime}+1)}-\mathbb{E}^{(t^{\prime})}||_{h}^{2}\rightarrow 0$ and $\eta_{x}^{t^{\prime}}\rightarrow 0$ ( $\eta_{y}^{t^{\prime}}\rightarrow 0$ ) as $t\rightarrow\infty$ . Completely analogous to Cases 1 and 2, $\eta_{i,j}^{(t^{\prime}+1)}-\eta_{i,j}^{(t^{\prime})}\rightarrow 0$ as $t\rightarrow\infty$ holds. Since only finitely many updates occur between $\tau$ and $t$ , thus $\eta_{i,j}^{(t)}\rightarrow\eta_{i,j}^{(\tau+1)}=0$ as $t\rightarrow\infty$ . At this point, the convergence of the algorithm has been established. ∎
## 5 Numerical results
In this section, we conduct numerical experiments to verify the accuracy and effectiveness of the proposed algorithm. We focus on comparing the newly proposed HLR methods with the original single mesh method from multiple perspectives, including metrics like accuracy, number of iteration steps and computational time.
### 5.1 Accuracy and comparison for Poisson’s equation with exact solution
Consider $\Omega=(0,4)\times(0,4)$ with periodic boundary conditions. The Poisson’s equation has the analytical solution for inhomogeneous dielectric permittivity
$$
\phi(x,y)=\cos\left(\frac{\pi x}{2}\right)\sin\left(\frac{\pi y}{2}\right),\quad\text{and}\quad\varepsilon(x,y)=2+\cos\left(\frac{\pi x}{2}\right)\cos\left(\frac{\pi y}{2}\right),
$$
where the source term $\rho$ is analytically computed according to Poisson’s equation $-\nabla\cdot(\varepsilon\nabla\phi)=\rho$ such that the unique solution is given by the first equation of (5.1).
In our calculations, the Poisson’s equation is solved by local algorithms via three distinct relaxation methods over a series of computational grids. Based on the given charge density distribution, the electric field is first initialized to satisfy Gauss’s law by the same procedure. Upon convergence of the relaxation iteration, numerical $L^{\infty}$ errors of the electric fields are calculated. As shown in Table 1, both the forward and zigzag HLR methods achieve second-order convergence with respect to the number of grid points $N_{x}=N_{y}=N$ in terms of $L^{\infty}$ error, a convergence behavior consistent with that of the single mesh update method. Notably, the derived convergence characteristics align with the conclusions reported in [22].
| 32 64 128 | 8.157469E-3 2.051296E-3 5.135733E-4 | - 1.9916 1.9979 | 8.157469E-3 2.051296E-3 5.135728E-4 | - 1.9916 1.9979 | 8.157469E-3 2.051296E-3 5.135728E-4 | - 1.9916 1.9979 |
| --- | --- | --- | --- | --- | --- | --- |
| 256 | 1.284438E-4 | 1.9994 | 1.284400E-4 | 1.9994 | 1.284400E-4 | 1.9994 |
Table 1: The $L^{\infty}$ error and convergence order for three relaxation methods.
To elucidate the convergence behavior of the curl residual, we analyze its spatial profile along the $x$ -direction at the cross-section $y=0.5$ . The residual distributions at the initial stage and after 50 and 100 iterations for three relaxation methods with $N=128$ are compared in Fig. 3. While the single mesh relaxation rapidly suppresses high-frequency residuals, it converges slowly for low-frequency components. In contrast, the HLR methods effectively reduce residuals over a broad range of spatial scales, leading to faster decay of low-frequency error modes. Moreover, the zigzag HLR method consistently shows stronger attenuation of low-frequency modes than the forward HLR method.
*[Error downloading image: 2603.02763v1/figures/ncurl1.png]*
(a) (a) Single mesh
*[Error downloading image: 2603.02763v1/figures/ncurl2.png]*
(b) (b) Forward HLR
*[Error downloading image: 2603.02763v1/figures/ncurl3.png]*
(c) (c) Zigzag HLR
Figure 3: Distribution of the curl residual along the $x$ -direction at the cross-section $y=0.5$ for three relaxation methods. Panels (a)–(c) correspond to single mesh relaxation, forward HLR and zigzag HLR methods, respectively.
### 5.2 Time-dependent Poisson’s equation tests
In most electrolyte transport problems, charge concentrations evolve dynamically over time. This temporal variation necessitates the frequent solution of the Poisson’s equation at every time step to update the electric potential field. To focus exclusively on the numerical solution of the Poisson’s equation itself, we simplify the problem by neglecting the detailed charge transport dynamics and instead abstract the original time-dependent system into a sequence of independent Poisson’s equation solves. We then apply the three aforementioned relaxation schemes to carry out these numerical solutions and compare their performance.
Denote $\Omega=(0,4)\times(0,4)$ with periodic boundary conditions. Let the initial charge density distribution and the homogeneous dielectric permittivity be given by
$$
\rho(x,y)=0,\quad\text{and}\quad\varepsilon(x,y)=1.
$$
Denote the charge density at the $n$ th time step by $\rho^{(n)}$ , which is iteratively given by
$$
\rho^{(n+1)}=\rho^{(n)}+\frac{1}{M^{(n)}}\sum_{k=1}^{16}\left[a_{k}^{(n)}\cos\left(\frac{k\pi x}{2}\right)\sin\left(\frac{k\pi y}{2}\right)+b_{k}^{(n)}\sin\left(\frac{k\pi x}{2}\right)\cos\left(\frac{k\pi y}{2}\right)\right].
$$
for $M^{(n)}=64\sum_{k=1}^{16}(a_{k}^{(n)}+b_{k}^{(n)})$ , and $a_{k}^{(n)},\ b_{k}^{(n)}$ being random variables uniformly distributed on the interval $(0,1)$ .
Fig. 4 shows the average computational time per time step as a function of the grid points $N$ for three relaxation methods and an FFT-based solver, averaged over 100 time steps. At each time step, the iterations are terminated when the tolerance $\varepsilon_{tol}$ reaches $10^{-7}$ (left panel) and $10^{-9}$ (right panel). The results indicate that the HLR methods lead to a clear reduction in computational cost compared to the single mesh method. For $\varepsilon_{tol}=10^{-7}$ , HLR methods are still competitive with the FFT-based solver, requiring roughly eight times the computational time over the range of spatial step sizes considered. Owing to the weak dependence of the iteration count on grid refinement, together with the hierarchical structure of the algorithm, the overall computational complexity is approximately $\mathcal{O}(N\text{log}N)$ .
In the case of spatially varying dielectric permittivity, we set the initial charge density and permittivity are
$$
\rho(x,y)=0,\quad\text{and}\quad\varepsilon(x,y)=2+\cos(\frac{\pi x}{2})\cos(\frac{\pi y}{2}).
$$
The charge density evolves as time according to the same rule as given in (5.2). Fig. 5 presents the average computational time per time step as a function of the grid points $N$ for the three relaxation methods in the case of inhomogeneous dielectric permittivity. As FFT-based methods are not applicable to inhomogeneous media, the computational times for the homogeneous permittivity case are plotted as dashed lines for comparison. It can be seen that HLR methods continue to provide a substantial reduction in computational cost and maintain FFT-like computational complexity. Moreover, compared with the homogeneous dielectric permittivity case, no additional computational complexity is incurred.
*[Error downloading image: 2603.02763v1/figures/cht1.png]*
*[Error downloading image: 2603.02763v1/figures/cht2.png]*
Figure 4: Average iteration time per time step versus the number of grid points $N$ for three relaxation methods and an FFT-based method, averaged over 100 time steps for homogeneous dielectric permittivity case. The stopping tolerance at each time step is set to: $\varepsilon_{tol}=10^{-7}$ (left), $\varepsilon_{tol}=10^{-9}$ (right).
*[Error downloading image: 2603.02763v1/figures/cht3.png]*
*[Error downloading image: 2603.02763v1/figures/cht4.png]*
Figure 5: Average iteration time per time step versus the number of grid points $N$ for three relaxation methods averaged over 100 time steps for inhomogeneous dielectric permittivity case. The stopping tolerance at each time step is set to: $\varepsilon_{tol}=10^{-7}$ (left), $\varepsilon_{tol}=10^{-9}$ (right). The dash lines represent the FFT results of Figure 4.
### 5.3 Poisson–Boltzmann equation tests
The Poisson–Boltzmann (PB) equation governs electrostatic interactions at thermodynamic equilibrium [12, 16, 42, 45]. In the presence of fixed charge density and spatially varying dielectric coefficient, the dimensionless PB equation can be written as
$$
-\nabla\cdot\varepsilon\nabla\phi=\kappa^{2}\left[\rho^{f}+\sum_{\ell}z^{\ell}\frac{N^{\ell}\exp(-z^{\ell}\phi)}{\int_{\Omega}\exp(-z^{\ell}\phi)d\bm{x}}\right].
$$
where $\kappa$ is the dimensionless parameter, $\rho^{f}$ is fixed charge density, $N^{\ell}$ is the total number of the $\ell$ -th charged particle and $z^{\ell}$ is its valence. Notably, solving the PB equation can be reformulated as a constrained minimization problem. The associated energy functional admits minimization via a Maggs-type algorithm [1, 47], yielding a numerical scheme that is not only straightforward and easy to implement, but also unconditionally stable. For a system consisting of $M$ types of particles, it admits minimizing the following energy functional subject to the constraints of discrete Gauss’s law and mass conservation:
$$
\mathcal{F}_{PB}=\int_{\Omega}\left(\frac{\varepsilon}{2}|\bm{E}|^{2}+\sum_{\ell=1}^{M}c^{\ell}\ln c^{\ell}\right)d\bm{x}.
$$
The Maggs-type algorithm for solving this minimization problem consists of two core iterative steps: field updates applied to all finest cell in the grid, and concentration updates executed for all pairs of adjacent movable grids [1]. A key feature of this algorithm is that both update steps rigorously preserve the discrete Gauss’s law and mass conservation constraints, which are essential for the physical consistency of Poisson–Boltzmann (PB) equation solutions. The field update step is identical to that used for solving the standard Poisson’s equation in this work, the only additional component required for the PB equation is the concentration update step. The HLR methods proposed in this paper are designed to accelerate the convergence of the Maggs-type solver, thereby improving the computational efficiency of PB equation solutions.
This section presents three-dimensional ( $3\text{D}$ ) numerical examples to validate the proposed methods. In the numerical example, we denote $\lambda_{0}$ as the fundamental length scale, corresponding to the Lennard–Jones diameter used in the simulation. A spherical colloidal particle with radius $r=3\lambda_{0}$ is placed at the center of a cubic simulation box $(-\frac{L}{2},\frac{L}{2})\times(-\frac{L}{2},\frac{L}{2})\times(-\frac{L}{2},\frac{L}{2})$ with periodic boundary conditions, where $L=30\lambda_{0}$ . The colloidal particle carries a total charge $Ze=60e$ , which is uniformly distributed over its surface. The solvent contains only one ionic species, namely monovalent counterions with valence $z^{1}=-1$ . To maintain overall charge neutrality, counterions with a total charge of $-Ze$ are initially uniformly distributed in the region outside the colloidal particle. The system is discretized on three-dimensional grids with different spatial resolutions, $N_{x}=N_{y}=N_{z}=32,~64$ and $128$ . An iterative stopping criterion is imposed such that the energy change associated with the field update is below $\varepsilon_{tol}=10^{-7}$ . The electrostatic energy functional is then minimized using the single mesh relaxation, forward HLR, and zigzag HLR methods to obtain the equilibrium counterion distribution.
Fig. 6 presents plots of cross-sectional at the plane $z=0$ for the cunterion concentrations at equilibrium obtained using the HLR methods in different spatial resolutions. As expected, due to the positive surface charge on the colloidal sphere, counterion accumulate in its vicinity. Upon reaching equilibrium, the cunterion concentrations obtained from both the forward HLR and zigzag HLR are nearly identical to those obtained from the single mesh update in [1].
*[Error downloading image: 2603.02763v1/figures/PB32.png]*
(a) (a) $32\times 32\times 32$
*[Error downloading image: 2603.02763v1/figures/PB64.png]*
(b) (b) $64\times 64\times 64$
*[Error downloading image: 2603.02763v1/figures/PB128.png]*
(c) (c) $128\times 128\times 128$
Figure 6: Counterion concentration obtained using HLR methods in different discretized grid. Rows from top to bottom represent forward HLR and zigzag HLR method.
Table 2 presents the computational costs of the curl-free step and the total runtime (including concentration moves), for all three update methods upon reaching equilibrium. Compared with the single mesh method, both the forward HLR and zigzag HLR methods reduce the total computational time required to reach the same equilibrium distribution, with particularly noticeable improvements in the curl-free step. For low-resolution grids, the total number of iteration steps required for the zigzag HLR method to reach equilibrium is not significantly smaller than that for the forward HLR method. Moreover, a single local relaxation in the zigzag HLR method involves more computational operations than in the forward HLR method, resulting in a comparable overall computational time between the two methods. However, as the grid resolution increases, the zigzag HLR method requires noticeably fewer iteration steps to reach equilibrium than the forward HLR method. Consequently, despite the slightly higher cost of the curl-free relaxation in the zigzag scheme, the total computational time for the zigzag HLR method becomes lower than that of the forward HLR method.
| 32 | | 13.3s | 35.1s | | 6.7s | 21.4s | | 7.5s | 21.3s |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 64 | | 320.4s | 834.7s | | 171.9s | 530.4s | | 209.4s | 547.7s |
| 128 | | 6494.7s | 16822.0s | | 3767.3s | 11381.7s | | 4006.1s | 10262.0s |
Table 2: Computational time required by single mesh method and HLR methods for the total simulation and for the curl-free relaxation alone.
### 5.4 Poisson-Nernst-Planck equation simulations
The Poisson–Nernst–Planck (PNP) equations provide a time-dependent mean-field continuum framework for modeling ion transport in a wide range of systems [4, 5, 9], describing ionic transport in electrolyte solutions through the coupling of Poisson’s equation with the Nernst–Planck equations. However, the classical mean-field approximation neglects ion–ion correlation effects. To address this limitation, modified PNP models have been developed to incorporate both long-range Coulomb interactions and short-range hard-sphere (HS) correlations [24]. Furthermore, the PNP model can be reformulated within the Maxwell–Ampère Nernst–Planck (MANP) framework [34] to enable charge transport modeling under a local algorithmic formulation. Within the MANP framework, the charge dynamics of a system with $M$ ionic species can be described by
$$
\left\{\begin{aligned} &\frac{\partial c^{\ell}}{\partial t}=\nabla\cdot\gamma^{\ell}\left(\nabla c^{\ell}-z^{\ell}c^{\ell}\bm{E}\right),\quad\ell=1,2,\ldots,M,\\[6.0pt]
&\varepsilon\frac{\partial\bm{E}}{\partial t}=\sum_{\ell=1}^{M}z^{\ell}\gamma^{\ell}\left(\nabla c^{\ell}-z^{\ell}c^{\ell}{\bm{E}}\right)+\boldsymbol{\Theta},\\[6.0pt]
&\nabla\cdot\boldsymbol{\Theta}=0,\\[6.0pt]
&\nabla\times\bm{E}=\mathbf{0},\end{aligned}\right.
$$
where $\gamma^{\ell}>0$ is the diffusion coefficient, $c^{\ell}$ is the concentration of the $\ell$ -th particle and $z^{\ell}$ is its valence. $\boldsymbol{\Theta}$ is introduced as a degree of freedom to enforce $\nabla\times\bm{E}=0$ , which implies the existence of the electric potential satisfying the Poisson’s equation in a connected spatial domain.
Consider a representative numerical test case. All numerical experiments are performed using the same parameters and settings as in [34], with details provided in Appendix A. The evolution process of particle concentration and the changes in the electric field during the evolution process has been studied. The snapshots of the cation concentration $c^{1}$ and electric field $|\bm{E}|$ at two time instants (T = 0.05 and 2) for three different local curl-free relaxation methods are shown in Fig. 7. It can be observed that with all three approaches, the anions are repelled around the negative fixed charges and subsequently diffuse to other areas. The similar evolution patterns observed for the three methods confirm that the HLR methods reproduce the expected behavior and are numerically consistent.
*[Error downloading image: 2603.02763v1/figures/3MANP2.png]*
(a) (a)
*[Error downloading image: 2603.02763v1/figures/3MANP1.png]*
(b) (b)
Figure 7: Snapshots of the cation concentration $c^{1}$ (a) and the magnitude of the electric field $|\bm{E}|$ (b) for different time. From top to bottom, the rows correspond to single mesh relaxation, forward HLR and zigzag HLR, respectively.
The choice of $\boldsymbol{\Theta}$ and its impact on the efficiency of the local curl-free condition has been discussed in [34]. Here, we adopt three specific forms of $\boldsymbol{\Theta}$ , which is given by the following expressions,
| | $\displaystyle\boldsymbol{\Theta}_{1}^{n}=$ | $\displaystyle 0,$ | |
| --- | --- | --- | --- |
where $\mathbf{J}^{\ell,n}=-\gamma^{\ell}(\nabla c^{\ell,n+1}-z^{\ell}c^{\ell,n+1}\bm{E}^{n})$ . Let $N_{x}=N_{y}=256$ , Table. 3 presents the time of relaxation steps required to satisfy the stopping criterion $\varepsilon_{tol}=10^{-8}$ at each time step, for different choices of $\boldsymbol{\Theta}$ and update methods.
It can be observed that, for $\boldsymbol{\Theta}$ taking the value $\boldsymbol{\Theta}_{1}$ and $\boldsymbol{\Theta_{3}}$ , the local algorithm combined with HLR methods requires only a few update rounds per time step, resulting in a significantly lower computational cost compared to the single mesh method. In these cases, the zigzag HLR method performs nearly the same number of iterations as the forward HLR method, and its higher per-iteration cost leads to the largest overall computational time. When $\boldsymbol{\Theta}=\boldsymbol{\Theta}_{2}$ , the single mesh method requires multiple iterations only during the first few steps, after which only a single iteration per time step is needed. Similarly, the HLR methods require on average only one iteration per step. As a result, the total computational time of the single mesh method is comparable to that of the forward HLR method. Since the zigzag HLR method has a higher cost per local relaxation, it remains the most time consuming.
| $\boldsymbol{\Theta}_{1}$ | | 0.0516s | | 0.0116s | | 0.0151s |
| --- | --- | --- | --- | --- | --- | --- |
| $\boldsymbol{\Theta}_{2}$ | | 0.0024s | | 0.0023s | | 0.0036s |
| $\boldsymbol{\Theta}_{3}$ | | 0.0076s | | 0.0041s | | 0.0043s |
Table 3: Average computational time per time step for the forward HLR, zigzag HLR, and single mesh relaxation methods under different values of $\boldsymbol{\Theta}$ .
## 6 Conclusions
We develop a variational algorithm for solving Poisson’s equation, inspired by the local algorithm proposed by Maggs and its underlying dynamical properties. The algorithm consists of two stages. First, an initial electric field satisfying Gauss’s law is constructed, which may not be curl-free. Second, the electric field is iteratively updated through local relaxations to progressively reduce the electrostatic energy until convergence. The limited effectiveness of the single mesh relaxation in damping low-frequency curl modes motivates the development of two hierarchical schemes: the forward HLR and the zigzag HLR method. These hierarchical strategies significantly accelerate convergence by efficiently reducing low-frequency errors. We further establish a rigorous convergence proof for the proposed algorithm, and numerical experiments confirm the high efficiency of the multilevel relaxations. Beyond Poisson’s equation, the proposed multilevel curl-free relaxation schemes demonstrate improved convergence in a broader class of Maggs-based local algorithms, as evidenced by numerical tests on the Poisson–Boltzmann equation and by the application of the MANP framework to the Poisson–Nernst–Planck equation. Owing to the strictly local updates, the proposed algorithms are well suited for large-scale parallel implementations in particle-based simulations. Moreover, the variational and multilevel nature of the framework suggests promising extensions to coupled kinetic field systems, such as the Vlasov-Poisson equation, which will be explored in future work.
## Acknowledgment
Z. Xu and H. Zhou are financially supported by National Key R&D Program of China (grant No. 2024YFA1012403), Natural Science Foundation of China (grants No. 12325113 and 12426304) and SJTU Kunpeng & Ascend Center of Excellence. Q. Yin is partially supported by the Hong Kong Polytechnic University Postdoc Matching Fund Scheme 4-W425.
## Data Availability
Data will be made available on request.
## Declarations
Competing interests The authors have not disclosed any competing interests.
## Appendix A Numerical Parameters and Settings
The benchmark in Section 5.4 took the same setup as [34]. The system is composed of binary ionic species. The fixed charge distribution is defined by the following equation,
| | $\displaystyle\rho^{f}(x,y)=e^{-100\left[\left(x+\frac{1}{2}\right)^{2}+\left(y+\frac{1}{2}\right)^{2}\right]}-e^{-100\left[\left(x+\frac{1}{2}\right)^{2}+\left(y-\frac{1}{2}\right)^{2}\right]}\qquad$ | |
| --- | --- | --- |
The dielectric permittivity is given by
$$
\varepsilon(x,y)=\frac{\varepsilon_{\text{max}}-\varepsilon_{\text{min}}}{2}\Big[\tanh(100d-25)+1\Big]+\varepsilon_{\text{min}},
$$
where $d$ is
$$
d=\begin{cases}\sqrt{\left(x-\frac{1}{2}\right)^{2}+\left(y-\frac{1}{2}\right)^{2}},&0<x\leq 1,\;0<y\leq 1,\\[6.0pt]
\sqrt{\left(x-\frac{1}{2}\right)^{2}+\left(y+\frac{1}{2}\right)^{2}},&0<x\leq 1,\;-1\leq y\leq 0,\\[6.0pt]
\sqrt{\left(x+\frac{1}{2}\right)^{2}+\left(y-\frac{1}{2}\right)^{2}},&-1\leq x\leq 0,\;0<y\leq 1,\\[6.0pt]
\sqrt{\left(x+\frac{1}{2}\right)^{2}+\left(y+\frac{1}{2}\right)^{2}},&-1\leq x\leq 0,\;-1\leq y\leq 0.\end{cases}
$$
Let $\gamma^{\ell}=0.01,\ z^{1}=1,\ z^{2}=-1,\ \varepsilon_{\text{min}}=2\times 10^{-4},\ \varepsilon_{\text{max}}=1.56\times 10^{-2}$ and the initial ion concentration is:
$$
c^{\ell}(x,y,0)=0.1,~\ell=1,2.
$$
## References
- [1] M. Baptista, R. Schmitz, and B. Dünweg (2009) Simple and robust solver for the Poisson-Boltzmann equation. Phys. Rev. E 80, pp. 016705. Cited by: §1, §3.2, §5.3, §5.3, §5.3.
- [2] A. H. D. Cheng and D. T. Cheng (2005) Heritage and early history of the boundary element method. Eng. Anal. Bound. Elem. 29, pp. 268–302. Cited by: §1.
- [3] D. E. Clapham (2007) Calcium signaling. Cell 131, pp. 1047–1058. Cited by: §1.
- [4] B. Corry, S. Kuyucak, and S. Chung (2003) Dielectric self-energy in Poisson-Boltzmann and Poisson-Nernst-Planck models of ion channels. Biophys. J. 84, pp. 3594–3606. Cited by: §5.4.
- [5] H. Daiguji, Y. Oka, and K. Shirono (2005) Nanofluidic diode and bipolar transistor. Nano Lett. 5, pp. 2274–2280. Cited by: §5.4.
- [6] L. Desiderio, S. Falletta, M. Ferrari, and L. Scuderi (2022) CVEM-BEM coupling with decoupled orders for 2d exterior Poisson problems. J. Sci. Comput. 92, pp. 96. Cited by: §1.
- [7] G. Dhatt, E. Lefrançois, and G. Touzot (2012) Finite element method. John Wiley & Sons. Cited by: §1.
- [8] P. Duhamel and M. Vetterli (1990) Fast Fourier transforms: a tutorial review and a state of the art. Signal Process 19, pp. 259–299. Cited by: §1.
- [9] B. Eisenberg and W. Liu (2007) Poisson–Nernst–Planck systems for ion channels with permanent charges. SIAM J. Math. Anal. 38, pp. 1932–1966. Cited by: §5.4.
- [10] F. Fahrenberger, Z. Xu, and C. Holm (2014) Simulation of electric double layers around charged colloids in aqueous solution of variable permittivity. J. Chem. Phys. 141, pp. 064902. Cited by: §1.
- [11] H. Feng, G. Long, and S. Zhao (2021) FFT-based high order central difference schemes for Poisson’s equation with staggered boundaries. J. Sci. Comput. 86, pp. 7. Cited by: §1.
- [12] G. Gouy (1910) Constitution of the electric charge at the surface of an electrolyte. J. Phys. 9, pp. 457–468. Cited by: §5.3.
- [13] W. S. Hall (1994) Boundary element method. In The boundary element method, pp. 61–83. Cited by: §1.
- [14] P. W. Hemker (1981) Introduction to multigrid methods. Nieuw Arch. Wisk. 3, pp. 71–101. Cited by: §1, §3.3.
- [15] R. W. Hockney and J. W. Eastwood (2021) Computer simulation using particles. CRC Press. Cited by: §1.
- [16] B. Honig and A. Nicholls (1995) Classical electrostatics in biology and chemistry. Science 268, pp. 1144–1149. Cited by: §1, §5.3.
- [17] A. F. Izmailov and M. Solodov (2003) Karush-Kuhn-Tucker systems: regularity conditions, error bounds and a class of Newton-type methods. Mathematical Programming 95, pp. 631–650. Cited by: §2.
- [18] J. Jiang and Z. Wang (2018) Improved local lattice Monte Carlo simulation for charged systems. J. Chem. Phys. 148, pp. 114105. Cited by: §3.2.
- [19] H. Kim, L. Chacon, and G. Lapenta (2005) Fully implicit particle-in-cell algorithm.. Bull. Am. Phys. Soc.. Cited by: §1.
- [20] A. Kirsch and F. Hettlich (2015) The mathematical theory of Time-Harmonic Maxwell’s equations. Appl. Math. Sci. 190, pp. 20. Cited by: §2.
- [21] L. Levrel and A. C. Maggs (2005) Monte Carlo algorithms for charged lattice gases. Phys. Rev. E 72, pp. 016715. Cited by: §3.2.
- [22] B. Li, Q. Yin, and S. Zhou (2024) Finite-difference approximations and local algorithm for the Poisson and Poisson-Boltzmann electrostatics. arXiv preprint arXiv:2409.15796. Cited by: Theorem 2.1, §2, §3.2, §4, §4, §5.1.
- [23] J. Liang, P. Liu, and Z. Xu (2021) A high-accurate fast poisson solver based on harmonic surface mapping algorithm. Commun. Comput. Phys. 30, pp. 210–226. Cited by: Remark 1.
- [24] M. Ma, Z. Xu, and L. Zhang (2021) Modified Poisson–Nernst–Planck model with Coulomb and hard-sphere correlations. SIAM J. Appl. Math. 81, pp. 1645–1667. Cited by: §5.4.
- [25] A. C. Maggs and V. Rossetto (2002) Local simulation algorithms for Coulomb interactions. Phys. Rev. Lett. 88, pp. 196402. Cited by: §1, §3.2.
- [26] A. C. Maggs (2002) Dynamics of a local algorithm for simulating Coulomb interactions. J. Chem. Phys. 117, pp. 1975–1981. Cited by: §1, §1, §3.2.
- [27] B. Maury (2001) A fat boundary method for the Poisson problem in a domain with holes. J. Sci. Comput. 16, pp. 319–339. Cited by: §1.
- [28] S. F. McCormick (1987) Multigrid methods. SIAM. Cited by: §1, §3.3.
- [29] G. Mesa, E. Dobado-Fuentes, and J. J. Sáenz (1996) Image charge method for electrostatic calculations in field-emission diodes. J. Appl. Phys. 79, pp. 39–44. Cited by: §1.
- [30] J. Milaszewicz (1987) Improving Jacobi and Gauss-Seidel iterations. Linear Algebra Appl. 93, pp. 161–170. Cited by: §1.
- [31] P. Monk (2003) Finite element methods for maxwell’s equations. Oxford university press. Cited by: §2.
- [32] H. J. Nussbaumer (1982) The fast fourier transform. Springer. Cited by: §1.
- [33] I. Pasichnyk and B. Dünweg (2004) Coulomb interactions via local dynamics: A molecular-dynamics algorithm. J. Phys.: Condens. Matter 16, pp. S3999. Cited by: §1.
- [34] Z. Qiao, Z. Xu, Q. Yin, and S. Zhou (2023) A Maxwell–Ampère Nernst–Planck framework for modeling charge dynamics. SIAM J. Appl. Math. 83, pp. 374–393. Cited by: Appendix A, §1, §5.4, §5.4, §5.4.
- [35] Z. Qiao, Z. Xu, Q. Yin, and S. Zhou (2023) Structure-preserving numerical method for maxwell-ampère nernst-planck model. J. Comput. Phys. 475, pp. 111845. Cited by: §1.
- [36] Z. Qiao, Z. Xu, Q. Yin, and S. Zhou (2025) Intrinsic local Gauss’s law preserving PIC method: a self-consistent field-particle update scheme for plasma simulations. arXiv preprint arXiv:2506.02407. Cited by: §1.
- [37] J. Rottler and A. C. Maggs (2004) A continuum, O (N) Monte Carlo algorithm for charged particles. J. Chem. Phys. 120, pp. 3119–3129. Cited by: §1, §1.
- [38] B. Roux (2002) Theoretical and computational models of ion channels. Curr. Opin. Struct. Biol. 12, pp. 182–189. Cited by: §1.
- [39] J. W. Ruge and K. Stüben (1987) Algebraic multigrid. In Multigrid methods, pp. 73–130. Cited by: §1.
- [40] T. Sorokina and S. Zhang (2022) An interpolated Galerkin finite element method for the Poisson equation. J. Sci. Comput. 92, pp. 47. Cited by: §1.
- [41] K. Stüben (2001) A review of algebraic multigrid. Numerical Analysis: Historical Developments in the 20th Century, pp. 331–359. Cited by: §1.
- [42] Y. Tu, Q. Pang, H. Yang, and Z. Xu (2022) Linear-scaling selected inversion based on hierarchical interpolative factorization for self Green’s function for modified Poisson-Boltzmann equation in two dimensions. J. Comput. Phys. 461, pp. 110893. Cited by: §5.3.
- [43] F. H. van der Heyden, D. Stein, K. Besteman, S. G. Lemay, and C. Dekker (2006) Charge inversion at high ionic strength studied by streaming currents. Phys. Rev. Lett. 96, pp. 224502. Cited by: §1.
- [44] E. Waisman and J. L. Lebowitz (1972) Mean spherical model integral equation for charged hard spheres i. method of solution. J. Chem. Phys. 56, pp. 3086–3093. Cited by: §1.
- [45] Z. Xu and A. C. Maggs (2014) Solving fluctuation-enhanced Poisson–Boltzmann equations. J. Comput. Phys. 275, pp. 310–322. Cited by: §5.3.
- [46] Q. Zheng, D. Chen, and G. W. Wei (2011) Second-order Poisson–Nernst–Planck solver for ion transport. J. Comput. Phys. 230, pp. 5239–5262. Cited by: §1.
- [47] S. Zhou, Z. Wang, and B. Li (2011) Mean-field description of ionic size effects with non-uniform ionic sizes: A numerical approach. Phys. Rev. E 84, pp. 021901. Cited by: §1, §5.3.
- [48] O. C. Zienkiewicz, R. L. Taylor, P. Nithiarasu, and J. Z. Zhu (1977) The finite element method. Elsevier. Cited by: §1.