{'fraction_non_alphanumeric': 0.05007393198257603, 'fraction_numerical': 0.020021580146265434, 'mean_word_length': 4.073803730738037, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strict derivation of angular-averaged Ewald potential

In this work we strictly derive an angular-averaged Ewald potential suitable for numerical simulations of disordered Coulomb systems. The potential was first introduced by E. Yakub and C. Ronchi without strict mathematical justification. Two methods are used to find the coefficients of the series expansion of the potential: based on the Euler-Maclaurin and Poisson summation formulas. The expressions for each coefficient is represented as a finite series containing derivatives of Jacobi theta functions. We also demonstrate the formal equivalence of the Poisson and Euler-Maclaurin summation formulas in the three-dimensional case. The effectiveness of the angular-averaged Ewald potential is shown by the example of calculating the Madelung constant for a number of crystal lattices.

In this work we strictly derive an angular-averaged Ewald potential suitable for numerical simulations of disordered Coulomb systems. The potential was first introduced by E. Yakub and C. Ronchi without strict mathematical justification. Two methods are used to find the coefficients of the series expansion of the potential: based on the Euler-Maclaurin and Poisson summation formulas. The expressions for each coefficient is represented as a finite series containing derivatives of Jacobi theta functions. We also demonstrate the formal equivalence of the Poisson and Euler-Maclaurin summation formulas in the three-dimensional case. The effectiveness of the angular-averaged Ewald potential is shown by the example of calculating the Madelung constant for a number of crystal lattices.

I. INTRODUCTION
The Coulomb potential plays a fundamental role in numerous theoretical and applied problems [1]. Any system of charged particles is characterized by long-range electrostatic interaction which originates the main problem in the mathematical description of such systems. Electrostatic energy of an infinite electroneutral system of charged particles is a conditionally convergent series, and its sum depends on the summation order [2]. The solution to this problem for systems with a translational symmetry was proposed by Ewald [3]. By addition and subtraction of normally-distributed screening charges the original sum may be transformed into two rapidly converging sums. The correctness of the Ewald's summation technique is justified experimentally. On the other hand, for disordered systems of particles (ionic liquids, gases, plasma) the Ewald's technique is artificial and reduntant. Nevertheless, in atomistic simulations of plasma the Ewald's technique is widely used despite the fact that the computational effort scales as O(N 3/2 ) [4]. There are also other approaches concerning the problem of long-range potentials, including fast multipole [5] and particle-mesh-based methods [6]; however, they become efficient only for sufficiently large systems (at least 10 5 particles).
In 2003 E. Yakub and C. Ronchi [7] proposed an angular-averaging technique for the Ewald potential and showed numerically the consistency of the new potential. Their result was presented as a power series depending on the distance between particles and it was stated that all the coefficients of the series except for the first two were equal to 0. The potential was then used in many computational works due to its obvious efficiency [8][9][10][11]. However, no strict mathematical derivation of the angular-averaged Ewald potential was published in the literature. Therefore, the question remained whether the new potential is approximate?
In this work we discuss a strong derivation of the angular-averaged Ewald potential from the original Ewald potential. We investigate the coefficients of the power series for the potential including their dependence on the smearing parameter. We strictly show that all the coefficients except for the first two tend to 0 in the case of point charges. Two methods are used to find the coefficients of the series expansion of the potential: based on the Euler-Maclaurin and Poisson summation formulas. The expressions for each coefficient is represented as a finite series containing derivatives of Jacobi theta functions. We also demonstrate the formal equivalence of the Poisson and Euler-Maclaurin summation formulas in the three-dimensional case. The physical meaning of the potential is discussed including the fulfillment of the elecroneutrality condition. Finally, we demonstrate the convergence of the Madelung constant for a number of crystal lattices using the direct summation with the angular-averaged Ewald potential for up to 2 × 10 7 particles.

II. SUMMATION PROBLEM
Consider a cubic cell of a volume L 3 which contains N point particles. Each i-th particle has a charge Q i and position r i in the cube. The system is electroneutral: Periodic boundary conditions are assumed, so the cell repeats itself in three mutually perpendicular directions. It means that a particle with a position r i in the cell has infinite number of images with positions r i + nL. Here, n is an integer vector n = (n x , n y , n z ), n x , n y , n z ∈ Z. According to Coulomb's law, the total potential energy E of such an infinite system is (we use Gaussian units): The summation is performed over all integer vectors n; the prime means that the terms with n = 0 are omitted if i = j. So a particle i interacts with all its replica images, but not with itself. This sum is conditionally convergent; thus, to obtain the correct answer one has to use a special Ewald summation technique [3]. The main idea is to add and subtract a normally distributed screening charge with a standard deviation 2L/δ [2, p. 294]; δ > 0 is a dimensionless parameter. Below we assume the dependence of all values on δ; point charges correspond to the case δ → ∞. Ewald summation procedure results in [12, p. 346]: where erfc(x) is the complementary error function. The summation n =0 means that the term n = (0, 0, 0) is omitted for all i, j. Here, n = |n|, r ij = r i − r j , r ij = |r ij |. If δ 1, the terms of the order erfc(δn)/n can be omitted [13, p. 34]. Thus, formula (3) simplifies to: where We will call Eqs. (7)-(8) the Ewald potential and Eq. (6) the Ewald formula. The unary potential φ 1 does not depend on r ij ; the pair potential φ 2 (r ij ) defines the electrostatic interaction and is angular dependent. In disordered media, such as electrolyte or plasma, this dependence is confusing and results in additional complications. Thus, our goal is to somehow make the Ewald potential spherically symmetric.
To do it, we use the approach of E. Yakub and C. Ronchi [7]. Following them, we average Eq. (8) over all directions of n at a distance r ij , since all spatial orientations are equivalent: The only factor to average is the cosine (n · r ij = nr ij cos θ): 2π 0 cos (2πnr ij cos θ/L) dψ = L/r ij 2πn sin (2πnr ij /L) .
Thus, we get an averaged pair potential φ a 2 (r ij ): We expand it in the converging series of r ij , expanding erfc(δr ij /L) and sin(2πnr ij /L) into the Taylor series: with the coefficients Here, we introduced the notation One can easily find C 0 : which is two times larger than the unary potential φ 1 . To compute C k for k ≥ 1, we need to sum an infinite series over n. Further computations are made for k ≥ 1.
We are going to include the zero term n = 0 in the sum of Eq. (13). Since f k (0) = δ 1,k : Here, δ ij is the Kronecker delta. Now Eq. (13) transforms into the following expression: for k ≥ 1. Thus, we need to exactly calculate the following series over all integer vectors n: Below we compute such a series for any δ > 0. The most interesting case of δ → ∞ will be considered separately (see Sec. III C 2). Also, we formally investigate the case δ → 0 (see Sec. III C 3). We provide numerical computations of series (18) at different values of δ for 1 ≤ k ≤ 4 in App. A. The summation of (18) is the key result of our work: it proves that E. Yakub's and C. Ronchi's results [7, Eqs. (6)-(8)] published without a strict proof are correct in the limit δ → ∞.

III. SUMMATION
We consider two ways to calculate series (18) using the Euler-Maclaurin and Poisson summation formulas. Both of these formulas transform the sought series into an integral form with some residual terms. Next, we formulate these formulas in case of three dimensions as theorems.

A. Formulation of summation formulas
There are more general formulations and relations for the Euler-Maclaurin theorem. They can be found in [14, Sec. 3, Eq. (5)], [15,Eq. (9)], [16]. We use a more practical and simple form of the formula.
Theorem III.1 (Euler-Maclaurin). Let D ∈ R 3 be a regular region with continuously differentiable boundary surface ∂D. Let f (r) : R 3 → R be a twice continuously differentiable function inD = D ∪ ∂D and let h be the unit outward normal to ∂D. Then where Here, q = (q x , q y , q z ) is an integer vector (q x , q y , q z ∈ Z), q 2 = q 2 x + q 2 y + q 2 z , qr = q x r x + q y r y + q z r z . The summation q =0 means that the term q = (0, 0, 0) is omitted. The Poisson formula [17,Theorem 6.11] imposes stronger conditions on the function f (r). where is a Fourier transform of f (r).
The definition of a Schwartz function can be found in [17, Definition 5.1.]. Both (19) and (21) should give the same result in practice. Next, we show how the Euler-Maclaurin formula (19) can be formally reduced to the Poisson formula (21).

B. Transformation of the Euler-Maclaurin to the Poisson summation formula
Let us begin by proving the following statement.
Theorem III.3. Let D ∈ R 3 be a regular region with continuously differentiable boundary surface ∂D. Let f (r) : Proof. First, we use the following rule for integration by parts [18, Theorem 37.2]: where h is the unit outward normal to ∂D; u(r) is a scalar-valued function and w(r) is a vector-valued function.
Since f (r) and D satisfies the conditions of Theorem III.1, Eq. (19) can be used. Next, we integrate the last term in (19) by parts (u(r) = G(r), w(r) = ∇f (r)): Second, the first Green's identity [19, Chapter 7, Eq. (G1)] will be used: where v(r) is a scalar-valued function. Using (26), we obtain: Substituting Eqs. (25), (27) into the Euler-Maclaurin summation formula (19), we get: It is easy to find ∆G(r): Then, the last term in (28) is simplified: Now we can include the q = 0 term into sum (30): which leads us to the desired relationship (23).
Corollary III.3.1. The residual term in the Euler-Maclaurin formula can be written as follows: This form of the residual term in the Euler-Maclaurin formula is a more appropriate for practical calculations due to the absence of a surface integral. If one considers a sphere of radius R as a region D and take the limit R → ∞, formally the Poisson formula (21) is obtained (see App. B). However, the class of functions f (r), for which such a limit is valid, remains uncertain.
It is easier to sum series (18) using the Poisson formula (21), that we will perform next. Nevertheless, in App. B we provide a summation of (18) using the Euler-Maclaurin formula.

C. Summation using the Poisson formula
To calculate (18), we use the Poisson summation formula (21): where is a Fourier transform of f k (n). The summation is now performed over an integer vector q = (q x , q y , q z ), q x , q y , q z ∈ Z. M (a, b, x) is the confluent hypergeometric function defined by the series: where a (s) denotes the rising factorial:

General formula
By definition: Since (1 − k) (k) = 0, (1 − k) (k+1) = 0 and so on, the series is truncated: Substituting (38) to (34), we get series (18) in the following form: To perform the summation over q, we use the multinomial theorem: The summation over α 1 , α 2 , α 3 ≥ 0 is performed only if α 1 + α 2 + α 3 = s. In this way, we separated the variables so that all sums became one-dimensional. Each internal sum is related with a Jacobi theta function with zero argument: where ϑ 3 (0, x) is defined by: Thus, we get the final formula substituting (40) to (39): So the summation over unrestricted three-dimensional argument is transformed into a finite sum. We tested formula (43) for 1 ≤ k ≤ 4 using Wolfram Mathematica [20]: numerical summation of n f k (n) and symbolic calculation of the right part of (43) for 0.8 ≤ δ ≤ 3 gives the same results with a machine accuracy. Here, the limitation δ ≥ 0.8 appears because Mathematica fails to compute the final numerical result for small δ. Therefore, we hope, that formula (43) will be useful for numerical calculations of Jacobi theta function ϑ 3 (0, x) derivatives (see App. A).
The most interesting and practical result is produced in the limit δ → ∞. It corresponds to an infinitely small width of the normally distributed charge.
2. The limit of infinitely small width 2L/δ In the limit δ → ∞, theta function becomes constant: Thus, only zero-order derivatives in (43) gives a non-zero result; therefore the only term at α 1 = α 2 = α 3 = 0 = s contributes to (43). The final result for all k ≥ 1 is: This asymptotic behavior is valid even for relatively small values (see App. A, Fig. 5). Using Γ(k + 1/2) = (2k)! √ π 2 2k k! , we obtain all the coefficients for k ≥ 1: This key result was presented in [7] without any proof. Now the averaged pair potential takes a simple form: The full potential energy then reads: The first constant term in Eq. (48) is eliminated by C 0 in the second term due to the electroneutral condition (1): Total energy results in where r m = ( 3 4π ) 1/3 L is the radius of the sphere 4π 3 r 3 m = L 3 with equivalent volume L 3 . We will call Eq. (51) the averaged potential.
Eq. (51) can be used in practical calculations of two-component plasma energy and is independent of any external parameters. In Sec. IV, we describe the main properties of the averaged potential.
3. The limit δ → 0 In the limit δ → 0, theta function shows hyperbolic behavior: First, we calculate the derivative of theta function in the limit δ → 0: The sought (43) then become: We rewrite the sum over α 1 , α 2 , α 3 as follows and introduce the notation: One can easily find that S(1) = 1. We could not derive an explicit expression for S(k) for any k ≥ 2; nevertheless exact symbolic calculations using Mathematica [20] give S(k) = 0 for any 2 ≤ k ≤ 60. We suppose, that that results in n f k (n) = δ 1,k , δ → 0.
One can see the dependence of series (18) on δ for δ 1 in Fig. 5.

IV. ANALYSIS OF AVERAGED POTENTIAL
Remember that our system has periodic boundary conditions. Since that the interaction potential must have the following properties along some chosen direction: • be periodic; • has a minimum at some point; • be symmetrical relative to the point of minimum.
To illustrate this, let us consider a simple one-dimensional periodic system (see Fig. 1) with only one ion in the main cell. The solid line shows the potential energy interaction between some trial unit charge at position r and the system of periodically repeated ions. First, if the position of a trial charge is the same as one of the ions, the energy should be infinite. Second, if the position of a trial charge is equidistant from the two ions (at the cell edge), the energy should take a minimum value. These considerations show why a potential must have the properties described above.
The averaged potential (51) reaches a minimum value 3/(2r m ) > 0 at a point r = r m , but for r > r m , it behaves incorrectly. Therefore, we must consider expression (51) only up to a distance r m ; for r > r m we redefine ϕ(r) by zero: ϕ(r) = 0.
Then to calculate the potential at the point r i , one has to use the following algorithm: 1. Move to the reference point r i of the selected ion i;

Calculate the energy of its interaction with each
Thus, we proceed to consider a spherical cell, which we superimpose on a periodic cubic cell (see Fig. 2). Now the selected ion is affected not only by particles in the main cell but also particles in the sphere. The total number of particles in the sphere N s depends on the sphere center position. Now examine what charge density is created by each particle in the sphere. Consider just one particle in a sphere at the position r 1 . It creates, at some point r, a potential: We calculate the charge density ρ(r) at a point r: using the Poisson equation. The Laplacian of the averaged potential has the following form: Then the charge density: We see that this point particle is not a Coulomb one in the usual sense. In addition to the point density Q 1 δ(r − r 1 ), it creates a uniformly distributed charge of the opposite sign in the entire sphere; its magnitude is 3Q1 4πr 3 m . This particle can be treated as an ordinary Coulomb point particle + some additional charge around it. The interaction of this additional charge with some other particle is determined by an additional cubic term 1 2 (r/r m ) 3 in the averaged potential (51). Moreover, the charge density (62) is such that the entire sphere is electrically neutral: Thus, the averaged potential (51) describes the interaction of spheres with radius r m and zero charge; the spheres interact with each other only if the distance between their centers is less than r m . Finally, we compare the averaged potential with the Ewald potential along the three primary crystal directions (see Fig. 3) and with the pure Coulomb potential. We picture the Ewald potentials up to a minimal value since further behavior is trivial (see Fig. 1 and the reasoning above). The averaged potential is plotted for r ≤ r m . All of these potentials tend to the Coulomb one at small distances, in particular, the averaged one: But for large distances they behave differently. The curve of the averaged potential is situated between the curves of the Ewald potentials along [100] and [111] directions. Another difference is the position of the minimum. The smallest (r/L = 1/2) and the largest ( √ 3/2) positions are for the Ewald potential along the directions [100] and [111], respectively. The value r m /L ≈ 0.62 for ϕ(r) is between them. The fact that r m /L > 1/2, which makes difficulties during the numerical calculations, will be discussed further (see Sec. V A). Now we will show some applications and advantages of the obtained potential in practice calculations.

A. Computation algorithm
Since the averaged potential is truncated at r = r m , there is a discontinuity in the energy at r ≥ r m . It can lead to problems during numerical calculations and simulations [21, p. 302]. To avoid this, we shift the averaged potential to make it zero at r ≥ r m : whereφ The potentialφ(r) will be used in numerical calculations. We rewrite the last term in (65) in a more simple form due to the electroneutrality condition (1): Thus, we have the following formula for energy [7, Eqs. (7), (8)]: The range of r ij in the second sum of Eq. (68) is r m > L/2. It sets up another problem: ion i is affected not only by particles in the main cell but also by their images [7, Sec. III]. Formally, one can apply the same technique, as in traditional atomistic simulations, using the cut-off radius of a potential. So, every particle in the main cubic cell is surrounded by a sphere with radius r m ; all interactions of the central particle with other particles and images inside the sphere are summed (see Fig. 2). The details can be found in [10, see Sec. 3 and Fig. 2].

B. Madelung constant
In this section, we apply the averaged potential to compute the energy of a two-component system of charges. For simplicity, we consider several ordered systems of stationary charges, as serious problems arise in simulations of pure Coulomb dense systems of moving charges. We show that even for this case our calculations give accurate results. For ordered systems, the spherical cell is not electroneutral (in the sense of Eq. (1)). However, if the number of charges increase the total charge i Q i /N tends to zero in the spherical cell. Thus, the calculation accuracy improves with the number of charges.
We will examine several ordered systems and calculate their Madelung constants using the averaged potential [8]. For this case, the formula for the Madelung constant M is: where r 0 is the nearest neighbor distance and Z i = Q i /e is the charge of a particle with respect to the elementary charge. M is independent of i since all ions of the same sort are in equivalent positions. The results are presented in Secs. V B 1, V B 2, V B 3. All the lengths are given in the units of a, the length of a unit cubic cell. In all cases, increasing the number of ions N leads to a more accurate value of the Madelung constant. We see that the decrease in relative charge Z/N does not necessarily increase the accuracy of M , here Z = Ns i=1 Z i and N s is the total number of particles in the sphere around a chosen ion. Moreover, the absolute value of Z can increase with N . Nevertheless, the convergence for M is obviously observed.

NaCl
The unit cell of NaCl consists of 8 ions. Positions of Na + and Cl − ions are shown in Tab. I. The dependence of the Madelung constant (69) on the number of ions is given in Tab. II. Table I: Coordinates of ions in a unit cubic cell with L/a = 1 of rock salt, r0/a = 1/2.

Coordinates
Na

C. Scaling and advantages
It is of practical importance to compare the performance of computations between the exact Ewald formula and averaged potential. The Madelung constant takes the following form with the exact Ewald formula (3) for energy: For comparison, we make calculations using (69) and (70) for NaCl at different N . The following parameters for the Ewald potential are chosen: n x , n y , n z = −20, . . . , 20 and δ = 10. We perform sequential calculations on a CPU Intel Core i7-7700HQ 2.8 GHz. Calculation time for different number of ions is present in Fig. 4. The results were linearly fitted. Although both dependencies are linear, the calculation by the formula (69) is approximately 2600 times faster than by (70). In [10, Fig. 3] it was also shown that the averaged potential gives faster results compared to the other methods used in paper [10].
The difference between the exact result and the one obtained using the averaged potential (69) is about 0.08% for equivalent computation time (Fig. 4, horizontal black line). We hope, that for a disordered system this error will be even smaller.
The averaged potential has a significant advantage: it doesn't depend on any external parameters that affect the convergence. On the other hand, the influence of parameters on the results with the Ewald potential is significant [13,.
Thus, the averaged potential seems very perspective for the numerical simulations of Coulomb system. In addition, the averaged potential (51) has a fairly simple analytical form compared to the Ewald potential (8), that makes it attractive for analytical studies and obtaining exact results.

VI. CONCLUSION
A strict derivation of the angular-averaged Ewald potential is proposed using the Euler-Maclaurin and Poisson formulas. Thus, mathematically exact character of the potential is established. Additionally, the formal equivalence of the Euler-Maclaurin and Poisson formulas is proved. From a physical point of view, the averaged potential describes the interaction of two spheres with radius r m = (3/(4π) 1/3 )L, where L is the size of a cubic computational cell. Each sphere contains a point charge in the center and a compensating uniformly distributed in the entire sphere charge of the same value and opposite sign. The spheres interact with each other only if the distance between their centers is less than r m . Thus, the long-range Coulomb interaction in an disordered point system of charges is equivalent to the interaction of electrically neutral spheres with a finite-range potential; the range of the potential depends on the size of a cubic computational cell. Technically, the calculation of the interaction energy is straightforward: every particle in the main cubic cell is surrounded by a sphere with radius r m ; all interactions of the central particle with other particles and images inside the sphere are summed. Our computations of the Madelung constant for a number of crystal lattices show the efficiency of the angular-averaged potential for systems containing up to 2 × 10 7 particles. n f 4 (n) = − 3δ 9 e −3δ 2 8π 15/2 16δ 6 (∂ϑ 3 (0, x) /∂x) 3 x=e −δ 2 + 24 2δ 2 − 7 e δ 2 δ 4 ϑ 3 0, e −δ 2 (∂ϑ 3 (0, x) /∂x) 2 x=e −δ 2 +2δ 2 ϑ 3 0, e −δ 2 (∂ϑ 3 (0, x) /∂x)| x=e −δ 2 24δ 4 ∂ 2 ϑ 3 (0, x) /∂x 2 x=e −δ 2 + 4δ 4 − 42δ 2 + 105 e 2δ 2 ϑ 3 0, e −δ 2 +ϑ 2 3 0, e −δ 2 8δ 6 ∂ 3 ϑ 3 (0, x) /∂x 3 The dependence of series (43) on δ for 1 ≤ k ≤ 4 are shown in Fig. 5 by the solid lines; the dashed lines represent asymptotic behavior at δ → ∞ (45). The difference between the exact and asymptotic results is small even at δ ≥ 2.5; at small δ the series tends to 0 for k ≥ 2 and to 1 for k = 1 as it was predicted (58). We hope, that Eqs. (43), (A1)-(A4) will be usefull to numerically calculate derivatives of Jacobi theta function ϑ 3 (0, x).