text:
# 1 Introduction

## 1 Introduction

Recent analysis of the AGASA data shows anisotropy in arrival directions of cosmic rays with energies $`10^{17.9}10^{18.3}`$ eV, with excesses from two directions near the Galactic Centre (4.5$`\sigma `$) and the Cygnus region (3.9$`\sigma `$) (Hayashida et al. 1999). The existence of a point like excess at $`7.5^\mathrm{o}`$ from the Galactic Centre (GC) has been confirmed by the analysis of the SUGAR data (Bellido et al. 2001). Hayashida et al. suggested that such point like excesses might be caused by relativistic neutrons which are able to reach the Earth from distances as large as that to the GC. These neutrons could be produced in hardonic collisions of cosmic rays (Medina-Tanco & Watson 2001, Takahashi & Nagataki 2001) which have been accelerated: (1) by a massive black hole associated with the Sgr A (Levinson & Boldt 2002); (2) by very young pulsars (Blasi, Epstein & Olinto 2000, Giller & Lipski 2002, Bednarek 2002); or (3) by shock waves of supernovae exploding into their own stellar winds (Rhode, Enslin & Biermann 1998). Relativistic protons, likely responsible for injection of these neutrons, escape from the source and may reach the Earth after propagation in regular and turbulent galactic magnetic fields.

In this paper we concentrate on the details of propagation of protons with energies $`10^{18}`$ eV from a source located in the general direction of the GC (but at different distances from the Earth) applying two models of the galactic magnetic field in the Galactic Plane (GP) and halo. The existence of a large magnetic halo extending several kpc out of the GP in the perpendicular direction $`z`$ is suggested by the lack of cosmic ray deficit out of the GP in the AGASA data at $`10^{18}`$ eV (see Clay 2001) and by observations of nearby spiral galaxies, e.g. NGC 253 - $`z>10`$ kpc, $`B7\mu `$G (Beck et al. 1994), NGC 4631 - $`z>8`$ kpc, $`B=2\mu `$G (Golla & Hummel 1994), NGC 891 and NGC 4561 - $`z3`$ kpc, $`B=1\mu `$G (Sukumar & Allen 1991).

The propagation effects of charged particles with extremely high energies (EHE) through the galactic regular magnetic field has been already analysed by many groups (e.g. Giller et al. 1994, Stanev 1997, Zirakashvili et al. 1998; Biermann et al. 2000, Harari et al. 2000 Alvarez-Muniz et al. 2001, and references therein). In these papers the method of solving the ”inverse” problem has been applied i.e. calculating trajectories of oppositely charged particles emitted from the observation point (Karakuła et al. 1972). This method is useful when considering many sources in the Galaxy but it is not better than the straightforward following the particle trajectory from its source to the observer, when one point source is considered. The propagation of particles accelerated in the GC has been also studied recently by assuming that the field is axially symmetric allowing for the collection of all particles arriving in an annular stripe with radius equal to the distance of the Sun from the GC (e.g. Giller & Zielińska 2000, O’Neill, Olinto & Blasi 2001). Such method, however, can not be applied to the propagation calculations from a source located at an arbitrary site in the Galaxy because the whole system is not symmetric then. Therefore the propagation of particles with energies $`10^{18}`$ eV from a discrete source situated in an arbitrary place in the Galaxy requires calculating trajectories of particles injected by the source and registering those intersecting a vicinity of the Earth. Such studies of CR propagation from the source located in or close to the GC in the context of the above mentioned observational results were carried out by Clay (2000), Clay et al. (2000), Bednarek, Giller & Zielińska (2001).

The cosmic ray anisotropies obtained from the above studies, although strongly dependent on the magnetic field adopted, have shown that it is very unlikely that a source of heavy nuclei with $`E10^{18}`$ eV would produce a compact excess on the sky, unless very close to us. Therefore, we will consider here the possibility that the observed excess is due to cosmic ray protons and study their propagation from some point sources located in the direction of the GC. In particular, we have studied the point source image as a function of time. We have also studied its dependence on the model of the regular magnetic field in the Galaxy as well as on the proton energy. We have also considered whether the observed excess near the GC could be caused by particle emission from a single pulsar. In a future paper we shall discuss the propagation from sources located in the Cygnus region (second excess found by the AGASA group) and a source at the Galactic anticentre to find out if a relatively close single pulsar of the Crab type can contribute significantly to the observed cosmic rays at EeV energies. Note that a dominant contribution of a single source to the cosmic ray spectrum at lower energies $`10^{15}10^{16}`$ eV (the knee region) has been suggested by Erlykin & Wolfendale (1997).

## 2 Sources of EHE particles within the Galaxy

We assume that particles with different energies are ejected isotropically by a short lived source, most likely a very young pulsar. Such very young pulsars (with milisecond periods) are presumably formed during the supernova type Ib/c explosions. The precursors of these supernova types are probably low mass Wolf-Rayet or oxygen-carbon stars rotating very fast and having light envelopes after explosion. Recent observations of diffuse hot plasma emitting X-rays (Yamauchi et al. 1990) suggest that in the past $`10^5`$ years about $`10^3`$ supernovae have exploded in the GC. Let’s assume that at least one of them had parameters allowing acceleration of protons to energies above $`10^{18}`$ eV.

We follow the suggestion that pulsar winds are able to accelerate particles to energies $`E`$ corresponding to the full potential drop available across the polar cap region (Gunn & Ostriker 1969, Blasi, Epstein & Olinto 2000),

$`E={\displaystyle \frac{eBR^3\mathrm{\Omega }^2}{c^2}}{\displaystyle \frac{6.6\times 10^{19}B_{13}}{P_{\mathrm{ms}}^2}}\mathrm{eV},`$ (1)
where $`\mathrm{\Omega }=2\pi /P`$, $`P=10^3P_{\mathrm{ms}}`$ s is the pulsar period, $`R=10^6`$ cm is the radius of the neutron star, $`B=10^{13}B_{13}`$ G is its surface magnetic field, $`e`$ is the elementary charge, and $`c`$ is the velocity of light. Eq. (1) allows us to constrain the parameters of the pulsar able to accelerate protons to energies $`E10^{18}`$ eV. The following condition has to be fulfilled

$`P_{\mathrm{ms}}8B_{13}^{1/2}.`$ (2)
If the pulsar loses its rotational energy, $`E_{\mathrm{rot}}`$, only on emission of the dipole radiation with the power $`L`$, then its period at specific time $`t`$ is determined by the equation

$`\dot{E}_{\mathrm{rot}}=LI\mathrm{\Omega }\dot{\mathrm{\Omega }}=B^2R^6\mathrm{\Omega }^4/6c^3`$ (3)
where $`I=1.4\times 10^{45}`$ g cm<sup>-2</sup> is the neutron star moment of inertia. As a result of this energy losses the period of the pulsar changes in time according to

$`P_{\mathrm{ms}}^2(t)=P_{0,ms}^2+3.5\times 10^8B_{13}^2t,`$ (4)
where $`t`$ is in seconds, and $`P_{0,ms}`$ is the pulsar initial period. By reversing Eq. (4) and applying Eq. (2), we estimate the time elapsed from the pulsar formation, $`t_{\mathrm{acc}}`$, during which protons will be accelerated above $`10^{18}`$ eV (assuming that $`P_{0,ms}P_{\mathrm{ms}}(t_{\mathrm{acc}})`$). It is

$`t_{\mathrm{acc}}60B_{13}^1\mathrm{yr}.`$ (5)
Since this time is relatively short, when compared to the time scale of particle propagation, we can consider such injection of relativistic protons as instantaneous.

The hypothesis for the origin of particles with such energies in pulsars makes the model flexible since, in principle, the pulsar can be born at an arbitrary site in the Galaxy. Such proposition can give a better explanation of the observations (i.e. the second AGASA excess in the Cygnus region) and allows investigation of particle propagation from the source shifted from the GC itself (see Bellido et al. 2001), or in the direction of the GC but located at a different distance. Since the propagation of protons with EeV energies in the Galactic magnetic field takes tens to hundreds of thousand years longer than that of radiation, this pulsar may not necessarily be visible presently as a source of $`\gamma `$-rays, neutrinos or neutrons, unless it is immersed in a dense molecular cloud accumulating protons (a model recently considered by Bednarek 2002). Almost instantaneous injection of $`10^{18}`$ eV protons by such pulsar together with a particular structure of the regular magnetic field in the Galaxy may result in quite unusual images of the source at different times after injection.

## 3 The structure of the Galactic magnetic field

Protons ejected from a point-like source propagate in the Galactic magnetic field which consists of regular $`B_{\mathrm{reg}}`$ and irregular $`B_{\mathrm{irr}}`$ components. We have adopted two different models for the regular Galactic magnetic field. The first model (our model I) proposed and described in detail by Urbanik, Elstner & Beck (1997) bases on observations of near-by spiral galaxies. $`B_{\mathrm{reg}}`$ has a toroidal component, confined mainly to the disk, and a large scale poloidal component, extending up to $`z=10`$ kpc. The total $`B_{\mathrm{reg}}`$ does not exceed 2$`\mu `$G anywhere. For the detailed structure of the magnetic field in the halo see Fig. 2 in the Urbanik et al. paper.

As the second possibility (our model II) we adopt the bisymmetric field model with field reversals and odd parity (BSS-A) proposed by Han & Qiao (1994) and applied for the particle propagation purposes by e.g. Stanev (1997). This model incorporates the knowledge from experimental observations of our Galaxy as well as many other galaxies. It consists of the GP component in which the field strength at a point (r, $`\theta `$) is described by

$`B(r,\theta )=B_0(r)\mathrm{cos}[\theta \beta \mathrm{ln}(r/r_0)],`$ (6)
where $`r_0`$ is the Galactocentric distance of the location with maximum field strength at $`l=0^\mathrm{o}`$, $`\beta =5.67`$, and $`r_0=10.55`$ kpc. $`B_0(r)`$ is taken to be $`3R_{\mathrm{GC}}/r`$ $`\mu `$G above $`r=4`$ kpc and constant below this distance from the centre of the Galaxy, where $`R_{\mathrm{GC}}=8.5`$ kpc. The radial and azimuthal components of the magnetic field in the halo are described by

$`|B(r,\theta ,z)|=|B(r,\theta )|exp(|z|/z_0),`$ (7)
with two scale heights $`z_0=1`$ kpc for $`|z|<0.5`$ kpc and $`z_0=4`$ kpc for $`|z|>0.5`$ kpc with the field direction at the disk crossing unchanged (see Stanev 1997). The magnetic field component perpendicular to the GP, $`B_\mathrm{z}`$, is assumed constant with the value of 0.3$`\mu `$G and is always directed to the north.

In both models we describe the irregular field, $`B_{\mathrm{irr}}`$, as a sum of many plane waves with isotropically distributed wave vectors and amplitudes corresponding to the Kolmogorov power spectrum. Its mean value is 2 $`\mu `$G in the disk with $`z\pm 500`$ pc, and 0.5 $`\mu `$G in the spherical halo with radius 20 kpc. The irregularity scale is, however, different in the disk and the halo: the longest wave is 150 pc in the disk and 7 kpc in the halo.

## 4 Propagation of EHE protons

We calculate numerically the proton trajectories within the range of energies corresponding to those of the AGASA excess ($`10^{17.9}10^{18.3}`$ eV) and the SUGAR excess ($`10^{17.9}10^{18.5}`$ eV) from the direction of the GC. For $`1.2\times 10^6`$ protons ejected isotropically from a point-like source located at three different points in the Galaxy we record the parameters (numbers, directions) of particles intersecting a sphere with the radius of 250 pc centred on the Earth. These events are considered as observed by a detector on the Earth. Other nuclei with energy $`E`$ and charge $`Z`$, propagating in the magnetic field by a factor $`\alpha `$ stronger behave exactly as protons with energies $`E/(\alpha Z)`$.

We consider the sources of charged particles located: (a) exactly at the GC at the distance of 8.5 kpc from the Sun; (b) 2 kpc from the Sun towards the GC (for this distance the potential source of the SUGAR excess lays still within the region of the galactic disk); (c) 8.5 kpc from the Sun towards the direction of the SUGAR excess displaced from the GC direction by $`7.5^\mathrm{o}`$, and $`2.5^\mathrm{o}`$ below the GP (about $`400`$ pc from GP). Below, we discuss the calculation results for these three cases.

### 4.1 A source in the Galactic Centre

The GC seems to be one of the most likely site for particle acceleration to energies above $`10^{18}`$ eV. In Sect. 2 we argue that the best candidate source is a very young neutron star with a milisecond period and surface magnetic field $`B10^{13}`$ G. Since the particles can be acclerated by such a source for a short time, we assume the instantaneous injection of protons with energies $`1,2,3\times 10^{18}`$ eV.

First we discuss the results of calculations of proton trajectories from a point-like source exactly in the GC for the two galactic field models. The numbers of particles arriving at different times after injection (the time of flight along the straight line being subtracted) to the sphere centered on the Earth are displayed in the form of histograms in Figs. 1h, 2h, and 3h for model I, and in Figs. 1p, 2p, and 3p for model II. It becomes evident that the distribution of the arrival times of particles with energies by a factor 2-3 larger is completely different. The bulk of particles arrive to the observer within $`2.5\times 10^4`$ years (model I) and $`5\times 10^4`$ years (model II) for $`3\times 10^{18}`$ eV (Figs. 1h and 1p), up to $`10^5`$ years (model I) and $`10^6`$ years (model II) for $`10^{18}`$ eV (Figs. 3h and 3p).

In Figs. 1,2, and 3 from a) to f) (model I) and from i) to n) (model II) we show maps (in galactic coordinates, with longitude increasing to the left) with the arrival directions of protons intercepting the sphere around the Earth within consecutive time delay intervals chosen accordingly to the particle energy and magnetic field model (see figure captions). Maps summed up over time are in Figs. g) and o) showing the direction distribution in the case of a steady source.

Let’s first concentrate on the results for model I. The most interesting feature for protons with energies $`(23)\times 10^{18}`$ eV is the particle clustering in multiple images of the source. These images appear at different places on the sky at different times after injection. A large number of protons reach the Earth’s vicinity from directions close to the GC (shifted by about $`10^\mathrm{o}`$ towards positive longitudes) creating an extended source with the radius of about $`20^\mathrm{o}30^\mathrm{o}`$. Note that this location is consistent with the direction towards the excess of particles obtained by AGASA and SUGAR. Therefore, for the galactic magnetic field structure proposed by Urbanik et al. the source of these particles can be actually located exactly in the Galactic Centre !

For protons with energies $`10^{18}`$ eV the arrival directions become much more scattered (see Figs. 3a to 3h). In this case, protons arrive from a large part of the sky, almost independently of time after injection, apart from the peak for the first $`2\times 10^4`$ years. Increasing slightly the field strength will cause the discovered features shifting to higher energies and fitting better to the energy range where the actual excess of particles has been detected.

The arrival directions of protons and their time distributions are completely different for the field model II. There are no protons arriving directly from the actual position of the source at the GC (out of $`10^6`$ ejected) up to $`3\times 10^{18}`$ eV. For protons with $`3\times 10^{18}`$ eV, only a single image of the source is visible at a high negative latitude. It is created mainly by protons arriving with relatively small time delay with respect to the rectilinear propagation, i.e. within less than $`2\times 10^4`$ years (see Fig. 1p). For lower proton energies the image of the source is also centred on high galactic latitudes becoming broader and stronger. Particles arrive to the Earth much later than for model I i.e. after $`(24)\times 10^5`$ years for $`2\times 10^{18}`$ eV and $`(27)\times 10^5`$ years for $`3\times 10^{18}`$ eV. In spite of the instantaneous injection the anisotropy due to the source would be visible in the same directions for a long time.

By comparing our calculation results on the proton anisotropy for the two magnetic field models we conclude that the propagation of charged particles is very sensitive to their energies and to the structure and strength of the magnetic field. The two models, both based on experimental observations, give totally different predictions concerning the particle angular distribution on the sky, meaning that one should be very careful with drawing any conclusions based on one particular model of the Galactic magnetic field.

### 4.2 A source at 2 kpc towards the Galactic Centre

The distance to the source of particles responsible for the AGASA-SUGAR excess can not be constrained by the observations. Therefore it is reasonable to investigate the case of a source located much closer than the GC.

In Figs. 4,5 and 6 we show maps of the arrival directions of protons intercepting the sphere around the Earth for a source at 2 kpc. The maps summed up over time are shown in Figs. g) and o). For this relatively small distance in all considered cases (different proton energies 1,2, and $`3\times 10^{18}`$ eV and both structures of the galactic magnetic field) a large number of protons reach the Earth by moving almost along straight lines (see initial peaks in histograms showing the time distribution of the arriving particles). Therefore the excess of relativistic protons should be observed for a short time only.

In model I protons injected with $`3\times 10^{18}`$ eV can reach the Earth only sporadically, if the time delay after injection is large (see Figs. 4 b-f and h). However, the number of the delayed protons quickly increases with decreasing energy and at $`10^{18}`$ eV the particles are arriving to us for $`2\times 10^5`$ years. The intensity of these particles differ significantly at specific time intervals showing several, periodic, maxima, seen clearly in Figs. 5h and 6h. The magnetic field $`B_{\mathrm{reg}}`$ in the halo turns out to be directed from the Sun towards the source, so that this periodic arrival corresponds to the particle multiple gyro-orbits. During the maximum intensities particles arrive from a great circle on the sky, roughly perpendicular to the GP and crossing it at $`l90^\mathrm{o}`$ and $`270^\mathrm{o}`$.

The maps of arrival directions and their time structure differ significantly for model II also for this case when the source is relatively close to us (2 kpc). There is clear deflection of the instantaneuos image of the source towards the south, increasing for lower energies. For $`3\times 10^{18}`$ eV, there exists a secondary extended image shifted by a large angular distance towards negative galactic longitudes at the time interval $`(1.32)\times 10^5`$ years (see Figs 4k,l, and p). For lower energies (2$`\times 10^{18}`$ eV) a large fraction of protons are delayed by up to $`10^6`$ years after injection. They arrive mainly from the southern hemisphere due to the asymmetric form of the $`B_\mathrm{z}`$ component in the halo. An interesting thing is that there is an elongated image in the directions opposite to that of the source ($`l180^\mathrm{o}`$, $`\delta <0`$).

### 4.3 A source towards the AGASA-SUGAR excess

As we mentioned above the excess seen in the SUGAR data is consistent with a point-like source but displaced from the GC direction by $`7.5^\mathrm{o}`$. The maximum probability map of the AGASA excess is shifted from the GC direction even more. Therefore, it is reasonable to consider a source located at a certain distance from the exact position of the GC. We have calculated trajectories of isotropically injected protons assuming that the source is 8.5 kpc away but in the direction of the SUGAR excess. As before we collect protons with energies $`1,2,4,10\times 10^{18}`$ eV intersecting the sphere around the Earth for the two magnetic field models. Now the source is located about $`400`$ pc below the GP. The obtained maps (integrated over time) are presented in Fig. 7. However, there are no significant differences in comparison to the calculations from the source located exactly in the GC. Nevertheless, the image sizes for both cases are rather large at the considered energies and would be close to the AGASA excess size for $`E4\times 10^{18}`$ eV.

## 5 Single pulsar as a plausible source

Let us consider whether a single pulsar could be responsible for the observable excess from the GC. The power emitted by the pulsar accelerating protons can be expressed by the proton energy $`E`$ from Eqs. 1 and 3,

$`L=cE^2/6e^25.5\times 10^{40}E_{18}^2\mathrm{erg}\mathrm{s}^1,`$ (8)
where $`E=10^{18}E_{18}`$ eV. It is difficult to state at present if such powerful pulsar is present at the GC. One remnant of a young supernova with age of $`80`$ years (G0.570-0.018) has been recently reported by Senda et al. (2001). Also a strong $`\gamma `$-ray source ($`2\times 10^{37}`$ erg s<sup>-1</sup>) with a very hard spectrum is observed by the EGRET instrument from the GC (Mayer-Hasselwander et al. 1998). The nature of this source is at present unknown.

The rate, $`r_{\mathrm{obs}}`$, at which protons arrive to the Earth (i.e. intersect the sphere with the radius $`R_\mathrm{E}=250`$ pc) from the ’AGASA-SUGAR source’ can be estimated from the flux $`F_{\mathrm{obs}}=9\times 10^{14}`$ m<sup>-2</sup> s<sup>-1</sup> derived from the SUGAR observations. It is

$`r_{\mathrm{obs}}=\pi R_\mathrm{E}^2F_{\mathrm{obs}}1.7\times 10^{25}\mathrm{particles}s^1.`$ (9)
On the other side the number of particles $`\mathrm{\Delta }N`$ injected by the pulsar with energies between $`E_1`$ and $`E_2`$ would be

$`\mathrm{\Delta }N={\displaystyle \frac{c^2\eta I}{ZeR_{\mathrm{NS}}^3B}}\mathrm{ln}(E_2/E_1){\displaystyle \frac{2.6\times 10^{44}\eta }{B_{13}}}\mathrm{ln}{\displaystyle \frac{E_2}{E_1}}\mathrm{protons},`$ (10)
where $`Ze`$ is the particle charge, and $`\eta `$ is the efficiency of particle acceleration defined as the ratio of the number of injected particles to that of the maximum number possible given by the Goldreich & Julian (1969) density at the pulsar light cylinder. Our present calculations have shown that the source image changes significantly even if proton energy goes from 3 to $`4\times 10^{18}`$ ev, i.e. changes by 30%. Therefore, we think that it is quite likely that the observed excess is actually due to a narrower energy band than the quoted factor of 4 (unless the source is quite close to us) and we adopt that $`E_2/E_1=1.2`$ (with $`E3\times 10^{18}`$ eV). Then we have

$`\mathrm{\Delta }N{\displaystyle \frac{5\times 10^{43}\eta }{B_{13}}}\mathrm{protons},`$ (11)
For a pulsar located in the GC the fraction $`f`$ of these particles giving the image close to the source (for model I, Fig 1h) is $`5\times 10^4`$ (about 550 particles out of $`1.2\times 10^6`$) i.e. 2.5 times larger that that for a straightforward propagation. As these particles arrive to the Earth vicinity within $`\mathrm{\Delta }t5\times 10^3`$ years, we should expect for their average rate

$`r_1=f{\displaystyle \frac{\mathrm{\Delta }N}{\mathrm{\Delta }t}}2\times 10^{29}{\displaystyle \frac{\eta }{B_{13}}}\mathrm{protons}s^1.`$ (12)
Comparing this with $`r_{\mathrm{obs}}`$, we obtain that $`\eta /B_{13}10^4`$ for this case (a value probably not unreasonable !).

It may be, however, that the observed excess is due to a more delayed image of a source located somewhere else. Let’s assume (somewhat arbitrarily) that we can apply our calculation results to this case as well. For example, the secondary compact images (Fig. 1b,c,d,e) are visible for $`\mathrm{\Delta }t2\times 10^4`$ years and the fraction of the emitted particles producing them equals $`f7\times 10^5`$ ($`80`$ particles per one secondary image, out of $`1.2\times 10^6`$). Then their average rate is

$`r_25\times 10^{27}{\displaystyle \frac{\eta }{B_{13}}}\mathrm{protons}s^1,`$ (13)
and, to agree with $`r_{\mathrm{obs}}`$, $`\eta /B_{13}3\times 10^3`$. If the pulsar was closer then the efficincy $`\eta `$ would have to be smaller (Figs. 4,5,6).

We can not, however, exclude a possibility that the observed excess is being produced by even more delayed protons with a smaller intensity and a longer arrival time $`\mathrm{\Delta }t`$. If the total observed CR flux was not produced by pulsars (but by some other sources) then such a bunch of particles could produce an increase on the average CR flux.

For the bisymmetric field model proposed by Han & Qiao (1994) the image of the source located at the GC does not appear in the direction to the real source up to $`3\times 10^{18}`$ eV, but is shifted from its real position by a large angle. Therefore, in this case protons can not be responsible for the observed excess.

## 6 Conclusions

The main conclusions of our calculations of the propagation of protons with energies $`10^{18}`$ eV, injected by an isotropic point-like source in the direction towards the Galactic Centre are the following:

* Protons injected instantaneously by a point-like source at the Galactic Centre, arriving to the Earth after propagation in the galactic regular and irregular magnetic fields, can form multiple images at directions completely different from those towards the source, as well as images shifted only slightly from the position towards the source, and large scale anisotropies (north-south asymmetry) depending on the proton energies and the time elapsed from injection. These multiple images appear only at relatively narrow energy ranges of particles between $`(13)\times 10^{18}`$ eV (for our magnetic field models). The arrival directions of protons with lower energies are broadly scattered. Particles with higher energies create an extended image centered on the source.
* The results of particle propagation for the two considered magnetic field models give totally different predictions for the particle angular distribution on the sky, meaning that one should be very careful when drawing conclusions based on one particular model.
* The application of the Han & Qiao magnetic model produces strong north-south anisotropy of particles with energies $`(23)\times 10^{18}`$ eV. Most of the delayed particles arrive from the southern Galactic hemisphere, in contrast to the models of extragalactic origin of the highest energy cosmic rays favouring the northern hemisphere (e.g. Biermann et al. 2000).
* Our results do not depend on the particular distribution of the irregular magnetic component, providing that its magnitude is not larger than assumed here. A detailed experimental study of the $`10^{18}`$ eV CR anisotropy should give information about the validity of the pulsar model and/or that of the large scale galactic magnetic field.
* It is interesting to note that the maximum proton fluxes predicted for a pulsar model are larger than that observed. Thus, the model of isotropic and instantaneous particle injection by a pulsar, located close to the GC, could explain the observed flux of particles in the AGASA-SUGAR excess as due to prompt protons (travelling almost along straight lines) or as due to the delayed protons. In the last case, however, the source would be in some other direction.

## Acknowledgments.

This work is supported by the KBN grants No. 2P03C 006 18 and 5P03D 025 21.


format: text

text:
# Simple device for scanning image

## I General description

In this letter we suggest a simple device for scanning image, which has much in common with usual TV camera with an electron beam replaced by an optical one. The scheme of the device is presented in Fig. 1, where (1) is the transparent electrode; (2) is the semiconductor layer with the gap, wide enough to be transparent for the light from the image; (3) is the semiconductor layer with the gap, narrow enough to be not transparent for the light from the image, and (4) is the electrode, or the substrate with high conductivity. Suppose, that some image is formed in the plane of the layer (3) by means of some optical system. Though the image produces the profile of photoconductivity in the layer (3), there is no current in the circuit, because the layer (2) is still an isolator — note, that the layer (2) is transparent for the image light, which causes no photoconductivity in this layer. Now, let us switch on the narrow reading beam of a short enough wavelength, to produce the photoconductivity in the layer (2). This beam produces the narrow path (5) with high conductivity in the layer (2), which gives rise to the current in the circuit. The value of this current is proportional to the local illumination intensity of the image. For the sake of simplicity we suppose, that the reading beam has absorbing length much shorter, than the thickness of the layer (3), so it does not produce any noticeable photoconductivity in the layer (3). Now, if we move the reading beam along the image, the current in the circuit will drop down, when the reading beam is passing through the dark areas of the image, and will increase, when the beam is passing through the bright areas. This process is similar to that in conventional TV camera, with the optical beam instead of the electron one.

## II Experimental illustration

To test the aplicability of the above scheme, we construct the simple installation shown in Fig. 2. The exact scheme of the above device is presented in the inset. Transparent electrode (1) and electrode (4) (notations of Fig. 1) are the water layers. The 3 mkm AlGaAs layer is playing the role of the layer (2) (see Fig. 1) with the relatively wide gap ($``$1.9 eV). And finally, the 300 mkm of semiconducting GaAs substrate is playing the role of the layer (3) with the narrow gap ($``$1.5 eV). The beam of the HeNe laser (0.63 mkm, or 1.97 eV), which plays the role of a reading beam, passes through a simple deflecting system (consisting of the loudspeaker, the mirror (1) and the lens (see Fig. 2)), and gives a small spot moving along the linear segment in the plane of the device.

The Nd laser beam (1.06 mkm, or 1.17 eV) formes a small spot in the same plane, which plays the role of the image. Manipulating with the mirror (2) one can force the reading beam to pass through this spot. When this takes place, the short pulses of current occur in the circuit. It should be pointed out that despite GaAs is transparent for the 1.06 mkm radiation, the photoconductivity still occur, due to the presence of the impurities and intrinsic defects in the semiconducting GaAs substrate.

The main property of the described device is its simplicity. This device may be preferable in all cases, when high energy electron beam can damage the light-sensitive surface. Also, the device of this kind can perform some integral transformations, similar to the Fourier transformation. For this purpose one should substitute the narrow running reading beam with the sinusoidal profile, which gives the sinusoidal distribution of the conductivity in the layer (2).

The author is not sure, that such a simple device was not described, or used before, and will be gratefull for any references and coments.


text:
# Non-monotonic variation with salt concentration of the second virial coefficient in protein solutions

## I introduction

Interactions between proteins in aqueous solutions determine their collective behavior, in particular their aggregation, their complexation with other macromolecules, and ultimately their phase behavior, including phase separation, precipitation and crystallization. Any theoretical analysis of the properties of protein solutions must rely on a reliable understanding of their interactions. A good example is provided by the control of protein crystallization, which is an essential prerequisite for the determination of protein structure by X-ray diffraction. While at present protein crystallization is still mostly achieved experimentally by “trial and error”, and on the basis of a number of empirical rules, there is clearly a need for a more fundamental understanding of the mechanisms controlling protein crystallization, and this obviously requires a good knowledge of the forces between protein molecules in solution, and of their dependence on solution conditions, including pH and salt concentration.

Protein interactions have various origins, and one may conveniently distinguish between direct and induced (or effective) contributions. Direct interactions include short-range repulsive forces, which control steric excluded volume effects, reflecting the shape of the protein, van der Waals dispersion forces, and electrostatic forces associated with pH-dependent electric charges and higher electrostatic multipoles carried by the protein residues. Other, effective, interactions depend on the degree of coarse-graining in the statistical description and result from the tracing out of microscopic degrees of freedom associated with the solvent and added electrolyte, i.e. the water molecules and microions. Tracing out the solvent results in hydrophobic attraction and hydration forces, while integrating over microion degrees of freedom leads to screened electrostatic interactions between residues, the range of which is controlled by the Debye screening length, and hence by electrolyte concentration.

However, while coarse-graining through elimination of microscopic degrees of freedom, leading to state-dependent effective interactions is a priori a reasonable procedure to describe highly asymmetric colloidal systems, where particles have diameters of typically hundreds of nm and carry thousands of elementary charges, this is obviously less justified for the much smaller and less charged proteins. In particular the assumption of uniformly charged colloid surfaces, leading to spherically symmetric, screened interactions between the electric double layers around colloid particles, as epitomized by the classic DLVO (Derjaguin-Landau-Verwey-Overbeek) potential, ceases to be a reasonable approximation at the level of nanometric proteins carrying typically of order 10 elementary charges. The reason is that length scales which are widely separated in colloidal assemblies, become comparable in protein solutions, while the discreteness of charge distributions on proteins can no longer be ignored, since the distance between two charged residues on the protein surface is no longer negligible compared to the protein diameter. Thus, electrostatic, as well as other (e.g. hydrophobic) interactions are much more specific in proteins, and must be associated with several interaction sites, rather than merely with the centers of mass as is the case for (spherical) colloidal particles.

Another very important distinction between colloids and protein solutions is that the forces between the former may be measured directly, using optical means, while interactions between proteins can only be inferred indirectly, from measurements by static light scattering of the osmotic equation of state which, at sufficiently low concentrations, yields the second osmotic virial coefficient $`B_2`$, the main focus of the present paper. The variation of $`B_2`$ with solution conditions yields valuable information on the underlying effective pair interactions between proteins. Moreover it was shown empirically by George and Wilson that there is a strong correlation between the measured values of $`B_2`$ and the range of solution conditions which favor protein crystallization. Only if the measured value of $`B_2`$ falls within a well defined “slot” can crystallization be achieved. If $`B_2`$ is too large, repulsive interactions predominate, leading to slow crystallization rates. On the other hand if $`B_2`$ is highly negative, strong attractions lead to amorphous aggregation.

The correlation between $`B_2`$ and crystallization may be rationalized by noting that protein crystals generally coexist with a fairly dilute protein solution, the thermodynamic properties (and in particular the free energy) of which are essentially determined by $`B_2`$. Coexistence between a dense solid phase and a dilute fluid phase is generally a signature of a very short-ranged attraction between particles as compared to their diameter.

For such short-ranged attractive interactions, the phase-separation into dilute and concentrated proteins solutions expected on the basis of a mean-field, van der Waals theory, is in fact pre-empted by the freezing transition, i.e. the critical (or “cloud”) point lies below the freezing line. The critical fluctuations associated with this metastable cloud point may lead to a significant enhancement of the crystal nucleation rate, while the position of the cloud point in the concentration-temperature plane is strongly correlated with the virial coefficient $`B_2`$.

The present paper focuses on the variation of $`B_2`$ with ionic strength of added salt. This is a particularly important issue since “salting out “ of protein solutions is one of the standard methods used to induce crystallization. An increase in salt concentration reduces the screening length and hence the electrostatic repulsion, allowing short range attractive forces (e.g. of hydrophobic or van der Waals origin) to come into play which will ultimately trigger nucleation. Recent experiments point to a non-monotonic variation of $`B_2`$ with increasing ionic strength, or to a pronounced shoulder in the $`B_2`$ versus ionic strength curve in lysozyme solutions. Closely related findings are the observation of a non-monotonic cloud point, and of a minimum in the solubility of lysozyme with increasing salt concentration; the solubility is obviously related to the osmotic virial coefficient. Similarly, the attractive interaction parameter $`\lambda `$, which controls the variation of the measured protein diffusion coefficient $`D`$ with volume fraction, was found to exhibit a sharp minimum upon an increase of ionic strength of lysozyme solutions; again, this interaction parameter strongly correlates with $`B_2`$.

Traditional models for the protein–protein interaction cannot easily reproduce such non-monotonic behavior of $`B_2`$ or related quantities. The “colloidal” approach based on spherical particles interacting via the screened Coulomb DLVO potential can only predict a monotonic decrease of $`B_2`$ with ionic strength. The same is true of models accounting for short-range attractions via Baxter’s “adhesive sphere” representation. In these models, which assume central pair-wise interactions, $`B_2`$ reduces to a simple integral of the Mayer function associated with the spherically symmetric potential. More recent calculations account for the asymmetric shape of proteins, or include several “sticky” sites at the surface of the protein.

In these traditional calculations, electrostatic interactions between proteins and microions are routinely treated within mean-field Poisson-Boltzmann (PB) theory, generally in its linearized version (as is the case for the classic DLVO potential). However, as explained earlier, all relevant length scales (i.e. protein diameter, mean distance between charged sites on the protein surface, and between co and counterions, as well as the Debye screening length) are comparable in protein solutions, so that the discrete nature of both the interaction sites, and of the co and counterions, can no longer be ignored. Moreover, Coulomb correlations are expected to be enhanced on protein length scales and may lead to strong deviations from the predictions of PB theory, which have recently been shown to induce short-range attractions, even between much larger colloidal particles.

The present paper takes into account the discrete nature of the microions within a “primitive model” description of the electrolyte, and presents results of Molecular Dynamics calculations of the equilibrium distribution of co and counterions around two proteins and of the resulting osmotic virial coefficient $`B_2`$. Two models of the charge distribution on the surface of the spherical proteins will be considered. In the colloid-like model the charge is assumed to be uniformly distributed over the surface, while in the discrete charge model, the charges are attached to a small number of interaction sites. The latter model will be shown to lead to a distinctly non-monotonic variation of $`B_2`$ with ionic strength, as observed experimentally.

The paper is organized as follows: The model and key physical quantities are introduced in section II. Simulation details are described in section III. Results of the simulations are presented and discussed in section IV, while conclusions are summarized in section V. A preliminary account of parts of the results was briefly reported elsewhere.

## II Models, effective forces and second virial coefficient

The globular proteins under consideration are modeled as hard spheres of diameter $`\sigma _p`$, carrying a total (negative) charge $`Ze`$. Within a “primitive model” representation, the molecular granularity of the aqueous solvent is ignored, and replaced by a continuum of dielectric permittivity $`ϵ`$, while the monovalent counterions and salt ions are assumed to have equal diameters $`\sigma _s`$ and charges $`\pm e`$.

Two models are considered for the charge distribution on the surface of the protein. In the “smeared charge model” (SCM), the total charge $`Ze`$ is assumed to be uniformly distributed over the spherical surface, which is the standard model for charge-stabilized colloidal suspensions, involving highly charged objects. According to Gauss’ theorem, the SCM is equivalent to the assumption that the total charge $`Ze`$ is placed at the center of the sphere. In the “discrete charge model” (DCM), point charges $`(e)`$ are distributed over a sphere of diameter $`\sigma _d=\alpha \sigma _p`$ (with $`\alpha <1`$, i.e. slightly inside the protein surface), in such a way as to minimize the electrostatic energy of the distribution. The resulting charge pattern, well known from the classic Thompson problem (see  for a recent review), is kept fixed throughout. Such Thompson patterns do not correspond to the true charge distribution on any specific protein (see  where a simple toy model of lysozyme with different charge ditribution corresponding to solutions of different pH is constructed) but do provide a well defined discrete model for any value of $`Z`$. Note that the discrete distributions are characterized by non-vanishing multi-pole moments, depending on the symmetry of the distribution for any specific value of $`Z`$, while the SCM implies vanishing multipoles of all orders.

At this stage the SCM and DCM models involve only excluded volume and bare Coulomb interactions (reduced by a factor $`1/ϵ`$ to account for the solvent) between all particles, proteins as well as microions.

The following physical quantities were systematically computed in the course of the MD simulations, to be described in the following section.

a) the density profiles of co and counterions around a single globular protein

$$\rho _\pm (r)=<\underset{j}{}\delta (\stackrel{}{r}_j^\pm \stackrel{}{r})>$$
(1)
Here $`\stackrel{}{r}_j^\pm `$ is the position of the $`j^{th}`$ microion of species $`\pm `$ relative to the protein center, while the angular bracket denotes a canonical average over the microion configurations. For an isolated SCM protein these profiles are spherically symmetric, and depend only on the radial distance $`r=|\stackrel{}{r}|`$. For isolated DCM proteins the profiles are no longer spherically symmetric, and may be expanded in spherical harmonics, as discussed in the Appendix. The anisotropy turns out to be weak, and only the spherically symmetric component (corresponding to averaging $`\rho _\pm (\stackrel{}{r})`$ over protein orientations) will be shown in the following.

b) The second quantity, which will be the key input in the calculation of $`B_2`$, is the microion averaged total force $`\stackrel{}{F}_1=\stackrel{}{F}_2`$ acting on the center of two proteins, placed at a relative position $`\stackrel{}{r}=\stackrel{}{r}_1\stackrel{}{r}_2`$; the force $`\stackrel{}{F}_1`$ is a function of $`\stackrel{}{r}`$. Its statistical definition was discussed earlier in the context of charged colloids, and it involves three contributions:

$$\stackrel{}{F}_1=\stackrel{}{F}_1^{(1)}+\stackrel{}{F}_1^{(2)}+\stackrel{}{F}_1^{(3)}.$$
(2)
$`\stackrel{}{F}_1^{(1)}`$ is the direct Coulomb repulsion between the charge distributions on the two proteins; $`\stackrel{}{F}_1^{(2)}`$ is the microion induced electrostatic force, while $`\stackrel{}{F}_1^{(3)}`$ is the depletion force which may be traced back to the in-balance of the osmotic pressure of the microions acting on the opposite sides of protein $`1`$ due to the presence of protein $`2`$. $`\stackrel{}{F}_1^{(3)}`$ is directly expressible as the integral of the microion contact density over the surface of the protein.

In the case of the SCM, the microion averaged force depends only on the distance $`r=|\stackrel{}{r}_{12}|`$ between the two proteins. For the DCM, on the other hand, $`\stackrel{}{F}_1`$ is a function of the relative orientations of the two proteins, as characterized by the sets of Euler angles $`\stackrel{}{\mathrm{\Omega }}_1`$ and $`\stackrel{}{\mathrm{\Omega }}_2`$, i.e. $`\stackrel{}{F}_1=\stackrel{}{F}_1(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)`$.

c) Once the force $`\stackrel{}{F}_1`$ has been determined as a function of $`\stackrel{}{r}`$, $`\stackrel{}{\mathrm{\Omega }}_1`$ and $`\stackrel{}{\mathrm{\Omega }}_2`$, one may then calculate an orientationally averaged, but distance resolved, effective protein-protein pair potential according to

$$V(r)=_r^{\mathrm{}}𝑑r^{}\frac{\stackrel{}{r}^{}}{|\stackrel{}{r}|}\stackrel{}{F}_1(\stackrel{}{r}^{},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)_{\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2}.$$
(3)
where the angular brackets $`<\mathrm{}>_{\stackrel{}{\mathrm{\Omega }}_1\stackrel{}{\mathrm{\Omega }}_2}`$ refer to a canonical statistical average over mutual orientations of the two proteins weighted by the Boltzmann factor of the effective potential $`V_{eff}(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)`$ such that $`V_{eff}(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)/\stackrel{}{r}=\stackrel{}{F}_1(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)`$. Explicitly, for any quantity $`A(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)`$,

$$<A>_{\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2}=\frac{𝑑\stackrel{}{\mathrm{\Omega }}_1𝑑\stackrel{}{\mathrm{\Omega }}_2A(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)\mathrm{exp}\{(V_{eff}(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)/k_BT)\}}{𝑑\stackrel{}{\mathrm{\Omega }}_1𝑑\stackrel{}{\mathrm{\Omega }}_2\mathrm{exp}\{V_{eff}(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)/k_BT\}}$$
(4)
d) The second virial coefficient $`B_2`$ finally follows from the expression

$$B_2=\frac{1}{2}𝑑\stackrel{}{r}\left[1b(r)\right]$$
(5)
where

$$b(r)=(\frac{1}{8\pi ^2})^2𝑑\stackrel{}{\mathrm{\Omega }}_1𝑑\stackrel{}{\mathrm{\Omega }}_2\mathrm{exp}(V_{eff}(\stackrel{}{r},\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2)/k_BT).$$
(6)
The angular integrations are trivial in the case of the SCM, where $`V_{eff}`$ depends on $`r`$. In the case of the DCM, one may use the identity

$$b(r)=\mathrm{exp}\left[_r^{\mathrm{}}𝑑r^{}\frac{d}{dr^{}}[\mathrm{ln}b(r^{})]\right],$$
(7)
to show that $`B_2`$ may be cast in a form similar to that appropriate for the SCM, namely

$$B_2=\frac{1}{2}𝑑\stackrel{}{r}\left[1\mathrm{exp}\{V(r)/k_BT\}\right]$$
(8)
where $`V(r)`$ is the potential of the orientationally averaged projected force, as defined in Eq. (3). As pointed out earlier, $`B_2`$ is directly accessible experimentally by extrapolating light scattering data to small wavevectors or by taking derivatives of osmotic pressure data with respect to concentration. Results will be presented in the form of the reduced second virial coefficient $`B_2^{}=B_2/B_2^{(HS)}`$ where $`B_2^{(HS)}=2\pi \sigma _p^3/3`$, i.e.

$$B_2^{}=1+\frac{3}{\sigma _p^3}_{\sigma _p}^{\mathrm{}}r^2𝑑r\left[1\mathrm{exp}\left\{V(r)/k_BT\right\}\right].$$
(9)
## III Simulation details

We study a pair $`(N_p=2)`$ of spherical proteins with center-to-center separation $`r`$, confined in a cubic box of length $`L=4\sigma _p`$, which also contained monovalent co and counterions in numbers determined by their bulk concentrations and overall charge neutrality. There are $`ZN_p`$ counterions dissociated from the protein surface, and added $`N_s`$ salt ion pairs such that the screening of proteins is implemented via $`N_+=N_s`$ coions and $`N_{}=N_s+ZN_p`$ counterions in simulation box. A snapshot of a typical equilibrium microion configuration around two proteins is shown in Figure 1 for the protein charge number $`Z=15`$. The two proteins were placed symmetrically with respect to the center along the body diagonal of a cubic simulation cell; periodic boundary conditions in three dimensions were adopted. $`L`$ was chosen such that the box length is much larger than the range of the total (effective) protein-protein interaction, so that the results are independent of $`L`$ for non-zero salt concentration. The long-range electrostatic interactions between two charged particles in the simulation box with periodic boundary conditions were modified using the Lekner summation method of images . For our model to be a rough representation of lysozyme, we chose $`\sigma _p=4nm`$, and three different protein charges $`Z=6,10`$ and $`15`$, corresponding to three different values of the solution pH. The microion diameter is fixed to be $`\sigma _c=\sigma _p/15=0.267nm`$.

For both the SCM and the DCM, the contact coupling parameter between a protein and a microion, namely $`\mathrm{\Gamma }=2e^2/[ϵk_BT(\sigma _p\sigma _d+\sigma _c)]`$ for the DCM, and $`\mathrm{\Gamma }=2Ze^2/[ϵk_BT(\sigma _p+\sigma _c)]`$ for the SCM, are comparable, and of the order of $`\mathrm{\Gamma }3`$ at room temperature. We fixed the dielectric constant of water to be $`ϵ=81`$ and the system temperature to be $`T=298K`$. Varying salt concentration for fixed protein charge $`Z`$ corresponds to a fixed solution pH.

Details of the runs corresponding to different salt concentrations are summarized in Table I. Note that the Debye screening length $`r_D`$, defined by

$$r_D=\sqrt{\frac{ϵk_BTV^{^{}}}{8\pi (N_s+Z)q_se^2}},$$
(10)
is less than $`10\AA `$ for salt concentration beyond 0.1M. Here

$$V^{^{}}=V\frac{\pi N_p}{6}(\sigma _p^3+\sigma _c^3),$$
(11)
is the accessible volume for salt ions such that the salt concentration $`C_s`$ is $`N_s/V^{^{}}`$. Thus, the point charges on the protein surface are effectively screened from each other . For each of the runs indicated in Table I, the distance-resolved effective forces and interaction potentials are calculated according to Eqs. (2,3). The statistical averages over microion configurations leading to $`\stackrel{}{F}_1^{(2)}`$ and $`\stackrel{}{F}_1^{(3)}`$ were evaluated from time averages in the Molecular Dynamics (MD) simulations.

## IV Microion distributions around a single protein

First, as a reference, consider a single protein ($`N_p=1`$) placed at the center of the simulation box. We calculated spherically averaged, radial microion density profiles $`\rho (r)=\rho _+(r)+\rho _{}(r)`$ in the immediate vicinity of the protein surface. For a single neutral sphere in a salted solution, results for $`\rho (r)`$ are drawn in Figure 2. There is a marked depletion in the microion density, signaled by a minimum of $`\rho (r)`$, well below the asymptotic bulk value. The depletion is enhanced upon increasing the salt concentration. At sufficiently high salt concentrations, this minimum is followed by a weak, but detectable, ion layer (see corresponding lines for runs 7 and 9 in Figure 2). The formation of a depletion zone is not a consequence of the direct (hard core) interaction between salt ions and the protein surface, since the position of the observed layer is significantly further away from the protein surface than one ion diameter. A rough estimate for the distance between layer and neutral sphere gives a value of $`2.5\sigma _s`$, or equivalently $`0.17\sigma _p`$. For runs 7 and 9, where the ion layer emerges, this distance is of the order of an average ion separation $`r_s`$ in the system and twice the Debye screening length $`r_D`$ as well (see Table I). Obviously, it is the small ion correlations which lead to the peak formation in the salt density profiles. As an intuitive argument, we posit that the lack of mutual polarization in dense salt solution near neutral surfaces causes the ion depletion. Qualitatively similar depleted density profiles were observed in Lennard-Jones system confined between neutral planes  and in Yukawa mixtures . Furthermore, an effective force that pushes a single ion toward regions of higher salinity is predicted within Debye-Hückel theory for interfacial geometries .

Next we consider a protein sphere with charge number $`Z=10`$. The density profiles of small ions are shown in Figure 3 for both the SCM and DCM. Figure 4 represents the corresponding total salt densities, as a sum of co- and counterion densities from Figure 3. At the lower salt concentration (up to run 5) the SCM and DCM models both yield an accumulation of the microion density near ion-protein contact, in semi-quantitative agreement with the prediction of standard Poisson-Boltzmann (PB) theory. For rising ionic strength the total microion density gets depleted near the protein surfaces, as in the previously considered case of a neutral sphere. Remarkably, this depletion occurs both with the SCM and DCM and contradicts the PB prediction. The intuitive picture is that a microscopic layer of counterions is formed around the proteins. An additional salt pair now profits more from the bulk polarization than from the protein surface polarization and is thus excluded from this layer. By normalizing the profiles to the total bulk density, this effect becomes visible as a depletion zone in Figure 4, where a noticeable difference between the SCM and DCM profiles also emerges. Whereas the DCM predicts a contact value $`\rho _c(r=(\sigma _p+\sigma _c)/2)`$ larger than the bulk value, SCM predicts a much stronger microion depletion near contact. Together with this, the contact value of the DCM model is always larger than that of the SCM model for the same salt concentrations. This finding illustrates the sensitivity of correlation effects to the assumed charge pattern at the surface of a protein. This correlation effect is, of course, absent in the (non-linear) PB and DLVO theories, which always predict a monotonically decreasing density profile $`\rho (r)`$ (see Figure 5). Within DLVO theory the density of plus and minus salt ions near protein surface in the SCM model are defined as

$`\rho _+(r)={\displaystyle \frac{Z_{DLVO}}{q_s}}{\displaystyle \frac{k_+^2}{4\pi }}{\displaystyle \frac{e^{k_Dr}}{r}}+{\displaystyle \frac{2\overline{\rho }_+\overline{\rho }_{}}{\overline{\rho }_++\overline{\rho }_{}}},`$ (12)
$`\rho _{}(r)={\displaystyle \frac{Z_{DLVO}}{q_s}}{\displaystyle \frac{k_+^2}{4\pi }}{\displaystyle \frac{e^{k_Dr}}{r}}+{\displaystyle \frac{2\overline{\rho }_+\overline{\rho }_{}}{\overline{\rho }_++\overline{\rho }_{}}}.`$ (13)
where

$`\overline{\rho }_+={\displaystyle \frac{N_s+2Z}{V^{^{}}}},\overline{\rho }_{}={\displaystyle \frac{N_s}{V^{^{}}}},Z_{DLVO}=Z{\displaystyle \frac{\mathrm{exp}(k_D\sigma _p/2)}{1+k_D\sigma _p/2}},`$ (14)
$`k_D^2=k_+^2+k_{}^2,k_{}^2=4\pi q_s^2e^2\overline{\rho }_{}/ϵk_BT,k_+^2=4\pi q_s^2e^2\overline{\rho }_+/ϵk_BT.`$ (15)
Nonlinear PB equations were solved via iterations of the potential of a homogeneously charged sphere placed at the center of a spherical cell. The cell radius was determined by the given protein concentration.

A direct comparison between SCM, DCM models and non-linear PB theory is shown in Fig. 6 for two of the higher salt concentrations from Fig. 4. For the intermediate salt concentration $`C_s=0.206`$ Mol/l (run 4) both simulation and theory predict a monotonic decrease of salt density away from the protein surface, whereas, in the case of a dense salt, $`C_s=0.824`$ Mol/l (run 7), simulation results strongly deviate from the PB prediction.

A multipole expansion of the total salt number density in the DCM, discussed in Appendix A, demonstrates that the higher order expansion coefficients are strongly damped and much weaker than the zero-order homogeneous term shown in Figs. 3, 4 and 6.

### A Effective force and $`B_2`$ for a protein pair

Next we calculate the angularly averaged effective interaction force $`F(r)=\frac{dV(r)}{dr}`$ and potential $`V(r)`$ between two proteins embedded in a sea of small salt ions. Simulation results for the simpler case of the SCM, are plotted in Figure 7 for $`Z=10`$ and compared to the DLVO predictions. There is a systematic deviation between the theoretical and simulation results. While the DLVO theory  potential

$$U^{(DLVO)}(r)=\frac{Z_{DLVO}^2e^2}{ϵr}\mathrm{exp}(r/r_D),$$
(16)
always results in repulsive forces, simulations indicate the possibility of an attraction between proteins for large salt concentrations. The force $`F(r)`$ at the higher salt concentrations $`C_s`$, shows a maximum at a distance $`r`$ nearly equal to the ion diameter. Note that, for the highest salt concentration considered, $`C_s=2.061`$ Mol/l (run 11), where the electrostatic interactions are almost completely screened out, the effective force $`F(r)`$ is dominated by entropic effects; it is reminiscent of entropic depletion force of hard sphere system. The corresponding potential is negative at short distances as shown in the inset of Figure 7, and is related to the depletion in the microion total density profiles $`\rho (r)`$ around an isolated protein, shown in Figure 4. We note that such an entropic attraction is not contained in DLVO theory. Its origin is also different from the salting out effect studied in  or the macroion overcharging effect studied in. In Figure 8 the salt dependence of the total interaction force $`F(r)`$ (Eq.(2)) is broken down into its components $`F^{(2)}`$ and $`F^{(3)}`$ for two values of $`r`$. This helps to show that at large salt concentrations, it is indeed the entropic component that causes the force to be attractive for run 11 in Figs 7 and 8a. Finally, we mention that the range of attraction observed here will depend on the electrolyte (salt ion) size. This feature of our model may hint at a cause for the salt specificity observed in salting out experiments on protein crystallization. In addition, as shown in , an entropic attraction from the electrolyte could lead even to phase separation in colloid-electrolyte mixtures. Including this term leads to an effective Hamaker constant – describing the dispersion interactions – that is lower, and in better accord with experimental findings .

The same calculations were carried out for the other two protein charges for the SCM model, $`Z=6`$ and $`Z=15`$, with qualitatively similar results to those obtained for $`Z=10`$.

It is clear that the effective forces and potentials between two proteins will no longer be spherically symmetric within the DCM model. For example, three distinguishable mutual orientations of the two proteins are schematically outlined in Figure 9a, corresponding to particular configurations of the Euler angles $`\stackrel{}{\mathrm{\Omega }}_1,\stackrel{}{\mathrm{\Omega }}_2`$ of the two proteins. Nevertheless, our simulation results, presented in Figure 9b, for the three orientations, show that the actual anisotropy of force is weak. However, this is no longer true for the Yukawa segment model , as shown in the inset of Figure 9b. Within this model, the total effective interaction potential between a pair of protein spheres is given by

$$U^{(YS)}(r)=\frac{1}{Z^2}\underset{k,n=1}{\overset{Z}{}}U^{(DLVO)}(\stackrel{}{r}_k\stackrel{}{r}_n),$$
(17)
where $`\stackrel{}{r}_k`$ and $`\stackrel{}{r}_n`$ represent the positions of the point unit charges of different proteins. We emphasize that the aelotopic (or nonisotropic) interactions incorporated in our DCM differ from those considered, for example, in Ref., where $`B_2`$ is calculated for a set of hydrophobic attractive patches on protein surface. Within our version of the DCM the third configuration in Fig. 9 (solid line), has the highest statistical weight of the three cases pictured (see Eqs.(3,4)). If, on the other hand, the point charges on the protein are replaced by attractive patches, then the configuration with two points nearly touching (dot-dashed line in Figure 9), is the statistically most favorable conformation. Similar arguments hold within a molecular model for site-specific short-range attractive protein-protein interactions,. Results for distance-resolved forces within the DCM model are shown in Figure 10a, for $`Z=10`$. When the salt concentration is less than $`C_s0.2`$ Mol/l, the results are similar to those of the SCM model: i.e. for low ionic strength, the force is repulsive, while for high ionic strength there is an attraction near contact followed by a repulsive barrier. The distinguishing property of the DCM is the nonmonotonicity of the force with the increase of ionic strength. This, in turn, gives rise to the nonmonotonic behavior of the spherically averaged interaction potential $`V(r)`$ shown in Figure 10b. This feature of $`V(r)`$ manifests itself in the following way in Figure 10b: the potential is first strongly reduced as $`C_s`$ is increased, then its amplitude and range increase very significantly at intermediate concentrations ($`C_s1`$ Mol/l), before it nearly vanishes at the highest salt concentrations. Note that $`V(r)`$ even becomes slightly attractive at contact ($`r=\sigma _p`$) for $`C_s2`$ Mol/l. Similar effects are also observed for $`Z=6`$ and $`Z=15`$ (see Figure 11), suggesting that the effect is generic for discrete charge distributions.

Once the effective potential $`V(r)`$ is known, it is straightforward to calculate the second osmotic virial coefficient using Eq.(8). In doing so, however, one should keep in mind that it is the total interaction that enters $`B_2`$. Real proteins also exhibit an additional short-range interaction, as seen, for example, in experimental studies of the osmotic pressure and structural data for lysozyme, or in fits to its phase-behavior. This attraction stems from hydration forces, van-der-Waals interactions, and other molecular interactions that are, to a first approximation, independent of salt concentration. Hence, we have taken the expected short-range attraction between proteins into account by adding to the effective Coulomb potential in Eq.(9), an additional “sticky” sphere potential of the Baxter form,

$$\frac{V_{SHS}(r)}{k_BT}=\{\begin{array}{cc}\mathrm{}\hfill & r\sigma _p\hfill \\ \mathrm{ln}\left[\frac{12\tau \delta }{\sigma _p+\delta }\right]\hfill & \sigma _p<r<\sigma _p+\delta \text{ ,}\hfill \\ 0\hfill & r\sigma _p+\delta \hfill \end{array}$$
(18)
with potential parameters $`\delta =0.02\sigma _p`$ and $`\tau =0.12`$, which yield reasonable osmotic data for lysozyme solutions in the high salt concentration regime. This square well potential is isotropic by nature and ignores the directionality in hydrophobic attraction between proteins . Short range attractions lead to “energetic fluid” behavior, where the crystallization is driven primarily by the details of the interactions, instead of being dominated by the usual entropic hard-core exclusions. This suggests that the directionality may be very important to details of the protein crystallization behavior. However, for the physically simpler behavior of the virial coefficient, the directionality can be ignored as a first approximation. For simplicity, we assume the parameter $`\tau `$ to be independent of electrolyte conditions, although a weak dependence based on experimental observations is reported in . The addition of $`V_{SHS}(r)`$ strongly magnifies the nonmonotonicity of $`B_2`$ stemming from the nonmonotonic behavior of $`V(r)`$ near contact.

Results for $`B_2^{}`$ as a function of salt concentration are shown in Figure 12 for three different protein charges. There is a considerable qualitative difference between the predictions of the SCM and the DCM models for the variation of $`B_2^{}`$ with monovalent salt concentration $`C_s`$ for each protein charge $`Z`$. Whereas the SCM (dashed curves in Figure 12) predicts a monotonic decay of $`B_2^{}`$ with $`C_s`$, the DCM leads to a markedly non-monotonic variation, involving an initial decay toward a minimum (salting-out) followed by a subsequent increase to a maximum (salting-in) and a final decrease at high $`C_s`$ values (salting-out). The location of the local minima shifts to higher/lower values of $`C_s`$ for larger/smaller protein charges $`Z`$. Thus for larger protein charge one needs a higher salt concentration to achieve the “salting out” conditions conducive to protein crystallization . Even though the effective Coulomb potential between proteins is small, with an amplitude only a few percent of the thermal energy $`k_BT`$, its effect on $`B_2`$ is dramatically enhanced by the presence of the strong short-range attractive Baxter potential.

We have also compared the effective potentials shown in figure 12 with a recently proposed scaling collapse of protein osmotic virial coefficients . This scaling effect, observed for a number of experimental conditions, can be explained with simple arguments based on Donnan equilibrium. To lowest order, the effects of salt concentration and protein charge on $`B_2`$ are to subsumed in the following approximate scaling relation:

$$B_2^{(0)}=B_2Z^2/4C_s,$$
(19)
where $`B_2^{(0)}`$ is the bare virial coefficient, independent of charge effects. As shown, for example in Figure 1 of reference , this simple relation hold remarkably well above a salt concentration of $`C_s0.25M`$ for a wide range of experimental measurements of $`B_2`$ for lysozyme, which all tend to a plateau value of $`B_2^0/B_2^{(HS)}(2.7\pm 0.2)`$. One implication of this observed scaling is that the attractive interactions that govern $`B_2^{(0)}`$ are indeed roughly independent of salt concentrations above $`C_s0.25M`$. When we applied the same scaling procedure to our $`B_2`$ curves, a similar plateau develops for both the DCM and the SCM models, albeit with $`B_2^{(0)}`$ slightly less negative than that found in the experiments. One could, of course, very easily match to experiments by adjusting the value of $`\tau `$, but to keep contact with our earlier work, we don’t do so. Clearly the scaling does bring the DCM and SCM $`B_2`$’s close together for a given $`Z`$, but for different $`Z`$ (related to solution pH), the scaling collapse is not as good as that seen in experiments, since we observe a larger $`C_s`$ before it sets in. Nevertheless, considering the high density of co and counter-ions in the simulation, it is remarkable that a simple Donnan theory based on ideal gas terms performs so well.

The origin of the non-monotonic variation of $`B_2^{}`$ with $`C_s`$ can be traced back to the subtle correlation effects which cause an enhancement of the effective Coulomb repulsion at intermediate salt concentrations in the DCM. These effects cannot be rationalized in terms of simple mean-field screening arguments . The protein-microion correlations are of a different nature to those in the SCM, where they lead to a much more conventional, monotonic decay of $`B_2`$ with $`C_s`$, similar to that expected from a simple screening picture.

In order to gain further insight into the physical mechanism responsible for the unusual variation of the effective Coulomb potential and of $`B_2`$ with salt concentration in the DCM, we consider the influence of a second near-by protein on the microion distribution near protein-ion contact. We have computed the difference between “inner” and “outer” shell microion contact densities for $`Z=10`$, as schematically illustrated in the inset to Figure 14. The local microion density is no longer spherically symmetric, due to the interference of the electric double-layers associated with the two proteins. The difference $`\mathrm{\Delta }\rho =\rho _{in}\rho _{out}`$ between the mean number of microions within a fraction of a spherical shell of radius $`R=0.6\sigma _p`$ subtended by opposite $`60^{}`$ cones, is plotted in Figure 14 versus salt concentration. $`\mathrm{\Delta }\rho `$ is always positive, indicating that microions ( mainly counterions) tend to cluster in the region between the proteins, rather than on the opposite sides. This follows because they can lower the total electrostatic energy by being shared between two proteins. However, there is a very significant difference in the variation of $`\mathrm{\Delta }\rho `$ with salt concentration $`C_s`$, between the SCM and the DCM models. Both exhibit similar behavior for lower salt concentrations $`C_s0.5`$ Mol/l; for example, both show a small maximum around $`0.2`$ Mol/l. But for salt concentrations above $`0.5`$ Mol/l, the SCM predicts a monotonic decrease of $`\mathrm{\Delta }\rho `$, while the DCM leads to a sharp peak in $`\mathrm{\Delta }\rho `$ for $`C_s1`$ Mol/l. This highly non-monotonic behavior clearly correlates with the non-monotonicity observed in Figs. 1011 and 12. The basic mechanism can be summarized as follows: For the DCM, the excess number of microions between the two proteins leads to an excess entropic pressure or force, as demonstrated in Figure 15, which is the origin of the increased repulsion between proteins around $`C_s=1`$ Mol/l. The enhanced microion density arises from subtle crowded charge correlation effects that cannot easily be understood at a mean-field level.

## V conclusion

In conclusion, we have calculated the effective interactions and the second osmotic virial coefficient $`B_2`$ of protein solutions incorporating the electrostatics within the “primitive” model of electrolytes. In this way we include nonlinear screening, overscreening, and correlation effects missed within the standard Poisson-Boltzmann description. For discrete charge distributions, the interactions and related $`B_2`$ vary in a non-monotonic fashion for increasing ion strength while for the smeared charge model, a standard workhorse of colloidal physics, this effect is absent. These correlation induced effects are missed within non-linear PB theory, and similar coarse-graining techniques taken from the theory of colloids. In addition to this, our simulations indicate the necessity of taking entropic forces into account when treating systems on the nanoscale. These forces are believed to be essential in the salting-out effect  and could lead to an attraction even between neutral globular proteins .

Our MD calculations can easily be extended to the more complex (pH dependent) charge patterns of realistic proteins. In fact, in some cases it may be easier to do a full MD simulation than to solve the non-linear PB equations in a very complicated geometry. We expect similar mechanisms to those found for the DCM to be active there, leading, for example, to an enhanced protein-protein repulsion at intermediate salt concentration. Since the second osmotic virial coefficient determines much of the excess (non-ideal) part of the chemical potential of semi-dilute protein solutions, we expect the non-monotonicity of $`B_2`$ to have a significant influence on protein crystallization from such solutions in the course of a “salting-out” process. The non-monotonic behavior also suggests the possibility of an inverse, “salting-in” effect, whereby a reduction of salt concentration may bring $`B_2`$ into the “crystallization slot”. The sensitivity of $`B_2`$ to ion-correlation effects may help explain the salt specificity of the Hofmeister series . Finally, we stress that our non-monotonicity is qualitatively different from that observed for added non-adsorbing  and absorbing  polymers.

###### Acknowledgements.

The authors are grateful to R. Piazza, I.L. Alberts, P.G. Bolhuis, G. Bricogne, J. Clarke, S. Egelhaaf, J.F. Joanny, D. Rowan, R.Blaak and W.C.K. Poon for useful discussions, and to Schlumberger Cambridge Research and the Isaac Newton Trust for financial support.

## A Multipole expansion of counterion density around protein in DCM

We expand the numerically calculated counterion density around a protein with charge $`Z`$ in a Laplace series of spherical harmonics:

$$\rho (r,\theta ,\phi )=\underset{n,m}{}C_{nm}(r)P_n^m(\mathrm{cos}(\theta ))\mathrm{cos}(m\phi )+S_{nm}(r)P_n^m(\mathrm{cos}(\theta ))\mathrm{sin}(m\phi )$$
(A1)
The multipole spherical expansion coefficients (MSEC) are calculated during an MD simulation via

$`C_{00}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i),`$ (A2)
$`C_{10}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)\mathrm{cos}(\theta _i),`$ (A3)
$`C_{11}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)\mathrm{sin}(\theta _i)\mathrm{cos}(\phi _i),`$ (A4)
$`S_{11}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)\mathrm{sin}(\theta _i)\mathrm{sin}(\phi _i),`$ (A5)
$`C_{20}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i){\displaystyle \frac{3\mathrm{cos}^2(\theta _i)1}{2}},`$ (A6)
$`C_{21}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)3\mathrm{sin}(\theta _i)\mathrm{cos}(\theta _i)\mathrm{cos}(\phi _i),`$ (A7)
$`S_{21}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)3\mathrm{sin}(\theta _i)\mathrm{cos}(\theta _i)\mathrm{sin}(\phi _i),`$ (A8)
$`C_{22}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)3\mathrm{sin}^2(\theta _i)\mathrm{cos}(2\phi _i),`$ (A9)
$`S_{22}`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (\stackrel{}{r}\stackrel{}{r}_i)3\mathrm{sin}^2(\theta _i)\mathrm{sin}(2\phi _i),`$ (A10)
where $`i`$ runs over counterions.

The protein charge is chosen to be $`Z=12`$, so that the rotation symmetry axis through two surface charges has a fivefold symmetry. The $`xz`$ and $`xy`$ plane projections of the protein charge pattern are shown in Figure 16.

The variation of the nonzero MSEC versus distance from the protein center are shown in Figure 17.