text:
# Metal abundances and ionization conditions in a possibly dust-free damped Ly𝛼 system at 𝑧=2.3Based on observations collected at the European Southern Observatory, Chile (Programme ID: 65.O-0411(A))

## 1 Introduction

Damped Ly$`\alpha `$ systems (DLAs), defined as those QSO absorption systems with the highest column densities, $`N(H\text{i})>2\times 10^{20}`$ cm-2, are widely believed to occur in the ISM of the high redshift progenitors of present-day galaxies (Wolfe Wolfe (1993)). It is fortunate that such large column densities permit very accurate metal abundance measurements because, since DLAs contributed most of the neutral gas mass when the Universe was 10% of its present age (e.g., Storrie-Lombardi et a. Storrie (2000)), they provide us with a powerful means to trace the chemical evolution of star forming galaxies over a considerable fraction of the Hubble time. Since the absorbing gas in DLAs is expected to be mostly neutral, the bulk of metals should be in neutral or singly ionized form (depending on whether or not the first ionization potential IP $`>13.6`$ eV). Thus, the many transitions of these “low ions” that are redshifted to the optical range can be used to compute \[X+/H0\] $``$ \[X/H\] $`\mathrm{log}(\mathrm{X}/\mathrm{H})\mathrm{log}(\mathrm{X}/\mathrm{H})_{\mathrm{}}`$.

With the advent of 8-10m-class telescopes high resolution spectroscopy of distant QSOs and detailed abundance studies in DLAs have become possible with relatively high accuracy ($`10`$ %). The picture that has emerged is that of $`Z1/101/300Z_{\mathrm{}}`$ absorbers with abundance patterns that, except for the Fe-peak ratios, resemble that of Milky Way metal-poor halo stars (Lu et al. Lu (1996); Pettini et al. Pettini (1997); Prochaska & Wolfe Prochaska (1999); Boisse et al. Boisse (1998)). These observations seem to support the underlying hypothesis that at least part of the DLAs have undergone a similar chemical evolution history as our own Galaxy.

On the other hand, the observed sub-solar abundances \[X/H\] observed in the Galactic ISM are interpreted as due to a partial incorporation of the atoms from the gas phase into dust grains (e.g., Savage & Sembach Savage (1996)). Since dust is formed through condensation from the gas phase, any measure of the absorbing column density of a refractory element will underestimate the true number of atoms. Moreover, observations of diffuse clouds in the Galactic ISM have shown the clear trend of different elements to condense into grains in different proportions, depending on condensation temperature (Savage & Sembach Savage (1996)). Fe, Cr and Ni are among the most depleted elements in the ISM, while O, P, Ar, S, and Zn have generally very small depletion factors. Since DLAs are likely to sample the ISM of galaxies at earlier evolutionary stages, it seems reasonable to raise the question as to whether the measured abundances at high redshift have also been modified by differential dust depletion (e.g., Lu et al. Lu (1996); Pettini et al. Pettini (1997)). Comparing abundances of refractory and non-refractory elements in high redshift DLAs can lead to estimates of the dust-to-gas ratio in those environments (Pettini et al. Pettini (1997); Pettini et al. Pettini1 (1999)); assuming particular depletion patterns can in turn be used to correct the measured gas-phase abundances (Vladilo et al. Vladilo1 (1998)).

However, any assessment of the relative abundances in DLAs is difficult as long as the interplay between dust and nucleosynthesis effects is not well understood, and may not be equal in high redshift DLAs and the local ISM.

The most direct way to disentangle these two contributions to the observed abundance pattern is to find DLAs with as small a level of depletion as possible. Such systems, though, are quite scarce and only few cases are reported in the literature (Molaro et al. Molaro1 (2000); Pettini et al. Pettini2 (2000); Prochaska & Wolfe Prochaska3 (1999)). In this paper we present VLT UVES observations of the $`z=3.01`$ QSO HE 2243$``$6031 (discovered by the Hamburg/ESO QSO Survey; Reimers & Wisotzki Reimers1 (1997)), and argue that the DLA observed at $`z=2.33`$ shows negligible amounts of dust. We have been able to accurately measure abundances for a set of elements of different nucleosynthetic origins and with differing susceptibility to incorporation into dust grains, including the detection of the “rare” elements P and Ar.

The paper is organized in three parts. First we address the question of the physical conditions in the single velocity components and perform a cloud-by-cloud analysis of the relative abundances (§ 4); next we look at the uncertainties introduced by possible ionization effects (§ 5); and the final part (§ 6) is devoted to studying dust-depletion effects and relative abundances. The conclusions are outlined in § 7.

## 2 Observations and data reduction

HE 2243$``$6031 was observed in service mode with the UVES instrument at the ESO Kueyen telescope in August 2000 under good seeing conditions ($`0.8\mathrm{}`$). An extra 1-hour exposure had been previously taken during science verification of the instrument. Using two dichroic modes the whole optical range was covered with a total exposure time of $`\mathrm{28\hspace{0.17em}800}`$ sec (Table 1).

After bias-subtracting and flat-fielding of the individual CCD frames, the echelle orders were extracted and reduced with self-implemented routines running under MIDAS. The algorithm attempts to reduce the statistical noise to a minimum by fitting the seeing profile with Gaussians. Flux values are assigned with a variance according to the Poisson statistics and the read-out noise, while cosmic-ray hits are assigned with infinite variances. The extraction proceeds in 3 steps: (1) A Gaussian with unconstrained parameters is fitted to the sky profile at each wavelength using the Levenberg-Marquardt method (Press et al. Press (1986)); (2) the variation of width and position with respect to the orders is then fitted along the dispersion direction with low-order polynomials; (3) step (1) is repeated, this time with width and position fixed at the values given by the fit solution, so only the amplitudes can vary.

The extracted orders were wavelength calibrated using as reference Th-Ar spectra taken after each science exposure. All wavelength solutions were accurate to better than typically $`1/10`$ pixel. The wavelength values were converted to vacuum heliocentric values and each order of a given instrumental configuration was binned onto a common linear wavelength scale of $`0.043`$ Å pixel-1. The reduced orders were then added with a weight according to the inverse of the flux variances. Finally, the flux values were normalized by a continuum that was defined using cubic splines over featureless spectral regions. The spectral resolution is FWHM $`6.5`$ km s-1. The typical signal-to-noise ratio varies between S/N$`80`$ at 3700–5700 Å, and S/N$`50`$ at 5700–8000 Å (becoming worse at both ends of the spectrum).

In addition to the UVES data, a low resolution (FWHM $`4`$ Å) spectrum of HE 2243$``$6031 was obtained on October 12 2001 using the Boller and Chivens (B&C) Spectrograph on the Baade 6.5m Telescope at Las Campanas Observatory. This spectrum was used to better define the quasar continuum in the spectral region around the damped Ly$`\alpha `$ line.

## 3 Line parameters

Column densities $`N`$ and Doppler parameters $`b`$ were obtained by fitting Voigt profiles to each velocity component. The fits were performed using the minimum-$`\chi ^2`$ code FITLYMAN (Fontana & Ballester Fontana (1995)). Since most of the species under study present 2 or more transitions, the results of the fits were in general stable and quite independent of the fitting region chosen.

As an alternative approach to obtain column densities in DLAs, the apparent optical depth method (Savage & Sembach Savage1 (1991)) has been often used to get integrated column densities. The method is straightforward and gives accurate results for a resolved and non-saturated line (Lu et al. Lu (1996); Prochaska & Wolfe Prochaska (1999)). In our case, although most of the ions have at least one transition likely not to be saturated, we use instead the fitting approach, because it provides us with better column density estimates of partially blended components, and of those lines having evident departures from the linear part of the curve-of-growth. In addition, unlike previous studies, in which redshift and Doppler parameter of a given velocity component are held fixed for all ions, in our fitting procedure all fit parameters were allowed to vary.111Traditionally, line profile fitting of DLA systems has been performed by assuming the same Doppler parameter for all ions at given velocity component, i.e., that bulk motion dominates over thermal motion. This is a reasonable choice as long as DLAs occur in cold gas and the ions under study are massive; however, this does not prevent us a priori to chose the free-parameter approach described above.

The fit results are summarized in Table 2, while the column labeled $`N_{\mathrm{fit}}`$ in Table 3 lists the overall column densities which resulted from summing the individual component column densities. All line parameters come from the Voigt profile fits, with the exception of $`N(Zn\text{ii})`$ and $`N(Cr\text{ii})`$, which were obtained using the apparent optical depth method. Following is a description of the fits to individual species.

### 3.1 Description of the fits

#### Fe ii.—

We begin by fitting Fe ii because this is the ion with the most complete set of transitions redward of the Ly$`\alpha `$ forest. We fitted simultaneously 5 Fe ii transitions. They are shown in Fig. 1, plotted in velocity space with respect to $`z=2.330000`$, and along with the fitted profiles (smooth line). The $`\lambda 2344`$ and $`\lambda 2382`$ transitions provided a fitting interval only at the position of the weaker components. We found that 14 velocity components were necessary to model the whole profiles. Their positions are indicated by the vertical dashed lines. We used these line positions to identify – but not to tie – velocity components in the remaining ions and throughout the paper we refer to them as to “C1”, “C2”, etc. We were not able to fit the absorption feature between C11 and C12. There appears not to be enough information on the line profile for FITLYMAN to find the best line centroid. Instead of defining a fixed line position ‘by hand’ we preferred to leave the feature unfitted. An optical depth analysis shows that the contribution to the total Fe ii column density should be less than 1%.

#### Ar i.—

We have detected Ar i via the $`\lambda 1048`$ transition. The $`\lambda 1066`$ line unfortunately does not provide any useful constraint to the fits. Although the $`\lambda 1048`$ transition lies in the forest, the spectrum is clean at the expected position of its strongest components. Absorption by Ar i is detected in C3, C4 and C11, which correspond to the three strongest lines in Fe ii ($`50`$ % of the total Fe ii). The small line widths suggest these lines are very likely not due to weak Ly$`\alpha `$ interlopers. Moreover, the position of the lines match well those of Fe ii, with exception maybe of C4, due to blending with a Ly$`\alpha `$ forest line (which makes the fit an overestimate of $`N(Ar\text{i})`$ in this component). This is the second ever measurement of argon at high redshift after the $`z=3.39`$ detection by Molaro et al. (Molaro (2001)), and one of the few outside the Galaxy available (Vidal-Madjar et al. Vidal-Madjar (2000); Levshakov, Kegel & Agafonova Levshakov (2001)).

#### P ii.—

We have detected P ii via the $`\lambda 963`$ and $`\lambda 1152`$ transitions (Fig. 2). We are quite confident of the detection and fit results because they are based on two transitions that happen to fall in a relatively unabsorbed part of the forest. This is only the third detection of extragalactic phosphorus so far (Molaro et al. Molaro (2001), Outram et al. Outram (1999)). The fits give reasonably good results for C3, C4, C6, and C8. C11 is blended with a Ly$`\beta `$ line at $`z=2.743903`$ with $`\mathrm{log}N=14.57`$ and $`b=25.1`$. These parameters resulted from fitting simultaneously the P ii lines, the Ly$`\beta `$, and its corresponding Ly$`\alpha `$.

#### S ii, N i, Si ii, Ni ii, Al ii.—

We did not encounter big difficulties by fitting these ions because all of them but Al ii show more than one transition, although we were able to fit all 14 components only for Si ii. The Al ii $`\lambda 1670`$ line is saturated over most of the components, so column densities of the stronger components must be treated as lower limits. In S ii and N i there is obvious blending with forest lines. C1 and C2 in N i were treated as one single component.

#### O i, C ii.—

For these ions most transitions have saturated components or they lie in the forest, so only a small set of components could be fitted. C3 in C ii and C8 in O i are uncertain because of the small fitting interval available. C1 and C2 in O i were treated as one single component.

#### Mn ii.—

Mn ii is not detected (only the weak $`\lambda 2606`$ is available, the other transitions falling in the red cross-disperser gap at $`8600`$ Å). We derive a $`3\sigma `$ upper limit of $`\mathrm{log}N<11.6`$ for a single line. If this value is scaled according to the strongest Fe ii component, it gives a total Mn ii $`3\sigma `$ limit of $`\mathrm{log}N(Mn\text{ii})<12.43`$.

#### Zn ii and Cr ii.—

To estimate column densities for Zn ii and Cr ii we used the Zn ii $`\lambda 2026`$ and Cr ii $`\lambda 2056`$ transitions. These are the strongest transitions available and should not be blended with unrelated lines. The $`\lambda 2062`$ line is partially lost within telluric lines.

Unfortunately, the $`\lambda 2026`$ and $`\lambda 2056`$ lines lie in the noisiest part of our spectrum, where spurious features introduced by traces of fringing could not be removed properly. In addition, the Zn ii line happens to fall in the overlapping region of two echelle orders, where the continuum definition is inaccurate. This was reflected in the impossibility of profile fitting all components. For this reason we used the apparent optical depth method to constrain the respective column densities. Integrating between $`v=180`$ and $`v=+100`$ km s-1 leads to $`\mathrm{log}N(Zn\text{ii})=12.22\pm 0.03`$ and $`\mathrm{log}N(Cr\text{ii})=13.24\pm 0.02`$. We have excluded the velocity region \[+100,+180\], due to spurious features that introduce extra noise (excluding C13 and C14 introduces an uncertainty no larger than $`0.5`$%, though). We stress, however, that due to the weakness of the lines these values are subject to continuum uncertainties (which are not considered in the formal errors).

#### Doubly ionized species.—

Fig. 7 displays the expected position in velocity space of all doubly ionized species covered by our spectrum. N ii is also included in the figure. It is of great importance to estimate or put useful limits on column densities for these ions, as they will constrain the ionization level of the system.

For Al iii we have at our disposal two transitions with apparently unsaturated lines, so the fit results can be considered reliable. Note that C13 and C14 are not detected.

For the rest of the ions we have attempted to put upper limits on column densities using different approaches. For S iii we use the transition at $`\lambda =1012`$ Å. Although the line lies in the forest, there is a relatively good match of the absorption profile with that of Fe ii, in line position as well as in line strengths. This suggests the S iii is present in the DLA gas and tracks the low-ion profiles. However, since the lines stand out only poorly we did not attempt a fit; instead, we used the total apparent column density, which gives $`N(S\text{iii})=14.05\pm 0.01`$ for $`v=`$ \[$`130,+80`$\].

To constrain Fe iii we use the $`\lambda 1122`$ transition. This line is also in the forest and appears so severly contaminated that no resemblance with Al iii or low-ion profiles is recognized. An upper limit was estimated by assuming a Fe iii line lies at C3 and has $`b=b_{Fe\text{ii}}=8.0`$ km s-1. The maximum column density allowed by the data then becomes $`\mathrm{log}N=13.4`$, which, weighted with C3 from the Al iii fits yields a total column density of $`N(Fe\text{iii})<14.18`$. Note, however, that the Fe iii profile does not exclude $`N(Fe\text{iii})=0`$.

For P iii we performed the same exercise using the expected position of the $`\lambda 998`$ transition, which also shows no black absorption. At C3 we get $`\mathrm{log}N(P\text{iii})<13.2`$ for $`b=8.4`$ km s-1, or $`\mathrm{log}N(P\text{iii})<13.98`$ cm-2 for the whole profile. Again, since the line is in the forest and we do not recognize an Al iii-like absorption pattern, we cannot exclude $`N(P\text{iii})=0`$.

Finally, the $`\lambda 1083`$ profile, despite forest contamination, also resembles the low-ion pattern. This is interesting given that nitrogen is expected to be mainly neutral due to its high IP of $`14.5`$ eV. Assuming that the absorption feature at C11 is a N ii $`\lambda 1083`$ line and weighting with C11 from the Al iii fits we estimate that $`N(N\text{ii})<14.55`$.

### 3.2 Neutral hydrogen column density

Despite its high S/N and coverage of higher order Lyman series lines, obtaining $`N(H\text{i})`$ from the UVES spectrum alone posed some difficulties owing to the non-Voigt shape of the damped Ly$`\alpha `$ line. To estimate $`N(H\text{i})`$ we used instead the B&C spectrum, where the fit to Ly$`\alpha `$ was much better constrained than in the high resolution spectrum. This shows the difficulties of defining a continuum around an absorption feature that extends over several echelle orders.

We fitted the damped Ly$`\alpha `$ line with 14 individual components whose positions were fixed at the Fe ii velocity components. The H i column density was varied as to keep the $`N(Fe\text{ii})/N(H\text{i})`$ ratio constant for all components, and the $`b`$-parameter was assumed the same for all components. The best fit we obtained implies an overall column density of $`\mathrm{log}N(H\text{i})=20.67\pm 0.02`$ cm-2. Since the same value is obtained if a single fit component is assumed, we conclude that the assumption of constant Fe ii/H i is unimportant in so far as the $`N(H\text{i})`$ value is concerned. The fitted profiles are shown in Fig. 8.

We note that the fit solution still leaves underabsorbed “voids” in the blue wings of the Ly$`\beta `$ line. These are produced by low column-density Ly$`\alpha `$ systems that do not show the low-ionization species studied in this paper.

Of course the underlying assumption in our fitting procedure that Fe ii/H i remains constant over all velocity components, although expected if the gas is well mixed, must not hold if Fe ii is for instance inhomogeneously dust depleted. Next section deals with such possible variations.

## 4 Inferences from cloud-to-cloud variations

In this section we investigate cloud-to-cloud variations of the relative abundances. There are 3 possible sources of abundance variations among clouds: (1) differential depletion due to differing composition and destruction histories of dust grains; (2) different ionization conditions due to density variations and/or local changes in the ionization field (§ 5); and (3) simply differing chemical compositions.

### 4.1 No abundance gradients

Fig. 9 plots $`\mathrm{log}N_i(\mathrm{X})/N_i(Fe\text{ii})`$ versus $`\mathrm{log}N_i(Fe\text{ii})`$ for all velocity components for which column densities and associated errors are available. In each panel the values have been normalized by $`N_{\mathrm{corr}}(\mathrm{X})/N_{\mathrm{corr}}(Fe\text{ii})`$, where $`N_{\mathrm{corr}}`$ is the total column density listed for each element in Table 3. Since we were not able to obtain individual component values of $`N(Zn\text{ii})`$ we use instead $`N(Fe\text{ii})`$ as reference, bearing in mind that iron is possibly subject to differential dust depletion. C1 and C2 have been integrated and treated as a single component due to possible blending that could have affected the fits; the new value corresponds to the point at $`\mathrm{log}N(Fe\text{ii})=13.02`$.

From Fig. 9 we see that the $`N(\mathrm{X})/N(Fe\text{ii})`$ ratio remains, within measurement errors, fairly constant over almost two and half orders of magnitude in $`N(Fe\text{ii})`$. This is specially evident in Si ii, for which we have the most robust column density estimates and coverage of all Fe ii components. The same trend holds even for those ions for which we were able to obtain $`N`$ only for the strongest components (Ni ii, N i, P ii, S ii, Ar i). The major discrepancies with respect to the assumed abundances are found in O i and C ii, but these are easily explained by the saturated C3 in both transitions used for C ii, and by the overfitted C8 in O i. Al ii also shows discrepant values between weak and stronger components; this is most likely due to saturation effects affecting the fits to high column densities components (which underestimate the true $`N(Al\text{ii})`$).

From the data on Si ii we deduce that the single column densities diverge from the assumed value $`\mathrm{log}N(Si\text{ii})=\mathrm{log}N(Fe\text{ii})+0.44`$ on average by only $`0.15`$ dex. The constancy of $`N(Si\text{ii})`$ over such wide range of Fe ii column densities has important implications for this DLA. First, silicon is an $`\alpha `$-element. If the delayed injection of processed material from SN Type I and II had undergone differing histories from cloud to cloud, one would expect such variations to be reflected in the individual Si ii/Fe ii ratios; indeed, we would expect to observe as large a scatter in Si ii/Fe ii as the dispersion observed among different DLAs of the same metallicity (e.g., Lu et al. Lu (1996))). Secondly, silicon is a refractory element. Locally, the depletion of Si and Fe differ by much larger factors than what is observed here, over all ISM environments from the cool disk to the warm halo (Savage & Sembach Savage (1996)). In addition, since the gas-phase abundances of refractory elements decrease with $`\mathrm{log}H\text{i}`$ (e.g., Wakker & Mathis Wakker (2000)) – and thus with $`\mathrm{log}N(Fe\text{ii})`$ in this DLA– we would expect departures from the assumed abundances at higher values of $`\mathrm{log}N(Fe\text{ii})`$. Altogether, the small variations we observe in Fig. 9 mean both the chemical evolution in single clouds has been the same, and dust is either absent or has different properties from the ISM. Such small variations have been pointed out just for one other DLA (Prochaska & Wolfe Prochaska2 (1996)) toward Q0201+365 at $`z=2.5`$. Those authors also suggested the absence of dust could explain their observations. The novelty of our approach is that we present the first cloud-by-cloud analysis where that trend is observed (instead of integrating the column densities over particular velocity windows). From an instrumental point of view, the homogeneous behavior of these abundances imply the $`N(Fe\text{ii})`$-weighted corrections that we will apply to some of the low-ion column densities are physically acceptable.

The behavior we see in this DLA has not been observed in the local ISM, which shows a range of values of roughly $`1`$ dex (Spitzer & Fitzpatrick Spitzer (1995); Welty et al. Welty1 (1999)), nor in the LMC (Welty et al. Welty1 (1999)). It is well known that such lines of sight probe regions in different environments; therefore, it seems conceivable that this DLA encounters an homogeneous population of clouds.

Under the assumption that Fe ii tracks H i this line of sight probes the interval $`\mathrm{log}N(H\text{i})1820`$ cm-2, which, if in our Galaxy, would mean probing a mixture of halo and disk clouds (Wakker & Mathis Wakker (2000); and references therein). Wakker & Mathis (Wakker (2000)) find a remarkable anticorrelation of the metal abundance with H i in the local ISM, which they interpret as the environmental dependence expected if both disk (with low gas-phase abundances) and halo (with high gas-phase abundances) are plotted together as a function of H i. Our observations do not fit into this scenario, as – again assuming $`N(H\text{i})N(Fe\text{ii})`$ – they show no dependence of abundances with respect to H i column density (Fig. 10). Two possible explanations may be responsible for this discrepancy: (1) Fe ii simply does not track H i but is an increasing function of H i, in which case, a higher density cloud would have a higher metallicity; or (2) this particular line of sight probes just one type of environment. The above discussion seems to support the second alternative.

## 5 Ionization

If high redshift DLAs occur in dense protogalactic clumps, then one can expect a soft ionizing radiation field which is not able to penetrate gas with $`N(H\text{i})>10^{20}`$ cm-2. As a result, all elements but those having their first ionization potential $`>13.6`$ eV: N i, O i, Ar i, should be primarily in singly ionized form, with little gas in ionized form (e.g., Viegas Viegas (1995)). More recently, several authors have investigated the systematic effects of ionized gas on the determination of metal abundances from the low-ionization species (e.g., Howk & Sembach Howk (1999); Vladilo et al. Vladilo (2001)).

In the present DLA the moderately low overall $`N(H\text{i})`$, and the detection of some doubly ionized ions such as Al iii having similar velocity profiles as their singly ionized counterparts, both suggest ionization may play a role in bringing part of the low ions to a higher, in principle undetectable ionization level. In this section we argue against ionization affecting considerably our measurements of the metal abundances even for the single clouds.

### 5.1 Self-shielding

Our observations suggest that H i is distributed over all of the single components. Under the assumption $`\mathrm{log}N(H\text{i})=\mathrm{log}N(Fe\text{ii})+5.75`$ for each of the velocity components, the median value of an H i cloud is derived to be $`\mathrm{log}N(H\text{i})=19.26`$, i.e., at least half of the clouds would be considerably ionized, were they part of a classical Lyman-limit system (e.g., Lopez et al. Lopez (1999)). Given the complexity of the present system, however, it is also conceivable that self-shielding does keep the single clouds in a more neutral environment. Evidence for this comes from the constancy of the low-ion/Fe ii ratio down to the lowest measurable Fe ii column densities. Fig. 10 shows the ratio $`\mathrm{log}(N(Si\text{ii})/N(Fe\text{ii}))`$ as a function of $`\mathrm{log}N(Fe\text{ii})`$. The Si ii/Fe ii ratio of integrated column densities in DLAs should be only mildly influenced by ionization effects (Vladilo et al. Vladilo (2001); Howk & Sembach Howk (1999)). We interpret the datapoints in Fig. 10 as an extension of that property down to much lower densities, thus suggesting self-shielding in the present system.

Next, to address the question of how ionization may influence the single-cloud abundance ratios, we compare the fit results for the single components of Fe ii and Al iii, which are good proxies for the low and medium-ionization gas, respectively. As for Fe ii, Al iii column densities also fulfill the requirement of coming from unsaturated and resolved lines. Ideally, we would have liked to have used Al ii instead of Fe ii to disentangle the effect of possible variations in the relative abundances (despite the results shown in § 4), but the former ion has only one transition ($`\lambda 1670`$), and some of the components are saturated. Triply ionized species such as C iv or Si iv were not used as their profiles indicate they occur in distinct volumes. Fig. 11 shows the ratio $`\mathrm{log}(N(Fe\text{ii})/N(Al\text{iii}))`$ as a function of $`\mathrm{log}N(Fe\text{ii})`$ for all 12 velocity components for which column densities and associated errors are available. The two lower limits correspond to the non-detection of Al iii at the positions of C13 and C14, and the dashed line mark the Fe ii/Al iii calculated from the overall column densities. We use the same y-scale as in Fig. 10 to permit a comparison of the scatters. Contrary to the constant behavior of the N(Si ii)/N(Fe ii) ratio, there is a larger scatter in $`N(Fe\text{ii})/N(Al\text{iii})`$ of about 1 dex over the same interval in $`N(Fe\text{ii})`$. Given the constancy in the low-ion abundances, such scatter could be attributed to cloud-to-cloud variations in ionization.

Our $`\mathrm{log}N(Al\text{iii})/N(Al\text{ii})0.6`$ is typical of DLAs with similar $`N(H\text{i})`$, according to the compilation by Vladilo et al. (Vladilo (2001)). These authors find an empirical anti-correlation between $`N(H\text{i})`$ and the Al iii/Al ii ratio, suggesting that H i indirectly traces the level of ionization in the gas. A similar trend is not seen in our Fe ii/Al iii plot. We find that the two points representing the most ionized clouds (i.e., lowest Fe ii/Al iii ratio; corresponding to C1+C2 and C7) contribute only 17 % (6%) of the total Al iii (Fe ii) column density. Ignoring those two datapoints we find that Fe ii/Al iii remains in this DLA quite independent of Fe ii (and thus of H i in the single clouds). A linear regression analysis gives a slope of only $`m=0.16\pm 0.32`$, which, should Fe ii indeed track H i over these interval, is small when compared with the $`m=0.81\pm \pm 0.015`$ found by Vladilo et al. (Vladilo (2001)) for larger H i column densities. The conclusion again is that this system, although made of several individual absorbing clouds, can be treated as one single cloud exposed to the ionizing radiation. This is basically the same conclusion reached by Prochaska & Wolfe (Prochaska2 (1996)) for the $`z=2.5`$ DLA. Let us also note that a constancy of the ionization conditions is similarly observed in the local ISM through the small variability of ionized helium (Wolff, Koester & Lallement Wolff (1999)).

### 5.2 Spatial distribution of Al iii

Fig. 11 also provides hints on the spatial distribution of Al iii. In view of our finding that the low-ionization species in this DLA track each other independently of gas density, if Al iii is co-spatial with the low ions then we expect that they anticorrelate, because a decrease in H i (low-ion) column density implies an increase in ionization. Since little dependence of Al iii with Fe ii is observed, we conclude that a scenario of partially ionized clouds (i.e., with Al iii occuring in partially ionized interfaces bordering the neutral gas) is more plausible than co-spatiallity.

In such a scenario one would expect from the line widths to observe no correlation between Al iii and Al ii $`b`$-values. We obtain a dispersion of $`\sigma \sqrt{(\mathrm{\Delta }b)^2/11}=1.94`$ km s-1 for Al ii–Al iii, while $`\sigma =1.1`$ km s-1 for Al ii–Fe ii, and $`0.7`$ km s-1 for Fe ii–Si ii (including C13 and C14). If real, these differing dispersions are indicative of absorption occuring in different regions.222 We stress, however, the potentially large systematic errors that may be present in our measurements of $`b`$, especially due to saturation of some Al ii components.

### 5.3 Photoionization models

We performed photoionization models for this DLA and compared the simulations with the observed Si ii, Fe ii, and Al iii in order to establish the level of ionization of each single cloud. The simulation outputs were also used to quantify possible (if any) ionization corrections to the observed relative abundances.

We used the code CLOUDY (Ferland Ferland (1993)) to model the ionization state for a range of H i column densities and ionization parameters ($`U`$, the ratio of ionizing photons to the total hydrogen number density). CLOUDY gives column densities for a given set of the input parameters: gas metallicity $`Z`$, gas number density $`n_\mathrm{H}`$, and ionizing spectrum $`J_\nu `$. Of these, shape and intensity of the ionizing spectrum introduce the major uncertainties. We adopted the universal background ionizing field due to QSOs and galaxies calculated in Haardt & Madau (Haardt (1996); an electronic, updated version was kindly provided by the authors). This is the appropriate choice for our treatment of the individual clouds as isolated objects embedded in an intergalactic field. Stars as sources of photoionization are not likely to influence significantly these models: again referring to the work by Prochaska & Wolfe (Prochaska2 (1996)) the exponential fall off of a blackbody ionizing flux has a similar shape as the He ii Lyman-break in the Haardt & Madau spectrum for the energy range of interest. The intensity of the field at the Lyman-limit was taken to be $`J_{912}=4.4\times 10^{22}`$ ergs s-1cm-2Hz-1sr-1 at $`z=2.3`$ (consistent with results obtained using the proximity effect; e.g., Giallongo Giallongo (1996)). The geometry was plane parallel and illuminated from one side. The metallicity was fixed at \[Fe/H\] $`=1.0`$, and the relative abundances of interest \[Si/Fe\] and \[Fe/Al\] were fixed at the values listed in Table 3. For fixed values of $`U`$ we ran a grid of models by varying $`N(H\text{i})`$ in steps of $`0.2`$ dex. Each integration was stopped when the required H i value was reached. (i.e., the cloud size scales with H i in order to keep $`U`$ constant). The smoothed curves in figures 10 and 11 show the results of the simulations for $`\mathrm{log}U=2`$, $`3`$, and $`4`$.

Fig. 10 shows that in our model even at the very lowest ionization parameters a mild but significant increase of the Si ii/Fe ii ratio with decreasing Fe ii is expected. Compared to the model, the datapoints confirm that the behavior of the observed $`Si\text{ii}/Fe\text{ii}`$ ratio is not consistent with isolated clouds, and that some degree of self-shielding must exist in order to keep this ratio constant over the observed range of Fe ii column densities. We warn that the H i value at which the iterations are stopped to obtain a given (Fe ii,Si ii/Fe ii)-point is systematically lower than our assumed $`\mathrm{log}N(H\text{i})=\mathrm{log}N(Fe\text{ii})+5.75`$. In other words, Fe ii is not tied to H i in this model for $`N(Fe\text{ii})<13.5`$ cm-2. The simulated Fe ii values are indeed obtained for systematically higher H i column densities than the assumed from the observed Fe ii values. However, this fact and the use in the simulations of slabs illuminated just from one side both support the notion that i) the individual clouds are embedded in a more neutral environment than if they were isolated, and ii) the ionization is homogeneous over the DLA components. In consequence, treating the system as one single cloud with $`N(H\text{i})=20.67`$ should not introduce major systematic errors.

Concerning Al iii, note that all but 2 points in Fig. 11 lie within the contours $`\mathrm{log}U=4`$ and $`3`$, while all Si ii measurements are explained with $`\mathrm{log}U<4`$. If different ionization conditions explain the data, then Al iii and Si ii (and so Al ii) may occur in spatially distinct regions.

Finally, Fig. 12 considers the same photoionization models but with $`N(H\text{i})=20.67`$. The vertical axis now shows the ratio $`N(\mathrm{X}^i)/N(\mathrm{X}^{i+1})`$ for the variety of observed ions as a function of U. The thick sections of the curves indicate those regions allowed by the $`3\sigma `$ upper limits on the doubly ionized species Al iii, Fe iii and S iii. The softer limit is set by Al ii/Al iii, but recall that Al ii is probably underestimated due to mild saturation. Assuming $`\mathrm{log}U3.8`$, this plot shows the approximation \[X/H\] $``$ \[X+/H0\] should be accurate for the present DLA to within $`5`$ %, with the exception of Ar, which is still primarily ionized in this regime.

## 6 Metal abundances

### 6.1 Corrections to column densities

For all ions but Fe ii and Si ii, we were not able to determine column densities of each of the 14 velocity components. For weak transitions in noisy parts of the spectrum only the stronger components could be fitted, while for strong transitions only the weakest components where accessible due to blending of the saturated components. For those incomplete fits we applied a correction by scaling the summed column densities according to the results for Fe ii. We have shown in the previous section that this approach is acceptable in this DLA. The corrected values are listed in the column labeled $`N_{\mathrm{corr}}`$ of Table 3. Since for weak transitions the strongest components dominate the total column densities, the corrections are small and should not introduce large uncertainties; on the contrary, intrinsically strong transitions with only weak components fitted (O i and C ii), should suffer from systematic errors, and must thus be taken cautiously.

### 6.2 Metallicity of the DLA gas

Taking the non-refractory element Zn as reference, the gas metallicity for this DLA turns out to be \[Zn/H\] $`=1.10\pm 0.05`$ (or $`Z_{\mathrm{DLA}}=1/12Z_{\mathrm{}}`$). A slightly lower metallicity is obtained if Fe is used, \[Fe/H\] $`=1.26\pm 0.02`$. For $`z=2`$, this Fe abundance is quite at the upper end of current samples of DLA systems (Prochaska & Wolfe Prochaska1 (2000)), while the Zn value is much more representative of all measurements at $`z22.5`$ ((Pettini et al. Pettini (1997); Vladilo et al. Vladilo2 (2001)). This is because \[Zn/Fe\] $`>0`$ for the majority of the DLAs.

### 6.3 Dust content

The key property of the present DLA system is its apparent lack of dust. We measure \[Zn/Cr\] $`=0.01\pm 0.05`$, a value which is consistent with no dust-depletion, since Cr, unlike Zn, is depleted in the local ISM. Another check for the presence of dust is provided by the refractory elements Fe and Ni. The Zn/Fe and Zn/Ni ratios are expected to be supersolar if dust is present, because Fe and Ni are two of the most heavily depleted elements in the ISM. In addition, \[Zn/Fe\] $`0`$ for metal-poor stars of a wide range of metallicities (Gratton & Sneden Gratton (1991); Sneden, Gratton & Crocker Sneden (1991); but see the discussion in Prochaska et al. Prochaska4 (2000)), so any departure from this value could be attributed to dust depletion. We observe \[Zn/Fe\] $`=+0.16\pm 0.03`$ and \[Zn/Ni\] $`=0.18\pm 0.06`$. While the Zn/Fe ratio implies a certain degree of dust depletion, the Zn/Ni abundance gives further evidence in support of dust-free gas.  Altogether these 3 ratios indicate that the fraction of atoms that are missing from the gas phase is insignificant.

Nevertheless, these inferences all rely on the goodness of our zinc abundance, which was not determined by profile fitting of the velocity components. An independent argument against significant dust depletion comes from the ratio \[S/Si\] $`0.06\pm 0.03`$. Since (a) these two elements have a common nucleosynthetic origin, and (b) silicon is normally depleted in the ISM while sulfur is not (Savage & Sembach Savage (1996)), one expects this ratio to be supersolar if dust were present in this DLA, which is not observed.

In conclusion, it seems that we have found a new DLA lacking dust, the second such instance at high redshift. Other cases of DLAs with small dust content are reported in Molaro et al. (Molaro1 (2000)), in Pettini et al. (Pettini2 (2000)), and in Prochaska & Wolfe (Prochaska3 (1999)). From the Pettini et al. compilation the DLA with lowest dust content has \[Zn/Cr\] $`=0.13\pm 0.20`$ at $`z=0.86`$ toward Q0454+039, and from Prochaska & Wolfe, \[Zn/Cr\] $`=0.20\pm 0.12`$ in the $`z=1.920`$ DLA toward Q2206-199. In both cases, there appears to be some degree of dust depletion, at least stronger than in the present DLA. A more similar case as ours was reported by Molaro et al. (Molaro1 (2000)) who measured \[Zn/Cr\] $`=0.06\pm 0.08`$ at $`z=3.39`$ toward Q0000-2620. This shows that the present system is uncommon among the few documented low-dust DLAs.

### 6.4 $`\alpha `$/Fe-peak ratios

A diagnostic tool to probe the evolutionary stage in DLAs is the abundance ratio of $`\alpha `$ to Fe-peak elements. Due to their different origins, $`\alpha `$-capture elements being produced mainly in SNe Type II, and Fe-peak elements being delivered to the ISM by SNe Type Ia on a longer timescale, the $`\alpha `$ to Fe-peak ratio is expected to rapidly attain supersolar values in a chemically less mature system. Moreover, this ratio should decrease as cosmic time proceeds, unless possible environmental factors might influence it (Ellison & Lopez Ellison (2001)). In our Galaxy the values range from $`<[\alpha /\mathrm{Fe}]>`$ $`+0.5`$ in metal-poor stars down to $`0`$, as metallicity increases.

All $`\alpha `$ elements we measure in the DLA toward HE 2243$``$6031 follow each other quite closely (with the exception of Ar; see below). We obtain \[Si/Zn\] $`=+0.24\pm 0.05`$, and \[S/Zn\] $`=+0.18\pm 0.04`$ (recall that silicon is normally depleted in the local ISM; Savage & Sembach Savage (1996)). Even the $`\alpha `$-element oxygen, despite the larger measurement errors, fits well into the trend with \[O/Zn\] $`=+0.29\pm 0.22`$. On average, we obtain for Si and S $`<[\alpha /\mathrm{Zn}]>=+0.21\pm 0.06`$ (simple mean). This value is in good agreement with measurements of metal-poor halo stars (e.g., Ryan et al. Ryan (1996) and references therein) or the thick disk (Prochaska et al. Prochaska4 (2000)) at \[Fe/H\] $`=1.2`$, but is high when compared with the current sample of DLAs. The conclusion is: like in Galactic metal-poor stars, these ratios probably reflect the temporal delay between SNe of Type I and II. They show that, despite its high metallicity, this DLA gas has not undergone significant metal pollution from SN Type I yields. We re-emphasize that these ratios are likely unbiased from dust-depletion effects. Centurion et al. (Centurion (2000)) have used the \[S/Zn\] ratio to investigate a possible evolution of the $`\alpha `$/Fe-peak abundance ratio that is unbiased from dust effects. Although their sample is small (6 measurements) those authors do find a decrease in that ratio with increasing metallicity, which is what one would expect if metallicity keeps track of chemical evolution in DLAs and if the objects conform a chemically homogeneous sample. Our measure of \[S/Zn\] does not fit in that trend, being too high for \[Zn/H\] $`=1.1`$, but it much better conforms to the enhancement observed in metal-poor stars.

### 6.5 The odd-even effect

The odd-even effect, that is, the underabundance of odd-$`Z`$ elements relative to even-$`Z`$ elements of the same nucleosynthetic origin, is another known property of Halo stars. Since Fe is more prone to dust depletion than Mn in a variety of ISM environments, this ratio can be used as a discriminant between dust depletion and pure SN Type II enrichment. In fact, \[Mn/Fe\] averages $`0.3`$ in Galactic Halo stars at \[Fe/H\] $`=1`$ (Ryan et al. Ryan (1996)), whilst \[Mn/Fe\] $`+0.3`$ in the local ISM (Savage & Sembach Savage (1996)).

The Fe ii-corrected value of our $`3\sigma `$ upper limit for the non-detection of Mn ii yields \[Mn/Fe\] $`<0.51`$, which confirms the odd-even effect, is in accordance with the sample of Halo metal-poor stars in Ryan et al. (Ryan (1996)), but is quite low for that of thick disk stars in Prochaska et al. (Prochaska4 (2000)). The low value also conforms with previous DLA measurements (Lu et al. Lu (1996)).

We do not observe the same effect for P, the other odd-element observed in this DLA. We obtain \[P/Si\] $`=+0.10\pm 0.09`$ at \[P/H\] $`=0.8`$. This value is high if compared with \[P/Si\] $`=0.40`$ at \[P/H\] $`=2.3`$ (Molaro et al. Molaro (2001)) and \[P/Si\] $`=0.30\pm 0.09`$ at \[P/H\] $`=1.2`$ (Outram et al. Outram (1999)).

### 6.6 Underabundant Ar i

Our column density estimate for Ar i implies a significant underabundance of Ar with respect to other $`\alpha `$-chain elements. This corresponds to $`0.5`$ dex of what is expected if argon is overabundant by $`+0.2`$ dex (assuming argon tracks $`\alpha `$ elements \[Timmes, Woosley & Weaver Timmes (1995)\]). In the local ISM, Ar i is rarely found or is significantly below its cosmic abundance relative to H i (Sofia & Jenkins Sofia (1998); Jenkins et al. Jenkins (2000)). As suggested by Sofia & Jenkins in order to explain their finding of a low Ar i abundance in the local ISM, the underabundance can be explained if a significant part of the argon is ionized. In this DLA that fraction is $`60`$ %. Moreover, since Ar i has a much larger ionization cross section than H i, our measurement can also be explained in terms of moderately ionized gas, which brings Ar i into Ar ii in higher proportions than for the rest of the low-ions. The alternative explanation, a low Ar i abundance caused by dust depletion, is in our DLA ruled out given (1) the scarce propensity of Ar to dust-grain incorporation (Meyer et al. Meyer (1991); Sofia & Jenkins Sofia (1998)) and (2) the negligible amount of dust (6.3). Thus, these conditions resemble very much those in the local ISM described by Sofia & Jenkins where argon is systematically found underabundant, and confirm their suggestion that ionization, and not dust, is the responsible agent.

### 6.7 Nitrogen

Chemical evolution models predict the production of nitrogen in the CNO cycle in stellar interiors to have two different origins. While primary N is produced in intermediate-mass stars (3 – 8 $`M_{\mathrm{}}`$) without dependence on the initial metal content, secondary nitrogen is also produced in stars of any mass but with an initial supply of heavy elements. The net result is that the N/O ratio should be independent of O/H for primary N while N/O $``$ O/H for secondary N. In addition, since intermediate-mass stars evolve on much larger scales than the short-lived progenitors of SN Type II – responsible for the bulk of $`\alpha `$-elements like O, Si, S, etc. –, there should be a delay of $`5\times 10^8`$ yr between the release of primary N and that of O (Lu et al. Lu1 (1998) and references therein).

We measure \[N/O\] $`=1.03`$ at \[O/H\] $`=0.81`$. However, despite the accurate N abundance we were able to obtain, our oxygen estimate is rather uncertain.333Fig. 9 shows C8 departures from the trend followed by the other 4 components. This might be due to an overestimate of C8 (which dominates the total column density). Excluding C8 and re-scaling, however, implies lowering the total column density by $`0.7`$ dex. Although such value conforms better with the rest of the components, it translates into \[O/H\] $`=1.5`$, so the underabundance of N over O would become only \[N/O\] $`=0.2`$. Relying rather on the better estimates for Si and S, we find \[N/Si\] $`=0.98`$ at \[Si/H\] $`=0.86`$ (and basically the same values are found if S  is used). Referring to Fig. 1 in Lu et al. (Lu1 (1998)), this value is at the very low end of DLA measurements at this Si abundance; moreover, according to the fits to Galactic measurements, it is better explained as secondary nitrogen444But let us note that N ii could contribute with $`0.20`$ dex to the nitrogen abundance, thus bringing \[N/Si\] more in line with the Lu et al. values.. This is in agreement with the high $`\alpha `$/Fe-peak ratios we observe because both effects can be explained by an early phase in this protogalaxy’s chemical evolution, observed at a time when the release of Fe-peak elements and primary nitrogen has not yet begun.

## 7 Summary of the conclusions

1. We have studied in detail a $`z=2.33`$ DLA observed at high resolution toward the bright QSO HE 2243$``$6031. We have been able to detect and obtain reliable column densities for H i, N i, Si ii, P ii, S ii, Ar i, Cr ii, Fe ii, Ni ii, Zn ii, and Al iii; less accurate determinations for C ii, Al ii and O i; and to put upper limits on Mn ii, S iii, Fe iii, P iii, and N ii. The metallicity of this system is \[Zn/H\] $`=1.1`$.
2. By comparing column densities of single velocity components we have performed the first cloud-by-cloud analysis in a high-$`z`$ DLA so far. We find a tight correlation of all low-ionization species with respect to Fe ii extending over 2.5 orders of magnitude in $`N(Fe\text{ii})`$. We interpret this trend as a lack of significant changes in the dust, ionization, and intrinsic chemical-enrichment effects. If Fe ii tracks H i, such homogeneity differs from what is observed in the local ISM, meaning that our line of sight probably samples clouds located in one single environment.
3. We find a probable scenario for the occurrence of Al iii is ionized shells surrounding neutral cores of Al ii. Photoionization models show that the ratio of Al iii to low ions can be modeled with a wider range of higher ionization parameters than the ratio among low ions. Independently of this, the disparate line width differences between Al ii and Al iii also support the notion of separate Al ii and Al iii regions.
4. Using different abundance ratios between ISM refractory and non-refractory elements we have shown that the observations are compatible with negligible amounts of dust in this protogalaxy, although the measured \[Zn/Fe\] $`=+0.16\pm 0.03`$ may indicate a small degree of dust depletion.
5. We observe – unbiased by dust depletion effects – an enhancement of the $`\alpha `$/Fe-peak ratios of $`+0.2`$ dex for various elements, a marked odd-even effect in Mn, and a strong underabundance of N relative to Si and S, \[N/Si,S\] $`=1`$ at \[Si/H\] $`=0.86`$. All of these ratios support an environment that is in an early evolutionary stage, where the onset star formation has begun shortly before the DLA was observed. The absence of dust is in agreement with this notion, if there has not been time enough to form dust grains.
6. The observed underabundance of Ar i by $`0.5`$ dex is easily explained by a significant fraction of argon residing in a higher ionization phase. This situation resembles the conditions in the local ISM. Our ionization simulations demonstrate, however, such an ionized phase does not contain a considerable part of the other elements, so our abundance measurements are not affected by undetected ionized atoms.

### 7.1 Final comment

The observed uniformity of metal abundances across the line of sight to HE 2243$``$6031 has important implications for current models of galaxy formation. As mentioned before, it suggests the individual components are embedded in the same environment, which, in addition, must have undergone significant mixing. If DLAs are formed from the accretion of protogalactic small clumps (e.g., Haehnelt et al. Haehnelt (1998)) that have formed independently, then one expects each protogalaxy to contribute with gas of different dust and nucleosynthetic compositions. The homogeneity of metal abundances observed in this DLA is difficult to explain in this scenario, unless the mixing of metals has occurred in short enough timescales to counteract the enhancement of any original differences (e.g., through the metal release by SN Type I; see the discussion in Prochaska & Wolfe Prochaska5 (2002)). The chemically young gas we apparently observe might be a consequence of such a rapid process.

Although the statistical weight of our result on single-cloud abundances may not be sufficient to draw any general conclusions on the physical nature of DLAs, the observed uniformity in HE 2243$``$6031 seems to be a quite common property among DLAs (Prochaska et al. Prochaska6 (2001)).

On larger physical scales, also important for constraining the models is to investigate abundance variations among systems. There are few known ‘DLA groups’ along single lines of sight where the redshift separation is small enough to support a physical connection between the DLAs, but big enough to characterize the systems as separate entities (i.e., a few $`\mathrm{1\hspace{0.17em}000}`$ km s-1; Lopez et al. Lopez1 (2001); Ellison & Lopez Ellison (2001)). In the only such groups so far investigated at high resolution, Ellison & Lopez (Ellison (2001)) find very similar relative abundances among the DLA members, despite their different metallicities. While it is not clear whether such DLA groups later merge, environment may be an effective catalyzer of metal enrichment in DLAs at early evolutionary stages.

###### Acknowledgements.

We warmly thank Sara Ellison for her valuable comments on an earlier version of the paper; Ulysses Sofia for an informative discussion; Piero Madau and Francesco Haardt for having made available an updated version of their spectrum to us; Scott Burles for having taken the B&C spectrum of HE 2243$``$6031, and G. Vladilo, the referee, for a careful reading of the manuscript. JXP was supported by NASA through a Hubble Fellowship grant HF-01142.01-A awarded by STScI, and SL acknowledges financial support by FONDECYT grant N$`{}_{}{}^{\mathrm{o}}\mathrm{3\hspace{0.17em}000\hspace{0.17em}001}`$ and by the Deutsche Zentralstelle für Arbeitsvermittlung.


format: text

text:
# Electrostatics in Periodic Slab Geometries I

## I Introduction

The calculation of long range interactions due to Coulomb, gravitational, or dipolar particles is of broad interest from the astrophysics, the biophysics up to the solid state community. These interactions present a formidable challenge even to modern computers. Sophisticated methods such as fast multipole methods, tree code algorithms, Poisson grid solvers, or Ewald mesh methods, bring the complexity of N interacting particles down to an almost linear scaling for three-dimensional periodic systems. Often, however, one is interested in slab-like systems which are only periodic in two space dimensions and finite in the third, for example in problems involving electrolyte solutions between charged surfaces, proteins near charged membranes, thin films of ferrofluids, Wigner crystals, charged films, membranes, solid surfaces decorated with dipoles etc.

For such systems Ewald based formulas are only slowly convergent, have mostly $`𝒪(N^2)`$ scalings and no “a priori” error estimates exist widmann97a . Fast Ewald based methods have been recently put forward in kawata01b , and a non-Ewald method has been put forward in Ref.lekner91a  that is based on a resummation of the force sum. However, these methods are hampered due to non-controllable errors and an $`𝒪(N^2)`$ scaling respectively. Recently we proposed a new method called MMM2Darnold01b ; arnold01c  which has an $`𝒪(N^{5/3})`$ complexity and full error control that is based on a convergence factor approach similar to MMMsperb98a . However, this will still only allow simulations including up to a few thousand charges. There have been early attempts to use a 3D Ewald sum for these slab problems. The main idea is to fill only parts of the simulation box with charges and to leave some space empty, in an attempt to decouple the interactions in the third dimension shelley96a ; spohr97a ; yeh99a . Since each image layer is globally neutral, one hopes that their interactions decay as they become more and more distant, i.e. as the size of the gap is increased. In this way one could make use of any advanced 3D Ewald implementation, see also Ref. minary02a  for a variant of this idea.

In this paper we will follow the last suggestion and derive a term, called electrostatic layer correction (ELC), which subtracts the interactions due to the unwanted layers. The combination of that term with any three dimensional summation method with slab–wise summation order will yield the exact electrostatic energy. Since the change in the summation order is done by adding a very simple term, any three dimensional summation method with the standard spherical summation order can be used. The new term can be evaluated easily in a time linear in the number of charges, hence the whole method scales like the underlying standard summation method. We develop also an error formula for the maximal pairwise error in the energy and forces of the layer correction term, hence the precision of this method can be tuned to any desired value, when used in conjunction with other error estimates for the standard summation methodkolafa92a ; deserno98b . In the first section we will recapitulate the way how to correct the summation order via a modified dipole term. In the second section we will derive the layer correction term, and develop in the following section error estimates for its value. The applicability of our method will be demonstrated by a numerical analysis in the following section, and we end with our conclusion.

## II Changing the summation order

We consider a system of $`N`$ particles with charges $`q_i`$ and positions $`p_i=(x_i,y_i,z_i)`$ that reside in a box of edges $`L\times L\times h`$, where $`h=\mathrm{max}_{i,j}|z_iz_j|`$ is the maximal $`z`$–distance of two particles. The basic idea is to expand this slab system in the non-periodic $`z`$–coordinate to a system with periodicity in all three dimensions. More precisely, the original box of size $`L\times L\times h`$ is placed inside a box of size $`L\times L\times L_z`$ where $`L_z>>h`$ sufficiently large. Then this box is replicated periodically in all three dimensions. The result is a three-dimensional periodic system with empty space regions (“gaps”) of height $`\delta :=L_zh`$ (see Fig. 1). $`\delta `$ will be called gap size in the following.

Since the electrostatic potential is only finite if the total system is charge neutral, the additional image layers (those layers above or below the original slab system) are charge neutral, too. Now let us consider the $`n^{th}`$ image layer which has an offset of $`nL_z`$ to the original layer. If $`nL_z`$ is large enough, each particle of charge $`q_j`$ at position $`(x_j,y_j,z_j+nL_z)`$ and its replicas in the $`x,y`$-plane can be viewed as constituting a homogeneous charged sheet of charge density $`\sigma _j=\frac{q_j}{L^2}`$. The potential of such a charged sheet at distance $`z`$ is $`2\pi \sigma _j|z|`$. Now we consider the contribution from a pair of image layers located at $`\pm nL_z`$, $`n>0`$ to the energy of a charge $`q_i`$ at position $`(x_i,y_i,z_i)`$ in the central layer. Since $`|z_jz_i|<nL_z`$, we have $`|z_jz_i+nL_z|=nL_z+z_jz_i`$ and $`|z_jz_inL_z|=nL_zz_j+z_i`$, and hence the interaction energy from those two image layers with the charge $`q_i`$ vanishes by charge neutrality:

$$2\pi q_i\underset{j=1}{\overset{N}{}}\sigma _j(|z_jz_i+nL_z|+|z_jz_inL_z|)=4\pi q_inL_z\underset{j=1}{\overset{N}{}}\sigma _j=0\text{.}$$
(1)
The only errors occurring are those coming from the approximation of assuming homogeneously charged, infinite sheets instead of discrete charges. This assumption should become better when increasing the distance $`nL_z`$ from the central layer.

However, in a naive implementation, even large gap sizes will result in large errors yeh99a . This is due to the order of summation for the three dimensional Coulomb sum, which is spherical by convention. This order implies that with increasing shell cutoff $`S`$ the number of image shells grows faster than the number of shells of the primary layer, namely $`𝒪(S^3)`$ versus $`𝒪(S^2)`$ (see Fig. 2(a)). In other words, we include the unwanted terms faster than the actually wanted terms. Also the image layers are not really infinite charged sheets but are truncated due to the cut-off. Yeh and Berkowitzyeh99a  already suggested that this problem can be solved by changing the order of summation. Smith has shown that by adding to the Coulomb energy the term

$$E_c=2\pi M_z^2\frac{2\pi M^2}{3}\text{,}$$
(2)
where $`M=q_ip_i`$ is the total dipole moment, one obtains the result of a slab–wise summation instead of the spherical limit smith81a . Slab–wise summation refers to the sum $`_{|n|0}E_l(n)`$, where $`E_l(n)`$ denotes the energy, calculated in spherical summation order, resulting from the image layer with shift $`nL_z`$ in the $`z`$–coordinate. Technically this is the order where we first treat the original layer and then add the image layers grouped in symmetrical pairs (see Fig. 2(b)). Obviously this summation order fits much better to the charged sheet argumentation given above. Although this is a major change in the summation order, the difference given by Eq. (2) is a very simple term. In fact, Smith shows that changes of the summation order always result in a difference that depends only on the total dipole moment.

Applying this slab–wise summation order, Yeh and Berkowitz showed that a gap size of at least $`h`$ is normally sufficient to obtain an moderately accurate result. Therefore the result of a standard three dimensional summation method plus the shape–dependent term given by Eq. (2), which we refer to as a *slab–wise method*, can be used to obtain a good approximation to the result for the slab geometry with the same computational effort as for the underlying three dimensional summation method (no matter if a simple or sophisticated method is used). One drawback is that no theoretical estimates exist for the error introduced by the image layers. Therefore one might be forced to use even larger gaps to assure that no artifacts are produced by the image layers. One simple deducible artifact is that the pairwise error will be position dependant. Particles in the middle of the slab will see no effect of the image layers due to symmetry, and particles near the surface will encounter for the same reason the largest errors, which is definitely an unwanted feature for studying surface effects. Therefore averaging error measures like the commonly used RMS force error should not be applied without additional checks for the particles near the surfaces.

The other drawback is that normally the box now will have a significantly larger $`L_z/L`$. But at least for Ewald type methods the computation time is proportional to this fraction. This is easy to see as the number of $`k`$–space vectors in the $`z`$ direction must be proportional to $`L_z`$ to maintain a fixed resolution and therefore error. It is verified experimentally that a gap of at least $`h`$ is needed. For a cubic system $`h=L`$ therefore the computation time at least doubles.

Nevertheless because of the bad scaling of the known methods for slab geometries like the one by Parry parry75a ; parry76a  ($`𝒪(N^2)`$) or MMM2D arnold01b ; arnold01c  ($`𝒪(N^{5/3})`$), for particle numbers above $`N1000`$ using slab–wise methods is a great improvement.

## III The electrostatic layer correction term

We will now derive a term that allows to calculate the *exact* contribution of the image layers very efficiently, which we will call the electrostatic layer correction (ELC) in the following. For the following analysis there is no special restriction on $`h`$ except for $`h<L_z`$, which is true even if the $`L\times L\times L_z`$–box is completely filled.

The method presented here is heavily based on parts of MMM2D arnold01b . We start with a formal definition of the Coulomb energy of the slab system

$$E=\frac{1}{2}\underset{S=0}{\overset{\mathrm{}}{}}\underset{\begin{array}{c}n^2\times \{0\}\\ n_x^2+n_y^2=S\end{array}}{}\underset{i,j=1}{\overset{N}{^{}}}\frac{q_iq_j}{|p_ip_j+\mathrm{\Lambda }n|}\text{.}$$
(3)
$`\mathrm{\Lambda }=\text{diag}(L,L,L_z)`$ is a diagonal matrix describing the shape of the box. The image boxes are denoted with the vector $`n=(n_x,n_y,n_z)`$, where $`n_z=0`$ for now. The prime on the inner summation indicates the omission of the self–interaction $`i=j`$ in the primary box $`n=(0,0,0)`$ (i. e. the singular case). For the surrounding dielectric medium we assume vacuum boundary conditions.

We now expand the system to a fully three-dimensional periodic system, where $`L_z`$ determines the period in the $`z`$-coordinate as in the previous section. We can rewrite the energy as

$$E=E_s+E_c+E_{lc}\text{,}$$
(4)
where

$$E_s=\frac{1}{2}\underset{S=0}{\overset{\mathrm{}}{}}\underset{\begin{array}{c}n^3\\ n^2=S\end{array}}{}\underset{i,j=1}{\overset{N}{^{}}}\frac{q_iq_j}{|p_ip_j+\mathrm{\Lambda }n|}\text{.}$$
(5)
denotes the standard three-dimensional Coulomb–sum with spherical limit. To evaluate this expression one can use any of the efficient algorithms, starting with the classical Ewald summation up to modern methods like fast multipole methods greengard87a  or mesh based algorithmsdeserno98a . $`E_c`$ again denotes the shape–dependent term given by Eq. (2) and finally

$$E_{lc}=\frac{1}{2}\underset{\begin{array}{c}T\\ T>0\end{array}}{}\underset{n_z=\pm T}{}\underset{S=0}{\overset{\mathrm{}}{}}\underset{\begin{array}{c}n^2\times \{n_z\}\\ n^2=S\end{array}}{}\underset{i,j=1}{\overset{N}{}}\frac{q_iq_j}{|p_ip_j+\mathrm{\Lambda }n|}\text{.}$$
(6)
denotes the contribution of the image layers, for which we are going to derive a new expression in the following.

We start with the expression for the energy induced by an image layer at $`z`$–offset $`n_z0`$:

$$E_l(n_z)=\frac{1}{2}\underset{S=0}{\overset{\mathrm{}}{}}\underset{\begin{array}{c}n^2\times \{n_z\}\\ n^2=S\end{array}}{}\underset{i,j=1}{\overset{N}{}}\frac{q_iq_j}{|p_ip_j+\mathrm{\Lambda }n|}\text{.}$$
(7)
It can be shown rigorously, although this is non–trivial, that

$$E_l(n_z)=\frac{1}{2}\underset{\beta 0}{lim}\underset{n^2\times \{n_z\}}{}\underset{i,j=1}{\overset{N}{}}\frac{q_iq_je^{\beta |p_ip_j+\mathrm{\Lambda }n|}}{|p_ip_j+\mathrm{\Lambda }n|}\text{.}$$
(8)
This is a convergence factor approach with a convergence factor of $`e^{\beta |p_ip_j+\mathrm{\Lambda }n|}`$. Note that this approach is exact only for two-dimensional systems, for three-dimensional systems Eqs. (7) and (8) differ by a multiple of the dipole moment deleeuw80a ; deleeuw80b .

In arnold01a ; arnold01b  one can find a proof for this equation and an efficient way of calculating $`E_l`$ for charge neutral systems. We do not want to go through the full derivation again; it consists of the application of Poisson’s summation formula along both periodic coordinates and performing the limit $`\beta 0`$ analytically. One obtains

$$E_{lc}(n_z)=\frac{1}{2}\underset{i,j=1}{\overset{N}{}}q_iq_j\varphi (p_ip_j+\mathrm{\Lambda }n)\text{,}$$
(9)
where $`\varphi `$ is given by

$$\begin{array}{c}\varphi (x,y,z)=\frac{4}{L}\underset{\begin{array}{c}k_{}^2\\ k_x,k_y>0\end{array}}{}\frac{e^{2\pi |k_{}||z|/L}}{|k_{}|}\mathrm{cos}(2\pi k_xx/L)\mathrm{cos}(2\pi k_yy/L)+\hfill \\ \hfill \frac{2}{L}\left(\underset{k_y>0}{}\frac{e^{2\pi k_y|z|/L}}{k_x}\mathrm{cos}(2\pi k_yy/L)+\underset{k_x>0}{}\frac{e^{2\pi k_x|z|/L}}{k_y}\mathrm{cos}(2\pi k_xx/L)\right)\frac{2\pi }{L^2}|z|\text{.}\end{array}$$
(10)
$`k_{}=(k_x,k_y)`$ is a Fourier variable with integer values. $`\varphi `$ is an artificial pairwise potential that yields the total Coulomb energy and its derivative produces the pairwise forces for the periodic system.

For now we only have a formula for the contribution of one image layer, so we still have to sum over all $`n_z`$. This task can be performed analytically. The terms $`2\pi |z|/L^2`$ can be omitted since they are exactly the homogeneous sheet potential and we have seen before that this cancels out for charge neutral systems (see Eq. (1)).

The summation over $`n_z`$ of the remaining sums over $`p`$ and $`q`$ is fairly easy to perform using the geometric series (as these sums are absolutely convergent, exchanging the summation over $`n_z`$ and the summations over $`(k_x,k_y)`$ is possible). Combining the terms for $`\pm n_z`$ again we obtain

$$E_{lc}=\underset{i,j=1}{\overset{N}{}}q_iq_j\psi (p_ip_j+\mathrm{\Lambda }n)\text{,}$$
(11)
where

$$\begin{array}{cc}\hfill \psi (x,y,z)=& \frac{4}{L}\underset{k_x,k_y>0}{}\frac{\mathrm{cosh}(2\pi k_{}z_{ij}/L)}{k_{}(e^{2\pi k_{}L_z/L}1)}\mathrm{cos}(2\pi k_xx_{ij}/L)\mathrm{cos}(2\pi k_yy_{ij}/L)+\hfill \\ & \frac{2}{L}\underset{k_x>0}{}\frac{\mathrm{cosh}(2\pi k_xz_{ij}/L)\mathrm{cos}(2\pi k_xx_{ij}/L)}{k_x(e^{2\pi k_xL_z/L}1)}+\hfill \\ & \frac{2}{L}\underset{k_y>0}{}\frac{\mathrm{cosh}(2\pi k_yz_{ij}/L)\mathrm{cos}(2\pi k_yy_{ij}/L)}{k_y(e^{2\pi k_yL_z/L}1)}\text{.}\hfill \end{array}$$
(12)
The forces can be obtained from that by simple differentiation since the sums are absolutely convergent. Although the form in Eq.(12) has a much better convergence than the original form in Eq.(6), its main advantage is a linear computation time with respect to the number of particles $`N`$. To see this, the equation has to be rewritten using the addition theorems for the cosine and the hyperbolic cosine. For each $`k_{}`$ one first calculates the sixteen terms

$$\begin{array}{cc}\hfill \chi _{(c/s,c/s,c/s)}& =\underset{i=1}{\overset{N}{}}q_i\mathrm{cosh}/\mathrm{sinh}(2\pi k_{}z_i/L)\mathrm{cos}/\mathrm{sin}(2\pi k_xx_i/L)\mathrm{cos}/\mathrm{sin}(2\pi k_yy_i/L)\text{,}\hfill \\ \hfill \chi _{(x,c/s,c/s)}& =\underset{i=1}{\overset{N}{}}q_i\mathrm{cosh}/\mathrm{sinh}(2\pi k_xz_i/L)\mathrm{cos}/\mathrm{sin}(2\pi k_xx_i/L)\text{,}\hfill \\ \hfill \chi _{(y,c/s,c/s)}& =\underset{i=1}{\overset{N}{}}q_i\mathrm{cosh}/\mathrm{sinh}(2\pi k_yz_i/L)\mathrm{cos}/\mathrm{sin}(2\pi k_yy_i/L)\text{,}\hfill \end{array}$$
(13)
where the indices in the obvious way determine which of the functions cosine (hyperbolicus) or sinus (hyperbolicus) are used. Then we evaluate

$$\begin{array}{cc}\hfill E_{lc}=& \frac{4}{L}\underset{k_x,k_y>0}{}\frac{1}{(e^{2\pi k_{}L_z/L}1)k_{}}(\chi _{(ccc)}^2+\chi _{(csc)}^2+\chi _{(ccs)}^2+\chi _{(css)}^2\hfill \\ & \chi _{(scc)}^2\chi _{(ssc)}^2\chi _{(scs)}^2\chi _{(sss)}^2)+\hfill \\ & \frac{2}{L}\underset{k_x>0}{}\frac{1}{(e^{2\pi k_xL_z/L}1)k_x}\left(\chi _{(xcc)}^2+\chi _{(xcs)}^2\chi _{(xsc)}^2\chi _{(xss)}^2\right)+\hfill \\ & \frac{2}{L}\underset{k_y>0}{}\frac{1}{(e^{2\pi k_yL_z/L}1)k_y}\left(\chi _{(ycc)}^2+\chi _{(ycs)}^2\chi _{(ysc)}^2\chi _{(yss)}^2\right)\hfill \end{array}$$
(14)
Similar expansions using the same sixteen terms can also be found for the forces. Obviously this has linear computation time with respect to the number of particles, as the only summations over all the particles occur in the $`\chi _{}`$. But up to now there is still the infinite summation over $`k_{}`$. So the next task is to derive an estimate for the error induced by the replacement of the infinite sum by a finite one.

## IV Error estimates

Since $`E_{lc}`$ is written as sum over an alternative potential $`\psi `$, it is reasonable to derive an upper bound for the error from the calculation of $`\psi `$ only with a finite cutoff. From this upper bound, crude estimates for other error measures such as the RMS (root-mean-square) force error can be derived. Again these error estimates are taken from MMM2D arnold01a ; arnold01b . As we will show later, the error distribution is not uniform along the $`z`$–axis. The error gets maximal for particles near the borders. Since these particles will normally be those of special interest, the maximal pairwise error should be some magnitudes smaller than the thermal noise.

While the error bounds for MMM2D were only used to tune the algorithm, the error estimates for $`E_{lc}`$ can also be used to obtain an error bound for the slab–wise method from Ref. yeh99a , and hence one can determine “a priori” the necessary gap size to reach a preset precision. Therefore we also have to deal with small cutoffs, especially the case when no terms of $`E_{lc}`$ are added.

First we choose a cutoff $`R1`$ and then evaluate $`E_{lc}`$ only over the area

$$\begin{array}{cc}\hfill \mathrm{\Gamma }_R=& \left\{(k_x,k_y)^2\right|k_x,k_y>0,(k_x1)^2+(k_y1)^2<R^2\}\hfill \\ & \left\{(k_x,0)\times \{0\}\right|k_x<R\}\hfill \\ & \left\{(0,k_y)\{0\}\times \right|k_y<R\}\text{.}\hfill \end{array}$$
(15)
The three sets correspond to the three sums in Eq.(12). Therefore we actually evaluate

$$\begin{array}{cc}\hfill E_{lc}=& \frac{4}{L}\underset{\begin{array}{c}k_x,k_y>0,\\ (k_x1)^2+(k_y1)^2<R^2\end{array}}{}\frac{\mathrm{cosh}(2\pi k_{}z_{ij}/L)}{k_{}(e^{2\pi k_{}L_z/L}1)}\mathrm{cos}(2\pi k_xx_{ij}/L)\mathrm{cos}(2\pi k_yy_{ij}/L)+\hfill \\ & \frac{2}{L}\underset{0<k_x<R}{}\frac{\mathrm{cosh}(2\pi k_xz_{ij}/L)\mathrm{cos}(2\pi k_xx_{ij}/L)}{k_x(e^{2\pi k_xL_z/L}1)}+\hfill \\ & \frac{2}{L}\underset{0<k_y<R}{}\frac{\mathrm{cosh}(2\pi k_yz_{ij}/L)\mathrm{cos}(2\pi k_yy_{ij}/L)}{k_y(e^{2\pi k_yL_z/L}1)}\text{.}\hfill \end{array}$$
(16)
$`\mathrm{\Gamma }_R`$ may look more complicated then necessary. But this form enables us to find a rigorous upper bound for the error. An upper bound for the absolute value of the summands is

$$\begin{array}{c}\left|\frac{\mathrm{cosh}(2\pi k_{}z_{ij}/L)}{k_{}(e^{2\pi k_{}L_z/L}1)}\mathrm{cos}(2\pi k_xx_{ij}/L)\mathrm{cos}(2\pi k_yy_{ij}/L)\right|\hfill \\ \hfill e^{2\pi k_{}Lz/L}\frac{\mathrm{cosh}(2\pi k_{}z_{ij}/L)}{(1e^{2\pi k_{}L_z/L})k_{}}e^{2\pi k_{}Lz/L}\frac{\mathrm{cosh}(2\pi k_{}h/L)}{(1e^{2\pi k_{}L_z/L})k_{}}\text{.}\end{array}$$
(17)
Of course because the cosine hyperbolicus is monotonous, one could use any larger value for $`h`$. This is for example necessary in a priori estimations. Using this we find the upper bound for the maximal pairwise error in the potential by a simple approximation of the sums by integrals as

$$\tau _E:=\frac{1/2+(\pi R)^1}{e^{2\pi RL_z/L}1}\left(\frac{\mathrm{exp}(2\pi Rh/L)}{Lzh}+\frac{\mathrm{exp}(2\pi Rh/L)}{Lz+h}\right)\text{.}$$

Details can again be found in arnold01a . By an analogous derivation one can also find an upper bound on the maximal pairwise error on each of the three force components as

$$\begin{array}{cc}\hfill \tau _F:=\frac{1}{2(e^{2\pi RL_z/L}1)}(& \left(\frac{2\pi R+4}{L}+\frac{1}{L_zh}\right)\frac{\mathrm{exp}(2\pi Rh/L)}{(L_zh)}+\hfill \\ & (\frac{2\pi R+4}{L}+\frac{1}{L_z+h})\frac{\mathrm{exp}(2\pi Rh/L)}{(L_z+h)})\text{.}\hfill \end{array}$$
(18)
This is also an weaker bound for the potential. In other words, the maximal pairwise error on the forces is larger than the error in the potential. For $`R=1`$ one obtains an overall estimate of the magnitude of the contribution of the image layers, i. e. an error estimate for the slab–wise methods.

Note that Eq.(17) shows that the error in the potential or the force for a single particle will be largest if it is located near the gap, since there $`|z_{ij}|`$ will be maximal. This effect will increase with increasing $`R`$. Therefore when using the layer correction one must apply non–averaging error estimates such as our maximal pairwise error. Averaging error estimates such as the RMS force error might be misleading about the error on the particles of interest.

Moreover the error will decrease exponentially with the distance from the gap. Since the particles near the gap (i.e. the surface) are normally of special interest in simulations with slab geometry, averaging error measures like the RMS force error might be misleading about the effect of these errors.

All our error estimates show that the error drops exponential both with $`R`$ and $`L_z/L`$. The decay in $`R`$ means that it is easy to achieve high accuracies with our layer correction formula, while the decay in $`L_z/L`$ shows that slab–wise methods can achieve good accuracies without increasing $`L_z/L`$ too much.

Although we do not encourage using the RMS error measure as we explained above, we still want to give an upper bound on the average error in $`E_{lc}`$ and the RMS force error. We assume that the pairwise potentials of the different particle pairs are independent identically distributed random variables, which is true for homogenous random systems and normally a good assumption otherwise. Of course the self interaction of the particles, i. e. $`q_i^2\psi (0,0,0)`$ has to be omitted. Let $`\sigma _E`$ be the variance of this random variable, then it is easy to see that $`\sigma _E\tau _E^2`$. Using this one can show that

$$E_{lc}\underset{i}{}q_i^2\psi (0,0,0)Q^2\sqrt{\sigma _E}Q^2\tau _E\text{,}$$
(19)
where $`Q:=_iq_i`$.

Similarly we define $`\sigma _F`$ as the variance of the forces measured in the Euclidian norm. Then because of the component–wise maximal error estimate for the force, we have $`\sigma _F^23\tau _F^2`$, and one obtains

$$\sqrt{1/N\underset{i}{}|\mathrm{\Delta }F_{lc}^i|}\sqrt{3}Q^2/\sqrt{N}\tau _F\text{,}$$
(20)
where $`\mathrm{\Delta }F_{lc}^i`$ denotes the error in the layer correction force on particle $`i`$. Note that both estimates (19) and (20) are much larger than the real error as one expects $`\sigma `$ to be *much* smaller the maximal error (about $`24`$ magnitudes).

## V Numerical demonstration

In this section we want to show results from our implementation of the layer correction term (ELC). First we want to show that our maximal pairwise error bounds are correct. To this aim we place two particles randomly in a box of size $`1\times 1\times 0.8`$, so that we leave a gap of $`\delta =0.2`$ in a box of dimensions $`1\times 1\times 1`$. Fig. 2 shows the maximal potential and force error that occurred during 10000 evaluations and our estimates $`\tau _E`$ and $`\tau _F`$. Moreover we have included the result for a particle pair with a relative position of $`(0,0,h)`$, the worst case position. For such a position the error estimate is exact up to the approximation of the sum by an integral. As “exact” force we used $`R=30`$.

One can see that the maximal error estimates are always above the measured deviations. Even after 10000 random evaluations the maximum error is considerably lower than for the special pair with relative position $`(0,0,h)`$, and the error is still not very smooth. This effects are due to the Fourier representation with exponentially decaying coefficients, which makes the worst case error extremely rare. But in a real simulation the particle distribution is not necessarily homogeneous and to be on the safe side one has to deal also with the rare worst case error. Nevertheless Fig. 2 shows that the error coming from the image layers can be strictly controlled.

Next we investigate the computation times of our implementation of the ELC term. For different gap sizes $`\delta =0.05,0.1`$ and $`0.2`$ we show in Fig. 3(a) the computation time $`T_{CPU}`$ for different numbers of particles $`N`$. The systems were again consisting of uniformly randomly distributed particles of charges $`\pm 1`$. The maximal pairwise error was fixed to be $`10^4`$. The times are averages over 10 runs on a Compaq XP1000. The computation time for the same system consisting of 1000 charges using P<sup>3</sup>M is $`330ms`$ for a typical RMS force error of $`10^4`$. Therefore even the small gap of $`0.05`$ gives just the same order of computation time. For the normally used gap sizes of $`0.2L_z`$ the computation time is negligible compared to P<sup>3</sup>M. In Fig. 3(b) we show the RMS force errors that occurred. One can see the predicted $`Q^2/\sqrt{N}`$ behavior. We also show the theoretical upper bound $`\sqrt{3}Q^2/\sqrt{N}\tau _F`$ for $`R=13`$ and $`\delta =0.2`$, which is considerably above as expected.

As the last important fact we demonstrate that the errors for the layer correction indeed are maximal near the gap (i.e. near the surface). For Fig. 4 100 particles were put 100 times randomly in a box with a gap of $`\delta =0.1`$. Then for every particle the magnitude of the layer correction for $`R=40`$, which is a good approximation to the full $`E_{lc}`$, and the difference in the layer correction between $`R=5`$ and $`R=40`$ was drawn against the $`z`$-coordinate. Clearly the error always is largest near the gap. This effect increases with increasing $`R`$, which is easy to understand from the error formula. Therefore the full RMS error of the system might be completely misleading about the effect the errors have on the particles near the gap as we mentioned before. Nevertheless the figure shows that $`E_{lc}`$ with $`R=5`$ reduces the error near the surfaces by a factor of $`100`$.

## VI Conclusion

We have derived a term called ELC to efficiently calculate the contribution of the image layers in three dimensionally periodically replicated slab systems. ELC scales as the number $`N`$ of particles and has a rigorous error bound. Moreover this error bound can be used to estimate the size of the image layer contribution and therefore gives a bound on the error introduced by slab–wise methods as proposed by Yeh and Berkowitz. We have found that the error for these methods decays exponentially in $`L_z/L`$. However, the errors are not uniformly distributed over the slab, namely they are worst at the surfaces of the slabs. This strongly suggests to restrict the maximal pairwise error instead of the usually assumed RMS-errors.

In a forthcoming paper dejoannis02a  we will focus on the application of ELC to the standard Ewald method and to P<sup>3</sup>M. We will show how the error formulas of Kolafa and Perram kolafa92a  have to be adapted to allow non–cubic simulation boxes which is essential for using ELC with $`R=1`$, i. e. a slab–wise method. For all combinations we present numerical results which allow easily to decide which method is optimal for use in a real simulation.

## Acknowledgments

Financial support from the DFG “Schwerpunkt Polyelektrolyte” is gratefully acknowledged.


text:
# On Euler characteristic of equivariant sheaves

## 1. The main conjecture

### 1.1. Notations

In what follows $`k`$ denotes an algebraically closed field of characteristic $`p>0`$. Let $`\mathrm{}`$ be a prime number different from $`p`$. For an algebraic variety $`X`$ over $`k`$ we denote by $`𝒟(X)`$ the derived category of $`\mathrm{}`$-adic sheaves on $`X`$. Also we denote by $`\mathrm{Perv}(X)𝒟(X)`$ the subcategory of perverse sheaves. For any $`𝒟(X)`$ we denote by $`\chi ()`$ the Euler characteristic of $``$, i.e.

$$\chi ()=(1)^idimH_c^i(X,)=(1)^idimH^i(X,)$$
(1.1)
(cf.  for the proof of the equality).

Let $`T`$ be an algebraic torus over $`k`$. The following theorem is proved in  (cf. also  for a partial analogue in characteristic zero):

###### Theorem 1.2.

Let $`\mathrm{Perv}(T)`$. Then $`\chi ()0`$.

Let $`G`$ be a connected reductive algebraic group over $`k`$. We shall denote by $`G_{rs}G`$ the open subspace of regular semi-simple elements. We denote by $`\mathrm{Perv}_G(G)`$ the subcategory of $`\mathrm{Perv}(G)`$ consisting of perverse sheaves which are equivariant with respect to the adjoint action. We propose the following generalization of Theorem 1.2.

###### Conjecture 1.3.

Let $`\mathrm{Perv}_G(G)`$. Then $`\chi ()0`$.

###### Theorem 1.4.

Assume that $`\mathrm{Perv}_G(G)`$ is equal to the Goresky-MacPherson extension of its restriction to $`G_{rs}`$. Then $`\chi ()0`$.

###### Corollary 1.5.

Conjecture 1.3 holds for $`G`$ of semi-simple rank 1.

###### Proof.

Let $`𝒩`$ be the cone of unipotent elements in $`G`$ and let $`Z`$ be the center of $`G`$. Clearly it is enough to prove Conjecture 1.3 only for irreducible sheaves. However if $`\mathrm{Perv}_G(G)`$ is irreducible then either $``$ is equal to the Goresky-MacPherson extension of its restriction to $`G_{rs}`$ or it is supported on $`Z𝒩`$. The former case is covered by Theorem 1.4; hence it is enough to deal with the latter.

First of all there exists a connected component $`Z^{}`$ of $`Z`$ such that $``$ is supported on $`Z^{}𝒩`$. Thus one of the following is true:

1) $``$ is supported on $`Z^{}`$.

2) $``$ is equal to $`^{}(\overline{}_l)_𝒩[2]`$ where $``$ is an irreducible perverse sheaf on $`Z^{}`$ and $`(\overline{}_l)_𝒩`$ is the constant sheaf on $`𝒩`$.

3) $`=^{}`$ where $`^{}`$ is an irreducible perverse sheaf on $`Z^{}`$ and $``$ is an irreducible perverse sheaf on $`𝒩`$ whose restriction to the open orbit is isomorphic to (the only) non-trivial equivariant irreducible local system on this orbit (this local system is of rank 1 and its square is isomorphic to the constant sheaf).

In cases 1 and 2 our result follows immediately from Theorem 1.2 (note that $`Z^{}`$ is a torus and that in case 2 we have $`\chi ()=\chi (^{})`$). In case 3 it is known that $`H^{}(𝒩,)=0`$, hence $`H^{}(G,)=0`$, hence $`\chi ()=0`$. ∎

Remark. Theorem 1.2 has an ”ideological” explanation: namely in  O. Gabber and F. Loeser construct a ”Mellin transform” functor $`:𝒟^b(T)𝒟_{coh}^b(\text{Loc}_T)`$ where $`\text{Loc}_T`$ denotes the moduli space of tame rank one local systems on $`T`$ and $`𝒟_{coh}^b(\text{Loc}_T)`$ is the derived category of quasi-coherent sheaves on $`\text{Loc}_T`$. Moreover, for every $`\mathrm{Perv}(T)`$ the complex $`()`$ is actually a sheaf and $`\chi ()`$ is equal to the generic rank of $`()`$ (hence $`\chi ()0`$). We don’t know whether a similar explanation of Conjecture 1.3 is possible.

We expect that Conjecture 1.3 holds also when $`k`$ has characteristic 0 and perverse sheaves are replaced by holonomic D-modules. When $`G`$ is a torus and the D-modules in question have regular singularities a beautiful geometric proof of Theorem 1.2 was given in  where the authors explain how to compute the corresponding Euler characteristics using certain non-compact generalization of Kashiwara’s index theorem. It would be interesting to generalize the proof of  to the case of arbitrary reductive group. Some results in this direction have recently been obtained by V.Kiritchenko.

## 2. Proof of Theorem 1.4

### 2.1. The space $`\stackrel{~}{G}`$

Let $``$ denote the flag variety of $`G`$, i.e. the variety of all Borel subgroups of $`G`$. Let $`\stackrel{~}{G}`$ denote the variety of pairs $`(B,gB)`$. We have the natural maps $`\pi :\stackrel{~}{G}G`$, $`\alpha :\stackrel{~}{G}T`$ and $`\beta :\stackrel{~}{G}`$. It is easy to see that $`\alpha `$ is smooth and $`\pi `$ is proper. Let $`d=dimGdimT=2dim`$. Set $`\stackrel{~}{G}_{rs}=\pi ^1(G_{rs})=\alpha ^1(T_{rs})`$.

###### Lemma 2.2.

1. Let $`𝒢\mathrm{Perv}(T)`$. Assume that $`𝒢`$ is equal to the Goresky-MacPherson extension of its restriction on $`T_{rs}`$. $`\pi _!\alpha ^{}𝒢[d](\frac{d}{2})`$ is a perverse sheaf which is equal to the Goresky-MacPherson extension of its restriction to $`G_{rs}`$.
2. Every $`W`$-equivariant structure on a sheaf $`𝒢\mathrm{Perv}(T)`$ as above gives rise to a $`W`$-action on $`\pi _!\alpha ^{}𝒢[d](\frac{d}{2})`$.
3. Let $`\mathrm{Perv}_G(G)`$. Assume that $``$ is equal to the Goresky-MacPherson extension of its restriction to $`G_{rs}`$. Then there exists a $`W`$-equivariant sheaf $`𝒢\mathrm{Perv}(T)`$ which is equal to the Goresky-MacPherson extension of its restriction to $`T_{rs}`$ such that

$$=(\pi _!\alpha ^{}𝒢[d](\frac{d}{2}))^W.$$
(2.1)

###### Proof.

The first statement of Lemma 2.2 is well-known (it follows from the smallness property of $`\alpha `$ – cf.  and references therein)

It follows from 1 that in order to prove 2 it is enough to construct an action of $`W`$ on the restriction of $`\pi _!\alpha ^{}𝒢[d](\frac{d}{2})`$ to $`G_{rs}`$. Note that $`W`$ acts freely on $`\stackrel{~}{G}_{rs}`$ in the natural way and

$$G_{rs}=\stackrel{~}{G}_{rs}/W.$$
(2.2)
Moreover, the restriction of $`\alpha `$ to $`\stackrel{~}{G}_{rs}`$ is $`W`$-equivariant. Hence $`\alpha ^{}𝒢[d](\frac{d}{2})|_{\stackrel{~}{G}_{rs}}`$ is $`W`$-equivariant and thus (2.2) implies that $`\pi _!\alpha ^{}𝒢[d](\frac{d}{2})|_{G_{rs}}`$ has a natural action of $`W`$.

Let us prove 3. Choose an embedding $`i:TG`$. Let $`i_{rs}:T_{rs}G`$ be its restriction to $`T_{rs}`$. Let $`𝒢_{rs}=i_{rs}^{}[d](\frac{d}{2})`$. It follows from the $`G`$-equivariance of $``$ that $`𝒢_{rs}`$ is perverse. We let $`𝒢`$ be its Goresky-MacPherson extension to $`T`$. We leave the verification of (2.1) to the reader. ∎

###### Proposition 2.3.

Let $``$ be as in Theorem 1.4 and let $`𝒢`$ be as in Lemma 2.2(3). Then $`\chi ()=\chi (𝒢)`$.

Clearly Proposition 2.3 together with Theorem 1.2 imply Theorem 1.4.

### 2.4. Proof of Proposition Proposition 2.3

Set $`=\alpha ^{}𝒢[d](\frac{d}{2})\mathrm{Perv}(\stackrel{~}{G})`$. Clearly, $`H_c^{}(G,)=H_c^{}(\stackrel{~}{G},)^W`$.

Consider $`\alpha _!()`$. By the projection formula it is isomorphic to $`𝒢\alpha _!(\overline{}_l)_{\stackrel{~}{G}}[d](\frac{d}{2})`$. It is easy to see that the complex $`\alpha _!(\overline{}_l)_{\stackrel{~}{G}}[d](\frac{d}{2})`$ is constant with fibers isomorphic to $`H^{}(,\overline{}_l)`$. Indeed, let $`\delta :\times TT`$ denote the projection to the first multiple. Then $`\alpha _!(\overline{}_l)_{\stackrel{~}{G}}=\delta _!(\beta \times \alpha )_!(\overline{}_l)_{\stackrel{~}{G}}`$. However, the map $`(\beta \times \alpha ):\stackrel{~}{G}\times T`$ is a locally trivial fibration (in Zariski topology) with fiber isomorphic to $`𝔸^d`$. Thus $`(\beta \times \alpha )_!(\overline{}_l)_{\stackrel{~}{G}}(\overline{}_l)_{\times T}[d](\frac{d}{2})`$ and hence

$$\alpha _!(\overline{}_l)_{\stackrel{~}{G}}(\overline{}_l)_TH^{}(,\overline{}_l)[d](\frac{d}{2}).$$
(2.3)
This clearly gives rise to the isomorphism

$$H_c^{}(\stackrel{~}{G},)H_c^{}(T,𝒢)H^{}(,\overline{}_l).$$
(2.4)
We can also characterize the isomorphism (2.3) in the following way. To determine it uniquely it is enough to construct it on $`T_{rs}`$. Choose as before an embedding $`i:TG`$. There is a canonical Borel subgroup $`B`$ containing $`i(T)`$ (recall that the abstract Cartan group $`T`$ comes with a root system with a preferred set of positive roots). Then for every $`tT_{rs}`$ we have canonical isomorphism

$$\alpha ^1(t)G/T$$
(2.5)
(which depends however on the above choice). Clearly $`H_c^{}(G/T,\overline{}_l)[d](\frac{d}{2})=H^{}(,\overline{}_l)`$. Hence we have constructed an isomorphism between $`\alpha _!(\overline{}_l)_{\stackrel{~}{G}}[d](\frac{d}{2})|_{T_{rs}}`$ and the constant complex with fiber $`H^{}(,\overline{}_l)`$. It is easy to see that this isomorphism does not depend on the choice of the embedding $`TG`$ and coincides with (2.3).

Let $`W`$ act on the left hand side by means of the identification $`H_c^{}(\stackrel{~}{G},)=H_c^{}(G,\pi _!)`$ (note that by Lemma 2.2 the group $`W`$ acts on $`\pi _!`$). Let $`W`$ also act on the right hand side by means of the tensor product of the $`W`$-action on $`H_c^{}(T,𝒢)`$ coming from the $`W`$-equivariant structure on $`𝒢`$ and the natural $`W`$-action on $`H^{}(,\overline{}_l)`$.

###### Lemma 2.5.

The isomorphism (2.4) is also an isomorphism of $`W`$-modules with respect to the above actions.

###### Proof.

Let $`Z=G\underset{T/W}{\times }T`$ be the image of $`\stackrel{~}{G}`$ in $`G\times T`$ under the map $`\gamma =\pi \times \alpha `$. This is a closed subvariety of $`G\times T`$. It is invariant with respect to the $`W`$ action on the second multiple. Let $`Z_{rs}`$ denote the set of ”regular semi-simple” elements of $`Z`$ (i.e. the set of all elements of $`Z`$ whose projection to $`Z`$ is regular semi-simple).

Let $`𝒦=\gamma _!()`$. Then $`𝒦`$ is equal to the Goresky-MacPherson extension of its restriction to $`Z_{rs}`$. Since $`\gamma `$ induces an isomorphism $`\stackrel{~}{G}_{rs}Z_{rs}`$ and since the restriction of $``$ to $`\stackrel{~}{G}_{rs}`$ is $`W`$-equivariant it follows that $`𝒦`$ also has a natural $`W`$-equivariant structure as a perverse sheaf on $`G\times T`$ (where the $`W`$-action is as before on the second multiple). Thus $`W`$ acts naturally on $`H_c^{}(G\times T,𝒦)`$. We have the natural isomorphism $`H_c^{}(\stackrel{~}{G},)H_c^{}(G\times T,𝒦)`$ and by the definition the action of $`W`$ on $`H_c^{}(\stackrel{~}{G},)`$ introduced before Lemma 2.5 corresponds to the action of $`W`$ on $`H_c^{}(G\times T,𝒦)`$ introduced above.

Consider, on the other hand, the complex $`p_!(𝒦)`$ where $`p:G\times TT`$ denotes the natural projection. Clearly, we have $`p_!(𝒦)=\alpha _!()𝒢H^{}(,\overline{}_l)`$. On the other hand, since $`p`$ commutes with $`W`$ it follows that $`p_!(𝒦)`$ admits a natural $`W`$-equivariant structure. To prove Lemma 2.5 it is enough to show that under the identification $`p_!(𝒦)𝒢H^{}(,\overline{}_l)[d](\frac{d}{2})`$ this structure is equal to the tensor product of the original $`W`$-equivariant structure on $`𝒢`$ with the natural $`W`$-action on $`H^{}(,\overline{}_l)`$. Since every perverse cohomology of $`p_!(𝒦)`$ is equal to the Goresky-MacPherson extension of its restriction to $`T_{rs}`$ it follows that it is enough to check this equality only on $`T_{rs}`$. This, however, follows immediately from the description of the $`W`$-equivariant structure on $`\alpha _!(\overline{}_l)_{\stackrel{~}{G}}[d](\frac{d}{2})|_{\stackrel{~}{G}_{rs}}`$ given in Section 2.4 and the following observation:

Choose as before an embedding $`TG`$. Let $`N(T)`$ denote the normalizer of $`T`$ in $`G`$. Then $`W=N(T)/T`$ acts on $`G/T`$ and hence on $`H_c^{}(G/T,\overline{}_l)`$. By identifying $`H_c^{}(G/T,\overline{}_l)[d](\frac{d}{2})`$ with $`H^{}(,\overline{}_l)`$ we get an action of $`W`$ on the latter. This is the standard $`W`$-action on $`H^{}(,\overline{}_l)`$. ∎

It follows from Lemma 2.5 that $`\chi ()`$ is equal to the Euler characteristic of $`(H_c^{}(𝒢)H^{}(,\overline{}_l))^W`$. However, it is known that $`H^{}(,\overline{}_l)`$ lives only in even degrees and when we forget the grading it is isomorphic to the regular representation $`\overline{}_l[W]`$ of $`W`$. Hence

$`\chi ()=\chi ((H_c^{}(T,𝒢)`$ $`H^{}(,\overline{}_l)[d])^W)=`$
$`\chi ((H_c^{}(T,𝒢)\overline{}_l[W])^W)=\chi (H_c^{}(T,𝒢))=\chi (𝒢)`$
which finishes the proof.

### 2.6. Acknowledgements

The author is grateful to F. Loeser for an illuminating discussion on the subject held in the summer of 2001, to V. Kiritichenko for making her results available to the author before publication and to M. Kapranov for explaining the contents of .