Numerical Solution to Poisson’s Equation for Estimating Electrostatic Properties Resulting from an Axially Symmetric Gaussian Charge Density Distribution: Charge in Free Space

: Poisson’s equation frequently emerges in many fields, yet its exact solution is rarely feasible, making the numerical approach notably valuable. This study aims to provide a tutorial-level guide to numerically solving Poisson’s equation, focusing on estimating the electrostatic field and potential resulting from an axially symmetric Gaussian charge distribution. The Finite Difference Method is utilized to discretize the desired spatial domain into a grid of points and approximate the derivatives using finite difference approximations. The resulting system of linear equations is then tackled using the Successive Over-Relaxation technique. Our results suggest that the potential obtained from the direct integration of the distance-weighted charge density is a reasonable choice for Dirichlet boundary conditions. We examine a scenario involving a charge in free space; the numerical electrostatic potential is estimated to be within a tolerable error range compared to the exact solution


Introduction
Poisson's equation stands as a fundamental partial differential equation that plays a crucial role in different fields of physics.In electrostatics, Poisson's equation establishes a connection between the electric field ( ⃗ E), the electrostatic potential (Φ), and the electric volume charge density (ρ) within a given region of space.Mathematically, Poisson's equation is ∇ 2 Φ = −ρ/ε 0 , where ε 0 is the permittivity of free space and the electric field emerges as the gradient of the electrostatic potential ( ⃗ E = − ⃗ ∇Φ).The distribution of the electric charge density is essential in defining the source term within Poisson's equation, thereby directly affecting its solution of the electrostatic potential and the resultant electric field within a designated spatial domain.Typically, regions with higher electric charge density exhibit stronger electric fields, resulting in more noticeable potential gradients and larger variations in potential.The distribution of electric charge density also defines the shape of equipotential surfaces within the region.Uniform charge densities typically yield concentric equipotential surfaces centered around the charge, while nonuniform distributions lead to asymmetric equipotential surfaces.Furthermore, the distribution of electric charge density is essential for understanding and determining the boundary value problems, where specific potentials are defined at the region's boundaries.It rules the behavior of the electrostatic potential and the electric field within the given region.
This introductory study aims to numerically solve Poisson's equation to estimate the electrostatic potential and the electric field.The focus will be on an axially symmetric source charge density that conforms to a Gaussian distribution described by the general formula ρ(r) = ρ 0 exp(−r 2 /2σ 2 ).Here, ρ 0 is the maximum of the charge density at the center of the distribution, and σ is a parameter that controls the rate at which the electric charge density decreases as the distance from the center increases.By employing the relation Q = ρdτ and applying Gaussian integrals [1][2][3][4], it can be shown that ρ 0 = Q/(2π) 3 2 σ 3 .Accordingly, the general form of Poisson's equation is In cases of axial symmetry, the charge density depends only on the radial distance from the center, and it is symmetric around the axis of symmetry.Consequently, despite the typical spherical symmetry associated with Gaussian charge distributions, the geometry of the desired spatial domain is more closely aligned with cylindrical coordinates rather than spherical coordinates.Our choice to focus on a Gaussian charge distribution in general and its axially symmetric form stems from several reasons.Firstly, the analytical solution for such charge distribution can be obtained using the error function, which will be used to validate our numerical approach.Secondly, such distribution is commonly utilized to interpret various natural phenomena.For instance, it is crucial to investigate atmospheric electricity, thunderstorm dynamics, and associated phenomena.While the charge distribution within thunderstorms may not perfectly conform to a Gaussian model, the Gaussian distribution offers a simplified approach for approximating charge density within particular regions of thunderstorms [5,6].This approximation assists researchers in unraveling the mysteries of thunderstorms, refining lightning forecasts, and advancing our understanding of the Earth's atmosphere.

Finite Differencing of Poisson's Equation
From a numerical standpoint, the choice of the coordinate system is contingent upon several factors, including the distribution of the charge density, the geometry of the spatial domain, and the specified boundary conditions.In the context of this study, where an axially symmetric charge density distribution is examined, it is noteworthy that such a distribution naturally exhibits cylindrical symmetry and remains unaltered by variations in the azimuthal angle (ϕ).Consequently, employing cylindrical coordinates allows us to disregard this angle, effectively reducing the problem's dimensions from three to two for computational efficiency.Thus, Poisson's equation can be reformulated as a twodimensional problem in cylindrical coordinates as follows: To numerically solve Equation (2), we employ the Finite Difference Method (FDM) for its accuracy and simplicity [7][8][9].The problem is tackled within a rectangular domain in the (sz)-plane, where both s and z span from their initial values s 0 and z 0 to their final values s f and z f , respectively.A uniform grid or mesh is applied with equal increments of ∆s in the s-direction and ∆z in the z-direction, as shown in Figure 1.This approach involves dividing the sz-plane into a grid of mesh points, each with coordinates (s, z).These coordinates are defined by the constant spacing of ∆s and ∆z in the respective s and z directions, where Here, n s and n z represent the total number of mesh elements along the s and z directions, respectively.They are calculated based on the increments ∆s = (s f − s 0 )/n s and ∆z = (z f − z 0 )/n z .
Figure 1.The 2D grid in the sz-plane for FDM.The sweep direction in SOR refers to the sequence in which grid points are updated during each iteration.While the specific order of the sweep is irrelevant, it is essential to maintain consistency once established [10].In this context, updates in each iteration originate from the bottom-left corner of the grid and proceed horizontally to the right, advancing row by row from the bottom to the top.
At any point (i, j) representing a specific location in the domain, the exact potential is Φ(s i , z j ).In contrast, the estimated value Φ i,j can be determined using the second-order central finite difference approximations [11][12][13] of the derivatives of Φ at the same point.The FDM converts the differential equation into a system of algebraic equations by replacing the derivatives in the differential equation with finite differences at the mesh points.The accuracy of the FDM is derived from the fact that as the mesh size decreases, the solution at the mesh points converges to the exact solution of the differential equation.Meanwhile, the simplicity of the FDM is evident in its easy implementation.The discretized form of Equation ( 2) for the mesh points within the solution region can be expressed as follows: The above equation can be restructured as follows: where A i,j = B i,j = (1/∆z 2 )/(α), C i,j = ((1/2s i ∆s) − (1/∆s 2 ))/(α), D i,j = ((1/2s i ∆s) + (1/∆s 2 ))/(α), S i,j = −(ρ i,j /ε 0 )/(α), and α = −((2/∆s 2 ) + (2/∆z 2 )).
Deriving Equation (3) applies to grid points within the solution region.However, a different approach is used for mesh points on the symmetry axis.On this symmetry axis, the derivative of the potential with respect to s, which is evaluated at s = 0 and represented as ∂Φ/∂s| s=0 , is zero.By employing L'Hospital's rule, this condition leads to the following expression: As a result, when Equation ( 2) is evaluated at s = 0, it transforms into In addition, when s = 0, it implies that Φ i−1,j = Φ i+1,j .Therefore, the discretized form of Equation ( 4) can be rewritten as where a i,j = b i,j = (1/∆z 2 )/(α 0 ), d i,j = (4/∆s 2 )/(α 0 ), S 0 i,j = −(ρ i,j /ε 0 )/(α 0 ), and α 0 = −((4/∆s 2 ) + (2/∆z 2 )).Equation ( 5) provides a more straightforward way to calculate the potential at discrete points on the symmetry axis in terms of the potential at the neighboring grid points.

Successive Over-Relaxation (SOR) Method
Equations ( 3) and ( 5) constitute a large system of linear equations.Such a system can be effectively tackled by employing the SOR method.This method initiates with an initial approximate guess of the solution that serves as the basis for generating subsequent solution estimates iteratively [12,14,15].The value of Φ i,j at the (k + 1) th iteration is determined by utilizing the known (already updated) values of Φ i,j at the k th iteration, following the sweep directions depicted in Figure 1.Here, the superscript k (k = 0, 1, 2, 3, ... ) denotes the iteration number.Another crucial aspect of the SOR method involves accelerating the convergence of the iterative process by introducing an additional parameter termed the relaxation factor (ω).This factor accommodates the difference between the current (at the (k + 1) th ) and the already updated (at the k th ) potential values, thereby regulating the convergence rate.The estimated current potential is given by [10]: The performance of Equation ( 6) depends on the selection of ω, which typically ranges between 1.0 and 2.0.However, an optimal choice can markedly expedite convergence.The most favorable optimum value of the relaxation factor (ω opt ) for the constant grid spacing numerical setup is provided by the following expression [10]: where ξ = [{cos(π/n s ) + β 2 cos(π/n z )}/1 + β 2 ] 2 , and β = △z/△s is the grid aspect ratio.
The potential Φ k+1 is indeed considered a reliable estimate of the potential if it converges to a stable solution.The convergence to a stable solution implies that the iterative process of updating the potential has reached a point where further iterations do not significantly change the potential.This is a desirable property in numerical methods as it indicates that the solution is not oscillating or diverging but instead settling into a consistent value.The convergence criterion applied in this study stipulates that the weighted root mean square (WRMS) should be less than 1 [16], where Here, the multiplicative weights W i,j = 1/(RT Φ k+1 i,j + AT) are calculated by using the relative tolerance (RT) and absolute tolerance (AT) that are specified by the user.The WRMS is a measure of the magnitude of the residuals (the differences between the current and updated values) and provides an indication of the accuracy of the potential Φ k+1 .If the WRMS is less than 1, it suggests that the residuals are small and the estimated potential is converging.

Boundary Conditions
The SOR method is a powerful technique for solving large systems of linear equations that often arise from the FDM.As previously discussed, selecting the appropriate ω is crucial for ensuring the convergence of the SOR method.However, it is essential to recognize that convergence alone does not guarantee accuracy.To achieve accurate solutions, it is also essential to enforce Dirichlet boundary conditions, which specify the values of the solution at the boundaries of the desired spatial domain.These boundary conditions are critical as they provide the necessary constraints that shape the solution within the computational domain.Therefore, the optimal choice of ω and the precise application of Dirichlet boundary conditions are essential for obtaining reliable and accurate results when using the SOR method in conjunction with the FDM.
The boundary conditions must be carefully chosen to reflect the physical problem being addressed accurately.Poorly chosen boundary conditions can result in solutions that, despite being mathematically valid, fail to represent the physical reality accurately.One common and straightforward technique for enforcing Dirichlet boundary conditions is to set a constant value, typically zero, for the potential at the boundaries.This approach is convenient and computationally efficient, making it a practical choice for initial approximations or situations where high accuracy is not critical.However, for problems where accuracy is paramount, more elaborate boundary conditions that closely emulate the physical context are vital.
This study determines the boundary conditions using the method proposed by Liu and Pasko [17].This method was initially designed to numerically model the dynamics of positive streamers produced by a space charge and developing at various air pressures.Although the spatial domain used in their study was relatively small, the proposed method can be extended to larger domains.The method begins by calculating the direct integral solution of the electric potential at the boundaries of the simulation domain.Subsequently, the potential at interior points is computed using the previously discussed SOR method, with the known boundary potentials serving as inputs.The potential at the boundaries is obtained by directly integrating the distance-weighted charge density, with the result of this integration then multiplied by the complete elliptic integral of the first kind as follows: where the subscript "b" refers to "at boundary", (1 − k 2 sin 2 α) −1 dα is the complete elliptic integral of the first kind [18][19][20].This integral approach ensures that the boundary conditions accurately reflect the physical distribution of charge, yielding more accurate solutions within the simulation domain when combined with the SOR method.By incorporating this approach, the model inherently accounts for the spatial distribution of the charge density in determining the boundary potentials, thus enhancing the realism and reliability of the electrostatic field simulations.This method provides a robust framework for establishing boundary conditions.It improves the accuracy and fidelity of the numerical solutions obtained for Poisson's equation, making it highly appropriate for complex physical scenarios.

The Numerical Electric Field ⃗ E
Once the numerical distribution of the electrostatic potential is estimated, the subsequent calculation of the electric field becomes more straightforward.This process is facilitated by using cylindrical coordinates and the inherent cylindrical symmetry.In this context, the partial derivative of the electrostatic potential with respect to the azimuthal angle ( ∂Φ ∂ϕ ) results in zero due to the symmetry.Consequently, the formula for the electric field can be expressed as follows: Numerically, the electric field at points within the spatial domain is computed using central difference formulas [21].These formulas offer a balanced and accurate method for estimating derivatives and gradients, making them well suited for analyzing the behavior of the electric field within the domain.Conversely, at the boundaries, one-sided formulas, which include both forward difference and backward difference formulations, are employed to compute the electric field [21].These one-sided formulas are suitable for handling boundary points, allowing for the electric field to be determined while considering the boundary conditions.The above-mentioned methods are delineated below:

Application and Results: Charge in Free Space
In this study, we examine a Gaussian charge density distribution with a total charge of Q = 1 C and σ = 100 m, centered along the z-axis at a height of 1 km.We aim to estimate the electrostatic potential and the electric field within a specified rectangular spatial domain.The boundaries of this domain are defined by a radial distance, s, ranging from 0 to 1 km, and a vertical height, z, varying from 0 to 2 km.Consequently, Poisson's equation, as denoted by Equation ( 1), takes the following form: The rectangular domain is divided into a uniform mesh, maintaining a consistent spacing of △s = △z = σ/10 = 10 m.This setup yields n s = 100 and n z = 200, leading to β = 1 and ξ = 0.99.Although non-uniform spacing in the radial and vertical directions (β ̸ = 1) is possible, opting for square cells with a grid aspect ratio of β = 1 is beneficial.This uniformity simplifies the implementation of the FDM because consistent finite difference approximations can be applied across the entire grid, reducing code complexity and minimizing potential discretization errors [14].Additionally, uniform grid spacing helps maintain consistent accuracy throughout the computational domain.When grid cells are equal in size, truncation errors (errors from approximating derivatives) are evenly distributed, enhancing the overall solution accuracy.On the other hand, non-uniform grids can lead to varying truncation errors that complicate error analysis and reduce accuracy [22].Furthermore, many numerical methods, including the FDM, have stability requirements dependent on grid spacing, and a uniform grid helps ensure this stability.Non-uniform grids can introduce interpolation or discretization errors, which are minimized with a square grid [23].Furthermore, ω opt is determined to be 1.95.To ensure convergence, RT is set at 10 −6 , while AT of 0.01 is established.
The charge density distribution being examined is depicted in Figure 2, with distances expressed in units of σ, a convention maintained in all subsequent figures within this study.Normalized to its maximum value (ρ/ρ 0 ), the charge density distribution peaks at the coordinates (0, 10σ).Subsequently, it gradually diminishes as radial and vertical distances extend, eventually fading near the boundaries.The electrostatic potential derived numerically from Equation ( 12) within the defined spatial domain is represented in Figure 3.The solution reflects the bell curve characteristic of the Gaussian charge density distribution.The electrostatic potential peaks slightly above 7 × 10 7 V, coinciding with the maximum charge density.As radial and vertical distances increase, the electrostatic potential gradually diminishes, following the exponential decay pattern inherent in the Gaussian distribution.The decrease in the electrostatic potential is evident in its values, which drop to less than 1 × 10 7 V at the boundaries, highlighting the diminishing influence of the charge density distribution.Furthermore, it is noticeable that the spatial extent of the charge density distribution plays a crucial role in shaping the equipotential surfaces.As shown in Figure 3, these surfaces adopt concentric configurations around the charge maxima, highlighting the symmetrical distribution nature.Also, moving away from the center of the charge distribution leads to forming more expansive equipotential surfaces that exhibit smoother gradients, indicating a more uniform distribution of the electrostatic potential across the spatial domain.
Following the numerical estimation of the electrostatic potential, the resulting electric field, computed using Equation (11), is presented in Figure 4.In particular, Figure 4a illustrates the radial component of the electric field (E s ).Interestingly, this component is finite and non-zero, maintaining a value of approximately 1 × 10 4 V/m at the position corresponding to the peak of the charge density distribution.The non-zero E s is attributed to the domain setup and the charge distribution, where the necessary symmetry for this component to cancel out is absent.Beyond the coordinates (0, 10σ), the magnitude of E s gradually intensifies, reaching a peak slightly over 19 × 10 4 V/m near the coordinates (1.5σ, 10σ).Subsequently, the magnitude of E s diminishes, nearly approaching zero as the boundaries are neared.In contrast, Figure 4b illustrates the vertical component of the electric field (E z ).Above the height of 10σ, E z directs upwards, while below the same height, it points downwards.Unlike E s , the magnitude of E z cancels out, resulting in a net zero at the position corresponding to the peak of the charge density distribution.However, the magnitude of E z gradually amplifies when moving away from the coordinates (0, 10σ).Its peaks are observed approximately at (0, 10 ± 1.5σ).Like E s , the magnitude of E z diminishes rapidly as the boundaries are approached.Figure 4c presents the total electric field E, computed as the combination of E s and E z components, with its magnitude determined by the formula E = (E 2 s + E 2 z ).The non-zero E s at the (0, 10σ) contributes to a finite value of approximately 1 × 10 4 V/m for E at this location.Moving away from (0, 10σ), E undergoes changes.It first increases rapidly, peaking slightly over 19 × 10 4 V/m at about 1.5σ away from (0, 10σ).This peak suggests a strong electric field in this region due to the high change in the electrostatic potential.Beyond the distance of 1.5σ, E decreases gradually.This reduction is especially evident in regions with smoother gradients of electrostatic potential, which usually indicate a decrease in charge density and, consequently, a reduction in the magnitude of the overall electric field.To ascertain the accuracy of our numerical approach, we compare the numerically estimated electrostatic potential with its exact analytical counterpart.This comparison entails calculating the relative error between the numerical and analytical quantities across all mesh points within the spatial domain under consideration.From an analytical perspective, the electrostatic potential that arises from the given Gaussian charge distribution at a certain point in space can be calculated by integrating the contributions (dΦ) from all infinitesimal charge elements (dQ).Upon performing the necessary mathematical operations, it can be demonstrated that the following equation gives the exact electrostatic potential: Here, erf is the error function.
As previously mentioned, applying Dirichlet's boundary conditions is essential in solving the large system of equations produced by the SOR method.These conditions are crucial in guiding the numerical approach toward a meaningful solution.In some scenarios, opting for Dirichlet's boundary conditions by setting the potential at the boundaries to zero may be advantageous.However, while this simplification may make the problem more tractable, it might not fully capture the accurate solution.To investigate this scenario further, we examine the resulting relative error under the conditions outlined in this study, as depicted in Figure 5a, when employing these zero-potential boundary conditions.The figure shows a non-uniform distribution of relative error across the domain, with the error reaching 100% at the boundaries.While this error is expected due to the assumption of zero potential at the boundaries in the numerical solution, the relative error is approximately 12% at the location corresponding to the peak of the electrostatic potential at (0, 10σ).This discrepancy suggests that the numerical solution underestimates the true peak electrostatic potential by roughly 12%.Conversely, Figure 5b demonstrates the relative error when applying Dirichlet's boundary conditions obtained from Equation ( 9).Even with the nonuniform distribution throughout the considered spatial domain, the relative error remains under 0.5% at the position corresponding to the peak of the electrostatic potential and peaks at less than 3% at the boundaries.Such a minimal relative error suggests that the Dirichlet's boundary conditions computed using the charge density integral method, as per Equation ( 9), can be considered an appropriate choice, leading to a satisfactory solution.
As depicted in Figure 6, our study suggests that while a relative error at the boundaries of approximately 3% is within acceptable limits, reducing the grid size can further decrease this error.However, this improvement comes at the cost of increased computational time.Reducing the grid size necessitates additional calculations, which, in turn, demand more processing power and time.This trade-off between accuracy and computational efficiency becomes particularly significant when dealing with larger spatial domains.As the analysis extends to larger areas, the reduction in grid size can markedly increase the average computational time, potentially becoming a limiting factor.Therefore, selecting the grid size should balance the desired accuracy with the available computational resources.In our study's example, using a 10-meter grid size over a spatial domain spanning 2 km vertically and 1 km horizontally, the average computational time remained modest, at approximately 1.2 s on an 8-core M1 CPU.Our findings suggest that, for instance, halving the grid size reduces the relative error by a factor of about 3, but also increases the average computational time by about a factor of 10, as shown in Figure 6.
Finally, it is essential to note that while this study focuses on a Gaussian charge density distribution, which is a simplification that may not accurately represent all physical scenarios, real-world charge distributions can be more complex and may not conform to a Gaussian profile.The model used in this study allows other axially symmetric distributions to be explored, such as those that fall off exponentially (ρ ∼ exp −r ) or follow a power law (ρ ∼ r −α ).Additionally, this model can be extended to fully three-dimensional simulations, allowing for studying more complex charge distributions and geometries.Although this extension would require significant computational resources, it could provide more comprehensive insights into electrostatic properties.

Conclusions
Poisson's equation is a fundamental component in physics.When an exact solution is unattainable, a numerical solution becomes essential.Due to its simplicity and accuracy, the FDM stands out in transforming Poisson's equation into a system of algebraic equations.Thanks to its effectiveness, the resulting system can be effectively addressed using the widely favored SOR iterative technique.To ensure a meaningful and convergent solution, it is crucial to apply appropriate boundary conditions.While simplifying the problem

Figure 2 .
Figure 2. A Gaussian charge density distribution with axial symmetry, characterized by a total charge of Q = 1 C and σ = 100 m, scaled to its maximum value (ρ/ρ 0 ).The distribution is positioned along the z-axis at a height of 10σ and spans a rectangular spatial domain measuring 20σ in height and 10σ in radial extent.

Figure 4 .
The distribution of the electric field.(a) The radial component of the electric field (E s ).(b) The vertical component of the electric field (E z ).(c) The total electric field (E).

Figure 5 .
The relative error between the numerical and analytical electrostatic potential values across all mesh points.(a) Zero-potential boundary conditions employed.(b) Equation(9)-derived boundary conditions employed.

Figure 6 .
Figure 6.The relation between the grid size (horizontal axis), the relative error at the boundaries (left vertical axis), and the average computational time on an 8-core M1 CPU processor (right vertical axis).