text:
# Electrostatics in Periodic Slab Geometries II

## I  Introduction

Long-range forces, namely gravitation and electromagnetic interactions, are difficult to treat exactly in many-particle systems. In condensed matter, biological and solid state physics, there is considerable interest in the simulation of charged and polar molecules embedded in a slab in which the particles are replicated periodically in two directions. The applications range widely to include soap bubbles, cell membranes and electrochemistry. Surprisingly Ewald and other rapidly convergent methods for computing long-range interactions in infinitely periodic systems are computationally more demanding when the system is periodic in only one or two of the three spatial dimensions because of the breaking of spatial symmetry.

Several “two-dimensional Ewald” (EW2D) methods have already been proposedwidmann97a . The most successful is that introduced by Parryparry75a ; parry76a  and othersheyes77a ; deleeuw82a . This method can be used to obtain accurate results but is limited to small systems (i.e. $`10^2`$ to $`10^3`$ charges). This springs from the fact that the pair separations, $`𝐫_{ij}`$, are no longer separable in the Fourier summation because of the necessary decoupling of the aperiodic $`z`$ component from the periodic $`x`$ and $`y`$ components. Therefore it appears that the complexity (scaling of time with number of charges) of this method can be no better than quadratic. Also the one dimensional Ewald methodporto00a  suffers from the same problem.

A promising alternative to Ewald summation is possible using a convergence-factor techniquesperb01a . The basic Coulomb pair interaction is multiplied by one with the factor $`lim_\beta \mathrm{}e^{\beta r_{ijn}}`$. It is crucial to prove that the limit can be taken outside of the summation and now the modified potential is more readily susceptible to analysis. In a recent study,arnold01b  Arnold and Holm derive a two-dimensional version of this method with complexity $`N^{5/3}`$ along with accurate error formulas. This method is particularly attractive when extremely high accuracy is desired. We will use this method, abbreviated as MMM2D and the EW2D algorithm mentioned above as reference methods.

One idea, used several times to study water interfaces, introduces a spatial constraint within the simulation cell. If the box has dimensions $`L_x\times L_y\times L_z`$, the particles may have $`z`$-coordinates on $`[0,h)`$, while the space within $`[h,L_z)`$ (the gap size) remains empty empty-space . The primary cell is, as usual, replicated in all three directions periodically. The effect of this is to create a “primary layer” and an empty layer followed by an infinite array of intercalated image and empty layers. Figure 1 summarizes the essential geometrical considerations in this problem.

While this idea seems trivially to be correct, it has one flaw. The flaw can be realized by attempting to reproduce the results of an EW2D calculation using a three dimensional Ewald summation (EW3D) as just described. One will find that the results always differ by a dipole-moment dependent constant. This fact, although at first thought to be merely a slow convergence issue, was noticed by Spohr spohr97a . The EW3D formula contains a shape dependent dipole-term whose origin was mathematically proven in the 1980’s deleeuw80a ; deleeuw80b , that reflects the naturally chosen spherical order of summation in the conditionally convergent Coulomb sum. For the case of layers, however, it is necessary to use a slab-wise summation order. This was realized by Yeh and Berkowitzyeh99a , who applied a theory due to E. R. Smith smith81a  for infinite crystals of various shapes in order to obtain the correct dipole term for a slab-wise summation order. Another complication arises, if the periodic supersystem is surrounded by some medium with a different dielectric constant. This medium has a polarizing effect on the particles in the simulation cell no matter how large the supersystem is. For the remainder of this paper we consider only the simple case of zero contrast for which the dielectric constant is the same inside and outside the supersytem. This is often referred to as the vacuum boundary condition.

One of the advantages of using this method from a practical point of view is that any standard Ewald program may be used with only minor modification. It is easy to show numerically that the use of this new dipole term and sufficient spacing ($`L_zh`$), yields the same result as two-dimensional “brute force”, MMM2D and EW2D summation methods. As explained in a previous paper, referred to here as “paper I”arnold01d , errors introduced due to the image layer effects decay exponentially with the size of the gap. This method is sometimes referred to as a corrected three-dimensional Ewald sum (EW3DC)yeh99a . Like the original Ewald method it has a complexity of $`N^{3/2}`$ and the equations can be discretized onto a mesh yielding an $`N(\mathrm{log}N)^{3/2}`$ complexity. In paper I arnold01d , we derived an analytic electrostatic layer correction term, called ELC-term, that subtracts the interactions of the unwanted slab replicas. This term can be evaluated linearly in $`N`$, and has full error control. The use of the ELC term and the change of summation order, shortly called the ELC method, enables us to adapt any 3D summation method, such as the non-Ewald convergence-factorsperb01a  and multipole expansiongreengard87a  methods, to slab geometries. The latter of these, though in possession of a slightly better scaling, has such a considerable amount of overhead that it is useful only for much larger systems than normally usedesselink94a . The abbreviations and complexities of some of the two-dimensional methods available at present are summarized in Table 1, where the ending LC always denotes the use of the ELC method. In this terminology the EW3D algorithm plus the ELC method is called EW3DLC and so on.

This paper is organized as follows. In Sec. II the equations required for implementation of EW3DLC are collected and discussed. In Sec. III the standard Ewald error formulas are extended for anisotropic systems and the error formula for the new layer correction term is described. The methodology is demonstrated on several test-systems in Sec. IV and the efficiency compared with other 2D methods is shown in Sec. V. Finally the results of this paper and the outlook for computational studies of slab geometries are summarized in Sec. VI. In the Appendix we give explicit formulas for the ELC-term contributions to the force and the diagonal part of the stress tensor.

## II EW3DLC in a Nutshell

The electrostatic energy of one section of a periodic slab follows from the classical formula for the Coulomb pair potential. Omitting the prefactor ($`1/4\pi ϵ_0`$),

$$E=\frac{1}{2}\underset{𝐧_{}}{\overset{}{}}\underset{i,j=1}{\overset{N}{}}\frac{q_iq_j}{|𝐫_{ij}+𝐧_{}|}.$$
(1)
This serves as a useful point of introduction to our notation: $`N`$ particles with charges $`q_i`$ and positions $`𝐫_i`$, reside in a simulation box of edges $`L_x\times L_y\times L_z`$. The image boxes are denoted using the vector $`𝐧_{}(L_x,L_y,0)`$ and $`𝐫_{ij}=𝐫_i𝐫_j`$. The prime on the inner summation indicates the omission of the primary box, $`𝐧_{}=(0,0,0)`$, when $`i=j`$ (the singular case). Strictly speaking the sum over $`𝐧_{}`$ is ill-defined without specification of the order of summation. This is a consequence of the fact that summations of the type $`|𝐧+c|^\xi `$ are conditionally convergent for $`0<\xi <d`$, where $`d`$ is the dimensional periodicity. The summation limit of relevance is an expanding cylinder or disc such that the summation terms appear in order of increasing $`|𝐧_{}|`$.

We convert this into a three-dimensional problem by nesting the $`L_x\times L_y\times h`$ particle volume within a larger $`L_x\times L_y\times L_z`$ box, thereby creating a screen for the ghost charges (the new image charges) as described in the introduction. The total electrostatic energy $`E`$ is, via combination of the celebrated Ewald identity and the recent convergence-parameter layer correction, rewritten in the computationally convenient form:

$$E=E^{(r)}+E^{(k)}+E^{(s)}+E^{(d)}+E^{(lc)},$$
(2)
where the individual terms are as follows,

$`E^{(r)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}{\displaystyle \underset{𝐧}{\overset{}{}}}q_iq_j{\displaystyle \frac{\text{erfc}(\alpha |𝐫_{ij}+𝐧|)}{|𝐫_{ij}+𝐧|}}`$ (3)
$`E^{(k)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{L_xL_yL_z}}{\displaystyle \underset{𝐤0}{}}{\displaystyle \frac{4\pi }{k^2}}e^{k^2/4\alpha ^2}\left|{\displaystyle \underset{j=1}{\overset{N}{}}}q_je^{i𝐤𝐫_j}\right|^2`$ (4)
$`E^{(s)}`$ $`=`$ $`{\displaystyle \frac{\alpha }{\sqrt{\pi }}}{\displaystyle \underset{i}{}}q_i^2`$ (5)
$`E^{(d)}`$ $`=`$ $`{\displaystyle \frac{2\pi }{L_xL_yL_z}}\left({\displaystyle \underset{i}{}}q_iz_i\right)^2`$ (6)
$`E^{(lc)}`$ $`=`$ $`{\displaystyle \frac{8\pi }{L_xL_y}}{\displaystyle \underset{k_x,k_y>0}{}}{\displaystyle \frac{1}{k_{}(e^{k_{}L_z}1)}}{\displaystyle \underset{p,q,r=1}{\overset{2}{}}}(1)^p\chi _{pqr}^2`$ (7)
$`{\displaystyle \frac{4\pi }{L_xL_y}}{\displaystyle \underset{k_x>0}{}}{\displaystyle \frac{1}{k_x(e^{k_xL_z}1)}}{\displaystyle \underset{p,q=1}{\overset{2}{}}}(1)^p\chi _{pq0}^2`$
$`{\displaystyle \frac{4\pi }{L_xL_y}}{\displaystyle \underset{k_x>0}{}}{\displaystyle \frac{1}{k_y(e^{k_yL_z}1)}}{\displaystyle \underset{p,q=1}{\overset{2}{}}}(1)^p\chi _{p0q}^2`$
where the summations run over $`𝐧(L_x,L_y,L_z)`$ and $`𝐤2\pi (/L_x,/L_y,/L_z)`$ and are order-independent. As usual, the symbol $`\alpha `$ denotes the Ewald parameter which: (1) determines the radius of the Gaussian charge distributions, (2) has units of inverse length and (3) can take any value without changing the net result for the energy, but has an optimum with respect to convergence of the summations. The first three terms (Eqs. 35) correspond to the usual Ewald real (r), Fourier (k) and self (s) contributions to the energy. The dipole term $`E^{(d)}`$ arises from the cylindrical or slab-wise summation limit as discussed in paper Iarnold01d  and at length in the introduction. It differs from the spherical version by a factor of $`3`$ and by the replacement of the total dipole moment with its $`z`$-component. Again, this dipole term applies only to vacuum boundary conditions. The final term $`E^{(lc)}`$ is the energy subtraction due to all image layers. This term was derived in paper I and is shown here in an expanded $`O(N)`$ form. These summations run over the variables $`k_x,k_y2\pi /L_x,2\pi /L_y`$ and the parallel combination of these is $`k_{}=\sqrt{k_x^2+k_y^2}`$. The component abbreviations are as follows,

$$\chi _{1/2,1/2,1/2}=\underset{i=1}{\overset{N}{}}q_i\mathrm{sinh}/\mathrm{cosh}(k_{}z_i)\mathrm{sin}/\mathrm{cos}(k_xx_i)\mathrm{sin}/\mathrm{cos}(k_yy_i)$$
(8)
and a subscript of zero stipulates omission of the corresponding trigonometric function. For example, $`\chi _{2,0,1}=q_i\mathrm{cosh}(k_{}z_i)\mathrm{sin}(k_yy_i)`$. The corresponding equations for the forces can be obtained by differentiation of Eqs. 37, i.e. $`𝐟_i=_iE`$.

The derivation of the layer correction energy is a crucial part of the method so we shall recapitulate it here in brief. The objective is to obtain a rapidly convergent linear expression for the total energy of interaction of the particles in the primary box with the ghost charges:

$$E^{(lc)}=\frac{1}{2}\underset{n_z0}{}\underset{𝐧_{}}{}\underset{i,j=1}{\overset{N}{}}\frac{q_iq_j}{|𝐫_{ij}+𝐧_{}|}.$$
(9)
It is quite similar in appearance to the original formula (Eq. 1), and can be treated analogously by convergence factors as in the study by Arnold and Holm arnold01b .

We have not here written explicitly the summation order but nevertheless this is an important consideration. sum-order  With the Poisson summation formula applied along each periodic direction a transformation of the in-plane variables, $`n_x`$ and $`n_y`$, into Fourier conjugates $`k_x`$ and $`k_y`$ is achieved. The derivation depends heavily on the requirement that the system is overall neutral and results finally in

$`E^{(lc)}={\displaystyle \underset{i,j=1}{\overset{N}{}}}q_iq_j(`$ $`{\displaystyle \frac{8\pi }{L_xL_y}}{\displaystyle \underset{k_x,k_y>0}{}}{\displaystyle \frac{\mathrm{cosh}(k_{}z_{ij})}{k_{}(e^{k_{}L_z}1)}}\mathrm{cos}(k_xx_{ij})\mathrm{cos}(k_yy_{ij})+`$
$`{\displaystyle \frac{4\pi }{L_xL_y}}{\displaystyle \underset{k_x>0}{}}{\displaystyle \frac{\mathrm{cosh}(k_xz_{ij})\mathrm{cos}(k_xx_{ij})}{k_x(e^{k_xL_z}1)}}+`$
$`{\displaystyle \frac{4\pi }{L_yL_x}}{\displaystyle \underset{k_y>0}{}}{\displaystyle \frac{\mathrm{cosh}(k_yz_{ij})\mathrm{cos}(k_yy_{ij})}{k_y(e^{k_yL_z}1)}}).`$
It is fortunate (in contrast to EW2D) that the coordinates appear as they do because after decomposition of the $`\mathrm{cos}`$ and $`\mathrm{cosh}`$ terms, the $`N^2`$ summation collapses, yielding the order $`N`$ Eq. 7. For reference purposes we have assembled in the appendix the full expressions for the force and the diagonal parts of the stress tensor due to the ELC-term.

## III Error Estimates

The most important measures of error in an electrostatic calculation are the root mean squared error in the particle forces,

$$\mathrm{\Delta }f:=\left(\frac{1}{N}\underset{i=1}{\overset{N}{}}\mathrm{\Delta }𝐟_i^2\right)^{1/2},$$
(11)
where $`\mathrm{\Delta }𝐟_i`$ is the difference between the computed and exact force on particle $`i`$, and the absolute value of the error in the total electrostatic energy, $`\mathrm{\Delta }E`$. We will discuss mainly the former, which is relevant to molecular dynamics simulations.

The total error can be separated into terms originating from the cutoffs used for the Ewald and layer correction,

$$\mathrm{\Delta }f\sqrt{\mathrm{\Delta }f_r^2+\mathrm{\Delta }f_k^2+\mathrm{\Delta }f_l^2},$$
(12)
where the individual components are determined in the same way as in Eq. 11. Analytical error estimates for the real- and reciprocal-space terms for a cubic simulation box have been previously derived.kolafa92a ; petersen95a  Shortly we will extend these results to obtain the formulas for arbitrary box sizes, but we write first the result,

$`\mathrm{\Delta }f_r(\alpha ,r_c)`$ $`=`$ $`{\displaystyle \frac{2}{\sqrt{L_xL_yL_z}}}{\displaystyle \frac{Q^2}{\sqrt{N}}}{\displaystyle \frac{\mathrm{exp}[\alpha ^2r_c^2]}{\sqrt{r_c}}},`$ (13)
$`\mathrm{\Delta }f_k(\alpha ,k_c)`$ $`=`$ $`{\displaystyle \frac{2\alpha }{\sqrt{\pi L_yL_z}}}{\displaystyle \frac{Q^2}{\sqrt{N}}}{\displaystyle \frac{\mathrm{exp}[(\pi k_c/\alpha L_x)^2]}{\sqrt{k_c}}},`$ (14)
where $`Q^2=_{i=1,N}q_i^2`$, $`r_c`$ is the real-space minimum image cutoff (typically but not necessarily $`\frac{1}{2}\mathrm{min}[L_x,L_y,L_z]`$) and $`k_c`$ is an integer defining the truncation of the reciprocal-space summation. These naturally reduce to the standard formulas when the box is cubic. The optimum set of convergence parameters is obtained by setting $`\mathrm{\Delta }f_r^2=\mathrm{\Delta }f_k^2`$. This reduces the number of independent variables to one yielding implicit relations $`\alpha (k_c),r_c(k_c)`$. Then one can easily scan $`k_c`$ to find the optimum (i.e. fastest) $`k_c^{}`$.

If the reciprocal space computation is instead handled by FFT methods such as P3M, inaccuracies arise due to the assignment and de-assignment of the charges onto a mesh. The error for P3M has been estimated (for a cubic system) by

$$\mathrm{\Delta }f_{P3M}(\alpha ,M,P)=\frac{Q^2}{L^2\sqrt{N}}\left(\frac{\alpha L}{M}\right)^P\sqrt{\alpha L\sqrt{2\pi }\underset{k=0}{\overset{P1}{}}a_k^{(P)}\left(\frac{\alpha L}{M}\right)^{2k}}$$
(15)
where $`M`$ is the number of mesh points per edge and $`P`$ is the interpolation order. The expansion coefficients, $`a_k^{(P)}`$ are constants given by Deserno and Holm.deserno98b  A more accurate numerical error estimate is also given by these authors but Eq. 15 is useful since it captures the function form of the error.

The Ewald terms have the advantage of well-behaved error distributions. This is not the case when turning to the layer correction term. Consider for example the probability distribution of $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }f_{lc}`$ over the ensemble of random particle configurations as shown in Fig. 2. The former distribution was generated under conditions of negligible layer correction error while the latter distribution is for layer correction error itself under conditions where it is significant (i.e. small gap). The distribution of the total EW3DLC error fits very well to a Gaussian curve. The error distribution of the layer correction term on the other hand is skewed to the left and fits well to a stretched exponential $`p(x)x\mathrm{exp}(ax)`$. The extended tail of this distribution means that exceptions to the average error will occur with non-negligible frequency. This feature makes it difficult to derive an accurate RMS error estimate for the layer correction. But what is possible and more useful is the derivation of a stricter error criterion, namely an upper bound on the RMS error. A rigorous upper bound was derived in the first paper yielding,

$`\mathrm{\Delta }f_{lc}(l_c,L_z)`$ $``$ $`{\displaystyle \frac{Q^2}{\sqrt{N}}}{\displaystyle \frac{\sqrt{3}}{2(e^{2\pi l_cL_z/L}1)}}(({\displaystyle \frac{2\pi l_c+4}{L}}+{\displaystyle \frac{1}{L_zh}}){\displaystyle \frac{e^{2\pi l_ch/L}}{L_zh}}`$ (16)
$`+({\displaystyle \frac{2\pi l_c+4}{L}}+{\displaystyle \frac{1}{L_z+h}}){\displaystyle \frac{e^{2\pi l_ch/L}}{L_z+h}}),`$
where $`l_c`$ is the integer-valued radial cutoff for truncation of the layer correction forces.

Before leaving this topic a reminder is due of the bias present in slab geometries. Both the Ewald and layer correction terms have an error bias with respect to the $`z`$-position of a particle. The error tends to be higher near the walls of the system as if it were “hotter” in those regions. In fact, one can imagine a simulation in which the overall RMS error in the forces is smaller than the RMS random force but due to this bias, the RMS force-error at the surfaces is larger than the thermal noise, amounting to true hot zones at the surfaces. Therefore it is recommended to exercise some caution in the treatment of these effects. The layer correction error in particular has a strong surface bias, and it is therefore quite appropriate to use an upper bound in place of an RMS error estimate. In other words, the layer correction should be carried out to higher accuracy than the real and reciprocal space parts.

### Anisotropic Real-Space Error

Let us write down the real-space force first,

$$𝐟_i^{(r)}=\underset{j=1}{\overset{N}{}}q_iq_j\underset{𝐧}{\overset{}{}}𝐫_{ijn}\left[\frac{2\alpha }{\sqrt{\pi }}\frac{e^{\alpha ^2r_{ijn}^2}}{r_{ijn}^2}+\frac{\text{erfc}\left(\alpha |𝐫_{ijn}|\right)}{|𝐫_{ijn}|r_{ijn}^2}\right].$$
(17)
This is always summed using a minimum image spherical cutoff $`|r_{mjn}|<r_c`$ regardless of the system geometry. Therefore the error in truncation of Eq. 17 for the force on particle $`i`$ is

$$\delta 𝐟_i=q_i\underset{j=1}{\overset{N}{}}\underset{𝐧:|𝐫_{ijn}|>r_c}{}q_j𝐫_{ijn}g(|𝐫_{ijn}|),$$
(18)
where $`g`$ is an abbreviation for the term in square brackets. We square this and obtain

$$(\delta 𝐟_i)^2=q_i^2\underset{jk}{}q_jq_k\underset{|𝐫_1|>c}{}\underset{|𝐫_2|>r_c}{}𝐫_1𝐫_2g(|𝐫_1|)g(|𝐫_2|)+q_i^2\underset{j=1}{\overset{N}{}}q_j^2\underset{𝐧:|𝐫_1|>r_c}{\overset{2}{}}𝐫_1^2g^2(|𝐫_1|),$$
(19)
where $`𝐫_1,𝐫_2`$ are the vectors pointing from particle $`i`$ to the $`n^{th}`$ images of $`j`$ and $`k`$ respectively. Since $`g`$ is an even function and $`q_1,q_2`$ will take negative and positive values the average of this quantity over all the $`i`$ particles in the system contains, in the limit of large, uncorrelated $`N`$, only the diagonal terms. Thus

$$(\delta 𝐟_i)^2=q_i^2\underset{j=1}{\overset{N}{}}q_j^2\underset{n:|r_1|>r_c}{}𝐫_1^2g^2(|𝐫_1|).$$
(20)
Moreover in the statistical limit the summation $`\mathrm{\Sigma }_n`$ is independent of particles $`i`$ and $`j`$ and we have

$$(\delta 𝐟_i)^2^{1/2}=|q_i|Q\mathrm{\Sigma }_n^{1/2},$$
(21)
where $`Q^2=q_j^2`$. We will proceed by treating the summation as a spherical integral. First we write the variable explicitly

$$|𝐫_1|=\sqrt{(x+sL_x)^2+(y+tL_y)^2+(z+uL_z)^2}=\rho $$
(22)
where $`s,t,u`$. In Cartesian coordinates we have the following integration

$$\mathrm{\Sigma }_n\frac{8}{L_xL_yL_z}_{\stackrel{~}{c_s}}^{\mathrm{}}𝑑\stackrel{~}{s}_{\stackrel{~}{c_t}}^{\mathrm{}}𝑑\stackrel{~}{t}_{\stackrel{~}{c_u}}^{\mathrm{}}𝑑\stackrel{~}{u}\rho ^2g^2(\rho )=\frac{4\pi }{LL_yL_z}_{r_c}^{\mathrm{}}𝑑\rho \rho ^4g^2(\rho )$$
(23)
where we have applied $`\stackrel{~}{s}=s/L_x`$ etc.

The final step involves the application of the asymptotic integral expansion,

$$_z^{\mathrm{}}f(x)e^{ax^2}𝑑x\frac{f(z)e^{az^2}}{2az}$$
(24)
where the validity of this relation is discussed by Kolafa & Perram. We can apply this formula twice to obtain our mean-squared particle error;

$$\text{erfc}(\alpha \rho )\frac{e^{\alpha ^2\rho ^2}}{2\alpha \rho }$$
(25)
and

$$_c^{\mathrm{}}𝑑\rho e^{2\alpha ^2\rho ^2}\left[\frac{2\alpha }{\sqrt{\pi }}+\frac{1}{2\alpha \rho ^2}\right]^2\left(\frac{1}{\sqrt{\pi r_c}}+\frac{1}{4\alpha ^2r_c^{5/2}}\right)^2e^{2\alpha ^2r_c^2}.$$
(26)
Finally we can obtain the root mean squared force error,

$$\mathrm{\Delta }f=\sqrt{N^1\underset{m=1}{\overset{N}{}}\delta f_m}\frac{Q}{\sqrt{N}}\left(\frac{(\delta f_m)^2}{q_m^2}\right)^{1/2}=\frac{2Q^2}{\sqrt{N}}\frac{1}{\sqrt{L_xL_yL_z}}\left(\frac{1}{\sqrt{r_c}}+\frac{\sqrt{\pi }}{4\alpha ^2r_c^{5/2}}\right)e^{\alpha ^2r_c^2}.$$
(27)
Neglecting the smaller $`\alpha ^2`$ term yields Eq. 13.

### Anisotropic Reciprocal-Space Error

First we write the generalized anisotropic $`k`$-space force on particle $`i`$,

$$𝐟_i^{(k)}=\frac{q_i}{L_xL_yL_z}\underset{j}{}q_j\underset{𝐤0}{}\frac{4\pi 𝐤}{k^2}\mathrm{exp}\left(\frac{k^2}{4\alpha ^2}\right)\mathrm{sin}(𝐤𝐫_{ij}),$$
(28)
where again $`𝐤2\pi (/L_x,/L_y,/L_z)`$.

The summation will be approximated using a cutoff $`\tau `$ and therefore only terms satisfying

$$𝐤^2\tau ^2$$
(29)
will be added. This equation defines the interior of an ellipsoid whose axes are $`L_x/2\pi `$, etc.

Now we compose the cutoff error for a single pair of particles, for which it will prove useful to write using complex notation. Using the fact the function in front of the $`\mathrm{sin}`$ is odd and that half of the vectors $`𝐤`$ are identical but opposite in sign to the other half we obtain

$$\delta 𝐟_{ij}=\frac{q_iq_j}{L_xL_yL_z}\underset{|𝐤|>\tau }{}\frac{4\pi 𝐤}{k^2}\mathrm{exp}\left(\frac{k^2}{4\alpha ^2}\right)i\mathrm{exp}(i𝐤𝐫_{ij}).$$
(30)
Taking the square of this yields

$$(\delta 𝐟_{ij})^2=\left(\frac{q_iq_j}{L_xL_yL_z}\right)^2\underset{|𝐤_1|>\tau }{}\underset{|𝐤_2|>\tau }{}\frac{16\pi ^2𝐤_1𝐤_2}{k_1^2k_2^2}e^{(k_1^2+k_2^2)/4\alpha ^2}\mathrm{exp}(i𝐤_1𝐫_{ij})\mathrm{exp}(i𝐤_2𝐫_{ij}).$$
(31)
At this point we remind ourselves that, in spite of how it looks, this expression is indeed real and positive-valued. We want to obtain the average of this function over all the particle pairs in the system. In the limit of many uncorrelated particles the following expression is exact,

$$\mathrm{exp}(i𝐤_1𝐫_{ij})\mathrm{exp}(i𝐤_2𝐫_{ij})_{ij}=\delta (𝐤_1,𝐤_2).$$
(32)
The details of this require the observation that $`<s_1c_2>=<s_1s_2>=<c_1c_2>=0`$ for $`𝐤_1𝐤_2,𝐤_2`$ and $`<s^2>=<c^2>`$ where the abbreviations $`s_1=\mathrm{sin}(𝐤_1𝐫_{ij})`$ and $`c_1=\mathrm{cos}(𝐤_1𝐫_{ij})`$ have been used. Finally we reach a greatly simplified version of the root-mean-squared pair-wise error:

$$(\delta 𝐟_{ij})^2^{1/2}=\frac{|q_i||q_j|}{L_xL_yL_z}\left(\underset{|𝐤|>\tau }{}\frac{16\pi ^2}{k^2}e^{k^2/2\alpha }\right)^{1/2}.$$
(33)
Now we will be somewhat explicit in treating the summation above. We rewrite the summation in terms of its fundamental variables and length scales,

$$\underset{|𝐤|>\tau }{}\frac{e^{k^2/2\alpha ^2}}{k^2}=\frac{L_x^2}{4\pi ^2}\underset{\rho ^2>k_c^2}{}\frac{\mathrm{exp}[a\rho ^2]}{\rho ^2}.$$
(34)
Here we are taking $`a=2\pi ^2/\stackrel{~}{\alpha }^2=2\pi ^2/(\alpha L)^2`$, $`\lambda =L_x/L_z`$, $`\gamma =L_x/L_y`$, $`\rho ^2=l^2+m^2\gamma ^2+n^2\lambda ^2`$, $`l,m,n`$ and $`k_c`$ is the integer cutoff which defines an ellipsoidal solid upon $`^3`$. With this notation, all of the units of the expression are carried by the factor $`L_x^2`$. This equation can be approximated by integration; in Cartesian coordinates

$$8\frac{L_x^2}{4\pi ^2}_{S1}𝑑l_{S2}𝑑m_{S3}𝑑n\frac{\mathrm{exp}[a(l^2+m^2\gamma ^2+n^2\lambda ^2)]}{l^2+m^2\gamma ^2+n^2\lambda ^2},$$
(35)
where the limits which we do not write here, will collapse into a much simpler form. Conversion to spherical coordinates requires only the substitutions $`\stackrel{~}{n}=\lambda n`$, $`\stackrel{~}{m}=\gamma m`$ and we have

$$\frac{L_x^2}{\pi }\gamma ^1\lambda ^1_{k_c}^{\mathrm{}}\rho ^2𝑑\rho \frac{\mathrm{exp}(a\rho ^2)}{\rho ^2}\frac{L_x^2}{\pi }\gamma ^1\lambda ^1\frac{\mathrm{exp}(ak_c^2)}{2ak_c}.$$
(36)
To summarize the results of this section, we have found

$$\underset{|𝐤|>\tau }{}\frac{e^{k^2/2\alpha ^2}}{k^2}L^2\frac{\stackrel{~}{\alpha }^2}{4\pi ^3\gamma \lambda l_c}\mathrm{exp}(2\pi ^2k_c^2/\stackrel{~}{\alpha }^2),$$
(37)
a result that was obtained by imposing a single-parameter ellipsoidal cutoff which makes simple spherical integration possible. That is to say, if we had not restricted the degrees of freedom of the cutoff object, then we would be faced with the daunting task of ellipsoidal integration. Two assumptions were made, namely that the function varies slowly enough across the grid points so that integration is valid, and secondly that the integral theorem mentioned by Kolafa & Perram is accurate.

We can write the RMS error, first we have the pair-wise error approximated as

$$(\delta 𝐟_{ij})^2^{1/2}=\frac{|q_i||q_j|}{L_yL_z}\frac{2}{\sqrt{\pi }}\frac{\stackrel{~}{\alpha }}{\sqrt{\gamma \lambda k_c}}\mathrm{exp}(\pi ^2k_c^2/\stackrel{~}{\alpha }^2);$$
(38)
which leads to

$$\mathrm{\Delta }f\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{}}\underset{ji}{}(\delta 𝐟_{ij})^2}\frac{Q^2}{\sqrt{N}}\frac{2}{\sqrt{\pi }}\frac{1}{L_yL_z}\frac{\stackrel{~}{\alpha }}{\sqrt{\gamma \lambda k_c}}\mathrm{exp}(\pi ^2k_c^2/\stackrel{~}{\alpha }^2).$$
(39)
## IV Demonstration of Accuracy

In this section we demonstrate the credibility of the EW3DLC technique by means of several standardized test cases that have been used in other literature. These fall into two categories: (1) two-particle systems and (2) a single charge elevated above a planar array of 25 charges. These cases are known to emphasize the possible errors arising from the non-periodicity in the $`z`$-direction. Secondly, several tests of the error formulas are carried out to see that they are indeed applicable to non-cubic, inhomogeneous systems. For these tests it is sufficient to use a standard three-dimensional Ewald summation.

In Fig. 3 the results from a two-particle system are given. In this system the first particle is fixed and the second particle is moving in the $`z`$ direction. For the remainder of this paper we will use $`L_x=L_y=L`$ and unless otherwise noted $`L`$ is set equal to unity. In principle there is no limitation on the accuracy of the EW3DLC or EW3DC methods. The only visible errors (Fig. 3a) occur for the EW3DC method as $`z1`$ since the particle begins to approach the nearest image layer. This error can be eliminated (at some computational cost) by increasing $`L_z`$, but we are interested in illustrating the point.

A richer picture is obtained by plotting the errors on a semi-log scale (Fig. 3b). The EW3DLC curves are complex but these quirks vanish already in 3 or 4 particle systems and are more a signature of the underlying features of the Ewald equations. The key to the explanation of these features is the spherical minimum image cutoff which we chose as $`r_c=0.5`$. When $`z`$ reaches $`0.48`$, the particles are no longer in view of each other, nor are any of the images in view. This leads on the right side of the plot to a significant underestimate of the real-space energy and force and is dominated by the exponential in Eq. 3. The left half of the EW3DLC curves are dominated by Fourier-space error which is trigonometric with a period of $`(k_c+1)L/L_z`$ due to the complex exponential in Eq. 4. Thus the final conclusion is that the basic two-particle curves exhibit nothing problematic in terms of the force and energy signatures.

The next test-system is the 26-particle one originally used by Widmann and Adolf widmann97a : a square array of 25 particles in the $`z=0`$ plane with coordinates $`(x,y,0)`$ where $`x,y(0.1,0.3,0.5,0.7,0.9)`$, and the 26th particle hovers above the center at $`(0.5,0.5,0.2)`$. The values reported previously are easily reproduced using a layer correction. For example, setting $`(\alpha ,r_c,k_c,L_z)=(15.0,0.40,15,0.80)`$ produces $`86.5655`$ and $`10.3642`$ for the total energy and the $`z`$-component force of the 26th particle respectively. For the layer correction the the maximal pairwise error was set at $`10^6`$. The layer correction contribution in this case is important even though the box edge is four times the layer thickness, amounting to $`0.178234`$ and $`0.476241`$ for the force and energy respectively.

With these basic properties of the proposed algorithm established we may begin to probe its applicability to realistic systems. To this end we generate a “random” system of 100 unit charges (overall electrically neutral). We have followed the same convention as Ref. deserno98a  so that these systems can be reproduced by other researchers. A set of 300 random numbers $`n_1,n_2,\mathrm{},n_{300}`$ between 0 and 1 and the coordinates of the charges are taken as $`𝐫_1=(n_1L,n_2L,n_3L_z)`$, $`𝐫_2=(n_4L,n_5L,n_6L_z)`$, … The system was checked to make sure that there are no problematically small spacings. We use this system to check the applicability of the anisotropic error formulas for varying shape and inhomogeneity and in the next section to compare the speed of various algorithms.

For the error formulas it is sufficient to consider the standard Ewald summation (EW3D). In Fig. 4a the results from a standard Ewald summation are compared to the error formulas Eqs. 1314. For the reference forces we also used Ewald summation, convergent to 8 digit accuracy in the total energy. The error formulas are satisfactory for both cigar- and pancake-shaped systems (i.e. $`L_z>1`$ and $`L_z<1`$ respectively).

The second concern is the validity of the error formulas when part of the system is empty. This is tested for 10, 50 and 100% filled systems in Fig. 4b. There is no dependence on $`h`$in the error formulas and therefore the error curves should all be the same. In fact there is good agreement and the most substantial deviations occur in the unimportant large-$`\alpha `$ region. These deviations are evidence of only a small effect of EW3DLC-type inhomogeneity on the error formulas.

To summarize, we have shown that the EW3DLC method produces accurate and well controlled results and can be used for a broad range of shape and inhomogeneity. This is an important advantage over the EW3DC method which is more restrictive on the selection of $`L_z`$ even if the error estimate that we developed in the first paper is used to control the error.

## V Computing Time

Of course the critical consideration, once issues of accuracy are established, is the speed with which a given algorithm can be implemented. This makes the P3MLC the most efficient choice for large $`N`$ but for smaller systems other algorithms such as EW3DLC or MMM2D may have superior speed due to a lower overhead. This section attempts to give some quantitative information on the break-even points between these various algorithms. As usual, the relative times vary slightly with the programming efficiency and the software and hardware of the computer. All times in this study were collected on a Compaq XP1000 workstation.

To optimize the time required for such a many-parameter computation, it is necessary to understand the scaling of time with the variables – most importantly with the number of charges, $`N`$. Here a brief, qualitative discussion of the minimization is made. This discussion tends to be more transparent (with the benefit of hindsight) yet less rigorous than more detailed versions available elsewhere.

The computing time of the real-space part will in general require $`0.5N\rho (4\pi r_c^3/3)`$ operations where $`\rho =N/L^2h`$ (large $`h`$) is the density. Similarly the reciprocal-space calculation scales as $`NL_zk_c^3`$ since increasing $`L_z`$ requires a proportionate increase of the $`z`$-component reciprocal space cutoff. The layer correction scales as $`Nl_c^2`$ since it is a cylindrical summation. Thus the overall time required scales as

$$t_{elc}=a_1N^2r_c^3+a_2NL_zk_c^3+a_3Nl_c^2$$
(40)
As is evident in plots like Fig. 4, the sharpness of the crossover between the real- and reciprocal-space accuracy regimes is such that the choice of $`\alpha `$, $`k_c`$ and $`r_c`$ will be near the optimum if $`\mathrm{\Delta }f_r\mathrm{\Delta }f_k`$. The layer correction cutoff can be determined independently from the others and moreover is only a weak function of $`L_z`$ for large $`(L_zh)`$. Eq. 40 is then interpreted as a function of two variables, for instance the reciprocal cutoff and the desired accuracy $`t_{elc}=f(k_c,\mathrm{\Delta }f)`$. In order to neutralize the change in accuracy in Eqs. 1314, the Ewald parameter should scale inversely with the real cutoff and linearly with the Fourier cutoff, i.e. $`\alpha 1/r_ck_c`$. Thus we can strike a compromise or balance between the first two terms of Eq. 40 by fixing

$$r_{max}N^{1/6},k_{max}N^{1/6}$$
(41)
which reduces the overall scaling of the Ewald calculation to the well-known $`O(N^{3/2})`$. Note that the inclusion of the layer correction does not influence the scaling since it will only add a term linear in $`N`$ to Eq. 40. Furthermore, it becomes a smaller part of the overall calculation as $`N`$ increases.

A similar analysis applies to mesh-based algorithms, of which we will deal with the P3M method. The major difference is that the number of operations required to perform a FFT is independent of the number of charges and only depends on the number of mesh points as $`M\mathrm{log}M`$ per direction in space. The other part of the mesh routine is the charge assignment via polynomial interpolation. The order of interpolation, $`P`$, represents the number of grid points in each direction that a given charge is distributed onto. Therefore it is much like a cutoff parameter and the number of operations per charge scales as $`(P+1)^3`$. The time required for a P3M computation is

$$t_{p3mlc}=a_1N^2r_c^3+a_4L_z(M\mathrm{log}M)^3+a_5N(P+1)^3+a_3Nl_c^2$$
(42)
An intuitive analysis of Eq. 42 is more difficult but analogous to the preceding reasoning (and benefits from foreknowledge of the answer). Again we attempt to balance the $`N`$-exponent in the first two terms but this time a substantially better result is possible due to the much smaller exponent in the second term ($`0`$). From Eq. 15, the choice of $`\alpha M`$ will roughly neutralize the effect of $`M`$ and thus the link between the mesh size and real cutoff is $`r_c1/M`$. The best choice for the scaling of these terms is

$$r_c\frac{(\mathrm{log}N)^{1/2}}{N^{1/3}},M\frac{N^{1/3}}{(\mathrm{log}N)^{1/2}}$$
(43)
resulting in an overall scaling for the time of $`O(N(\mathrm{log}N)^{3/2})`$.

The prefactors in the above equations represent a mixture of “primitive” times for various computing operations as well as, for brevity, some unimportant physical parameters. But in general we expect a performance crossover which favors the Ewald method for low $`N`$ because of a smaller overhead and eventually at large $`N`$, the P3M method becomes much faster.

By generating random systems for a broad range of $`N`$ it is possible to obtain an estimate for the crossover. We compare what are perhaps the three best algorithms for slab-like systems in Fig. 5 at a fixed error, $`\mathrm{\Delta }f=10^2`$ $`e^2/4\pi ϵL^2`$, including the Coulomb prefactor that has been until now omitted. This level of accuracy is as a value that is practical for Langevin and Brownian dynamics simulations. For instance, in a typical $`N=3600`$, $`\rho =4.210^3`$ polyelectrolyte simulation, the random forces are of order $`210^3`$ times larger than the above RMS electrostatic force error.

The recently introduced non-Ewald MMM2D method we recall has a theoretical complexity $`O(N^{5/3})`$. The break-even point between the MMM2D and P3MLC methods at this low accuracy is approximately. $`N^{MMM2DP3MLC}=50`$ and this would increase monotonically with increasing accuracy. For reference the dotted line in Fig. 5 indicates $`N(\mathrm{log}(N))^{3/2}`$ scaling as is expected for the P3MLC method. The EW3DLC curve plotted shows two break-even points at $`N^{MMM2DEW3DLC}=200`$ and $`N^{EW3DLCP3MLC}=17`$. The shape of these curves for the low $`N`$ region of the plot is dictated by the computational overhead which can be linear or constant depending on the algorithm. It is necessary to make Verlet lists outside of the timing loop to avoid the quadratic determination of the pair distances in the real-space calculation. Alternatively, these can be subtracted at the end if one accurately measure the primitive time of this step, as we have done. Note also that the layer correction part of the calculation used roughly 20% of the computation time in the scaling regime and the upper bound on its error was fixed at $`10^2`$.

Finally we make a comparison between the speeds of EW3DLC and EW3DC. Since these methods have the same scaling there is only a difference in the prefactor. But the larger gap size required in EW3DC mandates an increase of the reciprocal space cutoff in order to maintain the same accuracy. Thus there is a computational tradeoff between computing the extra layer correction or expanding the region of reciprocal space summation. We have used three systems with different shapes and $`N=1000`$ for the comparison. The first system has a cubic particle volume sized $`1\times 1\times 1`$ while the other two are pancake and cigar shaped with $`h=0.5`$ and $`2.0`$ respectively. The parameters $`(\alpha ,r_c,k_c,l_c,L_z)`$ were optimized with respect to fixed RMS errors of $`3.0,4.8`$ and $`1.910^2`$ for the pancake, cube and cigar respectively. Again the ELC error was chosen conservatively by setting the upper bound to $`10^2`$. These values were chosen in order to compensate for the different volumes of these systems so that the error at the same density is the same for all three shapes.

The results of this analysis are displayed in Table 2. Under these conditions, the EW3DLC method is only slightly (30%) faster. The reason is that much of the work in these calculations is performed in the real-space and therefore the extra amount of reciprocal-space required by EW3DC due to the enlarged gap is not a major issue. In fact it is almost compensated by the time used by layer correction term. For larger $`N`$, EW3DLC should become faster than EW3DC since the reciprocal-space time becomes much longer than the layer correction time. The same is true about higher accuracy conditions as is observable from the analytical error estimates (Eqs. 13-14). Moreover, in a mesh implementation, the same comparison should yield a more substantial advantage to the layer correction since more of the computational effort is expended in reciprocal space.

## VI Discussion

We have presented and tested a method with broad applicability and low computational cost for computing the electrostatic forces or energy of a cell of charges with two-dimensional periodicity. This method has a complexity that is virtually linear when applied to a mesh, so that only minor further improvements in computation time can be attempted through a reduction of computational overhead. The algorithm differs from EW3DC only by the effort to program the ELC-term, but once this is done the errors are easier to control and there is a greater degree of freedom in choosing the box edges. This makes it easy to apply to cubic versions of the P3M algorithm (which was one of our initial motivations). For larger systems and better accuracy as may be required in Monte Carlo and some molecular dynamics simulations the method should be considerably faster than EW3DC.

It is interesting to estimate the largest system tractable with modern computing. The P3M algorithm can be optimized to be at least twice as fast as the version presented here, so an $`N=10000`$ calculation could be carried out in approximately. 2 seconds. A 100-node parallel computer with a 90% scale factor could handle a system of 0.45 million particles in 1 second – which is sufficient for a many molecular dynamics and Monte Carlo simulations.

Finally we can draw attention to several issues suitable for further research. Most importantly we have not treated the case of dissimilar dielectric materials at the slab surfaces. Because of the symmetry of the consequent charge reflections there could be a way to handle the problem efficiently within the current framework. Another topic is the clarification of the dipole term for slab-wise summation for different dielectric constants at the interface between the periodic supersystem and the surroundings.

###### Acknowledgements.

The authors thank Zuowei Wang and Florian Müller-Plathe for helpful discussions. Financial support from the DFG “Schwerpunkt Polyelektrolyte” is gratefully acknowledged.

## ELC-term Forces and Pressure

For quick reference the force and pressure formulas for the layer correction are given here. The layer correction force on particle $`i`$ is

$`𝐟_i^{(lc)}`$ $`=`$ (44)
$`{\displaystyle \frac{8\pi }{L_xL_y}}({\displaystyle \underset{j}{}}q_iq_j[\mathrm{\hspace{0.33em}2}{\displaystyle \underset{k_x,k_y>0}{}}{\displaystyle \frac{k_x\mathrm{cosh}(k_{}z_{ij})}{k_{}(e^{k_{}L_z}1)}}\mathrm{sin}(k_xx_{ij})\mathrm{cos}(k_yy_{ij})+`$
$`{\displaystyle \underset{k_x>0}{}}{\displaystyle \frac{\mathrm{cosh}(k_xz_{ij})}{(e^{k_xL_z}1)}}\mathrm{sin}(k_xx_{ij})],`$
$`{\displaystyle \underset{j}{}}q_iq_j[\mathrm{\hspace{0.33em}2}{\displaystyle \underset{k_x,k_y>0}{}}{\displaystyle \frac{k_y\mathrm{cosh}(k_{}z_{ij})}{k_{}(e^{k_{}L_z}1)}}\mathrm{cos}(k_xx_{ij})\mathrm{sin}(k_yy_{ij})+`$
$`{\displaystyle \underset{k_y>0}{}}{\displaystyle \frac{\mathrm{cosh}(k_yz_{ij})}{(e^{k_yL_z}1)}}\mathrm{sin}(k_yy_{ij})],`$
$`{\displaystyle \underset{j}{}}q_iq_j[\mathrm{\hspace{0.33em}2}{\displaystyle \underset{k_x,k_y>0}{}}{\displaystyle \frac{\mathrm{sinh}(k_{}z_{ij})}{(e^{k_{}L_z}1)}}\mathrm{cos}(k_xx_{ij})\mathrm{cos}(k_yy_{ij})+`$
$`{\displaystyle \underset{k_x>0}{}}{\displaystyle \frac{\mathrm{sinh}(k_xz_{ij})}{(e^{k_xL_z}1)}}\mathrm{cos}(k_xx_{ij})+{\displaystyle \underset{k_y>0}{}}{\displaystyle \frac{\mathrm{sinh}(k_yz_{ij})}{(e^{k_yL_z}1)}}\mathrm{cos}(k_yy_{ij})]).`$
The diagonal components of the pressure tensor can be obtained from the partition function. For example the normal component is $`P_{}=P_{zz}=\frac{1}{L_xL_y}\frac{E}{L_z}`$. This yields

$`P_{}^{(lc)}={\displaystyle \frac{8\pi ^2}{L_x^2L_y^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}q_iq_j(\mathrm{\hspace{0.33em}2}{\displaystyle \underset{k_x,k_y,>0}{}}`$ $`\nu (k_{})\mathrm{cos}(k_xx_{ij})\mathrm{cos}(k_yy_{ij})+`$
$`{\displaystyle \underset{k_x>0}{}}`$ $`\nu (k_x)\mathrm{cos}(k_xx_{ij})+`$ (45)
$`{\displaystyle \underset{k_y>0}{}}`$ $`\nu (k_y)\mathrm{cos}(k_yy_{ij}))`$ (46)
where

$$\nu (u):=\frac{\left(\mathrm{sinh}uz_{ij}\right)(z_{ij}/L_z)e^{uL_z}\left(\mathrm{sinh}uz_{ij}\right)(z_{ij}/L_z)\left(\mathrm{cosh}uz_{ij}\right)e^{uL_z}}{(e^{uL_z}1)^2}.$$
(47)
For one of the parallel components we find

$`P_{xx}^{(lc)}={\displaystyle \frac{8\pi ^2}{L_x^2L_y^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}q_iq_j(\mathrm{\hspace{0.33em}2}{\displaystyle \underset{k_x,k_y,>0}{}}`$ $`\left[{\displaystyle \frac{k_x^2}{k_{}^2}}\nu (k_{})+{\displaystyle \frac{k_y^2\mathrm{cosh}(k_{}z_{ij})}{k_{}^3L_z(e^{k_{}L_z}1)}}\right]\mathrm{cos}(k_xx_{ij})\mathrm{cos}(k_yy_{ij})+`$
$`{\displaystyle \underset{k_x>0}{}}`$ $`\nu (k_x)\mathrm{cos}(k_xx_{ij})+`$
$`{\displaystyle \underset{k_y>0}{}}`$ $`{\displaystyle \frac{1}{k_yL_z}}{\displaystyle \frac{\mathrm{cosh}(k_yz_{ij})}{(e^{k_yL_z}1)}}\mathrm{cos}(k_yy_{ij})).`$ (48)
Obviously this is cumbersome since after decomposing the trigonometric and hyperbolic terms to obtain a formula that can be evaluated linearly in $`N`$ the expanded version will contain over 80 terms. If rapid pressure calculations are not demanded one can omit this term and increase $`L_z`$. One third of the trace of this tensor will as usual yield the scalar pressure, $`E^{(lc)}/3V`$.

Another consideration in obtaining the pressures is the treatment of the dipole term. Fortunately, a shear or distortion on the system does not result in a distortion in the shape of the supersystem. This is in contrast to the treatment of standard Ewald summation because the spherical summation order requires a complex analysis of the consequences of a distorted sphere.smith94a


format: text

text:
# Calculation Of Secondary Particles In Atmosphere And Hadronic Interactions

## 1 Introduction

Reliable calculations of flux of secondary particles in atmosphere, produced by the interaction of primary cosmic rays, are essential for the correct interpretation of the large amount of experimental data produced by experiments in the field of astroparticle physics. The increasing accuracy of modern experiments demands also an improved quality of the calculation tools. Different ingredients are required to produce a useful calculation model. Essentially they can be reduced to three important classes: the primary cosmic ray spectrum, the modelization of the environment (atmosphere, geomagnetic field, etc.) and a model of particle production in the hadronic shower following the collisions of primary c.r.’s with the atmosphere nuclei. The uncertainty on the primary spectra is dominated by the systematics of the experiments devoted to their measurement: in the light of recent measurements by AMS and BESS, such an uncertainty is about $`\pm `$5% below 100 GeV/nucleon, increasing to $`\pm `$10% at 10 TeV/nucleon. As far as the environmental description is concerned, the large amount of geophysical data now available allows, in principle, to achieve a high level of accuracy. On the other hand, the knowledge of the features of particle production in hadronic interactions is still affected by important uncertainties ($`\stackrel{>}{_{}}\pm `$15%). Since we have not yet a calculable theory for the non-perturbative QCD regime, we remain with many different attempts of building interaction models, which are tuned by comparison to existing experimental data. Sometimes, these models can give satisfactory outputs only in restricted fields of applications. In this work we discuss in some more detail the situation of these attempts to describe hadronic interactions, trying to evidentiate the advantages of the so called “microscopic” models, i.e. those which try to embed as much as possible of the current theoretical ideas in terms of elementary constituents and of their fundamental interactions. We shall make explicit reference to the set of models contained in the FLUKA MonteCarlo code.

## 2 Requirements for Interaction Models

Cosmic ray physics is particularly demanding from the point of view of particle production models, since it is in general necessary to consider a wide range of primary energy and projectiles. If we take for example the case of atmospheric neutrinos, even if we limit ourselves to the class of “contained” or “partially contained” events in Super–Kamiokande, namely events with $`E_\nu `$ in the range from 0.2 to few tens of GeV’s, primary cosmic rays from 1 GeV up to at least 1 TeV per particle must be considered. Furthermore, primary cosmic rays are composed of protons and nuclei, from He to Fe (higher mass can be neglected with good approximation). This is quite a different situation with respect to the standard case of particle physics, where, in general, almost mono-energetic beams are considered. In addition, while in particle physics the attention is mainly devoted to energy deposition, here instead the details of single interactions are fundamental to obtain a flux prediction. As previously stated, it is not possible to rely on a unique model capable of giving the same quality of results at all energies, and for all kinematics regimes. Therefore particular attention is necessary in order to assure the right continuity across the transition region between the different regimes.

Last, but not least, there is the necessity of dealing with different nuclear species. So far, this problem has been mainly solved by recurring to the so called “superposition” model, where a nucleus of mass number A and energy E0 is considered to be equivalent to A nucleons, each one having energy E0/A. The question if this is a totally acceptable approximation is still an argument of discussion.

## 3 Different approaches: parametrized vs “microscopic” models

In building a suitable model for particle production, we can identify two main different attitudes: parametrized codes and “microscopic” models. In the first case, one relies upon analytical formulas derived from some general phenomenological features of particle production, and with parameters that are obtained from fits to experimental data. An example of these fundamental properties of hadronic interactions, is Feymann scaling, which is known to be a rather good (although not completely exact) approximation, especially in the forward (“fragmentation”) kinematic region which is dominant in secondary production by cosmic rays. This property can be expressed as follows. The number of pions of energy $`E_\pi `$ produced in an interaction by a primary proton of energy $`E_0`$ is well represented expression:

$$\frac{dn_\pi }{dE_\pi }(E_\pi ,E_0)\frac{1}{E_\pi }F\left(\frac{E_\pi }{E_0}\right)$$
(1)
where the function $`F(x)`$, $`x`$ = $`E_\pi /E_0x_{Feynman}`$, is approximately independent from the primary energy $`E_0`$, and decreases monotonically from a finite value for $`x0`$, to zero for $`x1`$. The shape of the curve can be easily expressed by means of a combination of elementary functions with just a few parameters that can be extrated from experimental data sets. A noticeable example of this kind of approach is at the basis of the work which has been carried on by the Bartol group for many years, producing many valuable results, and in particular the prediction for atmospheric neutrino fluxes. For this purpose, they constructed the TARGET numerical model, a module which can be easily inserted in any cascade program. TARGET considers hadron interactions on light nuclei, like Oxygen and Nitrogen, subdividing the available energy between leading nucleons and other produced hadrons on the basis of an assumed elasticity function. Pions, kaons, etc. are then produced according to parametric formulas reproducing the scaling properties described above. Experimental data at different energies fix the parameters and guide the evolution of multiplicities as a function of energy. At low and intermediate energies, resonance production is considered. Care has been taken to assure event by event energy conservation. The advantage of this kind of approach, mainly used in the framework of a 1-dimensional description111an upgrade of the TARGET model, also in view of 3–D applications, has been presented in , is that it can lead to the comprehension of some important and general properties of particle production in terms of analytical expressions. The price to pay is the lack of generality and the spoiling of correlations among reaction products.

The second line of approach is instead the use of models which try to describe interactions in terms of the properties of elementary constituents. In principle one would like to derive all features of “soft” interactions (low-$`p_T`$ interactions) from the QCD Lagrangian, as it is done for hard processes. Unfortunately the large value taken by the running coupling constant prevents the use of perturbation theory. Indeed, in QCD, the color field acting among quarks is carried by the vector bosons of the strong interaction, the gluons, which are “colored” themselves. Therefore the characteristic feature of gluons (and QCD) is their strong self-interaction. If we imagine that quarks are held together by color lines of force, the gluon-gluon interaction will pull them together into the form of a tube or a string. Since quarks are confined, the energy required to “stretch” such a string is increasingly large until it suffices to materialize a quark-antiquark couple from the vacuum and the string breaks into two shorter ones, with still quarks at both ends. Therefore it is not unnatural that because of quark confinement, theories based on interacting strings emerged as a powerful tool in understanding QCD at the soft hadronic scale (the non-perturbative regime). Different implementations of this idea exist, having obtained remarkable success in describing the features of hadronic interactions. Some of these codes have already found applications in astroparticle physics. We can quote a few major examples:

1. the atmospheric neutrino flux by M. Honda et al., which has been obtained by a combinations of microscopic codes embedded into an original shower code;
2. the CORSIKA shower code, which offer the choice among different microscopic models.
3. the already mentioned general purpose MonteCarlo code FLUKA, which is now applied also to cosmic ray physics by us and other authors, and that will be later described in more detail.

In our opinion, in the microscopic approach each step has sound physical basis and allows to reach a deep understanding of the phenomena and a high reliability of predictions. The performances are optimized comparing with particle production data at single interaction level. The final predictions are obtained with a minimal set of free parameters, fixed for all energies and target/projectile combinations. Results in complex cases as well as scaling laws and properties come out naturally from the underlying physical models. The basic conservation laws are fulfilled “a priori”. A microscopic model can reach a very high level of detail, at least in principle, and therefore is a good choice when aiming at precision calculations. The price to pay is the loss of simplicity and flexibility: there are no more simple analytical guidelines which allow to understand the basic properties. Furthermore microscopic codes are more demanding than parametrizations in terms of computing power. Parametrized models (if parametrizations are performed at the level of single interactions) are instead useful as a first, fast and flexible approach.

Instead, models tuned on “integral data”, like calorimeter resolutions, thick target yields etc., can be very inaccurate at the level of single interactions, as shown in ref.  for the case of GEANT-GHEISHA: such a model cannot be used to obtain a reliable calculation of particle fluxes. In our opinion this might be a problem for the low energy ($`<`$80 GeV) calculations performed with CORSIKA, when GHEISHA is selcted under that energy threshold.

In the following section we shall concentrate on the example of the FLUKA MonteCarlo code.

## 4 The FLUKA model

The modern FLUKA is an interaction and transport MonteCarlo code able to treat with a high degree of detail the following problems:

* Hadron-hadron and hadron-nucleus interactions 0-100 TeV
* Electromagnetic and $`\mu `$ interactions 1 keV-100 TeV
* Charged particle transport - ionization energy loss
* Neutron multigroup transport and interactions 0-20 MeV
* Nucleus-nucleus and hadron-nucleus interactions 0-10000 TeV/n: under development

Here we shall review the two hadronic models which are used inside FLUKA to describe nonelastic interactions:

* The “low-intermediate” energy one, PEANUT, which covers the energy range up to 5 GeV
* The high energy one which can be used up to several tens of TeV, based on the color strings concepts sketched in the previous section.

The nuclear physics embedded in the two models is very much the same. The main differences are a coarser nuclear description (and no preequilibrium stage) and the Gribov-Glauber cascade for the high energy one.

### 4.1 The PEANUT Model

Hadron-nucleus non-elastic interactions are often described in the framework of the IntraNuclear Cascade (INC) models. This kind of model was developed at the very beginning of the history of energetic nuclear interaction modelling, but it is still valid and in some energy range it is the only available choice. Classical INC codes were based on a more or less accurate treatment of hadron multiple collision processes in nuclei, the target being assumed to be a cold Fermi gas of nucleons in their potential well. The hadron-nucleon cross sections used in the calculations are free hadron–nucleon cross sections. Usually, the only quantum mechanical concept incorporated was the Pauli principle. Possible hadrons were often limited to pions and nucleons, pions being also produce or absorbed via isobar (mainly $`\mathrm{\Delta }_{33}`$) formation, decay, and capture. Most of the historical weaknesses of INC codes have been mitigated or even completely solved in some of the most recent developments , thanks to the inclusion of a so called “preequilibrium” stage, and to further quantistic effects including coherence and multibody effects.

All these improvements are considered in the PEANUT (PreEquilibrium Approach to NUclear Thermalization) model of FLUKA. Here the reaction mechanism is modelled in by explicit intranuclear cascade smoothly joined to statistical (exciton) preequilibrium emission  and followed by evaporation (or fission or Fermi break-up) and gamma deexcitation. In both stages, INC and exciton, the nucleus is modelled as a sphere with density given by a symmetrized Woods-Saxon  shape with parameters according to the droplet model  for A$`>`$16, and by a harmonic oscillator shell model for light isotopes (see ). The effects of the nuclear and Coulomb potentials outside the nuclear boundary are included. Proton and neutron densities are generally different. Binding Energies are obtained from mass tables. Relativistic kinematics is applied at all stages, with accurate conservation of energy and momentum including those of the residual nucleus. Further details and validations can be found in .

For energies in excess of few hundreds MeV the inelastic channels (pion production channels) start to play a major role. The isobar model easily accommodates multiple pion production, for example allowing the presence of more than one resonance in the intermediate state (double pion production opens already at 600 MeV in nucleon-nucleon reactions, and at about 350 MeV in pion-nucleon ones). Resonances which appear in the intermediate states can be treated as real particles, that is, they can be transported and then transformed into secondaries according to their lifetimes and decay branching ratios.

### 4.2 The Dual Parton Model for high energy

A theory of interacting strings can be managed by means of the Reggeon-Pomeron calculus in the framework of perturbative Reggeon Field Theory, an expansion already developed before the establishment of QCD. Regge theory makes use explicitly of the constraints of analyticity and duality. On the basis of these concepts, calculable models can be constructed and one of the most successful attempts in this field is the so called “Dual Parton Model” (DPM), originally developed in Orsay in 1979 . It provides the theoretical framework to describe hadron-nucleon interaction from several GeV onwards. In DPM a hadron is a low-lying excitation of an open string with quarks, antiquarks or diquarks sitting at its ends. In particular mesons are described as strings with their valence quark and antiquark at the ends. (Anti)baryons are treated like open strings with a (anti)quark and a (anti)diquark at the ends, made up with their valence quarks.

At sufficiently high energies, the leading term in high energy scattering corresponds to a “Pomeron” ($`IP`$) exchange (a closed string exchange with the quantum numbers of vacuum), which has a cylinder topology. By means of the optical theorem, connecting the forward elastic scattering amplitude to the total inelastic cross section, it can be shown that from the Pomeron topology it follows that two hadronic chains are left as the sources of particle production (unitarity cut of the Pomeron). While the partons (quarks or diquarks) out of which chains are stretched carry a net color, the chains themselves are built in such a way to carry no net color, or to be more exact to constitute color singlets like all naturally occuring hadrons. In practice, as a consequence of color exchange in the interaction, each colliding hadron splits into two colored system, one carrying color charge $`c`$ and the other $`\overline{c}`$. These two systems carry together the whole momentum of the hadron. The system with color charge $`c`$ ($`\overline{c}`$) of one hadron combines with the system of complementary color of the other hadron, in such a way to form two color neutral chains. These chains appear as two back-to-back jets in their own centre-of-mass systems. The exact way of building up these chains depends on the nature of the projectile-target combination (baryon-baryon, meson-baryon, antibaryon-baryon, meson-meson). Let us take as example the case of nucleon-nucleon (baryon-baryon) scattering. In this case, indicating with $`q_p^v`$ the valence quarks of the projectile, and with $`q_t^v`$ those of the target, and assuming that the quarks sitting at one end of the baryon strings carry momentum fraction $`x_p^v`$ and $`x_t^v`$ respectively, the resulting chains are $`q_t^vq_p^vq_p^v`$ and $`q_p^vq_t^vq_t^v`$, as shown in fig. 1.

Energy and momentum in the centre-of-mass system of the collision, as well as the invariant mass squared of the two chains, can be obtained from:

$`E_{ch1}^{}`$ $``$ $`{\displaystyle \frac{\sqrt{s}}{2}}(1x_p^v+x_t^v)`$
$`E_{ch2}^{}`$ $``$ $`{\displaystyle \frac{\sqrt{s}}{2}}(1x_t^v+x_p^v)`$
$`p_{ch1}^{}`$ $``$ $`{\displaystyle \frac{\sqrt{s}}{2}}(1x_p^vx_t^v)=p_{ch2}^{}`$ (2)
$`s_{ch1}`$ $``$ $`s(1x_p^v)x_t^v`$
$`s_{ch2}`$ $``$ $`s(1x_t^v)x_p^v`$
The single Pomeron exchange diagram is the dominant contribution, however higher order contributions with multi-Pomeron exchanges become important at energies in excess of 1 TeV in the laboratory. They correspond to more complicated topologies, and DPM provides a way for evaluating the weight of each, keeping into account the unitarity constraint. Every extra Pomeron exchanged gives rise to two extra chains which are built using two $`q\overline{q}`$ couples excited from the projectile and target hadron sea respectively. The inclusion of these higher order diagrams is usually referred to as multiple soft collisions.

Two more ingredients are required to completely settle the problem. The former is the momentum distribution for the $`x`$ variables of valence and sea quarks. Despite the exact form of the momentum distribution function, $`P(x_1,..,x_n)`$, is not known, general considerations based on Regge arguments allow to predict the asymptotic behavior of this distribution whenever each of its arguments goes to zero. The behavior turns out to be singular in all cases, but for the diquarks. A reasonable assumption, always made in practice, is therefore to approximate the true unknown distribution function with the product of all these asymptotic behaviors, treating all the rest as a normalization constant.

The latter ingredient is a hadronization model, which must take care of transforming each chain into a sequence of physical hadrons, stable ones or resonances. The basic assumption is that of chain universality, which assumes that once the chain ends and the invariant mass of the chain are given, the hadronization properties are the same regardless of the physical process which originated the chain. Therefore the knowledge coming from hard processes and $`e^+e^{}`$ collisions about hadronization can be used to fulfill this task. There are many more or less phenomenological models which have been developed to describe hadronization (examples can be found in ). In principle hadronization properties too can be derived from Regge formalism .

It is possible to extend DPM to hadron-nucleus collisions too , making use of the so called Glauber-Gribov approach. Furthermore DPM provides a theoretical framework for describing hadron diffractive scattering both in hadron-hadron and hadron-nucleus collisions. General informations on diffraction in DPM can be found in  and details as well as practical implementations in the DPM framework in .

At very high energies, those of interest for high energy cosmic ray studies (10–105 TeV in the lab), hard processes cannot be longer ignored. They are calculable by means of perturbative QCD and can be included in DPM through proper unitarization schemes which consistently treat soft and hard processes together. The interested reader can find more informations as well as practical implementations and results in .

DPM exhibited remarkable successes in predicting experimental observables. The quoted references include a vast amount of material showing the capabilities of the model when compared with experimental data. However, it must be stressed that other models are available, but most of them share an approach based on string formation and decay. For example, the Quark Gluon String Model  has been developed more or less in parallel with DPM. This model shares most of the basic features of DPM, while differing for some details in the way chains are created and in the momentum distribution functions.

## 5 Benchmarks of the FLUKA Model

The predictions of FLUKA have been checked with a large set of experimental data collected in accelerator experiments. Here we shall limit ourselves to show only a few examples, among the most important in view of the application of the code to cosmic ray applications.

Two sets of data are of particular relevance to check the quality of a model to be used for the calculation of atmospheric neutrino fluxes. These concern p-Be collisions and are reported in fig. 2: in ref. the central rapidity region has been mainly explored, while in ref. the forward region has been investigated. In both cases the agreement of FLUKA predictions is quite good.

Measurements of $`\pi ^\pm `$ and $`K^\pm `$ production rates by 400 GeV/c protons on Be targets were performed by Atherton et al.  for secondary particle momenta above 60 GeV/c and up to 500 MeV/c of transverse momentum. Recently the NA56/SPY (Secondary Particle Yields) experiment  was devoted to directly measure these yields in the momentum region below 60 GeV/c. The SPY experiment measured the production at different angles $`\theta `$ and momenta $`P135`$ GeV/c down to 7 GeV/c for pions, kaons, protons and their antiparticles, using a 450 GeV/c proton beam impinging on Be targets. These data were extremely valuable to improve the hadronization model of FLUKA so to arrive at the present version. FLUKA is in agreement with the Atherton and the SPY measurements at the level of $`20\%`$ in the whole momentum range of all secondaries, with the exception of a few points mostly for negative kaons. The case of pions is reported in fig.3. Also the $`\theta `$ dependence of the measured yields is reasonably described by FLUKA. The measured $`K^\pm /\pi ^\pm `$ ratios are reproduced to better than $`20\%`$ below 120 GeV/c.

This example was of particular relevance, since other attempts, using for instance the hadronic interfaces of GEANT (and in particular GEANT-GHEISHA) yielded a much worser agreement, as shown in ref..

## 6 Example of calculations of particles in atmosphere

In the last years, the FLUKA interaction models has been used to produce new predictions for the atmospheric neutrino fluxes within a full 3D calculation. These fluxes have been also considered by the Super–Kamiokande experiment. In the framework of the same group, a new calculation has been developed, choosing again another microscopic code based on the Dual Parton Model: DPMJET-III. It gives results close to those obtained with FLUKA.

In the last two years a considerable amount of work has been devoted to cross check the validity of the calculation model. As far as the FLUKA approach is concerned, at least two remarkable results can be quoted:

1. The reproduction of the features of primary proton flux as a function of geomagnetic latitude as measured by AMS, thus showing that the geomagnetic effects and the overall geometrical description of the 3–D setup are well under control. In addition, the same work shows that also the fluxes of secondary $`e^+e^{}`$ measured at high altitude (eventually the last stage of the chain decay of produced mesons) are reproduced.
2. The good reproduction of the data on muons in atmosphere as measured by the CAPRICE experiment, both at ground level and at different floating altitudes, when starting from the same primary flux (Bartol fit) used to generate atmospheric neutrinos. See the quoted reference for relevant plots and numbers.

The fluxes of atmospheric muons are strictly related to the neutrino ones, because almost all $`\nu `$’s are produced either in association, with, or in the decay of $`\mu ^\pm `$. Therefore it is possible to conclude that, for that choice of primary spectrum, the $`\nu `$ fluxes predicted by FLUKA are probably in the right range. To a large extent the agreement between the original HKKM and Bartol calculations of the $`\nu `$ fluxes, despite they started from different estimates of the primary flux and different hadronic interaction models, is not casual, but the result of the $`\mu `$ constraint. Furthermore, the agreement exhibited by the FLUKA simulation for muons of both charges gives confidence on the predictions of FLUKA for the parent mesons of muons (mostly pions).

The shower simulations in atmosphere have been compared also to the most recent hadron spectra at different latitudes and altitudes, obtaining remarkable agreement. As an example, in fig. 4 we compare MonteCarlo results to the hadron flux measured with the KASKADE experiment.

## 7 Conclusions

The phenomenological study of hadronic interactions is still a fundamental issue for astroparticle physics. Now that, at least for primary energy lower than 100 GeV, the uncertainties on primary spectrum have been substantially reduced by the quality of new experiments like AMS and BESS, the need for a better quality model of hadronic interactions is even more necessary if accuracy of predictions has to be pursued. For this goal the “microscopic” codes are mostly recommended, thanks to their predictive power in a very large kinematic region, constrained by a limited number of parameters. Uncertainties on the modelling of hadronic interactions will remain a fundamental issue, and probably only new data, if experimental systematics can be kept under reasonable control, will help model builders. The HARP experiment at CERN is aiming at this goal. This kind of activity is beneficial not only for particle physics and astrophysics, but also for applied science, since these calculations are necessary to understand radiation fluxes in the Earth’s atmosphere, and this is of great interest for civil aviation and for the design of satellite activities. The FLUKA MonteCarlo model is already being used for this purpose: doses to commercial flight are the subject of a work in progress, together with the development of a specific model for heavy ion transport and interaction: this will be of the utmost importance for dose and damage calculation in space aircrafts.