text:
# Towards an anisotropic atom-atom model for the crystalline phases of the molecular S8 compound

## I Introduction

Our knowledge of the complex phase diagram of elemental sulfur, as well as of the statical, dynamical and thermodynamical properties of all its condensed phases, is far from complete steudel ; meyer1 . Its complexity is due to a rich variety of molecular allotropes (open and cyclic chains of different lengths) that correlate with a large variety of chemical and physical properties in the crystalline and liquid phases steudel ; meyer1 .

The most stable, and the most studied, sulfur allotrope at ambient temperature and pressure (STP) is S8, a crown-shaped cyclic molecule. Three crystalline phases, $`\alpha `$-, $`\beta `$\- and $`\gamma `$-S8, have been reported at ambient pressure and T$`{}_{}{}^{<}{}_{}{}^{}`$400 K. Figs. 1a, b and c show their structures. The orthorhombic $`\alpha `$-S8 can be measured alfa1 ; alfa2  at temperatures up to 385.8 K, its metastable melting point. At 369 K there is a solid-solid $`\alpha `$-$`\beta `$ phase transition, but it is most often observed when the samples are studied from high to low temperatures. Monoclinic $`\beta `$-S8 beta2  is the most stable phase in the 369-393 K temperature range. $`\beta `$-S8 crystals are usually obtained from the melt steudel  and they show an orientational order-disorder phase transition at 198 K beta1 . A third crystalline allotrope $`\gamma `$-S8 gam1  (Fig. 1c) has been observed, with a density, at STP, 5.8 % higher than that of $`\alpha `$-S8.

Most theoretical studies of the S8 crystalline phases, up to very recently, have been limited to calculations of the intramolecular mode frequencies for the isolated S8 molecule and their splitting when they are embedded in a $`\alpha `$-S8 matrix, including calculations of their infrared and Raman intensities alfafreq . Constant-volume molecular dynamics (MD) simulations with the same aim have been also performed cardini1 .

We recently reported claudio1 ; claudio2  a series of constant pressure-constant temperature classical MD simulations of $`\alpha ,`$ $`\beta `$ and $`\gamma S_8`$ crystals using an extremely simple and flexible molecular model, a slight modification of the one proposed in the constant volume MD simulations of $`\alpha `$-S8 cardini1 . The simple intermolecular potential was of the Lennard-Jones (LJ) atom-atom type, except that its parameters were fitted in ref. claudio1,  to the experimental data on the STP crystalline structure of $`\alpha `$-S8 and an estimated value of its configurational energy, calculated as in refs. claudio1, ; claudio2,  and based on the experimental STP value of the heat of sublimation of $`\alpha `$-S8 hsubalfa  (it is also explained in the Calculations section of this paper). These LJ parameters ($`\epsilon `$iso= 1.70 kJ/mol and $`\sigma `$iso=3.39 Å) were afterwards used throughout the calculations of all samples of $`\alpha ,`$ $`\beta `$ and $`\gamma S_8`$ crystals at different temperatures and zero kbar pressure claudio1 ; claudio2 . The intramolecular potential included a harmonic well term for the SSS angle and a torsional potential of the double well type for the sulfur chain, different from hydrocarbon chains that show a triple well torsional potential. These intramolecular coordinates fully describe all the low frequency intramolecular modes that mix with lattice modes alfafreq ; cardini1  and can, therefore, be relevant for the lattice stability or in the onset of structural phase transitions.

Our calculations showed that the simple and flexible molecule model succesfully reproduced most experimental data available on the three crystals claudio1 ; claudio2 , including structural and dynamical data of $`\alpha `$\- and $`\gamma `$-S8 at STP, the high temperature orientationally disordered phase of $`\beta `$-S8 at 375K and its low temperature ordered phase. Nevertheless, below 200 K this simple isotropic LJ model did not reproduce the 100 K structure of $`\alpha `$-S8, as determined by a neutron and X-ray measurements alfa3 . Although no structural phase transition was measured for $`\alpha `$-S8 crystals alfa3 , when decreasing the temperature of the sample below 200 K, our orthorhombic sample turned out to be unstable using this isotropic atom-atom potential, and transformed to a monoclinic structure that we labeled $`\alpha `$’-S8 claudio1 ; claudio2  (Fig. 1).

Assuming that the experimental data of ref. alfa3,  do not correspond to a metastable state, the problem can be correlated to the extremely simple model molecule and in this work we propose to take into account the anisotropy of the sulfur atom due to its 3s2 3p4 electronic configuration. This idea is not new, it gave excellent results in similar problems. For example, in the case of alkane hydrocarbon chains it is known that an unique isotropic atom representing the CH2 units is able to reproduce many facts of the high temperature phases alkane2 , but agreement with the measured crystalline structures is attained only when hydrogen atoms are explicitly included alkane1 . An intermediate model was proposed in ref. alkane3, , consisting in displacing an unique isotropic interaction site from the C position to the center of mass of the CH2 unit, this improvement was enough to study the liquid phase of alkanes but not their crystalline phases alkane4 . The ”all atoms” model has afterwards been succesfully applied to many studies, like the conformations of flexible molecules chains , the equation of state of alkanes alkane4  and the condensed phases of glycerol glycerol .

Also in ref. williams,  the parameters of an aspherical atom-atom model were fitted to reproduce 11 crystalline structures of elemental sulfur, the location of the sites that simulate the lone-pair electronic distribution was refined. The calculations were performed in the harmonic approximation, for rigid molecules, and no attempt to include temperature, intramolecular soft modes, lattice stability or vibrational frequencies was performed.

In the following sections we present the models of the anisotropic intermolecular potential that are analized in the present calculations, the method of calculation, the obtained results for $`\alpha `$-, $`\alpha `$’-, $`\beta `$ and $`\gamma S_8`$ crystals and our conclusions.

## II The anisotropic atom-atom potential

For the present calculations we analize two simple anisotropic atom-atom LJ intermolecular potential models, based on the experimental charge distribution of $`\alpha `$-S8 alfa3 , where the location of the sulfur lone pair orbitals were approximately determined. Our first anisotropic case molecule model model of refs. claudio1, ; claudio2,  replaces each isotropic LJ atomic site by four sites at a distance $`d_{ani}`$ from the S atom, two of them located along the SS bonds and the other two are generated with the S4 element of symmetry located at the atomic S site. This multi-site model (that we called ANI-S4) is a simple way to represent the atomic anisotropy.
For the distance $`d_{ani}`$ we considered two values. The first one is $`d_{ani}`$=0.6 Å, the experimental value determined at 100K alfa3 . The second is a reduced value $`d_{ani}`$=0.1 Å, that was included after taking into account that our four sites anisotropic model is extremelly simple and do not include a site at the S atom, making artificially large the anisotropy of the S atom.

The new LJ parameters of these sites were adjusted to fit the experimental structure and the estimated value of the potential energy of $`\alpha `$-S8 at 300K, obtaining ($`\epsilon `$ani=0.225 kJ/mol, $`\sigma `$ani=2.87 Å) for $`d_{ani}`$=0.6 Å, and ($`\epsilon `$ani=0.115kJ/mol, $`\sigma `$ani=3.35 Å) for $`d_{ani}`$=0.1 Å. It has to be taken into account that, should we have taken $`d_{ani}`$=0Å, the LJ parameter’s values should have been: $`\epsilon `$ani= $`\epsilon `$iso/16=0.107kJ/mol and $`\sigma `$ani=$`\sigma `$iso=3.39 Å.

As we show in the following sections, due to the results obtained with the ANI-S4 potential model, for both $`d_{ani}`$ distances, a further an extremely anisotropic atom model was tested, for the sake of completeness. The last model locates the two external sites almost perpendicular to the plane determined by the corresponding S atom and the two nearest S neighbors, at $`d_{ani}`$=0.6 Å from the S atom (ANI model). This approximated location is also suggested by the experimental data of ref. alfa3,  and the MO calculations that they include, on the sulfur containing simple molecule H2S2. As we explain in the next section, the most stable numerical simulations were obtained when the two sites are generated by simplely displacing the two sites at a distance $`d_{ani}`$=0.6 Å from the location of the corresponding atom Sj (j=1,8) and orienting them along the vectors $`{}_{}{}^{}{}_{}{}^{+}`$\[(Sj+2 -Sj+1)+(Sj-2-Sj-1)\]. That is, the two external sites are located in a plane parallel to that formed by the next two nearest SS bonds. This direction is approximately perpendicular to the plane formed by the Sj-1SjSj+1 atoms and this model is considered a case study that represents an extremely large atomic anisotropy. The LJ parameters of the ANI model were also adjusted to fit the experimental structure and estimated potential energy of $`\alpha `$-S8 at STP, obtaining in this case $`\epsilon `$ani=0.15 kJ/mol and $`\sigma `$ani=3.17 Å.

The LJ parameters of each multisite model, all of them adjusted to fit the experimental STP data of $`\alpha `$-S8, were afterwards used to perform the series of MD runs of the above mentioned four crystalline phases, at zero pressure and in the 50K - 450K temperature range.

The intramolecular terms of the potential model, for all molecular models, are essentially equal to those of refs. claudio1 ; claudio2 , with the only difference that the force constants take now a value about 10 % lower, in order to improve the fit of the intramolecular frequencies to the experimental data. The intramolecular potential for the S-S-S bending angles consists of a harmonic well

$$V(\beta )=\frac{1}{2}C_\beta (\beta \beta _0)^2,$$

with a force constant of Cβ/kB=22000 K/rad2 and $`\beta _0`$= 108 deg., kB is the Boltzmann constant. The intramolecular potential for torsion ($`\tau `$) angles is a double well:

$$V(\tau )=A_\tau +B_\tau \mathrm{cos}(\tau )+C_\tau \mathrm{cos}^2(\tau )+D_\tau \mathrm{cos}^3(\tau )\text{,}$$

with $`A_\tau `$/kB=51.473 K, $`B_\tau `$/kB=664.573 K, $`C_\tau `$/kB=2068.092 K and $`D_\tau `$/kB=476.453 K. These parameters describe a double well with minima at $`\tau `$=$`{}_{}{}^{}{}_{}{}^{+}`$ 98.8 deg., and a barrier height of about 8 kJ/mol at $`\tau `$=180 deg. The barrier height at $`\tau `$=0 deg. is of 27 kJ/mol, out of the range of energies explored in these simulations.
The bond lengths S-S are held constant at their mean experimental value of 2.06 Å.

## III Calculations

The method of calculation, as well as the performed series of numerical simulations are equal to those included in refs. claudio1, ; claudio2, , the only change is in the potential model used.

Briefly, the phase diagram and dynamical bulk properties of S8 crystals, as given by these simple molecular models, are calculated at 0 kbar pressure in the (N,P,T) ensemble, by a series of classical MD simulations. The external pressure is considered isotropic and the MD algorithm allows volume and shape fluctuations of the sample, in order to balance the internal stresses with an external isotropic pressure algor2 . The algorithm controls pressure $`via`$ an extended system which includes as extra variables the MD box parameters, periodic boundary conditions are applied. The temperature control of the sample follows the approach of Nosé algor4 ; algor5 , which also includes an external variable. The equations of motion of these flexible molecules are integrated using the Verlet algorithm for the atomic displacements and the Shake algorithm for the constant bond length constraints on each molecule algor1 . The MD algorithm is based in that used in a study of black Newton films bubbles1 , afterwards extended to constant-pressure constant-temperature simulations via an extended Lagrangian newsimul1 ; newsimul2 . The final version is similar to that used in ref. zully-c4f8, . In the present work, as the interaction sites of our anisotropic atom model are not coincident with the S atomic sites, the algorithm employed to translate the forces from massless to massive sites is based on that written for rigid moleculesalgor6 , extended to the case of flexible molecules.

The runs in the (N,P,T) ensemble were performed every 25 or 50 K in the range 50-450 K. At each point of the phase diagram the samples are equilibrated for 20000 to 30000 time steps (each one of 0.01ps.) and measured in the following 20 ps. Near the phase transitions the stabilization times were increased several times.
The exception are the series of runs corresponding to the ANI-S4 model with $`d_{ani}`$=0.6Å, that were performed with a time step of 0.005ps. Nevertheless, these runs, after equilibration of the sample at the desired T and P, and in a free trajectory in the phase space, give a large deviation of the total energy of the system (including the aditional coordinates of the extended system), implying that the numerical algorithm should be improved or that the molecular model do not calculate the sample as a stable one. For this model we only include, for the sake of comparison, the calculated structure and configurational energy that is measured (after the initial transient) in the equilibration trajectory, an algorithm that following ref. algor3,  simulates the system in the canonical ensemble.

The samples of $`\alpha `$-S8 and $`\alpha `$’-S8 crystals consisted of 3\*3\*2 orthorhombic cells (288 molecules). The sample of $`\beta `$-S8 crystals consisted of 4\*4\*4 cells (384 molecules) and that of $`\gamma `$-S8 crystals of 5\*3\*5 cells (300 molecules).

The free energy of all samples is estimated, at each point of the phase diagram, taking into account the contribution of the calculated vibrational modes to the entropy of the sample bornhuang . This simple method is widely used in lattice dynamics calculations: In a first step the calculated vibrational density of states VDS($`\nu `$) is normalized to the total number of degrees of freedom of the MD sample. This normalized function VDS($`\nu `$) gives the number of oscillators within $`\nu `$$`{}_{}{}^{}{}_{}{}^{+}`$$`\mathrm{\Delta }\nu `$ and in a first approximation, by considering a system of harmonic oscillators, the calculation of their contribution to the entropy is straightforward bornhuang . In the same way the internal energy and enthalpy of the system can be derived, the enthalpy being directly comparable to thermodynamic measurements.

## IV Results

Using the anisotropic LJ model ANI-S4, with distances $`d_{ani}`$=0.1 and 0.6Å, the orthorhombic $`\alpha `$-S8 sample does not spontaneously show the structural phase transition to the monoclinic $`\alpha `$’-S8, as was calculated with the isotropic LJ model when the temperature of the orthorhombic sample was decreased below 200 K claudio1 ; claudio2 . Nevertheless, if a new series of simulations are performed on a sample that initially has the structure of the monoclinic $`\alpha `$’-S8, this phase is calculated stable, in the time scale of our simulations (about 200 ps.), with a potential energy lower than that of $`\alpha `$-S8. Fig. 2 shows the calculated configurational energy (the time and sample average of the potential interaction energy per molecule) and volume per molecule as a function of temperature, for the four crystalline samples, and using the ANI-S4 model with $`d_{ani}`$=0.1 Å. Table I includes our results for the ANI-S4 model with $`d_{ani}`$=0.1 and 0.6 Å at selected points of the phase diagram.

A dramatic change is observed when the ANI potential model is used. These results are included in Fig. 3 and they clearly show a large diference with those of Fig. 2. The extremely anisotropic case model ANI gives results that, in general, compare better with some of the experimental data than those included in Fig. 2. In particular, Fig. 3 shows that the low temperature orthorhombic $`\alpha `$-S8 phase is now calculated with a lower configurational energy (and volume per molecule) than the monoclinic $`\alpha `$’-S8 phase suggested by the isotropic and ANI-S4 models. Nevertheless, this phase is calculated as stable with all tested models and this fact implies that $`\alpha `$’-S8 might be a metastable phase of elemental sulfur.
As we have also discussed in ref. claudio1, , the uncertainty in our comparison with the scarce experimental data on $`\beta `$\- and $`\gamma `$-S8 phases is mainly due to contradictory measurements of $`\gamma `$-S8 crystals. Density measurements at 300K imply that their density is 5.8% higher than that of $`\alpha `$-S8 at STP meyer1 ; gam1 . On the other hand the X-ray measurements gam1  determine a crystalline cell with a larger volume per molecule than that of $`\alpha `$-S8, implying a 2.3% lower density. Further measurements of this phase are needed.

The orientationally disordered phase of $`\beta `$-S8 is reproduced by the isotropic claudio1 , the ANI-S4 and ANI anisotropic models. At high temperatures, the molecules at (0,0,0) and (1/2,1/2,1/2) are dynamically disordered. The inertial axis perpendicular to the molecular plane librates around its equilibrium orientation, whilst the molecule reorientates within the plane. We calculate a characteristic reorientational time (within molecular plane) of 9 ps at 375 K with ANI-S4 ($`d_{ani}`$=0.6Å) and of about 50 ps at 375 K with ANI. Although no sharp transition is found for this crystal when decreasing the temperature of the sample, below 150 K both molecules are orientationally ordered and related by a C2 symmetry axes. No measurements on the molecular reorientational times and lattices parameters at diferent temperatures are available for comparison.

It has to be pointed out that, when analyzing the relative stability of two phases we have estimated the Gibbs free energy of our samples by including the contribution of the vibrational modes to the entropy. Fig. 4 compares the data of $`\alpha `$-, $`\alpha `$’-, $`\beta `$\- and $`\gamma `$-S8. At low temperatures the relative stability of the phases is equal to that obtained from the potential energy data. At high temperatures the merge of all phases is observed. No crossings (which identify phase transitions) are found.

Table I compares the experimental and calculated (using the isotropic and anisotropic models) structural data at the starting point of the phase diagram of our series of simulations. The calculated average pressure, not included in the Table is 0 kbar in all cases. For example, for $`\alpha `$-S8 at 300K is P=0.0(5) kbar.

As in refs. claudio1, ; claudio2, , the $`\alpha `$\- $`\beta `$ phase transition was not found in our calculations, in accord with the experimental fact that the transitions are promoted by the disorder existing in the real samples meyer1 , if not, metastable states can be maintained for days. As calculated in refs. claudio1, ; claudio2,  with the isotropic model, also for the anisotropic atom-atom model the solid-liquid phase transition is located at temperatures near the experimental value ($``$400 K), but only for a disordered sample of S8 molecules. This sample is analized with a MD algorithm that allows changes in the MD box volume, but not in its angles, as a function of pressure and temperature. The disorder is generated by initially increasing the temperature of a cubic sample to 500 K and afterwards following its evolution by increasing and decreasing the temperature in the 350-450 K range. The result is totally similar to Fig. 4 of our ref. claudio1,  and it is not included here. This last calculation is only of theoretical interest, and we did not intend to simulate the real liquid phase of elemental sulfur. The experimental data show that the disorder, in the case of the real solid-liquid phase transitions, is generated because the S8 molecules start to dissociate when the temperature of the sample is near the melting point and our calculations should be valid very near the transition, when the fraction of broken molecules is small.

## V Conclusions

It has to be pointed out that, as we are dealing with molecular crystals, as long as no chemical reactions or electronic excitations are involved, the classical approximation is good enough for the calculation of structural and dynamical properties, provided that the inter- and intramolecular potential gives a good description of the real interactions. Moreover, simple model molecules can be very helpul to elucidate the complex phase diagram of the different crystalline allotropes of elemental sulfur, some of them with a large number of molecules per cell and others with molecular reorientational disorder.

In refs. claudio1, ; claudio2,  we verified that, surprisingly, a large amount of crystalline properties can be very nearly reproduced by a simple and isotropic atom - atom LJ model. The most discordant point being the finding of a low temperature monoclinic sample $`\alpha `$’-S8, that has not been observed alfa2  but, nevertheless, might be a metastable phase of this elemental compound.

Assuming that this $`\alpha `$’-S8 phase is, probably, an artifact of the simple isotropic intermolecular potential model, in this paper we analized two anisotropic case molecular models.

Our calculations show that an anisotropic ANI-S4 model with local S4 symmetry calculates the same stable structures and dynamical properties than the isotropic model, in the time scale of our simulations (about 200 ps.), being the isotropic model more sensitive to changes in the thermodynamic parametes P and T than the anisotropic ANI-S4 one.

The results obtained with the ANI case model are, instead, qualitatively very different from the other models. The four crystalline phases are calculated as polymorphic structures, but their relative cohesive energy is changed. With this model molecule the $`\alpha `$-S8 phase has a lower configurational energy than $`\alpha `$’-S8 in all the studied range of temperatures, in agreement with available measurements. The calculated density of $`\gamma `$-S8 crystals is in accord with X-ray measurements, but not with the density data of ref. meyer1, ; gam1, .

In conclusion, because of the amount of experimental data that these two case model molecules nearly reproduce, and due to the fact that the $`\alpha `$’-S8 allotrope is calculated stable (for about 200ps.) with three model potentials, we propose a thorough measurement of the crystalline phases of S8, including different techniques employed to grow the crystals at different temperatures. Additional measurements of the $`\beta `$\- and $`\gamma `$-S8 crystals are also needed, in particular measurements of molecular reorientational times for $`\beta `$-S8 and lattice parameters as a function of temperature and pressure. Accurate measurements of the density of $`\gamma `$-S8 are still lacking.

Detailed calculations of the lone pair orientation in elemental sulfur molecules, including open chains and cyclic molecules of different size, are presently being carried out using the GAMESS (version 2001) gamess  program. Lastly, a further and simple anisotropic model will be developed with which we plan to extend our calculations to the crystalline phases of the molecular allotropes and the liquid phase of elemental sulfur.

###### Acknowledgements.

This work was partially supported by the grant PIP 0859/98 of CONICET and Fundación J. Balseiro-2001. C.P. thanks partial support from FOSDIC.

Table I: Experimental and calculated lattice parameters: a, b, c (Å) and monoclinic angle $`\beta `$(deg.). $`iso`$ denotes the calculations performed with the isotropic atom-atom modelclaudio1 ; claudio2  and ANI those performed with the corresponding anisotropic model. Z is the number of molecules in the cell, V (Å3) is the volume per molecule. See the text for a discussion about the $`\gamma `$-S8 data.

| Phases |  | Z | a(Å) | b(Å) | c(Å) | $`\beta `$(deg.) | V(Å3) |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $`\alpha `$-S8 | Ref. alfa1 | 16 | 10.4646(1) | 12.8660(1) | 24.4860(3) | - | 206.046(98) |
| (T=300 K) | $`iso`$ |  | 10.34(10) | 13.20(3) | 24.32(7) | - | 207.3(1.0) |
| (T=300 K) | ANI-S4 ($`d_{ani}`$=0.1Å) |  | 10.48(15) | 13.01(3) | 24.29(30) | - | 207.2(2) |
| (T=300 K) | ANI-S4 ($`d_{ani}`$=0.6Å) |  | 10.41(3) | 13.08(4) | 24.24(12) | - | 206.8(8) |
| (T=300 K) | ANI ($`d_{ani}`$=0.6Å) |  | 10.24(4) | 12.78(4) | 25.07(3) | - | 205.3(3) |
| $`\beta `$-S8 | Ref. meyer1 | 6 | 10.778 | 10.844 | 10.924 | 95.8 | 211.70 |
| (T=370 K) | $`iso`$ |  | 11.01(11) | 10.75(3) | 10.79(10) | 97.48(7) | 211.3(7) |
| (T=375 K) | ANI-S4 ($`d_{ani}`$=0.1Å) |  | 11.00(1) | 10.80(6) | 10.80(9) | 96.89(4) | 212.3(9) |
| (T=375 K) | ANI-S4 ($`d_{ani}`$=0.6Å) |  | 10.86(1) | 11.07(3) | 10.59(10) | 93.4(6) | 211.0(3) |
| (T=375 K) | ANI ($`d_{ani}`$=0.6Å) |  | 10.61(2) | 10.92(2) | 10.94(2) | 92.7(1) | 212.2(1) |
| $`\gamma `$-S8 | Ref. gam1 | 4 | 8.442(30) | 13.025(10) | 9.356(50) | 124.98(30) | 210.8(1.3) |
| (T=300 K) | $`iso`$ |  | 8.03(13) | 13.16(14) | 8.88(11) | 120(3) | 202.5(2.5) |
| (T=300 K) | ANI-S4 ($`d_{ani}`$=0.1Å) |  | 8.03(10) | 13.16(4) | 8.89(8) | 121(2) | 201.2(9) |
| (T=300 K) | ANI-S4 ($`d_{ani}`$=0.6Å) |  | 8.04(9) | 13.30(13) | 8.96(8) | 121.3(8) | 204.9(3) |
| (T=300 K) | ANI ($`d_{ani}`$=0.6Å) |  | 8.21(4) | 13.05(5) | 9.11(5) | 122.9(4) | 205.2(6) |
| $`\alpha `$’-S8 |  | 16 |  |  |  |  |  |
| (T=100 K) | $`iso`$ |  | 13.18(7) | 9.31(5) | 26.7(1) | 76.1(5) | 198.9(9) |
| (T=100 K) | ANI-S4 ($`d_{ani}`$=0.1Å) |  | 13.08(7) | 9.49(5) | 26.0(1) | 77.1(6) | 196.5(6) |
| (T=100 K) | ANI-S4 ($`d_{ani}`$=0.6Å) |  | 13.29(7) | 9.39(5) | 26.3(3) | 76.4(4) | 199.3(6) |
| (T=100 K) | ANI ($`d_{ani}`$=0.6Å) |  | 13.07(12) | 9.38(10) | 27.0(1) | 78.2(6) | 202.5(9) |

Figure:

Fig. 1: Crystalline phases of the molecular S8 compound. From left to right: a) orthorhombic $`\alpha `$-S8 at 300 K, b) monoclinic $`\beta `$-S8 at 150 K (with the high temperature orientationally disordered molecules in light colour), c) pseudo-hexagonal $`\gamma `$-S8 at 300 K, d) other view of the orthorhombic $`\alpha `$-S8 at 300 K, e) calculated monoclinic $`\alpha `$’-S8, at 100 K, with the isotropic LJ model claudio1 ; claudio2 .

Fig. 2: Configurational energy Uconf (kJ/mol) and volume per molecule (Å3) at 0 kbar and as a function of temperature T (K), as given by the anisotropic ANI-S4 model, with $`d_{ani}`$=0.1 Å. The lines are a guide to the eyes.

Fig. 3: Configurational energy Uconf (kJ/mol) and volume per molecule (Å3) at 0 kbar and as a function of temperature T (K), as given by the anisotropic ANI model, with $`d_{ani}`$=0.6 Å. The lines are a guide to the eyes.

Fig. 4: Gibbs free energy (kJ/mol) at 0 kbar and as a function of temperature T (K), as given by the anisotropic ANI model, with $`d_{ani}`$=0.6 Å. The lines are a guide to the eyes.


format: text

text:
# Molecular dynamics of flows in the Knudsen regime

In a landmark paper , Maxwell applied his kinetic theory  to the study of the slip length associated with fluid flow close to a solid surface. He postulated a linear combination of two kinds of behaviour when a molecule was incident on a wall – it is either specularly reflected or it undergoes enough collisions with the wall molecules so that it loses memory of the incident velocity and emerges with a velocity characteristic of the wall temperature. We revisit this problem with molecular dynamics simulations – as the properties of the wall are varied, we find both kinds of behaviour suggested by Maxwell, but the crossover from one to the other is not as simple as envisioned by him. In addition, rich behaviour is found in the fluid flow properties on varying the fluid density from the dense liquid regime to the Knudsen regime (large mean free path compared to the system size). Our studies of flows in the subcontinuum regime provide a unique window on the nature of boundary conditions at a fluid - solid interface and are relevant to the behaviour of nanoscale mechanical devices, high altitude aerodynamics, aerosol dynamics, isotope separation schemes, flow cooling techniques, and energy transfer in molecular collisions .

Studies of flows at nanoscales involve questions pertaining to a sub-continuum regime  for which experiments are often unable to provide the requisite answers. The behaviour of materials at the nanoscale is often quite different from that in the bulk . Numerical simulations based on molecular dynamics (MD) offer a natural tool to handle the interacting many-body nanoscale physics. A detailed understanding of the behaviour of materials at very small scales is necessary for the fabrication of miniature devices and for the manipulation of materials at the molecular scale. Such microscopic studies are useful for assessing the range of validity of traditional continuum theories. In addition, atomistic simulations are an invaluable tool for deducing constitutive relationships, which may then be fed into continuum calculations.

Here, we report on MD studies of the boundary conditions in single-fluid flows in the Knudsen regime  in which the mean free path, $`\lambda `$, becomes comparable to the characteristic device dimension, $`L`$. Such flows behave very differently from bulk, continuum expectations. There is considerable work on the modified boundary conditions based on either semi-heuristic kinetic theory arguments or using an analysis of the Boltzmann equation, but inevitably severe simplifications are made concerning the dynamics of the system . Furthermore, in such treatments, the solid boundaries are not handled in a realistic manner as having atomic structure on length scales comparable to that of the fluid constitutents .

We study Lennard-Jones fluids undergoing two kinds of flow, Couette and Poiseuille, between two solid walls with an atomic structure. Three kinds of walls, two attractive and one repulsive, are discussed. We consider a range of fluid densities to study the crossover from flows described by continuum fluid mechanics to dilute gas flows in the Knudsen regime.

In a Lennard-Jones fluid, two atoms separated by a distance $`r`$ interact with the potential

$$V_{LJ}(r)=\mathrm{\hspace{0.33em}4}ϵ\left[(\frac{r}{\sigma })^{12}(\frac{r}{\sigma })^6\right],$$
(1)
where $`\sigma `$ is the size of the repulsive core \[the potential is truncated at $`2.2\sigma `$\]. The flows take place between two solid planar walls parallel to the $`xy`$ plane with periodic boundary conditions imposed along the $`x`$ and $`y`$ directions. Couette flow was generated by moving the walls at constant velocity in opposite directions along the $`x`$-axis whereas Poiseuille flow was induced by an application of uniform acceleration $`g`$ on all the fluid molecules in the $`x`$-direction.

Following Thompson and Robbins , the wall of atoms form two  planes of an fcc lattice with each wall atom tethered to a fixed lattice site by a harmonic spring of spring constant $`k`$. Our studies where done for $`k`$=400 so that the mean-square displacement about the lattice sites did not exceed the Lindemann criterion for melting but was close to the value set by this criterion. The wall fluid interactions were modelled by another Lennard-Jones potential,

$$V_{wf}(r)=\mathrm{\hspace{0.33em}16}ϵ\left[(\frac{r}{\sigma })^{12}A(\frac{r}{\sigma })^6\right],$$
(2)
with the parameter $`A`$ varying between 1 and 0, corresponding to attractive and repulsive walls respectively. For the narrowest channel that we studied, the fluid atoms were confined to a space measuring 13.6$`\sigma `$ by 5.1$`\sigma `$ in the $`xy`$ plane and 12.75$`\sigma `$ between the walls. The “basic” walls in our studies contain 96 wall atoms each (in each wall plane, the atomic periodicity is 0.85$`\sigma `$). We have also considered the case of walls which are about 2.5 times denser and consist of 176 wall atoms each – we shall refer to these as “dense” walls.

The equations of motion were integrated using a fifth-order Gear predictor-corrector algorithm  with a step size of $`0.005\tau `$, where $`\tau =(m\sigma ^2/ϵ)^{1/2}`$ and $`m`$ is the atomic (wall and fluid) mass. The fluid and wall temperatures were fixed at $`k_BT/ϵ`$=1.1. We present results obtained with the use of Langevin thermostat  but the results obtained with the use of the Nose-Hoover thermostat  were virtually identical. The noise was applied in all directions during equilibration and only in the $`y`$-direction during the data-taking phase. The averages were obtained over runs which lasted at least 5000 $`\tau `$, and in the extreme dilution case, for up to 400000 $`\tau `$. An initial time of at least 400$`\tau `$ was spent for equilibration. The spatial averaging was carried out in slabs of width $`\sigma /4`$ along the $`z`$-axis.

In the dilute gas limit, the Knudsen number, $`Kn`$, is defined as $`\lambda /L`$, i.e., as $`(\sqrt{2}\pi \rho \sigma ^2L)^1`$, where $`\lambda `$ is the mean free path, $`\rho `$ is the fluid number density and $`L`$ is the characteristic width of the channel. A study of the density profile shows qualitatively different behaviours of the density profiles for the attractive and repulsive walls. In the former case, two adsorbed layers form near both walls whereas there is merely a depletion zone, whose thickness is of order $`\sigma `$, alongside the repulsive walls. The first layer at the attractive walls is fairly periodic, as measured by the static structure factor, and the resulting local order corresponds to either a square or hexagonal lattice depending on whether the wall is basic or dense. $`L`$ is defined as an effective channel width in which the flow takes place (in the attractive case this is the distance between centers of the second layers). In the center of the channel, the fluid density is essentially constant and this value is chosen as $`\rho `$ in the definition of $`Kn`$.

Figure 1 shows the Poiseuille flow velocity profiles for three different values of $`\rho `$ (or $`Kn`$) for the attractive wall case. We have confirmed with studies of channels of three different widths and the two wall densities that the crossover value of $`\rho `$ above which continuum fluid mechanics holds, $`\rho _K`$, is somewhat less than 0.19, corresponding to a $`Kn`$ of order 0.1 for the channel of the smallest width. In this regime, the profiles are clearly parabolic and essentially with no slip at the walls. In contrast, the velocity profile at the highest $`Kn`$ that we were able to study (corresponding to just 4 particles in the interior of the channel) is almost that of plug flow with a huge slip length. For repulsive walls, the plug flow regime extends to higher densities because such walls provide little coupling to the fluid and correspondingly, $`\rho _K`$ becomes of order 0.3. Furthermore, in the continuum regime, the velocity field exhibits substantial slip. A consistent behaviour is found for Couette flow with higher slip for the repulsive walls (compared to the attractive wall case) for $`\rho >\rho _K`$ and essentially zero fluid velocity at the walls at the highest $`Kn`$.

We now proceed to study the effect of interactions between the fluid particles undergoing Poiseuille flow on increasing $`\rho `$ within the channel. In Poiseuille flow, the velocity field is highest in the middle of the channel and is denoted by $`v_{\mathrm{max}}`$. At the lowest densities, one expects that $`v_{\mathrm{max}}`$ ought to scale as $`L`$. This is because the driving force is effective in accelerating the fluid particles only during the time spent traversing the channel between collisions with the walls (and the associated adsorbed layers). This time simply scales as the channel width. On the other hand, at large $`\rho `$, corresponding to usual fluid behaviour, $`v_{\mathrm{max}}`$ ought to scale as $`L^2`$ as given by the Navier-Stokes equation with no-slip boundary conditions. Figure 2 shows a plot of how $`v_{\mathrm{max}}`$ interpolates between the low and high Knudsen number behaviours as the fluid density is varied. Strikingly, $`v_{\mathrm{max}}`$ shows a maximum around $`\rho _K`$ with the peak being more pronounced for larger $`L`$’s. We have confirmed that this qualitative behaviour is essentially independent of the nature of the walls. Physically, $`\rho _K`$ corresponds to the situation when the fluid intermolecular interactions begin to dominate leading to a diminished influence of the wall in the interior of the fluid.

Figure 3 shows a vivid illustration of the role played by the intermolecular interactions in leading to an enhancement of $`v_{\mathrm{max}}`$ even as the fluid density is moderately increased from the high $`Kn`$ limit. For a repulsive wall (not shown), the particle is reflected approximately specularly and moves out of the vicinity of the wall immediately. The bottom panel shows qualitatively different behaviour in the (basic) attractive wall case. Here, again as suggested by Maxwell, a particle either thermalizes in the vicinity of the wall for a relatively long time and then escapes from it or occasionally still undergoes specular reflection. On increasing the fluid density, one gets the situation illustrated in the middle panel in which intermolecular collisions keep the particles longer within the channel thus reducing the time of localization near the walls and effectively increasing the velocity field in the middle of the channel (Figure 2). At still higher densities ($`\rho _K`$), the intermolecular interactions become sufficiently large so that liquid-like behaviour sets in and $`v_{\mathrm{max}}`$ begins to decrease.

A key quantity that characterizes the fluid solid interface is the collision kernel  – the probability density that a molecule striking the surface with a given velocity reemerges with a specific velocity after a given residence time. For the repulsive wall, we find that Maxwell’s specular reflection scenario holds exceedingly well and the residence time is negligible. On the other hand, for the attractive walls the distribution of residence times develops a long tail. Strikingly, the velocity distribution of a particle after the collision with the wall is substantially independent of the incoming velocity. As an illustration, Figure 4 shows the case of particles with unit speed incident normal to the wall. The inset shows the probability distribution of residence times whereas the main figure shows the distributions of the normal ($`z`$) and tangential ($`x`$) velocities of the scattered molecules. The results are in excellent accord with Maxwell’s thermal wall scenario in which the probability densities corresponding to $`v_x`$ and $`v_z`$ are

$$\varphi _x(v_x)=\sqrt{\frac{m}{2\pi k_BT}}exp(\frac{mv_x^2}{2k_BT})$$
(3)
and

$$\varphi _z(v_z)=\frac{m}{k_BT}v_zexp(\frac{mv_z^2}{2k_BT})$$
(4)
respectively. As expected for a thermal wall, the residence time is found to be uncorrelated with the outgoing velocity.

How does the behaviour crossover from that of specular reflection to that of a thermal wall? Following Maxwell1, a common hypothesis is a linear combination of both behaviours with weights $`(1f)`$ and $`f`$ respectively. Such an intermediate case is realized, e.g., on considering the scattering from a wall with the wall-fluid potential given by Eq. 2 with $`A`$=1/4. The velocity distributions of the scattered particles clearly deviate from pure thermal behaviour and the data might be interpreted as containing a specular component. However, the resulting distributions are more complex than the simple linear combination form . The behaviour of the collision kernel is further enriched on considering the role of the interactions between the fluid molecules as the Knudsen number decreases to a value of 0.27 corresponding to a density of $`\rho `$ = 0.066 (figure 5). The influence of other fluid molecules leads to the distributions more closely approaching the thermal ones. The nature of the boundary conditions at a fluid-solid interface is thus complex but akin to those envisioned by Maxwell.

We are indebted to Mark Robbins and Riina Tehver for many valuable discussions. This work was supported by KBN (Grant No. 2P03B-025-13), NASA, and the Petroleum Research Fund administered by the American Chemical Society.

Figure captions

FIG. 1. Velocity profiles of Poiseuille flow with $`g=0.01\sigma /\tau ^2`$ for the smallest width channel with basic attractive walls. The values of the interior fluid density and the corresponding Knudsen number are indicated. The filled circles denote the highest $`Kn`$ case. The inner wall layers are at the edges of the figure. Effectively, one observes no slip boundary conditions for the two higher density cases (the slip length is around 1.7$`\sigma `$ measured with respect to the wall location, but note the presence of the adsorbed layers at the wall).

FIG. 2. Plot of the fluid velocity in the middle of the channel in Poiseuille flow with $`g=0.01\tau /\sigma ^2`$ as a function of fluid density (and the Knudsen number in the inset). The plot is for three different channel widths. The $`L_0`$ in the figure denotes the narrowest channel. For this case, both attractive walls were considered – the basic (hexagons) and higher (asterisks) wall densities. $`2L_0`$ and $`4L_0`$ denote channels that are effectively twice and four times as wide as the smallest channel.

FIG. 3. Plot of the time dependence of the $`z`$-coordinate of a typical fluid particle undergoing Poiseuille flow in the narrowest channel. The horizontal dotted lines indicate the $`z`$-coordinates of the wall molecules. The solid horizontal line refers to the centers of the adsorbed first layers.

FIG. 4. The probability density distributions for the $`z`$ (normal) and $`x`$ components of the outgoing velocity of a fluid molecule after collision with the wall and the associated adsorbed layer. The solid histograms are obtained based on 3000 scattering events. The dotted lines are the analytic predictions of Maxwell, eqs. 3 and 4, discretized with a bin size of 0.1 $`\sigma /\tau _0`$ as in the simulations. The inset shows the distribution of residence times. The residence time is measured as the time interval between a fluid molecule leaving a force-free region (near the wall) and its reemergence. The peak of the curve (at a non-zero value of $`\tau `$) corresponds to a virtually instantaneous bounce back of the fluid molecule and accounts for the time needed to approach the wall and return.


text:
# P-Vortices, Gauge Copies, and Lattice Size

## 1 Introduction

In the last few years there has been renewed interest, within the lattice gauge theory community, in the center vortex theory of confinement . The revival of this old idea is due to a number of numerical studies, which all indicate that center vortices are ubiquitous in the QCD vacuum and give rise to the linear confining potential . Recently, however, a paper by Bornyakov, Komarov, Polikarpov, and Veselov (BKPV)  appeared which questions the validity of a procedure, known as center projection in maximal center gauge, which was essential in many of these studies. The claim is that when large numbers of random gauge copies are used in fixing the gauge (in an effort to minimize the Gribov copy problem), the center-projected string tension underestimates the full string tension by as much as $`30\%`$. Our purpose in the present article is to show that, while BKPV have certainly raised an important issue, their actual data was strongly affected by the rather small lattice volumes used in the numerical simulations. We will show that the conclusions drawn from this data are not sustained by simulations on larger lattices, which generally agree with our previously reported results .

In the direct version of maximal center gauge , the procedure is to maximize

$$R=\underset{\mu }{}\underset{x}{}\text{Tr}_A[U_\mu (x)]$$
(1)
by an iterative over-relaxation procedure, where $`\text{Tr}_A[U]`$ is the trace of $`U`$ in the adjoint representation. Let $`R_n`$ denote the value of $`R`$ after $`n`$ over-relaxation sweeps. When $`R_n`$ is judged to have converged, e.g. according to a criterion of the form

$$\frac{R_nR_{n50}}{R_n}<\delta $$
(2)
then the link variables $`U_\mu (x)`$ are projected onto the nearest center element $`Z`$. In SU(2) lattice gauge theory, the projection is simply

$$Z_\mu (x)=\mathrm{signTr}[U_\mu (x)]$$
(3)
Center-projected Wilson loops, Creutz ratios, etc. are observables computed from the center-projected link variables. It was shown in a number of studies (see, in particular, ref. ), that

* thin vortex excitations of the projected lattice, known as “P-vortices”, are located roughly in the middle of thick center vortices on the unprojected lattice;
* projected Creutz ratios $`\chi _{cp}(I,I)`$ are close to the asymptotic string tension on the unprojected lattice (“center dominance”);
* the density of P-vortices, from $`\beta =2.3`$ onward, scales according to asymptotic freedom;
* removing center vortices (located via the projected lattice) from unprojected lattices also removes confinement and chiral symmetry breaking, and brings the topological charge on the lattice to zero .

Because these numerical results have such strong implications for the QCD confinement mechanism, it is important that the center projection procedure be examined critically. The first obvious question $``$ why should this procedure work at all? $``$ was addressed in ref. . There it was shown that in the absence of Gribov copies (i.e. if the gauge can be fixed to a global maximum of $`R`$), then center projection in maximal center gauge will always locate a thin vortex inserted anywhere on the lattice. This was dubbed the “vortex-finding property” of maximal center gauge, and is certainly a necessary condition for its success. However, the vortices found in the QCD vacuum are not thin vortices, but are necessarily of finite thickness in physical units. Moreover, maximal center gauge is plagued with Gribov copies, since the over-relaxation scheme converges only to a local maximum of $`R`$, which will be slightly different for every gauge copy of a given lattice configuration.

In view of the Gribov copy problem, and the finite thickness of vortices, one must rely on empirical checks of vortex-finding via center projection. In this article we will study the sensitivity of center-projected Creutz ratios and the P-vortex density with respect to: (i) the number $`N_{copy}`$ of random gauge copies used for maximal center gauge fixing; (ii) the lattice size; and (iii) the convergence parameter $`\delta `$ in eq. (2). We will here only discuss center projection in direct maximal center gauge, since this is the case treated in ref. , and it is also the center gauge with which we have the most experience. It should be noted, however, that alternatives to maximal center gauge do exist; these include Laplacian center gauge  (which is free of the Gribov copy problem), as well as two other recent proposals .

## 2 $`𝑵_{𝒄𝒐𝒑𝒚}`$ and Lattice Size Dependence

We begin with the $`N_{copy}`$ dependence, whose importance was recently emphasized by BKPV . As noted above, when the over-relaxation gauge-fixing procedure is applied to different gauge copies of a given lattice configuration, slightly different values of $`R`$ are obtained. One way to minimize the gauge copy dependence is to carry out the over-relaxation procedure on a number $`N_{copy}`$ of random gauge copies, perform center projection on the copy with the largest value of $`R`$, and evaluate observables. Data obtained in this way, over a range of $`N_{copy}`$ values, can then be extrapolated to the $`N_{copy}\mathrm{}`$ limit.111P-vortex positions extracted from different Gribov copies are closely correlated, as one would expect if they locate physical objects . It is possible, however, by choosing a very special starting configuration with links gauge-fixed to have mainly positive trace, to destroy the vortex-finding property of center projection even for thin vortices . This is one aspect of the Gribov copy problem in maximal center gauge. Since the vortex-finding property is essential to center projection, it is evident that any unusual variation of the iterative gauge-fixing procedure, which destroys this property, must be avoided. In the original simulations of ref. , only three gauge copies were used, and there was no attempt to extrapolate to $`N_{copy}\mathrm{}`$. Projected Creutz ratios $`\chi _{cp}(I,I)`$ were found to be close to the asymptotic string tensions reported by Bali et al. , for all $`I2`$. BKPV, however, calculate Creutz ratios in the range $`N_{copy}[1,20]`$, and extrapolate to $`N_{copy}\mathrm{}`$ by fitting their data to the functional form

$$\chi _{cp}^{N_{copy}}(I,I)=\chi _{cp}(I,I)+\frac{c(I,I)}{\sqrt{N_{copy}}}$$
(4)
The result reported by BKPV is that projected string tensions, at $`\beta =2.4,2.5`$, underestimate the full string tension by about $`20\%`$ at $`N_{copy}=20`$, and by as much as $`30\%`$ in the extrapolation to $`N_{copy}\mathrm{}`$. This result suggests that the center dominance previously reported in maximal center gauge was a numerical accident, a result of using too few gauge copies.

However, apart from the number of gauge copies, there is one other notable difference between the BKPV simulations and previous work. This is the matter of lattice size. Data reported in ref.  was obtained at $`\beta =2.3`$ and $`\beta =2.4`$ on $`16^4`$ lattices, and at $`\beta =2.5`$ on a $`22^4`$ lattice. BKPV, on the other hand, used only $`12^4`$ lattices at $`\beta =2.3`$ and $`\beta =2.4`$, and lattice size $`16^4`$ at $`\beta =2.5`$. This raises the question of whether the BKPV results, obtained on the smaller lattices, were seriously contaminated by finite-size effects. To find out, we have repeated the center-projection calculation at $`\beta =2.3`$ and $`\beta =2.5`$ on a variety of lattice sizes, for $`N_{copy}[1,20]`$. For the convergence parameter in eq. (2), we have used $`\delta =2\times 10^7`$.

In Fig. 1 we display results for the Creutz ratio $`\chi _{cp}^{N_{copy}}(4,4)`$ vs. $`N_{copy}`$ at $`\beta =2.5`$, for lattice sizes ranging from $`8^4`$ to $`28^4`$. Two features of this data are immediately apparent. First, there is indeed a slow downward trend in the Creutz ratio as $`N_{copy}`$ increases, as noted by BKPV, but this effect is much more pronounced on smaller lattices than on larger lattices. Second, although Creutz ratios on the smaller lattices grossly underestimate the full string tension, the data appears to steadily increase towards the full asymptotic string tension, reported by Bali et al. , as the lattice size increases. These trends in the data are by no means unique to the particular Creutz ratio $`\chi _{cp}(4,4)`$ at $`\beta =2.5`$, but are typical of all of our results. For completeness we display, in Figs. 2 and 3, some other projected Creutz ratios $`\chi _{cp}^{N_{copy}}(I,I)`$ for $`I`$ in the range $`I=25`$, at couplings $`\beta =2.3`$ and $`\beta =2.5`$. In Figs. 2-6, solid (dashed) lines indicate the value (errorbar) of the asymptotic string tension on the unprojected lattice, reported in ref. .

In Fig. 4 we show the projected Creutz ratios $`\chi _{cp}^{N_{copy}}(I,I)`$ for $`N_{copy}=5,10,15,20`$ on the largest lattices we have used: $`20^4`$ at $`\beta =2.3`$ and $`\beta =2.4`$, and $`28^4`$ at $`\beta =2.5`$. We also show the values of these Creutz ratios extrapolated to the $`N_{copy}\mathrm{}`$ limit, using the fitting function (4) suggested by BKPV. As usual in maximal center gauge, all the $`\chi _{cp}(I,I)`$ for $`I2`$ are close to the asymptotic string tension, and these latest results are not far from our earlier results reported in ref. . As another way of showing lattice size dependence, we take the average of the projected ($`N_{copy}\mathrm{}`$) Creutz ratios $`\chi _{cp}(I,I)`$ in the range $`I=25`$

$$\chi _{av}\frac{1}{4}\underset{I=2}{\overset{5}{}}\chi _{cp}(I,I)$$
(5)
Figure 5 shows the lattice size dependence of $`\chi _{av}`$. Note the approach of $`\chi _{av}`$ to the full asymptotic string tension ($`\sigma _{SU(2)}=0.035`$ at $`\beta =2.5`$) as lattice size increases.

### 2.1 Gauge-Fixing Convergence Criterion

In addition to lattice volume and $`N_{copy}`$ dependence, it is also worthwhile to check that the numerical results are stable when the gauge-fixing convergence criterion is strengthened, i.e. when the constant $`\delta `$ in eq. (2) is reduced. When $`\delta `$ is chosen too large, the center projected Creutz ratios come out significantly too high. This is illustrated in Fig. 6, which shows results for projected Creutz ratios with convergence criteria $`\delta =10^2,10^3,10^4,2\times 10^7`$ at $`\beta =2.5`$ ($`24^4`$ lattice). The weakest convergence criterion, corresponding to $`\delta =10^2`$, is clearly insufficient for accurate results, but Creutz ratios obtained with the two smallest values of $`\delta `$ are fairly consistent, indicating that these numbers are not far from the $`\delta 0`$ limit.

## 3 Vortex Density

Next we consider the scaling of the vortex density. The lattice P-vortex density $`p`$ is the total number of P-vortex plaquettes (i.e. plaquettes on the projected lattice with $`ZZZZ=1`$), divided by the total number of plaquettes on the lattice. This quantity is proportional to the average area taken up by P-vortices per unit lattice volume, and is determined from the center-projected plaquette expectation value

$$p=\frac{1}{2}(1W_{cp}[1,1])$$
(6)
If this quantity scales as predicted by asymptotic freedom, then we would have $`p=p_{af}`$, where

$$p_{af}=\frac{\rho }{6\mathrm{\Lambda }^2}F(\beta )$$
(7)
Here $`\rho `$ is the vortex density (average vortex area per unit volume) in physical units of inverse area, and

$$F(\beta )=\left(\frac{6\pi ^2}{11}\beta \right)^{102/121}\mathrm{exp}\left[\frac{6\pi ^2}{11}\beta \right]$$
(8)
In Fig. 7 we plot the lattice vortex density rescaled by the asymptotic freedom expression

$$\stackrel{~}{p}\frac{p}{F(\beta )}$$
(9)
which should be constant in the large $`\beta `$ limit, if $`p`$ scales according to asymptotic freedom. There appears to be good evidence for this kind of scaling, already for $`\beta 2.2`$, in agreement with previous results . In Fig. 7 and Table 1, the vortex densities at $`\beta =2.3,2.4,2.5`$ are taken from the largest lattices ($`20^4`$, $`20^4`$ and $`28^4`$, respectively), and extrapolated to the $`N_{copy}\mathrm{}`$ limit by a fit to

$$p^{N_{copy}}=p+\frac{c}{\sqrt{N_{copy}}}$$
(10)
The vortex densities at other values of $`\beta `$ (with $`N_{copy}=3`$) are just taken from our previous work. For comparison, in Table 1, we display our values for $`\stackrel{~}{p}`$, and the values of the rescaled asymptotic string tension

$$\stackrel{~}{\sigma }\sigma _{SU(2)}/F(\beta )$$
(11)
where $`\sigma _{SU(2)}`$ is the string tension (in lattice units) on the unprojected lattice, at the given $`\beta `$ value. It is interesting to note that the scaling of vortex density $`p`$, in the range $`\beta =2.32.5`$, is substantially better than the the scaling of the full asymptotic string tension $`\sigma _{SU(2)}`$ in this range.

## 4 Vortex Thickness

Why are such relatively large lattices required, in order to avoid large finite size effects in center projected Creutz ratios and vortex densities? We believe that the relevant length scale here is associated with the vortex thickness. Center vortices are surface-like objects, in $`D=4`$ dimensions, which have a finite thickness in physical units. The thick vortex surface (or “core”) bounds a Dirac 3-volume, which represents the region of discontinuity of a singular gauge transformation associated with the vortex. In ref.  we have explained the vortex-finding property of maximal center gauge in terms of the global properties of this singular gauge transformation. On a small lattice, however, with an extension comparable to the vortex thickness, these global aspects of the vortex field may be almost absent, and the minimal Dirac 3-volume could be quite small. In that case, the argument of ref.  breaks down. For this reason, we expect center projection to be less effective at finding vortices on small lattices, leading to underestimates of both the vortex density and the projected Creutz ratios.

Assuming there is some truth in this explanation, center projection can only be accurate for lattices whose extension is large compared to the vortex thickness. There are three independent ways of estimating the thickness of center vortices, which can be deduced from either

1. the ratio of “vortex-limited” Wilson loops ;
2. the vortex free energy as a function of lattice size ;
3. the adjoint string-breaking length .

Vortex-limited Wilson loops are defined in the following way: $`W_n(C)`$ is a Wilson loop evaluated on a sub-ensemble of unprojected configurations, selected so that precisely $`n`$ P-vortices, in the corresponding center-projected configurations, pierce the minimal area of the loop. We can further make the restriction, for $`W_1(C)`$, that the negative P-vortex plaquette lies at (or touches) the center point of the rectangular loop. It is then expected that

$$\frac{W_1(C)}{W_0(C)}1$$
(12)
in the limit where the vortex core is entirely contained within the loop (cf. ref.  for a more extended discussion). In Fig. 8 we show the data for $`W_1/W_0`$ vs. loop area at $`\beta =2.3`$, taken from our previous work in ref. . Judging from this figure, the vortex appears to almost fit inside a $`5\times 5`$ loop, which leads to a rough estimate of the vortex radius, as it pierces a plane, of about 3 lattice spacings. At $`\beta =2.3`$ we have $`\sigma a^2=0.135`$, and taking $`\sigma =5\text{fm}^2`$, the lattice spacing is $`a=0.164`$ fm. A diameter of 6 lattice spacings then corresponds to a vortex thickness of about one fermi.

A second estimate is obtained from the recent calculation of vortex free energy vs. lattice size, carried out numerically by Kovács and Tomboulis . The vortex free energy is close to zero when the lattice extension is greater than the vortex thickness, and this again gives an estimate for the vortex thickness of a little over one fermi. Finally, if confinement is due to center vortices, then an $`R\times T`$ Wilson loop in the adjoint representation must change from a (Casimir scaling) area-law falloff to a (color-screening) perimeter-law falloff for charge separation $`R`$ greater than the vortex thickness . The adjoint string-breaking distance has been measured, by de Forcrand and Philipsen , to be $`1.25`$ fm, and this distance provides us with a third estimate of the vortex thickness, which is roughly consistent with the other two.

At $`\beta =2.5`$, one fermi corresponds to 12 lattice spacings. The lattice used by BKPV at $`\beta =2.5`$ was only 16 lattice spacings across, and this may simply be inadequate for center projection to reliably identify vortices, in view of the above estimates for the vortex thickness. In fact, the P-vortex density is quite low on small lattices, increasing sharply up to $`L=16`$ lattice spacings at $`\beta =2.5`$, where it begins to level off. This is illustrated in Fig. 9; the data points are the extrapolated values for $`p`$ at $`N_{copy}\mathrm{}`$.

In order to estimate the average distance between vortices, we need to know the vortex density. The P-vortex density, discussed in the previous section, is in fact an overestimate of the actual center vortex density . The reason is that P-vortices fluctuate within the thick vortex core, and P-vortex locations, while correlated among random Gribov copies, do vary somewhat from one random copy to another . Thus, although P-vortices are certainly near the middle of center vortices, they are unlikely to be exactly in the middle. A more accurate estimate of the center vortex density is arrived at by either “smoothing” the P-vortex surfaces, as explained in ref. , or else by taking the naive estimate of the number of vortices piercing a plane (per unit area in lattice units), extracted directly from the string tension (see, e.g., ref. )

$$f=\frac{1}{2}(1e^\sigma )$$
(13)
The two estimates agree fairly well , and for $`\beta =2.3`$ we find $`f=0.063`$. This implies an average distance of $`f^{1/2}4`$ lattice spacings between the centers of vortex cores piercing a plane. Since we have already estimated the vortex thickness at $`\beta =2.3`$ to be about 6 lattice spacings, its clear that there must be a substantial overlap between vortex cores, even on a very large lattice. There is nothing in principle wrong with that; vortex cores are not impenetrable objects, and their long-range effects are associated with Dirac 3-volumes, rather than the detailed structure of the core. Our findings for vortex thickness and separation simply indicate, in accordance with some old ideas of the Copenhagen group (Nielsen, Olesen, and Ambjørn in ref. ), that the QCD vacuum is more like a liquid of vortices than a dilute gas. Perhaps this is natural for objects which are both condensed and disordered.

## 5 Conclusions

We have found that center-projected lattices are more sensitive to finite size effects than are unprojected lattices, and precision results for center-projected Creutz ratios require lattice sizes which are large compared to one fermi. If accurate numbers (rather than just qualitative results) are required, then center-projected data from numerical simulations must be tested for convergence with respect to: (i) increasing the lattice size; (ii) increasing the number $`N_{copy}`$ of gauge copies used for center gauge fixing, and (iii) strengthening the gauge-fixing convergence criterion.

Our work was stimulated by the recent findings of BKPV , who have demonstrated a significant $`N_{copy}`$-dependence for center-projected Creutz ratios obtained on relatively small lattice volumes. In the large $`N_{copy}`$ limit, a large deviation ($`30\%`$) was found between the string tension on unprojected and projected lattices. Our result in the present article is that this $`N_{copy}`$ dependence is greatly reduced as lattice size increases, and center projected Creutz ratios $`\chi _{cp}(I,I)`$ appear to converge to values which are quite close to the asymptotic string tension obtained on unprojected lattices.

It is difficult to say whether the full and center-projected string tensions are in precise agreement in the infinite volume limit. To check such agreement, what is really required is an extrapolation to a triple limit: $`\text{volume}\mathrm{}`$, $`N_{copy}\mathrm{}`$, and $`\delta 0`$; our data is not yet adequate to extrapolate to this triple limit systematically. But we emphasize again that on the largest lattices we have used, and with results extrapolated (following BKPV) to the $`N_{copy}\mathrm{}`$ limit, the center projected Creutz ratios lie quite close to the asymptotic string tension, as illustrated in Fig. 4 above. Moreover, on the same large lattices and $`N_{copy}\mathrm{}`$ extrapolation, there is good evidence for asymptotic scaling of the vortex density at couplings $`\beta >2.2`$. The scaling of the vortex density, at the $`\beta `$ values studied, is in fact substantially better than the scaling of string tension itself.

Center projection in maximal center gauge is intended as a method for locating center vortices, on unprojected lattices, from the position of P-vortices on projected lattices. The correspondence of P-vortices with physical objects depends on

1. Scaling of the P-vortex density;
2. Correlation of P-vortex locations on the projected lattice with gauge-invariant observables on the unprojected lattice.

In particular, it is the correlation of P-vortices with unprojected Wilson loops, and specifically the ratios of “vortex-limited” Wilson loops $`W_n(C)/W_0(C)`$ on the unprojected lattice, which indicate that P-vortices correspond to center vortices (rather than to some other type of object) on the unprojected lattice . Moreover, if P-vortices locate center vortices, and if center vortices produce the full asymptotic string tension, then we must find

1. Center dominance; i.e. the equality of the string tension on projected and unprojected lattices.

The numerical simulations reported here were concerned with the scaling of the P-vortex density, and with center dominance. From the results of those simulations, we conclude that center projection in maximal center gauge remains a viable method of locating center vortices on full, unprojected lattices, and that these vortices probably account for the entire asymptotic string tension.

###### Acknowledgments.

Our research is supported in part by Fonds zur Förderung der Wissenschaftlichen Forschung P13997-PHY (M.F.), the U.S. Department of Energy under Grant No. DE-FG03-92ER40711 (J.G.), the “Action Austria-Slovakia: Cooperation in Science and Education” (Project No. 30s12) and the Slovak Grant Agency for Science, Grant No. 2/7119/2000 (Š.O.). J.G. also acknowledges the support of MaPhySto, Centre for Mathematical Physics and Stochastics, funded by the Danish National Research Foundation.


text:
# Saturation of low-energy antiproton annihilation on nuclei

## I Introduction

Experimental results for antiproton annihilation cross sections at very low energies ($`p_L<100`$ MeV/c), below the $`\overline{p}p\overline{n}n`$ charge-exchange threshold, have recently been reported for light nuclei . At these energies the total $`\overline{p}`$ reaction cross section consists only of $`\overline{p}`$ annihilation. Bianconi et al.  reported the first ever measurements of $`\overline{p}`$ annihilation on Ne in the momentum range of 53-63 MeV/c. Comparing with other data, they showed that whereas at relatively higher energies ($`p_L200600`$ MeV/c) ratios of $`\overline{p}`$ annihilation cross sections on different nuclei exhibit the well-known $`A^{2/3}`$ strong-absorption dependence, such ratios at very low energies defy any simple, obvious regularity. It has been demonstrated that the ‘expected’ $`ZA^{1/3}`$ dependence of these cross sections on the atomic charge $`Z`$ and mass number $`A`$ is badly violated . Antiproton annihilation cross sections at very low energies simply do not rise with $`A`$ as fast as is anticipated. For example, over the whole mass range studied so far (H, D, 4He, Ne), the Ne/H ratio of total annihilation cross sections is at least 6 times smaller than expected.

In the present Letter we study low-energy $`\overline{p}`$ annihilation on nuclei, using the optical model approach. Optical potentials have been very successful in describing strong-interaction effects in hadronic atoms , including $`\overline{p}`$ atoms . It has been noted, for pions, that the total reaction cross sections at low energies are directly related to the atomic-state widths, and that once a suitable optical potential is constructed by reasonably fitting it to the atomic level shifts and widths in the negative energy bound-state domain, these total reaction cross sections are reliably calculable . The recent publications  of experimental results of total cross sections for $`\overline{p}`$ annihilation on nuclei at very low energies raise the intriguing possibility of connecting these two energy regimes in a systematic way also for antiprotons. However, most of the data on annihilation cross sections for $`\overline{p}`$ are for very light nuclei, where the concept of a rather universal optical potential that depends on $`A`$ and $`Z`$ only through the nuclear densities is questionable. For this reason we use optical potentials, in the present work, mostly for crossing the $`E=0`$ borderline within the same atomic mass $`(A)`$ range, from bound states to scattering. This optical model approach is different from the allegedly model-independent scattering length approximation  which we have found to be less useful. We defer a discussion of the latter statement to a separate publication.

Among the targets used in the $`\overline{p}`$ annihilation measurements , the only ‘real’ nuclei are 4He and Ne, in the sense that the conventional optical model is not expected to be applied to the lighter targets of hydrogen or deuterium. In the present work we show that the very low-energy $`\overline{p}`$ total annihilation cross sections on 4He  and on Ne  are reproduced well by reasonable optical potentials fitted separately to the $`\overline{p}`$ atomic data at the relevant atomic mass range. These potentials are strongly absorptive, which leads to a remarkable saturation of the total reaction cross section with increasing $`A`$. Strong absorption has very recently been shown  to lead also to saturation of the widths of $`\overline{p}`$ atomic states and to the prediction of relatively narrow deeply bound $`\overline{p}`$ atomic states. In view of the close analogy between bound-state widths and total reaction cross sections (see Eqs. (1,2) below), it should come as no surprise that the underlying physics is the same. In fact, several $`\overline{p}`$D calculations  have found that the $`\overline{p}`$ absorptivity in deuterium, as suggested by the magnitude of the imaginary part of the $`\overline{p}`$D $`s`$-wave scattering length, is weaker than in hydrogen due partly to the repulsive effect of absorption. Lastly, we use $`\overline{p}`$ optical potentials fitted to a comprehensive set of atomic data across the periodic table to predict $`\overline{p}`$ total annihilation cross sections at very low energies for a wide range of $`A`$ values, with the exception of very light targets. We show such calculated cross sections at $`p_L=57`$ MeV/c and discuss their $`A`$ and $`Z`$ dependence, particularly with respect to the $`ZA^{1/3}`$ dependence arising from Coulomb focussing in the limit of very low energies. We also demonstrate the extent to which this dependence is violated at 57 MeV/c.

## II Optical potentials

In the present work we aim at connecting $`\overline{p}`$ atoms, at energies slightly below threshold, with $`\overline{p}`$ annihilation on nuclei at very low energies above threshold, using an optical potential $`V_{\mathrm{opt}}`$. Assuming for simplicity a Schrödinger-type equation, the width $`\mathrm{\Gamma }`$ of an atomic level is given (non-perturbatively) by:

$$\frac{\mathrm{\Gamma }}{2}=\frac{\mathrm{Im}V_{\mathrm{opt}}(r)|\psi (𝐫)|^2𝑑𝐫}{|\psi (𝐫)|^2𝑑𝐫}.$$
(1)
Here $`\psi (𝐫)`$ is the $`\overline{p}`$ full atomic wavefunction. The corresponding expression for the total reaction cross section at positive energies is

$$\sigma _R=\frac{2}{\mathrm{}v}|\chi (𝐫)|^2\mathrm{Im}V_{\mathrm{opt}}(r)𝑑𝐫.$$
(2)
Here $`\chi (𝐫)`$ is the $`\overline{p}`$ \- nucleus elastic scattering wavefunction and $`v`$ is the c.m. velocity.

There exist several optical model potentials  parameterized in terms of nuclear densities which are quite successful in reproducing the strong interaction effects in not too light nuclei, e.g. fitting all data for nuclei heavier than carbon. However, most of the good quality data on $`\overline{p}`$ annihilation at very low energies are for very light nuclei where a straightforward application of the potential deduced from heavier targets is not expected to yield good fits. We have therefore started the present analysis by studying $`\overline{p}`$ atoms of 3,4He.

The $`\overline{p}`$ nucleus optical potential used here is given by the ‘$`t\rho `$’ expression

$$2\mu V_{\mathrm{opt}}(r)=4\pi (1+\frac{\mu }{m})b_0\rho (r),$$
(3)
where $`m`$ is the nucleon mass, $`\mu `$ is the reduced $`\overline{p}`$-nucleus mass, $`b_0`$ is a complex parameter obtained from fits to the data and $`\rho (r)`$ is the nuclear density distribution normalized to $`A`$. Trying to fit level shift and width data in 3,4He  for the 2$`p`$ and 3$`d`$ states simultaneously, we always ended up with the calculated width of the 2$`p`$ state in 4He being too small, thus contributing an unacceptably large value to the total $`\chi ^2`$ of the fits. Handling each isotope separately did not change this situation and we then decided to fit only the data for the 2$`p`$ states for the two isotopes. This is somewhat unfortunate as one expects to be able to use $`l`$-independent potentials, as is indeed the case with heavier targets . However, since our main concern here was to use atomic potentials at (very low) positive energies where $`d`$ and higher partial waves contribute very little in light nuclei, we considered this procedure acceptable. Note that we are using the same potential for all the partial waves in the positive energy regime.

Very good fits to the data with reasonable values for the complex parameter $`b_0`$ could be obtained only when a ‘finite range’ was introduced in the form of a Gaussian folding of a $`\overline{p}N`$ interaction into the nuclear density distributions. This is in contrast to the case of heavier targets where such a procedure was not necessary. Nevertheless, we consider the resulting potential, with $`b_0`$=0.49+$`i`$3.0 fm and a Gaussian with a range parameter of 1.4 fm folded into the nuclear density distribution (hereafter referred to as potential (a)), as quite reasonable. It fits within the errors the measured strong interaction shifts and widths of the 2$`p`$ state in antiprotonic 3He and 4He.

## III Results and discussion

Turning to positive energies, we first calculate the total reaction (annihilation) cross section for 57 MeV/c $`\overline{p}`$ on 4He, in order to make it possible to compare results for He and Ne . Using the above potential (a), the calculated cross section is 901 mb, in excellent agreement with the interpolated value of 915$`\pm `$39 mb . For incident momenta of 47.0 and 70.4 MeV/c we calculate values of 1060 and 771 mb, respectively, where the experimental values  are 979$`\pm `$145 and 827$`\pm `$38 mb, respectively. Reasonably good agreement is therefore established between experiment and predictions made with a potential derived from fits to $`\overline{p}`$ atoms of 3,4He.

Figure 1 demonstrates the extreme strong-absorption conditions which are relevant to the $`\overline{p}`$ nucleus interaction at very low energies (and for $`\overline{p}`$ atoms). It shows calculated reaction cross sections for $`\overline{p}`$ at 57 MeV/c on 4He and Ne as function of the strength Im $`b_0`$ of the imaginary part of the potential (a) described above, with the rightmost edge corresponding to its nominal value Im $`b_0`$=3.0 fm. It is seen that as long as the absorptivity (Im $`b_0`$) is very weak, less than 1% of its nominal value, $`\sigma _R`$ is approximately linear in Im $`b_0`$, which according to Eq. (2) means that the $`\overline{p}`$ wavefunction depends weakly on Im $`b_0`$. However, already at below 5% of the nominal value of Im $`b_0`$ the reaction cross sections begin to saturate, much the same as for the widths of deeply bound $`\overline{p}`$ atomic states . The mechanism is the same in both cases, namely exclusion of the wavefunction from the nucleus due to the absorption, which reduces dramatically the overlap with the imaginary potential and consequently reduces the integrals in Eqs. (1,2). The onset of saturation is determined approximately by the strength parameter 2$`\mu `$(Im $`V_{\mathrm{opt}}`$)$`R^2`$, where $`R`$ is the radius of the nucleus. Thus, saturation of $`\sigma _R`$ in Ne starts at a smaller absorptivity than its onset in 4He. The Ne/4He ratio of $`\sigma _R`$ values changes, due to this effect, from about 15 in the perturbative regime of very weak absorptivity to about 3 in the strong-absorption regime. However, the precise, detailed pattern of the change also depends on the real part of the optical potential which may increase or reduce the exclusion of the $`\overline{p}`$ wavefunction from the nucleus.

Figure 2 shows calculated $`\overline{p}`$ reaction cross sections at 57 MeV/c across the periodic table. The dot-dashed line is for the above mentioned potential (a) which is expected to be valid only in the immediate vicinity of He. The dashed line (b) is for the first potential from Table 7 of Ref.  which fits $`\overline{p}`$ atom data over the whole periodic table, starting with carbon. This potential is not expected to fit data for very light nuclei, and indeed it does not fit the 4He annihilation cross section. However, it is noteworthy that for $`A>`$ 20 the two potentials predict almost the same cross sections and certainly display a very smooth and similar dependence on $`A`$. A smooth $`A`$ dependence of the calculated cross sections is generally expected once several partial waves contribute, as we elaborate below. Also shown in the figure is the recent experimental result  for Ne, with very limited accuracy. Furthermore, for $`A>`$ 20 the points along the solid line are the calculated $`\overline{n}`$ \- nucleus total reaction cross sections, obtained from potential (b) by switching off the Coulomb interaction. The solid line is a fit to an $`A^{1/3}`$ power law which appears to be appropriate to strong absorption of uncharged particles at these very low energies, in contrast to the prediction of approximate $`A`$ independence made in Ref.  for $`\overline{n}`$ annihilation at ultra low energies where only $`s`$ waves contribute. We checked that by reducing the absorptivity by 3 - 4 orders of magnitude, a proportionality of $`\sigma _R`$ to $`A`$ as expected for weak absorption, is indeed observed. Comparing the dashed line for negatively charged particles with the solid line for uncharged particles, it is clear that the $`\sigma _R`$ values obtained by including the Coulomb interaction are considerably enhanced with respect to the ones obtained without it. Attempting to fit the enhancement by a $`Z^\alpha `$ power law, we find that $`\alpha `$ varies roughly between 1/3 for Ne to 1/2 for Pb. Attempts to identify more precisely this enhancement across the periodic table have been recently addressed in Ref. .

Finally, we discuss the $`ZA^{1/3}`$ scaling considered in Refs.  and in references quoted therein. We note that an attractive Coulomb potential causes focussing of partial-wave trajectories onto the nucleus, an effect which at very low energies may be evaluated semiclassically (see Ref.  for the applicability of the semiclassical approximation to Coulomb scattering). The maximal orbital angular momentum $`l_{\mathrm{max}}`$, for which the $`\overline{p}`$-nucleus distance of closest approach is smaller than the nuclear radius $`R`$, at sufficiently low energy is given by

$$(l_{\mathrm{max}}+1/2)^22\eta kR,$$
(4)
where $`k`$ is the c.m. momentum and

$$\eta =2Z/k_La_\mathrm{B}^{(p)}$$
(5)
is the corresponding Coulomb parameter. Here $`k_L=p_L/\mathrm{}`$, where $`p_L`$ is the $`\overline{p}`$ laboratory momentum, and $`a_\mathrm{B}^{(p)}`$ is the Bohr radius of the $`\overline{p}p`$ atom. Since 2$`\eta >>kR`$ at very low energies and for high values of $`Z`$, the value of $`l_{\mathrm{max}}`$ greatly exceeds the value $`kR`$ which is appropriate to uncharged particles ($`\eta =0`$). For example, on Pb at $`p_L=57`$ MeV/c, $`l_{\mathrm{max}}`$ is about 6 for $`\overline{p}`$, but is less than 2 for $`\overline{n}`$. Therefore, for strong absorption, where all the lower partial waves are totally absorbed, one obtains

$$\sigma _R\frac{\pi }{k^2}\underset{l=0}{\overset{l_{\mathrm{max}}}{}}(2l+1)\frac{2\eta }{kR}\pi R^2>>\pi R^2,$$
(6)
which we rewrite as

$$\sigma _R(4\pi r_0/a_\mathrm{B}^{(p)})(ZA^{1/3}/kk_L),$$
(7)
with $`r_0`$ defined by $`R=r_0A^{1/3}`$. A similar approach was used by Fäldt and Pilkuhn  to calculate Coulomb corrections to $`\pi ^\pm `$-nucleus total cross sections near the 3,3 resonance. In the present case of very low energies their correction becomes the dominant term, and the resulting reaction cross section on the r.h.s. of Eq. (6) is substantially higher than the black-disk value $`\sigma _R=\pi R^2`$ which for uncharged projectiles may be perceived as the unitarity limit of $`\sigma _R`$. Indeed, the calculated saturation plateau values of $`\sigma _R`$ in Fig. 1 are several times larger than the appropriate black-disk values.

In Figure 3 we plot, for potential (b), ratios of calculated total reaction cross sections at 57 MeV/c to the scaling parameter $`S=ZA^{1/3}/kk_L`$. It is clear that $`\sigma _R`$ does not scale with $`S`$ across the periodic table. Also plotted is the asymptotic value

$$\sigma _R/S4\pi r_0/a_\mathrm{B}^{(p)}0.273,$$
(8)
with $`r_0=1.25`$ fm and $`a_\mathrm{B}^{(p)}=57.6`$ fm. For a given $`\overline{p}`$ beam energy this asymptotic value, arising from Coulomb focussing, is approached very slowly upon increasing $`Z`$.

Before concluding we briefly mention other derivations of the $`ZA^{1/3}`$ dependence of $`\sigma _R`$. The proportionality to $`Z`$ is usually derived from the well-known Coulomb enhancement Gamow factor which at low energies is well approximated by $`2\pi \eta `$ . However, at these very low energies, the Coulomb wavefunctions are far from being constant over the nuclear volume, so that the applicability of the Gamow factor appears questionable. The origin of a possible proportionality of $`\sigma _R`$ to $`A^{1/3}`$ is even less clear: for $`s`$-wave dominated $`\overline{p}`$ annihilation, as considered in Refs. , it is due to the imaginary part of the $`\overline{p}`$ scattering length which is assumed to closely follow the nuclear radius $`R`$. However, for $`\overline{p}`$ strongly absorptive potentials this imaginary part has invariably been found to be about 1 fm across the periodic table  (see also Ref. ), considerably less than $`R`$. It is inconceivable that a proportionality to $`R`$ can arise unless several partial waves contribute significantly.

## IV Conclusions

In conclusion, we have shown that the recently reported annihilation cross sections for $`\overline{p}`$ on 4He and Ne at 57 MeV/c are reproduced very well by optical potentials which fit $`\overline{p}`$ atomic data in the respective mass ranges (Fig. 2). For these strongly absorptive potentials, the asymptotic $`ZA^{1/3}`$ dependence of the total reaction cross section was derived from the Coulomb focussing effect at very low energies, but at 57 MeV/c it was found not to be satisfied across the periodic table except for unrealistically high values of $`A`$ (Fig. 3). The correct dependence is determined by the saturation property of total reaction cross sections at very low energies (Fig. 1). The apparent suppression of $`\overline{p}`$ total annihilation cross sections measured at $`p_L<100`$ MeV/c is a direct consequence of the strong absorption, and in particular its saturation property. Predictions have been given for $`\overline{p}`$ total annihilation cross sections at 57 MeV/c across the periodic table which would complement and enrich the information deduced from the measurements already available for $`\overline{p}`$ atoms.

A more extended discussion of optical model calculations and their implications for $`\overline{p}`$ atoms and $`\overline{p}`$ annihilation on nuclei, including a connection to $`\overline{p}p`$ annihilation at $`p_L<200`$ MeV/c, is included in a forthcoming publication .

CJB wishes to thank the Hebrew University for support for a visit during which this work was started.
This research was partially supported by the Israel Science Foundation.