A High-Accuracy Algorithm for Solving Problems of Electrostatics in a Nonhomogeneous Spatially Periodic Dielectric Medium

A high-accuracy economical iterative method is proposed for calculating the potential and the strength of the electric field in a three-dimensional inhomogeneous spatially periodic dielectric placed in an initially uniform electric field. The idea underlying the algorithm is that the potential is represented as a sum of a linear function and a spatially periodic correction, which can be expressed as an expansion in eigenfunctions of the Laplace operator that satisfy the appropriate periodicity conditions. The fast Fourier transform is used for an efficient numerical implementation of the proposed algorithm.

1. Problems of electrostatics in spatially periodic media naturally arise in physical chemistry and material science when the state of a system is controlled by applying an electric field [1]. The numerical solution of electrostatic problems is based mainly on three groups of methods, namely: (a) finite difference methods; (b) projective methods, such as variants of the Galerkin method; and (c) integral equation methods [2][3][4][5]. Along with them, spectral methods are used in a number of specific problems, in which case the sought solution is represented as an expansion in terms of eigenfunctions of the corresponding differential operator [4]. The spectral approach is effective for problems with periodic boundary conditions in a homogeneous medium, since the analytical solution can then be represented in the form of an expansion in terms of trigonometric functions. However, the situation is different in the case of inhomogeneous media, since the trigonometric functions are generally not the eigenfunctions of the differential operator of the problem in that case. Nevertheless, in spatially periodic problems, the trigonometric basis remains attractive for solution representation in the form of a linear combination of trigonometric functions when the characteristics of the medium are smoothly varying functions. In this paper, an analogue of the pseudospectral approach applied to problems in the theory of polymers [6] is proposed for solving electrostatic problems in spatially periodic dielectric media with smoothly varying permittivity. Based on this approach and the fast Fourier transform, a high-accuracy low-cost iterative method is developed for computing the electric field potential and strength.
2. The electric field potential Φ(r) in a medium with a sufficiently smooth positive permittivity ε(r), where r = (x, y, z) is the radius vector, satisfies the elliptic partial differential equation [1] (1) In the case of an inhomogeneous dielectric placed in an initially uniform electric field, it is convenient to represent the sought potential Φ(r) in the form (2) where E 0 is the electric field strength in a homogeneous dielectric medium with permittivity equal to the mean value of the permittivity ε(r) and ϕ(r) is an additional potential. Substituting the right-hand side of (2) into Eq. (1) leads to the equation The desired additional potential ϕ(r) on the righthand side of (2) satisfies the following periodic boundary conditions in the spatially periodic dielectric medium: (4) where f is either ε(r) or ϕ(r), Z is the set of integers, and {a k } is a basis in E 3 .
Note that the Poisson equation for the function ϕ(r), i.e., (5)  where ϕ(r) and ρ(r) satisfy the periodic conditions (4), determines the solution up to a constant component. Therefore, Eq. (5) has a unique solution in the subspace of potentials ϕ(r) with a zero mean. For such a solution of problem (5) with a zero constant component and boundary conditions (4), we use the notation (6) Equation (3) also determines the additional potential ϕ(r) up to a constant component, whose value can be set, without loss of physical generality, to zero. Substituting the right-hand side of Eq. (3) for the function ρ(r) in Eq. (6), we obtain another form of Eq. (3): In contrast to Eq. (3), Eq. (7) implies that the constant component of the potential ϕ(r) is zero.
Equation (7) with periodic boundary conditions can be solved by applying the Picard iteration method with a positive parameter τ: (8) or by applying other, more rapidly convergent iterative methods, including multilayered ones [7,8].
3. Consider an efficient method for computing the operator F[ϕ(r); ε(r)] with the help of the fast Fourier transform (FFT). For this purpose, we introduce a reciprocal (inverse) basis {b j } determined by the conditions (9) where δ kj is the Kronecker delta. Every sufficiently smooth periodic real-valued function f(r) can be expanded in a uniformly convergent trigonometric Fourier series, which can be termwise differentiated (i.e., the Fourier series for the derivatives are also uniformly convergent): (10) where m = (m 1 , m 2 , m 3 ) is an integer vector, Z is the set of integers, and the coefficients of series (10) are given by the equalities (11) The following differential operators (gradient and Laplacian) as applied to a sufficiently smooth function f(r) can also be expressed in the form of uniformly convergent Fourier series: The inverse of the Laplacian on the set of functions with a zero mean is given by the formula (13) Relations (10)-(13) underlie high-accuracy computations of the gradient and the Laplacian of the function on a uniform grid with the help of the discrete Fourier transform when the high-frequency harmonics of the function can be neglected.
To compute the Fourier coefficients in (12), we use the rectangle rule on a uniform grid in the parallelepiped Π: where (15) setting (16) where N denotes a finite set of integer vectors m and n with nonnegative integer components described by relations (15). Note that the rectangle rule on a straight line interval for trigonometric polynomials is closely related to the Gaussian quadrature rules for Chebyshev polynomials [9], namely, the rectangle rule on a uniform grid with N nodes is exact for trigonometric polynomials of degree n, where n < N/2 thus, it plays the role of a high-accuracy Gaussian formula on the set of periodic functions.
Equality (16) is the discrete Fourier transform of the grid function f(r n ) up to a multiplicative constant.
The grid function f(r n ) is recovered from its coefficients (16) with the help of the inverse discrete Fourier transform: Relations (16) and (17) allow us to use the FFT algorithm for high-accuracy low-cost computation of the Fourier coefficients and for the recovery of the original grid function from the computed coefficients. The Fourier coefficients of the discrete analogues of the gradient f m , the Laplacian f m (see (12)), and the , . Note that the discrete analogues (18)-(21) of the differential operators have advantages over their widespread difference approximations [2]. They coincide with the original operators (12) and (13) at grid points (14), (15) on the set of finite linear combinations of trigonometric functions of form (10) if the index components of the harmonics satisfy the inequalities 2m k < N k (k = 1, 2, 3). Moreover, these linear combinations can be used to calculate the functions and the corresponding operators (12) and (13) not only at grid nodes, but also at arbitrary points of the parallelepiped Π.
Remark. The expansion coefficients in the trigonometric Fourier series of the function lnε(r) in Eq. (7) must tend to zero rather quickly. To simplify the understanding of the conditions imposed on the spectrum of the function lnε(r), we consider the onedimensional equation (7). If the step size of a uniform spatial grid is equal to h, then the spatial Nyquist frequency is (2h) -1 . Therefore, the step size h is chosen so that the partial trigonometric sum of the Fourier series of the vector function ∇lnε(r), with harmonics involving wave numbers q satisfying the inequality |q| ≤ πh -1 approximates the original series and its derivative with sufficient accuracy determined by the particular application. Note that the right-hand side of Eq. (7) involves the scalar product of the gradients of lnε(r) and the sought additional potential ϕ(r). If each of these functions contains harmonics with the wave number modulus πh -1 , then their scalar product can contain harmonics with the wave number modulus 2πh -1 , i.e., can include spatial frequencies exceeding the Nyquist frequency. Therefore, the constraint on the mesh size must be at least doubled: the significant harmonics of lnε(r) (i.e., ones that carry significant information on the function) must include only wave vectors with |q| ≤ π(2h) -1 . In the case of the multidimensional equation (7), by q we mean each of the independent components of the wave vector q m .
The procedure for computing the right-hand side of Eq. (7) can be briefly described as follows.
Step 1. The function f(r n ) is sequentially set equal to ϕ(r n ) and lnε(r n ), and each time the discrete Fourier transform (16) is applied to it for computing the Fourier coefficients f m . Here, the vector coefficients f m are found according to the left formula in (18) and the discrete analogue of the gradient f(r n ) is computed using to formula (17), in which the scalar coefficients f m are replaced by the vector coefficients f m .
Step 2. Compute the scalar product ρ(r n ) = (E 0ϕ(r n ), lnε(r n )), where the vectors ϕ(r n ) and lnε(r n ) were found at the preceding step, and set the function f(r n ) equal to ρ(r n ), applying to it the discrete Fourier transform (16) for computing the Fourier coefficients f m . Here, the quantities f m are found according to formula (19) and the discrete analogue of the Laplacian inverse f(r n ) = ρ(r n ) is calculated using formula (17), in which the coefficients f m are replaced by the coefficients f m .
The computational complexity of the direct (16) and inverse (17) discrete Fourier transforms based on FFT is estimated as O(MlogM) operations, where M = N 1 N 2 N 3 . Due to this estimate and the possibility of efficient parallelization of the FFT algorithms [10], three-dimensional electrostatic problems in periodic dielectric media can be solved in real time by iterative methods with a large number of grid nodes in the parallelepiped Π.
4. The algorithm proposed above was used to compute the electric field in one-and three-dimensional spatially periodic dielectric model media consisting of components A and B characterized by permittivity ε A and ε B , respectively.
First, we consider a one-dimensional medium with lamellar morphology. The volume fractions of the components A and B are specified as (22) respectively. The local permittivity of the medium is defined as The uniform electric field E 0 = 1 is assumed to be aligned with the x axis. In this case, Eq. (3) for the additional potential ϕ(x) becomes (24) with the periodic boundary condition ϕ(x + 1) = ϕ(x).
To apply the numerical algorithm described in Steps 2 and 3, we introduce a uniform grid x n = n/N (n = 0, 1, …, N -1) on the interval [0, 1], where N = 2 k (k = 1, 2, …, 10). The trigonometric Fourier series expansion of the function (lnε(x)) x = ε x (x)/ε(x) (in terms of the functions exp(i2πnx)) has coefficients that converge rapidly to zero, so that, for |n| ≤ 20, the partial sum of this series approximates the function (lnε(x)) x up to 12-13 significant digits. A further increase in the number of series terms in the case of the standard mantissa length 16 for real numbers is not reasonable, since it does not improve the numerical accuracy, because harmonics with high indices contain only noise with characteristic amplitude values of ~10 -15 . The summation of such series is a well-known illposed problem [11]. An analysis of the harmonics of the numerical solution ϕ(x) and its first derivative ϕ x (x) shows that their maximum achieved numerical accuracy is 12-14 significant digits. Therefore, to compute the Fourier coefficients of these functions, it is sufficient to use grids with the number of nodes equal to N = 32, 64 and to use partial sums of their trigonometric series including harmonics with |n| ≤ 20. Grids with large numbers N of nodes do not improve the numerical accuracy. As an illustration, Fig. 1 presents the electric field strength (Fig. 1a) (25) and the energy density of the electric field (Fig. 1b) Now we compute the electric field in a threedimensional spatially periodic dielectric medium with diamond morphology [12]. The volume fractions of the components A and B used in the computations have the form (27) and the permittivity of the medium is specified as The numerical algorithm proposed above is used to compute the additional potential ϕ(r) satisfying Eq. (3) and the periodic conditions The numerical solution was determined on a uniform grid in a unit cube with 32 3 or 64 3 nodes. The accuracy of the computed electric field strength at any point of the cube was 6 and 10 significant digits, respectively. For numerical convergence of the iterations with relative accuracy of 10 -13 , we used 9 iterations of the Ng method [7] or 12 iterations of the successive approximation method (8) with τ = 1. Figure 2 presents level surface of the electric field energy density π + − + π − + + π− + + − π + + −∞ < < +∞ η = −η ( ) 0.5 0.044[cos(2 ( )) cos(2 ( )) cos(2 ( )) cos(2 ( ))], , , To conclude, we note once again that, for a fixed number of grid nodes, the accuracy of the numerical solution depends substantially on the rate at which the Fourier coefficients of the function lnε(r) decay with increasing indices of harmonics. The solution is computed more accurately when the rate of decay of the Fourier coefficients is higher.
The above-proposed method for solving Eq. (1) can also be used in the case of other boundary conditions that effectively reduce to periodic ones, such as homogeneous Dirichlet boundary conditions and reflection conditions.