Nuclear stopping power of antiprotons

The slowing down of energetic ions in materials is determined by the nuclear and electronic stopping powers. Both of these have been studied extensively for ordinary-matter ions. For antiprotons, however, there are numerous studies of the electronic stopping power, but none of the nuclear stopping power. Here, we use quantum-chemical methods to calculate interparticle potentials between antiprotons and different atoms, and derive from these the nuclear stopping power of antiprotons in solids. The results show that the antiproton nuclear stopping powers are much stronger than those of protons, and can also be stronger than the electronic stopping power at the lowest energies. The interparticle potentials are also implemented in a molecular dynamics ion range calculation code, which allows us to simulate antiproton transmission through degrader foil materials. Foil transmission simulations carried out at experimentally relevant conditions show that the choice of antiproton-atom interaction model has a large effect on the predicted yield of antiprotons slowed down to low (a few keV) energies.


I. INTRODUCTION
The production of stable antihydrogen atoms relies on the slowing down of antiprotons (p) with initial energies of the order of MeV's or keV's to thermal energies [1]. The first step for achieving this is to use thin degrader foils to reduce the kinetic energy to the order of a few keV's [2,3]. Even though this technology has been demonstrated to work, the antihydrogen yields are low, making it desirable to design higher yield degrading approaches [4,5]. Thus, it is important to understand the physics behind the retardation of antiprotons in materials.
The slowing down of ions in matter is conventionally described with the stopping power S, i.e., the energy loss of an energetic particle per path length travelled in the solid [6][7][8][9] (also called in some contexts the linear energy transfer (LET) [10]). The stopping power can be separated into nuclear (S n ), electronic (S e ) and nuclear reaction (S nr ) contributions, S = S n + S e + S nr [6,8,[11][12][13][14]. The nuclear reaction part is usually negligible at low energies, whereas the nuclear stopping becomes then more significant. Although the electronic stopping power of antiprotons in matter has been extensively studied [15][16][17][18][19][20][21], there have been no systematic theoretical studies The nuclear stopping power can be calculated by using classical scattering theory [6,7]. To get the spatial description of the pathway of an energetic ion in the material, one can use the binary collision approximation (BCA) [22][23][24][25] or molecular dynamics (MD) ion range simulations [26][27][28]. These computational approaches require the interaction potential between the atoms of the ion lattice as input data [29]. The interaction potentials of normal matter ions with atoms have been studied very extensively (see e.g., Refs. [7,[29][30][31][32]). However, while there are some calculations of antihydrogen interactions with hydrogen and helium [33][34][35][36][37][38], there are no systematic computational studies of the interaction energies of antiprotons with heavier atoms than He. Note that in simulations of the interaction process of energetic antiprotons with solid foils, one needs specifically the information about the interaction of a bare antiproton with atoms.
In this Article, we employ quantum chemical methods to calculate the interaction energy between antiprotons and atoms of several elements (Section II), in order to determine the nuclear stopping power of antiprotons (Section IV). We also use the determined potential for molecular dynamics simulations of antiproton transmission through energy degrading foils, under conditions relevant for near-future antihydrogen production experiments (Section V).

A. The interparticle potential of antiprotons
The full antiproton interaction potentials were calculated for H, Be, C, N, O, Al, Si, and Ti, because the degrading foil materials currently considered at CERN consist of these elements or their compounds. The interaction potential was also calculated for Ne in order to elucidate how antiprotons interact with closed-shell elements.
The interaction potential between the atoms and an antiproton was calculated at the second-order Møller-Plesset (MP2) [39] level using the resolutionof-the-identity approximation [40] and triple-ζ quality Gaussian-type basis sets augmented with polarization functions (def2-TZVP) [41]. The antiproton was modelled as a negative point charge. The hydrogen def2-TZVP basis set was employed for the antiproton. For Be, N, Ne, and Al, we also augmented the def2-TZVP basis set with diffuse functions to get a more accurate description of the low-energy range of the interaction potential [42][43][44]. The potential energy curves were calculated for the lowest electronic state of the atoms. Essentially identical potential-energy curves were obtained at the Hartree-Fock and MP2 levels of theory, showing that electron correlation does not play an important role in this context. The calculations were performed with the Turbomole code [45,46].  The potentials were calculated for selected interparticle separations in the range of [0.001,50] bohr. The results from the Turbomole calculations are given in Appendix A. The antiproton-atom interaction potentials V (r) used in the MD simulations were then obtained by subtracting the total energy calculated for a separation of 50 bohr from the total energy calculated for the rest of the interparticle distances. The calculated potentials for some of the elements are illustrated in the inset of Fig. 1 for low energies.
When the diffuse basis functions were included in the basis sets, none of the antiproton-element pairs showed any repulsion at any interparticle distance r. This shows that the strong Coulomb attraction between the antiproton and the nucleus of the atom dominates in the whole separation range. The high-energy part is almost unaffected by diffuse functions, and the difference in the results with and without them is always < ∼ 1 eV. Inaccuracies of the order of 1 eV in the long-ranged potential do not in any case significantly affect the simulated nuclear stopping power, because the antiproton beam energies are at least two orders of magnitude larger.
The second main feature of the calculated interaction potential is that the potential obtained by considering all electrons is clearly above the Coulomb interaction of the bare nuclear charges (Fig. 1), even though the potential is attractive in the whole range. Hence, the interaction between the antiproton and the atomic nucleus is strongly screened by the electron cloud. However, at larger separations there is a major difference due to chemical effects (as also evident from the inset of ( Fig. 1). Comparison of the screening functions in Fig. 2 also show that at separations of > ∼ 0.5Å the screening is completely different between the proton and antiproton.
A comparison of all the quantum chemically calculated potentials is shown in Fig. 3. The y axis energy E is plot-  ted in (a) as − log(−E) to allow for logarithmic plotting of negative values. The comparison in (a) shows that the bare Coulomb interaction (corresponding to a linear trend with slope of -2 in the log-log plot) is obtained only at distances below ∼0.1Å, i.e., at longer distances the screening effects are significant.

B. Comparison of polarizability values with experiments
The quality of the calculated interparticle potentials can be checked by comparison with experimental data on atomic polarizabilities. At large distances, r, the interaction between the antiproton and the atom is dominated by two terms. One is the charge-polarizability term (in atomic units) where α is the polarizability of the atom, and Z the charge of the approaching particle. When the atom is non-spherical having a quadrupole moment Q a , there will also be a quadrupole interaction where A is a factor that depends on the orientation of the atom. For filled-shell atoms only Eq. (1) survives. The polarizabilities from that fit compare favorably with standard values, see Table I. A fit to the MP2 points for neon and aluminum, as an example of an element with an atomic quadrupole moment, is shown in Fig. 4.
The reasonable agreement with experiments obtained here gives an indication the calculation results are reliable at large ( > ∼ 2Å) interparticle separations. On the other hand, at small separations ( < ∼ 0.1Å) all the potential data converge very accurately to the pure Coulomb potential of unscreened nuclei, as expected, see   [47], b Ref. [48] (exp), c Ref. [49] (exp). a). At intermediate separations there is no reference case that can be compared to, however, these two observations show the quantum chemical calculations can describe well the antiproton-atom interactions at least at small and large separations.
Finally, the interactions between atoms or molecules on one hand, and antiprotons on the other hand are sufficiently specific to be called antiproton chemistry.

C. Description of potentials with screening functions
To construct an analytical or an accurate numerical representation of the interparticle potentials from the quantum-chemically calculated interaction energies, we analyzed the potentials in terms of a screened Coulomb potential. This approach has been found to work well for repulsive interatomic potentials of normal matter [7,50,51]. The screening function φ(r) was obtained from the potential V (r) using where Z 1 = −1 is the charge of the antiproton, Z 2 > 0 the nuclear charge of the atom, and e the elemental charge.
The screening function illustrated in Fig. 2 for Si displays a nonlinear behaviour on the logarithmic scale. This is natural due to the shell structure of atoms. There is no known simple analytical function that can accurately describe the screening of repulsive pair potentials between atoms [7].
In this work, we consider two possible forms of the screening function. One is a sum of exponential functions, which follows common practise of repulsive potentials for normal matter [7] and hence can be easily implemented in existing codes. This form is presented in section II C 1. However, this form does not have the physically well motivated 1/r 3 or 1/r 4 form discussed in section II B. Hence we also constructed a new functional form for the screening, that has the correct limiting shape and also can be well fit to the data. This form is presented in section II C 2.

Exponential screening function φ exp
We found that a screening function consisting of a sum of a number of exponential functions can be fitted well to repulsion potentials that are deduced from quantum chemical calculations or from experimental data [7,52]. The calculated MP2 energies are here fitted to a screening function of the form The coefficient of the last term is b the resulting potential is a pure Coulomb potential of the bare nuclear charges at the smallest distances. The def2-TZVP basis set was used for all elements but Al. For Al, the data obtained with the def2-TZVP basis set augmented with the diffuse functions were fitted to Eq. (4) (without them, the obtained potential exhibited a spurious weak (about 0.1 eV) repulsive region at large separations). The fit is illustrated in Fig. 2 and the values of the fitted coefficients are listed in Table II. The fitting was performed with the Levenberg-Marquardt algorithm [53]. Fits with 3 or less exponential terms did generally not yield accurate fits to the calculated data at all distances. Typically, with 3 exponential terms, a good fit could not be obtained at the same time around distances of both ∼ 0.3 and ∼ 3Å. Very good fits were obtained with the 4 terms in Eq. (4), giving Pearson's χ 2 values about a factor of 10 smaller than the 3-term fits.
We note that in the H screening function, the first prefactor b 1 > 1, and hence the balancing prefactor b 4 is negative. This combination of factors was found to give the best fit to the data. We checked that, in spite of one negative prefactor, the sum of the exponentials is a monotonously decreasing function, as it should be for a purely attractive potential.
Comparison of thep-Si screening function with the standard universal repulsive Ziegler-Biersack-Littmark (ZBL) function [7] for H-Si in Fig. 2 also shows that the calculatedp-Si interaction potential is very different from the ZBL potential. Chemical interactions further enhance this difference at low energies, see inset in Fig. 1, which compares thep-Si interaction with the interaction energy of a regular proton with Si [29]. Thus, one should not use negative proton or repulsive hydrogen potentials to mimic antiproton interactions.

Screening function with limiting power law behaviour φ ep
The exponential fits were found to give very good fits to the data, especially at the smaller distances r that dominate the nuclear stopping. However, they do not have the limiting 1/r 3 or 1/r 4 form at large distances discussed in section II B. Hence we also constructed a new form of the screening function that has either the 1/r 3 or 1/r 4 behaviour when r → ∞. We found that combining two or three exponential terms with two simple Padé approximant-like [54] terms 3 (5) could be used to give good fits to all data. Here the fitted constants are d i and c i . Since the screening function is multiplied with the Coulomb term ∝ 1/r, the latter two terms 1/(1 + (d i r) g ), g = 2 or 3, lead to the desired 1/r 3 or 1/r 4 limit for large r values, respectively. Naturally when r → ∞, the g = 2 term will dominate, except if TABLE II. Coefficients X = bi or ai of the fitted antiproton screening function φ exp fit (r) in Eq. (4) obtained from the calculated potentials of the different elements. The last b term is always b4 = 1 − b1 − b2 − b3 to ensure that the potential is a pure Coulomb potential of the bare nuclear charges when The resulting fits are presented in Table III. In all cases, not all terms were needed to achieve a good fit, in which case d i = 0 for the unnecessary terms. In some cases, some of the terms have a negative prefactor. In these cases, we checked that the screening function is nevertheless monotonously decreasing at all distances (which includes the requirement that d 4 > 0, or in case d 4 = 0, that d 5 > 0).
A comparison of the screening function fits in the two different forms is shown in Fig. 5. The comparison shows that both functional forms give very good fits to the quantum chemical data. The very small differences are not expected to lead to significant differences in the nuclear stopping power or foil transmission simulation results, which is confirmed by results shown later in this Article. Unless otherwise mentioned, the results in the remainder of this article are obtained with the potentials constructed from the φ exp type screening functions.

III. ELECTRONIC STOPPING POWERS OF AL AND SI
To enable comparison with the nuclear stopping power, as well as molecular dynamics simulation of antiproton movement in materials, we constructed curves of the electronic stopping from the experimental data in Ref. 18 for the elements Al and Si of interest in the current study. An interpolation function was used to obtain a smooth  (5), obtained from the calculated potentials of the different elements. The last d term is always d5 = 1 − d1 − d2 − d3 − d4 to ensure that the potential is a pure Coulomb potential of the bare nuclear charges when r → 0. The fairly large number of decimals given for the d terms for oxygen are needed to ensure that the balance term d5 is obtained accurately. electronic stopping curve. The interpolation function was constructed such that it ensures that the stopping is linear with velocity at small velocities, to conform with stopping theory and experimental observations [19]. The data and fit are shown and compared to the electronic stopping of protons in Fig. 6.

A. Scattering calculations
Using the screening functions, one can construct continuous interparticle potentials for all distances The resulting continuous potentials are shown in Fig.  3b.
The continuous potentials were used as basis for classical scattering theory calculations to determine the nuclear stopping power of the antiprotons (S n ). This is obtained from the energy transfer of a single binary collision from the projectile with an initial energy E 0 to the sample atom as [7] where N is the atomic density of the target, σ the cross section, T is the energy loss in the binary collision, and b the impact parameter of the collision.
For atoms of normal matter, the interaction between the positive ion and the positive core of the sample atom is at high energies of > ∼ 10 eV purely repulsive [29,55]. Hence, the scattering trajectory has the form of hyperbolic orbits receding from each other, and Eq. (7) can be evaluated from the classical scattering integral in a straightforward way for any repulsive potential [7,55].
The situation is more complicated for antiprotons. Since our quantum-chemical calculations show that the interaction has the form of a screened, purely attractive potential, the antiproton approaches the nucleus. For small impact parameters and low kinetic energies, the antiproton may follow complex trajectories around the nucleus (this is possible since the potential is not a simple 1/r form [56]). It is hard to treat such trajectories using the binary scattering integral approach (although the binary collision approximation has been implemented for potentials with an attractive well [57], this approach cannot be directly used for a purely attractive potential).
To be able to deal with arbitrarily complex trajectories, we employed molecular dynamics simulations [31,32] to model the interactions of the antiproton with the atoms. We considered both the single-atom scattering and transmission through foils consisting of many atoms. We used the MDRANGE method [26] at the recoil interaction approximation [26,27] to simulate the retardation of antiprotons in materials. The MDRANGE method has previously been found to yield very good agreement between simulated and experimental ion range distributions [28,[58][59][60]. The nuclear stopping power was extracted from the scattering calculations by switching off the electronic interaction in MDRANGE simulations for a system consisting of one antiproton and a single atom. The adaptive time-step algorithm implemented in the MDRANGE code was found to be crucial for treating accurately the very high acceleration of the antiproton along trajectories bound in towards the nucleus [26]. Energy-conservation tests showed that constants of E t = 0.5 eV and k t = 0.01 A in the algorithm to select the time step ∆t [26] ∆t new = min yielded a good description of the trajectories down to interparticle distances of 3 fm. Here v max is the maximum velocity in the system and F max the maximum force that affects the ion from any other atom. If the antiproton came closer than 3 fm to the nucleus, a nuclear reaction was assumed to occur. The simulation was then stopped and the initial energy was added to the nuclear reaction stopping power S nr rather than S n . The analysis was carried out for initial antiproton energies (E 0 ) in the range of [2 eV, 1 MeV].
In addition to the tests of the accuracy of the calculations with the time step, we also tested the numerical reliability of the employed method by performing calculations for proton interactions described by the normal ZBL universal interatomic potential [7] using exactly the same MDRANGE approach. In the proton simulations, the method should give the same stopping power as obtained with the analytical ZBL stopping model.

B. Extraction of nuclear stopping
The antiproton and proton nuclear stopping power were obtained from the energy loss function T (E 0 , b) in the MDRANGE scattering calculations. The stopping power evaluated using Eq. (7) by numerical integration over b is shown in Fig. 8a for Si. The nuclear stopping power obtained for the proton agrees very well with the one calculated using the ZBL model, confirming that the numerical accuracy of the MDRANGE approach is good. The nuclear reaction stopping power is of the order of 10 −5 keV/nm, which is too small to be visible in the figures.
The results of the nuclear stopping power calculations in Fig. 8 show that the nuclear stopping power of antiprotons is in all cases clearly higher than that of protons at low energies (for clarity in the figure, proton data is shown only for Si). To analyze the reason to this, we compared the scattering of antiprotons with protons for identical initial energy E 0 and initial impact parameter b. Fig. 7 shows as an example case the scattering for the energy of 1 keV. The results show that for the same impact parameter b, the antiprotons scatter in larger angles than the proton, and lose more energy. The figure also shows that a pure Coulomb potential between charges of −1e and 14e gives very different results, highlighting the importance of obtaining a screened interparticle potential for antiproton-atom interactions.
The results in Fig. 8 also show a small bump around 10 eV in the antiproton S n in some elements. Analysis of individual ion trajectories showed that for these very low The solid curve is a fit of an interpolation curve constructed to ensure linear dependence with velocity at low velocities (the interpolation starts from S = 0 at v = 0, not visible in the log-log plot). The electronic stopping power for protons is the stopping curve obtained from the original ZBL model [7]. For the antiprotons, the nuclear stopping is shown based on both the exponential screening function fit φ exp and the exponen-tial+power law fit φ ep . The inset shows as a clarification Sp (from φ exp ) and Se for antiprotons only. b) Nuclear stopping Sn of antiprotons for all the solid-state elements considered in this study, compared to the electronic stopping Se for those elements for which it is available at low energies (Al and Si) energies, for impact parameters around 0.5Å, the antiproton has complex trajectories, e.g., doing a full circle around the atom. This leads to a higher energy transfer than the simpler hyperbolic paths observed at higher energies. In Fig. 8 b, we also compare the antiproton nuclear stopping with their electronic stopping power for those elements for which it is available from experimental data in Ref. 18 (Al and Si). The results show that at low energies, the antiproton nuclear stopping power for both Al and Si is almost an order of magnitude larger than the electronic ones. This is in sharp contrast to the case of protons, for which it is well established that the electronic stopping power dominates over the nuclear one at all energies and in all materials [7,61]. This results implies that when increasingly low energy antiprotons are used e.g., in antihydrogen experiments such as [2][3][4][5], it becomes crucial to include the antiproton nuclear stopping in the modelling of the process.
The results obtained with the two different screening fit functions for Si, φ exp and φ ep (Fig. 8 a), show a small variation at the lowest energies, but converge to be practically identical at the higher ones. A similar comparison for H and Ne (not shown) also showed a small difference in stopping powers below the maximum -smaller than the difference in Si -and essentially identical results around and above the stopping maximum. Since the fits are of comparable quality, it is not a priori clear which of the versions is more accurate. This indicates that in case very high accuracy nuclear stopping powers are needed for the very low ( < ∼ 100 eV) energies, the quantum chemical calculations would need to be expanded with additional data to enable higher-accuracy fits. However, we emphasize that the S n [φ exp ] -S n [φ ep ] difference is much smaller (of the order of 10 % at most) than the difference between the antiproton and proton stopping powers (of the order of a factor of 5), or the difference between the antiproton nuclear and electronic stopping powers. Thus the small uncertainty coming from the choice of the fitting functions does not affect the conclusions of this work. Fig. 9 presents the nuclear stopping powers S n [φ exp ] calculated for all the elements considered in this study, normalized with the densities. The electronic supplementary information package collected_stopping_powers.zip contains the original data used in the plot, as well as a plot of all stopping powers without density normalization. For the materials gaseous at room temperature, the same solid phase densities as those used in the ZBL model (SRIM) nuclear stopping power calculations were used. The density values are in all cases given in the figure legend. The data show that, as expected from basic stopping theory, the maximum in energy of the stopping increases with increasing target atomic number Z 2 . On the other hand, the normalization with atomic density shows that the density-scaled stopping power decreases with increasing Z 2 . This is because the energy transfer in binary collisions is the most efficient when the masses of the projectile and target are the same.

V. MOLECULAR DYNAMICS OF ANTIPROTON TRANSMISSION IN SI
The effect of the unexpectedly large antiproton stopping power on particle movement in solids and foil transmission experiments was tested by using the MDRANGE method [26] to simulate proton and antiproton transmission through Si foils with completely amorphous atomic structure [28]. The proton was in these simulation modelled with the ZBL repulsive interatomic potential [7], to be consistent with BCA codes. The ions were in all cases shot perpendicularly to the foil surface. The electronic stopping power was included in the simulations as a frictional force [26]. Unless otherwise mentioned, the electronic stopping was implemented including the Bohr model of the straggling of electronic stopping [6,62,63] in the simulations. The straggling of electronic stopping can have a significant effect on the range and energy distributions of low-energy transmitted protons [63,64]. The results in Figs. 10 and 12 show that in the current case it does have a weak, but statistically significant broadening effect on the range profiles and energy distributions of transmitted protons and antiprotons. For instance, the straggling (standard deviation) of the range distribution of antiprotons (Fig. 10) is 1781 ± 7Å with, and 1715 ± 12 A without, the straggling of the electronic stopping.

A. Comparison of proton and antiproton penetration depths
The large difference in the electronic stopping power and interparticle potentials between antiprotons and protons will naturally lead to a large difference in particle penetration depths. This difference is illustrated in Fig.  10, which compares the depth distribution of 100 keV antiprotons or protons implanted in bulk (infinitely thick) amorphous Si. Since the electronic stopping of antiprotons is weaker than that of protons (cf. Fig 6), the overall penetration depth of antiprotons is much deeper. For the use of energy degrading foils, this means that to slow down antiprotons to low energies, one needs to use clearly thicker foils than one would need for protons.
Since it would clearly be inappropriate to use proton electronic stopping for antiprotons, in the remainder of the paper we focus on simulations with the experimentally determined antiproton stopping, and examine how large effects the choice of the interparticle potential (nuclear stopping power) can have on the results.

B. Scattering in very thin foils
To illustrate the scattering effects, we first consider the case of antiprotons and protons passing a very thin 2 nm Si foil. The resulting trajectories are illustrated in Fig.  11. At low energies, the energy loss in the foil is much larger for antiprotons than for protons. Simulations were also performed for a hypothetical particle with the nuclear stopping power of a proton S n (p) and the electronic stopping power of an antiproton S e (p). Comparison of the antiproton results with those for the hypothetical particle shows that the nature of the interparticle potential leads to the most significant differences between antiproton and proton behaviour in the Si foil.
Inspection of individual trajectories showed that in foils, where energy loss can occur to several material atoms at the same time, the antiproton can be trapped in the potential well approaching the nucleus along helical trajectories, which is a possibility when the potential has a shape differing from that of a pure 1/r potential [56]. The antiproton trajectory can even form a stationary bound orbit [65]. This was observed to occur in a few cases, and was a major practical problem in the MD simulations, since a stationary antiproton-atom pair essentially takes the simulation into an infinite loop (the classical MD simulations do not include the possibility of quantum mechanical tunneling). To circumvent this, we stopped the simulation of a single antiproton case if it had been run for more than 100 million steps. Such cases were calculated into the nuclear reaction stopping power.
The larger scattering angle is demonstrated by the trajectories of the 10 keV particles in Fig. 11b, where the antiprotons are seen to be much more scattered than the protons. The large scattering angle also increases the fraction of backscattered antiprotons. The total energy loss of the particles was 66±1 eV for antiprotons and 43±1 eV for protons, i.e., the larger scattering leads to significant difference in the energy loss.

C. Antiproton transmission through degrading foils
As a case that is directly relevant to near-future antihydrogen experiments, we used MDRANGE to simulate the motion of 100 keV antiprotons that will be produced in the ELENA ring and degraded by foils to energies of a few keV. At these low energies, the antiprotons can be captured in a Penning-Ioffe trap, the immediate step prior to antihydrogen production. We used the simulations to determine the optimal film thickness with respect to maximum yield of transmitted protons in this energy range. We use here energies of 5 and 10 keV as examples of the limit for which transmitted antiprotons can be trapped. When the energy limit under experimental condition does not follow a sharp step function or it also depends on the emission direction, the analysis can be refined to account for that. The qualitative conclusions of this work are not expected to depend on minor changes in the energy limit.
To assess whether the nuclear stopping power of antiprotons matters for the experimental situation, we considered three cases: i) thep-Si interparticle potential determined in the current work, i.e. having the antiproton nuclear stopping power S n (p), ii) the hypothetical case of a particle with the electronic stopping power of an antiproton S e (p) and the proton-Si interatomic potential, i.e., the nuclear stopping power of a proton S n (p) and iii) the case of simulating antiprotons with no interparticle potential, i.e., zero nuclear stopping. For all cases, we simulated the antiproton passage until it was either captured by an atom in the foil, or transmitted through it. For the transmitted particles, statistics of their kinetic energy after passage through the foil E trans was collected.
The results on the energy distribution of transmitted particles is shown in Fig. 12 a) for one thickness (1550 nm) at which the production yield is appreciable with all three models. The simulations show that the fraction of transmitted antiprotons is affected by the choice of the nuclear stopping model. The two different screening functions φ exp and φ ep calculated for true antiprotons give identical results within the statistical uncertainty. However, the other nuclear stopping models (ii) and (iii) give clearly different results. For both true antiproton models (i), the useful transmission yield is 0.195±0.006, while for cases (ii) and (iii) it is 0.35±0.01 and 0.54±0.02, respectively, for the upper limit of 5 keV. Fig. 12 b) shows that the optimal film thickness for transmission up to 5 keV is 1500 nm. For an upper limit of 10 keV, the optimal thickness is found to be 1400 nm (Fig. 12 c), giving useful transmission yields by each model as (i) 0.65±0.02, (ii) 0.72±0.02 and (iii) 0.69±0.01.
The lower transmission yields of antiprotons are due to their larger scattering discussed above. The results in Fig. 12 b show that there are marked differences for all relevant film thicknesses.
The MD range calculation approach developed here will thus enable calculations for systematically predicting the antiproton transmission yield to low (keV) energies, and thus optimizing the film thickness and material.

D. Lateral range and emission angle distributions
As evident from Fig. 11, the antiprotons do not travel in straight paths through the material. The lateral straggling and distribution of angles after they are emitted may be of importance for the trapping of the antiprotons in the antihydrogen production apparatus.
To get insight on these effects, we analyzed the lateral (perpendicularly to the incoming beam direction) distribution of the antiprotons as well as the direction in which they move after emission. This analysis was done for the physically best motivated model of antiprotons that considers S e straggling. Fig. 13 shows the lateral range R ⊥ for 2000 nm and 1550 nm thick foils. In the 2000 nm case, all antiprotons stop in the foil (cf. Fig. 10), and the analysis was done on the final positions of the antiprotons. The lateral ranges R ⊥ (having a mean value of about 350 nm) are clearly smaller than the mean range of 1400 nm. This also explains why the chord range [6] distribution in Fig.  10 differs only slightly from the penetration depth distribution.
For the antiprotons travelling through 1550 nm foils, the analysis was done only for the final x and y positions of the antiprotons when they have exited the foil. These distributions are naturally narrower than those for the 2000 nm foil. Practically all the transmitted antiprotons exit the foil laterally within 500 nm from the initial impact point. Compared to the thickness of the foil, this shows the lateral spreading is within 1/3 of the foil thickness.
We also analyzed the distribution of outcoming emission directions of antiprotons that transmitted through the foils. The results in Fig. 14 show that for the foil thickness of 1550 nm, the distribution of emission angles is fairly wide, and that it does not depend strongly on the energy after transmission E trans . The higher-energy transmitted antiprotons have on average lower emission angles, since they have scattered less. For the same reason, the angular distribution of antiprotons exiting from a thinner foil is also peaked at lower angles.

E. Channeling effects
We also consider whether ion channeling [28,66] effects could enhance the transmission of antiprotons in the useful energy range. Antiproton channeling has previously been studied at higher energies [20]. Here we consider the case of 100 keV antiproton transmission for slowing down to the 0 -5 keV energy range. We assume perfect alignment of the antiproton beam with either a 100 or 110 crystal direction in crystalline Si. These simula- tions were run for the foil thickness of 1550 nm similarly as for a-Si. The results in Fig. 15 show that the fraction of ions transmitted in the energy range 0 -5 keV could indeed be somewhat enhanced by channeling. For the thickness of 1550 nm, the useful (0 -5 keV) transmission yield is 0.195±0.006 for amorphous Si, while for a 100 crystal direction it is 0.22±0.01 and for a 110 direction 0.21±0.01. However, the enhancement is small for all film thicknesses, and it is doubtful whether the minor improvement would justify the experimental efforts needed to use single-crystalline foils and orient them accurately in the beam direction.

F. Comparison with binary collision approximation
Finally, we analyze to which accuracy the computationally efficient BCA simulations can be used to simulate antiproton trajectories in thin films. For this purpose, we used the BCA code CASWIN that has been tailored to read in the same electronic stopping files as MDRANGE [25,67]. The good agreement between the MD and BCA results for protons in Fig. 10 confirms that the two codes are consistent with each other when the electronic stopping and interparticle potential are identical. Fig. 11b shows results of BCA simulations for case (ii). The actual antiproton-Si case (i) cannot be handled since BCA does not work for purely attractive potentials. The resulting good agreement with the MDRANGE case (ii) shows that conventional BCA simulations can be used to handle the electronic stopping power of antiprotons, but not the nuclear part. The minor MD-BCA difference for case (ii) appears because the straggling of electronic stopping is not considered in the BCA code.

VI. CONCLUSIONS
In conclusion, we have used ab initio quantum chemistry methods to calculate the antiproton-atom interaction potentials for a number of elements, and molecular dynamics simulations to determine the nuclear stopping power of antiprotons. The calculations show that the antiproton-atom interaction is purely attractive, and is screened at much further distances than the interactions of a proton with the same atom. This renders the nuclear stopping of antiprotons much stronger than that of protons, and also makes the antiproton nuclear stopping power stronger than the electronic one at low energies. The nuclear stopping power contributes significantly to the total deceleration of antiprotons in thin foils. Hence, the nuclear stopping power has to be considered when degrading foils are used for slowing down antiprotons to keV or lower energies.  The original Turbomole data for the antiproton-atom interaction energies is given in Tables IV -VI, as a function of the interparticle separation r. The lengths are given in units of bohrs and energies in units of hartrees (Ha). The large number of decimals is needed to avoid roundoff errors when the energy at infinity is subtracted out from the values to give an interatomic potential that ends at zero.