text:
# An Eternal Time Machine in 2+1 Dimensional anti-de Sitter Space

## I Introduction

Our time machine is produced by two particles orbiting a common center in circular orbits. In section two, we describe a useful coodinate system in adS. In section three we discuss the nature of the wedge identifications corresponding to the two particles. In section four, we construct the particle orbits and illustrate the resultant spacetime. In section five we examine the causal structure of the solution, and in the final section we compare with previous solutions.

2+1 dimensional anti-de Sitter space has been the subject of much recent investigation. Holst hol96  argued that two point masses on nearly radial orbits would create a Gott-type time machine as they passed each other near the center, providing that their relative velocities exceeded a certain mass-dependent threshold related to that for the original Gott configuration in Minkowski space gott . Matschull mat99  noted that two colliding particles would create a BTZ black hole ban92 . Holst & Matschull hm99  found that two light-like sources coming in from infinity and passing each other at a finite distance would create a rotating BTZ black hole with regions of closed timelike curves hidden inside the event horizon, and a “wormhole” connection to an additional universe. Many of the results of investigations persuing the connection between single point sources and the BTZ solution are summarized in Birmingham & Sen bs00 .

## II anti-de Sitter space

Anti-de Sitter space is a homogenous, isotropic space with a negative cosmological constant; 2+1 dimensional anti-de Sitter space can be represented as a three dimensional hyperboloid:

$$u^2v^2+x^2+y^2=r_0^2$$
(1)
embedded in a flat four dimensional space with the metric:

$$ds^2=du^2dv^2+dx^2+dy^2$$
(2)
where $`r_0^2=1/\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is the cosmological constant as usually defined in the Friedmann equations.

This space is not simply connected and contains closed timelike curves (i.e., the circles $`u^2+v^2=w^2r_0^2`$, with $`x`$ and $`y`$ constant and $`w^2>1`$.) These can be “unwrapped” to produce the covering space of anti-de Sitter space which we shall discuss in the remainder of this paper, and simply call anti-de Sitter space. This space has no closed timelike curves.

There are a number of ways to put coordinates on the hyperboloid; we shall be interested in a set that makes the conformal transformation clear. “Barrel coordinates,” discussed in Hawking & Ellis he73 , chapter $`5.2`$ and discussed for the 2+1 dimensional case by Holst & Matschull hm99 , are the most useful for our situation.

We define $`\tau `$, a time coordinate, $`\chi `$, a “radial” coordinate, and $`\varphi `$, an angular coordinate as follows:

$`u`$ $`=`$ $`r_0\mathrm{cos}\tau \mathrm{cosh}\chi `$
$`v`$ $`=`$ $`r_0\mathrm{sin}\tau \mathrm{cosh}\chi `$
$`x`$ $`=`$ $`r_0\mathrm{sinh}\chi \mathrm{cos}\varphi `$
$`y`$ $`=`$ $`r_0\mathrm{sinh}\chi \mathrm{sin}\varphi `$
the metric induced on the surface is then:

$$ds^2=r_0^2(d\chi ^2+\mathrm{sinh}^2\chi d\varphi ^2\mathrm{cosh}^2\chi d\tau ^2)$$
(3)
The time coordinate, $`\tau `$, runs from $`\mathrm{}`$ to $`\mathrm{}`$, $`\varphi `$, the angular coordinate is periodic, going from $`0`$ to $`2\pi `$, and $`\chi `$ goes from $`0`$ to $`\mathrm{}`$. We can set $`\mathrm{tan}\theta =\mathrm{sinh}\chi `$; since $`\chi `$ goes from $`0`$ to $`\mathrm{}`$, $`\theta `$ goes from $`0`$ to $`\pi /2`$. The new metric is then:

$$ds^2=r_0^2\left(\frac{d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2d\tau ^2}{\mathrm{cos}^2\theta }\right)$$
(4)
This demonstrates that anti-de Sitter space is conformally related to the product of a hemisphere and a real time axis. Performing the conformal transformation, and including the equator as a boundary at infinity, we find the “barrel” coordinates that we will use in this paper:

$$d\stackrel{~}{s}^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2d\tau ^2$$
(5)
which is related to the original metric by $`ds^2=r_0^2d\stackrel{~}{s}^2/\mathrm{cos}^2\theta `$. A property of adS is that null lines starting at the center can reach infinity and return in a finite coordinate time $`\mathrm{\Delta }\tau =2\pi `$.

We note that it is impossible to bring the timelike infinities $`\iota ^+`$ and $`\iota ^{}`$ in to a finite distance without collapsing the spacelike infinities to the origin; since we will only be dealing with processes that elapse over a finite coordinate time, this subtlety will not be an issue.

Each spacelike slice perpendicular to the time axis (i.e., $`\tau `$ constant) is thus mapped onto a hemisphere. A conformal stereographic projection (with $`\rho =\mathrm{tan}\theta /2`$) of the half sphere on to the plane:

$$ds^2=\frac{d\rho ^2+\rho ^2d\varphi ^2}{(1\rho ^2)^2}$$
(6)
will project this hemisphere onto a Poincaré disk of radius $`\rho =1`$.

## III Wedge identifications on the Poincaré disc

The effect of a point mass in $`2+1`$ dimensional Minkowski space is to remove a “wedge” with angle proportional to the mass of the particle. The two faces of the wedge are identified at constant times in the particle’s rest frame. This operation produces a conical singularity at the position of the particle, but leaves the rest of the space flat.

There is an analogous operation for introducing point masses in adS; a general formulation is described in Matschull mat99 . For our purposes, we need only know how to introduce a stationary point mass at the origin; our spacetimes will be constructed by a procedure of pasting together boosted solutions of this nature.

Our barrel metric gives us a notion of slices of constant time $`\tau `$ throughout the manifold; as we have seen above, each slice can be conformally mapped onto a Poincaré disk. The identification procedure is then simple: for a point mass of deficit angle $`\alpha `$ at the origin $`\chi =0`$, one identifies two geodesics in the slice of constant time $`\tau `$, emerging from the point particle, that are separated by $`\alpha `$ radians of rotation. This produces a conical singularity at the origin. In the limit as one approaches the origin ($`\chi <<1`$), the curvature term from the cosmological constant is negligible and the geometry near the origin approximates a piece of Minkowski space with a missing wedge of angular size $`\alpha `$ radians.

## IV The construction of the eternal time machine

Our time machine consists of two point masses orbiting a common center in circular orbits. We will construct it by finding the spacetime for a single point mass, boosting into the rest frame of what will be the common center, and putting together two such copies in such a way that the wedges do not intersect.

We begin with a description of the wedge identification lines emerging from single particle at the origin. At each slice of time $`\tau `$, there is a missing slice of spacetime of angle $`\alpha `$; the cross section at time $`\tau `$ looks like a pizza with a slice missing. The location of this missing slice in azimuth is arbitrary. Thus, we may rotate the azimuth of the location of the missing pizza slice at each epoch to produce a “rotating” missing wedge whose azimuth is equal to time $`\gamma `$ and makes one rotation in a period $`\gamma =2\pi `$. The two edges of the rotating missing wedge may then be described parametrically as follows:

$`u=\mathrm{cos}\gamma \mathrm{cosh}\lambda `$ $`u^{}=\mathrm{cos}\gamma \mathrm{cosh}\lambda `$
$`v=\mathrm{sin}\gamma \mathrm{cosh}\lambda `$ $`v^{}=\mathrm{sin}\gamma \mathrm{cosh}\lambda `$
$`x=\mathrm{cos}(\gamma \alpha /2)\mathrm{sinh}\lambda `$ $`x^{}=\mathrm{cos}(\gamma +\alpha /2)\mathrm{sinh}\lambda `$
$`y=\mathrm{sin}(\gamma \alpha /2)\mathrm{sinh}\lambda `$ $`y^{}=\mathrm{sin}(\gamma +\alpha /2)\mathrm{sinh}\lambda `$
where $`(u,v,x,y)`$ are the coordinates of an event on the trailing edge, and $`(u^{},v^{},x^{},y^{})`$ are the coordinates of the identified event on the leading edge, $`\gamma `$ is the time coordinate, equal to the proper time of the particle, $`\alpha `$ is the deficit angle, and $`\lambda `$, which runs from zero to infinity, parametrizes the distance along the edge, from the particle ($`\lambda =0`$) to infinity ($`\lambda =\mathrm{}`$). Identifications are to be made between points with the same $`\lambda `$ and $`\gamma `$ values. The $`\gamma `$ dependence in the $`x`$ and $`y`$ coordinates is what “rotates” the wedges.

We now boost this particle twice, by the same Lorentz boost, $`\psi `$, but in different directions; the first boost is in the $`xu`$ plane, the second in the $`yv`$ plane. In the new coordinate system, the equations of the two edges ($`\pm `$) are:

$`u`$ $`=`$ $`\mathrm{cos}\gamma \mathrm{cosh}\lambda \mathrm{cosh}\psi +\mathrm{cos}(\gamma \pm \alpha /2)\mathrm{sinh}\lambda \mathrm{sinh}\psi `$
$`v`$ $`=`$ $`\mathrm{sin}\gamma \mathrm{cosh}\lambda \mathrm{cosh}\psi +\mathrm{sin}(\gamma \pm \alpha /2)\mathrm{sinh}\lambda \mathrm{sinh}\psi `$
$`x`$ $`=`$ $`\mathrm{cos}(\gamma \pm \alpha /2)\mathrm{sinh}\lambda \mathrm{cosh}\psi +\mathrm{cos}\gamma \mathrm{cosh}\lambda \mathrm{sinh}\psi `$
$`y`$ $`=`$ $`\mathrm{sin}(\gamma \pm \alpha /2)\mathrm{sinh}\lambda \mathrm{cosh}\psi +\mathrm{sin}\gamma \mathrm{cosh}\lambda \mathrm{sinh}\psi `$
Again, points with the same $`\gamma `$ and $`\lambda `$ are to be identified. After the boosts, the $`\gamma `$ parameter remains the proper time.

To get a better intution for the behaviour of this spacetime, consider the line parametrized by $`\gamma `$ with $`\lambda =0`$ (i.e., the path of the particle itself.) The combination of these two boosts puts the particle into a circular orbit about the new origin of the coordinate system; the particle’s worldline is a helix. We note that the proper time of the particle, $`\gamma `$ is equal to the coordinate time $`\tau `$:

$`\varphi `$ $`=`$ $`\mathrm{arctan}(y/x)=\gamma `$
$`\tau `$ $`=`$ $`\mathrm{arctan}(v/u)=\gamma `$
The angular position of the particle as a function of $`\tau `$, in the barrel coordinates, is independent of the boost parameter $`\psi `$. To build our time machine, we shall introduce a second particle, $`\pi `$ radians out of phase with the first. The two particles will then helix around each other. The two particles orbit each other eternally without loss of energy through emission of gravity waves because there are no gravitational waves in $`2+1`$ dimensional spacetime.

Point particles in $`(2+1)`$ dimensional space exert no gravitational attraction for each other; the particles are kept in circular orbit by the overall gravitational attraction of the negative cosmological constant. (A net attraction is produced because the repulsion produced by the negative mass density of the cosmological constant is more than compensated by the gravitational attraction produced by its positive pressure in two dimensions, giving an overall attraction.)

\[We note in passing that, with this insight, we can write down a set of “corotating” coordinates, identical to the system described above but with the angular coordinate $`\varphi `$ replaced by $`\varphi ^{}=\varphi \tau `$. This gives the following metric (after performing the same conformal transformations as for the original “barrel” coordinates):

$$ds^2=\rho _0^2\left(\frac{d\theta ^2+\mathrm{sin}^2\theta (d\varphi ^2+2d\varphi ^{}d\tau )}{\mathrm{cos}^2\theta }d\tau ^2\right)$$
(7)
With this choice, there is a class of freely falling test particles, orbiting as do our orbiting point masses, that remain at the same radial and azimuthal positions for all times. While this system leads to needless complications for our purposes, it is interesting to note that such a coordinate choice leads one to view the adS manifold as a stationary, rotating universe. There are two such coordinate systems, the second being $`\varphi ^{}=\varphi +\tau `$, and rotating in the opposite sense. In these global coordinate systems, adS appears as a negatively curved, open, rotating – but not expanding – universe.\]

As shown in figure 2, the way in which we rotate the wedges in the particle rest frame means that we can do this without worrying about the wedges intersecting in a way that would complicate the simple wedge prescription of the one particle case. Our particluar choice of orientations for the wedges means that the “leading” and trailing” wedge faces are positioned symmetrically about the line emerging from the origin in the direction of each particle.

## V The causal Structure of the eternal time machine

We now wish to examine the causal structure of our system. Equation IV, along with our definitions of the barrel coordinates and our rule for identifying points has given us a full description of the spacetime. We can get a good handle on exactly what is happening by looking at the behaviour of the spacetime at spacelike infinity. In a manner similar to Matschull mat99 , we will examine the behaviour of null lines that remain always at infinity, and will then be able to make more general statements about the rest of the space.

We can find the identifications at infinity by taking $`\lambda `$ to infinity in the limit. In this limit, for a mass with deficit angle $`\alpha `$, the wedge identifications are given by:

$`\tau `$ $`=`$ $`\mathrm{arctan}\left({\displaystyle \frac{\mathrm{sin}\gamma +\mathrm{sin}(\gamma \pm \alpha /2)\mathrm{tanh}\psi }{\mathrm{cos}\gamma +\mathrm{cos}(\gamma \pm \alpha /2)\mathrm{tanh}\psi }}\right)`$ (8)
$`\varphi `$ $`=`$ $`\mathrm{arctan}\left({\displaystyle \frac{\mathrm{sin}\gamma \mathrm{tanh}\psi +\mathrm{sin}(\gamma \pm \alpha /2)}{\mathrm{cos}\gamma \mathrm{tanh}\psi +\mathrm{cos}(\gamma \pm \alpha /2)}}\right)`$ (9)
Where the point corresponding to a particular value of $`\alpha `$, $`\psi `$, and $`\gamma `$ is to be identified with the point specified by the same $`\psi `$ and $`\gamma `$, and with the $`\alpha `$ terms having opposite sign. For convenience, we have written these expressions using the arctangent; care must be exercised when dealing with points for which $`\varphi >\pi /2`$ or $`\varphi <\pi /2`$.

Following Holst & Matschull hm99 , we take the conformal boundary, and, cutting it along $`\varphi =0`$, unwrap it to form a plane. The wedge edges then appear as a parallel lines on this plane. In figure 3 we show this construction for $`\psi =0.25`$ and $`\alpha =\pi `$; the arrows indicate the identifications. We note that, at infinite distance from the origin, the rotating wedge edges move at a phase velocity equal to the speed of light in order to keep pace with the particle in the interior.

The heavy solid line in this figure shows the path of a light ray at infinity, rotating around the origin in the opposite sense to the two point masses (i.e., clockwise from the $`y`$ to the $`x`$ axis.) We see how, for sufficiently large deficit angles, the wedge identifications (depicted by dashed lines) result in the null ray closing on itself. For even larger deficit angles, when the identifications span a larger portion of $`(\tau ,\varphi )`$ space, these closed null lines become lines that propagate backwards in time, and one can now draw closed timelike curves on the surface at infinity.

The criteria for closed timelike lines can be found as follows. We require that $`\mathrm{\Delta }\varphi `$ and $`\mathrm{\Delta }\tau `$, the amount by which the light ray is identified along its direction of propagation in $`\varphi `$, and the amount by which the light ray is identified backward in $`\tau `$, respectively, for each wedge are sufficient to recover the time taken by the light ray to circle at infinity. In other words, the light ray can become closed by a combination of being retarded in time ($`\mathrm{\Delta }\tau `$) and by being advanced (clockwise) around the cylinder ($`\mathrm{\Delta }\varphi `$) across each wedge identification. If the two wedges taken together thus retard the light ray by a total of $`2\pi `$, the ray will close. More than that, and closed timelike curves at infinity result. Thus, to make closed timelike curves at infinity we require for each wedge:

$$\mathrm{\Delta }\tau +\mathrm{\Delta }\varphi >\pi $$
(10)
The expressions for $`\mathrm{\Delta }\tau `$ and $`\mathrm{\Delta }\varphi `$ would appear from equations 8 and 9 to be rather complicated. Surprisingly, however, they turn out to be very simple. We can examine the case where $`\gamma `$ is equal to zero (i.e., at the beginning of a cycle); by helical symmetry, our expressions for $`\mathrm{\Delta }\varphi `$ and $`\mathrm{\Delta }\tau `$ will be valid for all values of $`\gamma `$.

For this case, we have

$`\mathrm{\Delta }\tau `$ $`=`$ $`2\mathrm{arctan}\left({\displaystyle \frac{\mathrm{sin}(\alpha /2)\mathrm{tanh}\psi }{1+\mathrm{cos}(\alpha /2)\mathrm{tanh}\psi }}\right)`$
$`\mathrm{\Delta }\varphi `$ $`=`$ $`2\mathrm{arctan}\left({\displaystyle \frac{\mathrm{sin}(\alpha /2)}{\mathrm{tanh}\psi +\mathrm{cos}(\alpha /2)}}\right)`$
We remind the reader that, as mentioned before, in evaluating these expressions we must keep track of which quadrant we are in (using the signs of the numerator and denominator terms in parentheses) to recover the full range $`(0,2\pi )`$ of $`\mathrm{\Delta }\tau `$ and $`\mathrm{\Delta }\varphi `$, since the arctangent function has only the range $`(\pi /2,\pi /2)`$. We can find the expression for the sum $`\mathrm{\Delta }\tau +\mathrm{\Delta }\varphi `$ by judicious use of the tangent sum formula. We find, after much simplification, for each wedge, the result:

$$\mathrm{\Delta }\tau +\mathrm{\Delta }\varphi =\alpha $$
(11)
The dependence on $`\psi `$ has dropped out. Combining equation 11 with 10, we thus find that closed lightlike curves will be produced at infinity by two idenitical orbiting point masses when the sum of their deficit angles is equal to $`2\pi `$. If the sum exceeds $`2\pi `$, then closed timelike curves at infinity will result.

Now for point masses in 2+1 dimensional Minkowski space we expect that if the sum of the deficit angles of the two point masses exceeds $`2\pi `$ (more than a hemisphere of curvature), the spacetime closes like a “dunce cap” rather than extending to infinity. This induces additional mass points (at the “tip” of the dunce cap) so that the total mass in the closed space equals $`4\pi `$ as expected from the Gauss-Bonnet theorem in a closed space where the only curvature comes from the point masses. A spacelike section through anti-de Sitter space, however, has uniform negative curvature (it maps onto a Poincaré disk.) This negative curvature counteracts the positive curvature in the two point masses and allows the spacetime to remain open and extend to infinity even when the sum of the two point masses exceeds $`2\pi `$.

So we *are* allowed to have the sum of the point mass deficit angles that sum to greater than $`2\pi `$ while keeping a subset of points at infinity, by simply separating the two masses by sufficiently large boosts. In figure 4, we show this construction (produced with Mathematica) for $`\alpha =3\pi /2`$ for each point mass, and $`\psi =1`$. It is clear from this figure that the two wedges do not intersect, and that there are still two “stripes” of spacetime at infinity.

We now derive a condition on $`\psi `$ for the two wedges not to intersect or touch. To derive this condition, it will be convenient to go into a reference frame where one of the particles (the “stationary particle”) is at the origin. Then, the other particle (the “orbiting particle”) will orbit at a boost parameter of $`2\psi `$. The particle at the origin remains stationary, and its wedge, which we will call the “stationary particle wedge,” subtends an angle $`\alpha `$ with identifications at constant $`\tau `$.

Without loss of generality, we consider the case $`\gamma =0`$ (all other values of $`\gamma `$ can be obtained by rotating about the origin), and we remind the reader that our $`\varphi `$ coordinate increases going counterclockwise, starting from $`0`$ on the line going out from the origin through the orbiting particle at $`\gamma =0`$.

The lines of identification for the orbiting particle at $`\gamma =0`$ then start at $`\tau =0`$, $`\varphi =0`$ at $`\lambda =0`$, and run to some maximum $`\tau _{\text{max}}`$ and $`\varphi _{\text{max}}`$ as $`\lambda `$ goes to $`\mathrm{}`$ (we consider, again, without loss of generality because of the reflection symmetry, the leading wedge line only.)

The condition that the wedges never intersect becomes the following: the leading wedge identification line of the orbiting particle may not intersect any part of the stationary particle’s wedge in the range $`\tau =(0,\tau _{\text{max}})`$.

The stationary particle’s wedge planes, by construction, move around the origin with an angular velocity of unity. For the leading line on the leading edge of the orbiting particle, $`\tau (\lambda )`$ is always greater than $`0`$, and increasing monotonically moving out from the orbiting particle (i.e., as $`\lambda `$ increases.) We can thus write the angular position of the trailing stationary wedge plane as a function of $`\tau (\lambda )`$, the parameter along the leading line (if the leading line intersects the stationary particle’s wedges, it will do so first at the trailing wedge plane.)

We find:

$$\varphi _{\text{trail}}=(\pi \frac{\alpha }{2})+\tau (\lambda )$$
(12)
which shows the location of the trailing edge of the stationary particle at time $`\tau (\lambda )`$, where

$$\tau (\lambda )=\mathrm{arctan}\left(\frac{\mathrm{sin}(\alpha /2)}{\mathrm{coth}\lambda \mathrm{coth}\psi +\mathrm{cos}(\alpha /2)}\right)$$
(13)
Meanwhile, the $`\varphi `$ coordinate of the leading line on the leading wedge plane of the orbiting particle is written:

$$\varphi _{\text{lead}}=\mathrm{arctan}\left(\frac{\mathrm{sin}(\alpha /2)}{\mathrm{coth}\lambda \mathrm{tanh}\psi +\mathrm{cos}(\alpha /2)}\right)$$
(14)
Our condition for non-intersection of the wedges then reads:

$$\varphi _{\text{trail}}>\varphi _{\text{lead}}$$
(15)
for all values of $`\lambda `$ between 0 and $`\mathrm{}`$. In other words, the line parametrized by $`\lambda `$ for $`\gamma =0`$ on the leading edge of the orbiting particle should never hit the trailing wedge of the stationary particle.

We can see immediately that in the limit $`\psi \mathrm{}`$, the two arctangent terms become equal, and $`\varphi _{\text{trail}}\varphi _{\text{lead}}\pi \alpha /2`$. In other words, we can fit two particles with deficit angles of up to just below $`2\pi `$ each into our space, provided we separate them by a sufficient boost $`\psi `$.

For the case $`\alpha =\pi `$, the bounding condition on the deficit angles to produce a time travel region at infinity, we consider a small boost $`ϵ`$. Expanding equation 12 above to first order in $`ϵ`$, we find:

$$\varphi _{\text{trail}}=2ϵ\mathrm{tanh}\lambda $$
(16)
For small $`\lambda `$, $`\varphi _{\text{lead}}`$ goes to $`0`$ as long as $`ϵ>0`$, and we see that $`\varphi _{\text{trail}}\varphi _{\text{lead}}\pi /2`$, as expected. For large $`\lambda `$, $`\mathrm{coth}\lambda `$ is near unity, and we can again expand to first order in $`ϵ`$ to find:

$$\varphi _{\text{trail}}=2ϵ\mathrm{coth}\lambda $$
(17)
Thus, for large $`\lambda `$,

$$\varphi _{\text{trail}}\varphi _{\text{lead}}2ϵ(\mathrm{coth}\lambda +\mathrm{tanh}\lambda )$$
(18)
and thus a time machine where the masses have deficit angles of $`\pi `$ each requires only an infinitesimally small boost. This compares to the results of an earlier investigation of closed time-like curves in AdS for particles in near radial orbits mat99 ; there, in a situation similar to that of the original results for closed time-like curves from cosmic strings and point masses in a Minkowski background gott , it was found that two equal point masses passing each other on linear trajectories through the origin of AdS produce CTCs when their velocities relative to the center of mass exceed a certain amount:

$$\frac{v}{c}>\mathrm{cos}(\alpha /2)$$
(19)
where $`\alpha \pi `$.

So far, we have shown only the existence of closed lightlike curves on the surface at infinity. We now wish to show that this leads to the existence of closed timelike curves throughout the space.

We note that, without the missing wedges, test particles double boosted by $`\psi `$ orbit the center of the system with period $`2\pi `$, regardless of their radial distance from the center. Let us denote by $`f(\theta )`$ the amount, in $`\tau `$ and $`\varphi `$, that a double boosted particle, rotating at constant distance, is beamed back in time due to the missing wedges, where $`\theta `$ is the radial coordinate. For deficit angles that sum to greater than $`2\pi `$, $`f(\pi /2)`$ will be greater than $`2\pi `$. Since $`f`$ is a continuous function of $`\theta `$, this implies that there are closed timelike curves a finite distance away from the origin (i.e., $`f(\theta )>2\pi `$ for some $`\theta <\pi /2`$.) When the deficit angles sum to precisely $`2\pi `$, there exist only closed null lines at infinity.

Since our solution is stationary (in the rotating coordinate system described in equation 7 above), we can then draw a closed timelike curve through any event when the deficit angles sum to more than $`2\pi `$. We can draw a timelike curve between any two events; the timelike curve simply orbits the two point masses as many times as needed to “rewind” sufficiently backwards in $`\tau `$ to reach the event designated. (Recall that, from any finite distance you can reach the center in a time $`\mathrm{\Delta }\tau <\pi /2`$. So to go from any $`(\chi ,\varphi )`$ to any other $`(\chi ^{},\varphi ^{})`$ requires only a finite $`\mathrm{\Delta }\tau `$.)

## VI Discussion

We have constructed, in anti-de Sitter space, a system of two particles orbiting about a common center that leads naturally to the existence of closed time-like and null curves. We have derived analytic results on the constraints for the masses and relative separations of the two particles. Closed timelike curves fill the entire spacetime. Timelike curves may be drawn through any two events within the space; our “time machine” exists for all times.

Holst & Matschull hm99 , who also considered the effect of point particles in 2+1 dimensional adS, suggested a connection between their solution and the BTZ solution ban92 . The CTCs in their solution are hidden inside the event horizon of the rotating black hole. The CTCs in our solution, however, are not found in a single region bounded by Cauchy horizons, but penetrate the entire space. In their solution the point particles that generate the CTCs travel on lightlike paths and are “emitted” and “absorbed” at infinity at a time $`\mathrm{\Delta }\tau =\pi `$ apart. Regions before and after are connected with a BTZ wormhole with a time machine hidden inside. In our case, the particles simply orbit each other for all time, creating an eternal time machine.

SD acknowledges support from an NSF Graduate Research Fellowship and JRG acknowledges support from NSF grant AST-9900772.


format: text

text:
# Study of the Cabibbo-suppressed decay modes 𝐷⁰→𝜋⁻𝜋⁺ and 𝐷⁰→𝐾⁻𝐾⁺

## Abstract

Using data from the FOCUS (E831) experiment at Fermilab, we present a new measurement for the branching ratios of the Cabibbo-suppressed decay modes $`D^0\pi ^{}\pi ^+`$ and $`D^0K^{}K^+`$. We measured:

$`\mathrm{\Gamma }(D^0K^{}K^+)/\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)=2.81\pm 0.10\text{(stat)}\pm 0.06\text{(syst)}`$,

$`\mathrm{\Gamma }(D^0K^{}K^+)/\mathrm{\Gamma }(D^0K^{}\pi ^+)=0.0993\pm 0.0014\text{(stat)}\pm 0.0014\text{(syst)}`$, and

$`\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)/\mathrm{\Gamma }(D^0K^{}\pi ^+)=0.0353\pm 0.0012\text{(stat)}\pm 0.0006\text{(syst)}`$.

These values have been combined with other experimental data to extract the ratios of isospin amplitudes and the phase shifts for the $`DKK`$ and $`D\pi \pi `$ decay channels.

$`\text{The FOCUS Collaboration}^{}`$ thanks: See http://www-focus.fnal.gov/authors.html for additional author information. , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , ,

1. Introduction

In recent years, hadronic decays of charm mesons in many decay modes have been extensively studied. The theoretical models that have been developed mainly describe the 2-body decay modes and these models have led to several successful predictions. However the branching ratio $`\frac{\mathrm{\Gamma }(D^0K^{}K^+)}{\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)}`$ has, for a long time, been a puzzle of charm physics. The decay modes $`D^0K^{}K^+`$ and $`D^0\pi ^{}\pi ^+`$ are both Cabibbo-suppressed and in first order perturbative calculation both receive contributions from the same diagrams (external spectator and exchange). To first order in the SU(3) flavour symmetry limit, the above branching ratio should be one. However this ratio is reduced by a factor $`0.86`$ due to a phase space difference and increased by a factor of $`(f_K/f_\pi )^2=1.49`$ because of the different decay constants of the kaon and the pion. An overall ratio of $`1.29`$ is thus expected. Including SU(3) breaking effects, the expected ratio can increase to an upper limit of about 1.4 .

However, the measured ratio is close to $`2.5`$ . Both penguin diagrams and final state interactions (FSI) have been proposed as sources of such a high value. However the penguin interference  seems to be unable to explain this large value. Some theoretical models  propose the FSI as the solution of this puzzle.

In this paper we present a new measurement of the $`\frac{\mathrm{\Gamma }(D^0K^{}K^+)}{\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)}`$ branching ratio obtained using data from the FOCUS experiment as well as an isospin analysis of the $`DKK`$ and $`D\pi \pi `$ decay channels.

FOCUS is a charm photoproduction experiment  which collected data during the 1996–97 fixed target run at Fermilab. Electron and positron beams (with typically $`300\text{GeV}`$ endpoint energy) obtained from the $`800\text{GeV}`$ Tevatron proton beam produce, by means of bremsstrahlung, a photon beam which interacts with a segmented BeO target. The mean photon energy for triggered events is $`180\text{GeV}`$. A system of three multicell threshold Čerenkov counters performs the charged particle identification, separating kaons from pions up to $`60\text{GeV}/c`$ of momentum. Two systems of silicon microvertex detectors are used to track particles: the first system consists of 4 planes of microstrips interleaved with the experimental target  and the second system consists of 12 planes of microstrips located downstream of the target. These detectors provide high resolution in the transverse plane (approximately $`9\mu \text{m}`$), allowing the identification and separation of charm primary (production) and secondary (decay) vertices. The charged particle momentum is determined by measuring their deflections in two magnets of opposite polarity through five stations of multiwire proportional chambers.

2. Analysis of $`D^0\pi ^{}\pi ^+`$ and $`D^0K^{}K^+`$

The final states are selected using a candidate driven vertex algorithm . A secondary vertex is formed from the two candidate tracks. The momentum of the resultant $`D^0`$ candidate is used as a seed track to intersect the other reconstructed tracks and to search for a primary vertex. The confidence levels of both vertices are required to be greater than $`1\%`$. Two estimators of the relative isolation of these vertices are returned by the algorithm: the first estimator (*Iso1*) being the confidence level that tracks forming the secondary vertex might come from the primary vertex, while the second estimator (*Iso2*) is the confidence level that other tracks in the event might be associated with the secondary vertex. Once the production and decay vertices are determined, the distance $`L`$ between them and its error $`\sigma _L`$ are computed. The quantity $`L`$ / $`\sigma _L`$ is an unbiased measure of the significance of detachment between the primary and secondary vertices. These variables provide a good measure of the topological configuration of the event, so that appropriate cuts on them reject the combinatorial background effectively.

In addition to large combinatorial backgrounds, $`D^0\pi ^{}\pi ^+`$ and $`D^0K^{}K^+`$ decays are difficult to isolate because of a large reflection from the Cabibbo favored $`D^0K^{}\pi ^+`$ decays. Particle identification requirements for each decay mode have been chosen to optimize signal quality. To minimize systematic errors on the measurements of the branching ratios, we use identical vertex cuts on the signal and normalizing modes. The only difference in the selection criteria among different decay modes lies in the particle identification cuts.

In the $`D^0\pi ^{}\pi ^+`$ and $`D^0K^{}K^+`$ analysis, we require $`L`$ / $`\sigma _L`$ $`>`$ $`10`$, *Iso1* $`<10\%`$ and *Iso2* $`<`$ 0.5 $`\%`$. We also require the $`D^0`$ momentum to be in the range $`25250\text{GeV}/c`$ (a very loose cut) and the primary vertex to be formed with at least two reconstructed tracks in addition to the seed track. The Čerenkov identification cuts used in FOCUS are based on likelihood ratios between the various stable particle identification hypotheses. These likelihoods are computed for a given track from the observed firing response (on or off) of all the cells that are within the track’s ($`\beta =1`$) Čerenkov cone for each of our three Čerenkov counters. The product of all firing probabilities for all the cells within the three Čerenkov cones produces a $`\chi ^2`$-like variable $`W_i=2\mathrm{ln}(\mathrm{Likelihood})`$ where $`i`$ ranges over the electron, pion, kaon and proton hypotheses . All the kaon tracks are required to have $`\mathrm{\Delta }_K=W_\pi W_K`$ (kaonicity) greater than $`1`$, whereas all the pion tracks are required to have $`\mathrm{\Delta }_\pi =W_KW_\pi `$ (pionicity) exceeding $`1.5`$. Using the set of selection cuts just described, we obtain the invariant mass distributions for $`\pi ^{}\pi ^+`$, $`K^{}K^+`$ and $`K^{}\pi ^+`$ shown in Fig. 1.

In Fig. 1a the $`\pi ^{}\pi ^+`$ mass plot shows a broad peak to the left of the signal peak due to surviving contamination from $`D^0K^{}\pi ^+`$ events. The shape of the reflection peak has been determined by generating Monte Carlo $`D^0K^{}\pi ^+`$ events and reconstructing them as $`\pi ^{}\pi ^+`$. The mass plot is fit with a function that includes a Gaussian for the signal, a third-order polynomial for the combinatorial background and a shape for the reflection background obtained from the Monte Carlo simulation. The amplitude of the reflection peak is a fit parameter while its shape was fixed. The low mass region is excluded in the fit to avoid possible contamination due to other charm hadronic decays involving an additional $`\pi ^0`$. A least squares fit gives a signal of $`3453\pm 111`$$`\pi ^{}\pi ^+`$ events.

The $`K^{}K^+`$ mass plot, shown in Fig. 1b, is fit with a function similar to that for the $`\pi ^{}\pi ^+`$ fit. In the $`K^{}K^+`$ case, the reflection peak is on the right of the signal. The fit gives a signal of $`10830\pm 148`$$`K^{}K^+`$ events.

The large statistics $`K^{}\pi ^+`$ mass plot of Fig. 1c is fit with two Gaussians111With the lower statistics of the $`K^{}K^+`$ and $`\pi ^{}\pi ^+`$ signals a single Gaussian gives an adequate fit. with the same mean but different sigmas to take into account the different resolution in momentum of our spectrometer  plus a second-order polynomial. The fit gives a signal of $`105030\pm 372`$$`K^{}\pi ^+`$ events.

The fitted $`D^0`$ masses are in good agreement with the world average  and the widths are in good agreement with those of our Monte Carlo simulation.

3. Relative Branching Ratios

The evaluation of relative branching ratios requires yields from the fits to be corrected for detection efficiencies, which differ among the various decay modes because of differences in both spectrometer acceptance (due to different $`Q`$ values for the various decay modes) and Čerenkov identification efficiency.

From the Monte Carlo simulations, we compute the relative efficiencies to be: $`\frac{ϵ(D^0\pi \pi )}{ϵ(D^0KK)}=0.897\pm 0.003`$, $`\frac{ϵ(D^0K\pi )}{ϵ(D^0KK)}=0.963\pm 0.003`$ and $`\frac{ϵ(D^0K\pi )}{ϵ(D^0\pi \pi )}=1.074\pm 0.004`$. Using the previous results, we obtain the following values for the branching ratios: $`\frac{\mathrm{\Gamma }(D^0KK)}{\mathrm{\Gamma }(D^0\pi \pi )}=2.81\pm 0.10`$, $`\frac{\mathrm{\Gamma }(D^0KK)}{\mathrm{\Gamma }(D^0K\pi )}=0.0993\pm 0.0014`$, and $`\frac{\mathrm{\Gamma }(D^0\pi \pi )}{\mathrm{\Gamma }(D^0K\pi )}=0.0353\pm 0.0012`$.

Systematic uncertainties on branching ratio measurements can come from different sources. We determine three independent contributions to the systematic uncertainty: the *split sample* component, the *fit variant* component, and the limited statistics of the Monte Carlo.

The *split sample* component takes into account the systematics introduced by a residual difference between data and Monte Carlo, due to a possible mismatch in the reproduction of the $`D^0`$ momentum and the changing experimental conditions of the spectrometer during data collection. This component has been determined by splitting data into four independent subsamples, according to the $`D^0`$ momentum range (high and low momentum) and the configuration of the vertex detector, that is, before and after the insertion of an upstream silicon system. A technique, employed in FOCUS and in the predecessor experiment E687, modeled after the *S-factor method* from the Particle Data Group , was used to try to separate true systematic variations from statistical fluctuations. The branching ratio is evaluated for each of the $`4(=2^2)`$ statistically independent subsamples and a *scaled variance* $`\stackrel{~}{\sigma }`$ (that is, where the errors are boosted when $`\chi ^2/(N1)>1`$) is calculated. The *split sample* variance $`\sigma _{\text{split}}`$ is defined as the difference between the reported statistical variance and the scaled variance, if the scaled variance exceeds the statistical variance:

$`\sigma _{\text{split}}`$ $`=\sqrt{\stackrel{~}{\sigma }^2\sigma _{\text{stat}}^2}`$ $`if\stackrel{~}{\sigma }>\sigma _{\text{stat}}`$ (1)
$`\sigma _{\text{split}}`$ $`=0`$ $`if\stackrel{~}{\sigma }<\sigma _{\text{stat}}.`$
Another possible source of systematic uncertainty is the *fit variant*. This component is computed by varying, in a resonable manner, the fitting conditions on the whole data set. In our study, we changed the background parametrization (varying the degree of the polynomial), the fit function for the reflection peak (the reflection shape from the Monte Carlo was replaced by a Gaussian), and the use of two Gaussian for the fit of the peak of $`D^0\pi ^{}\pi ^+`$ and $`D^0K^{}K^+`$. The values obtained by the various fits are all a priori likely, therefore this uncertainty can be estimated by the simple average of the measures of the fit variants:

$`\sigma _{\text{fit}}=\sqrt{{\displaystyle \frac{_{i=1}^Nx_i^2Nx^2}{N1}}}.`$ (2)
Finally, there is a further contribution due to the limited statistics of the Monte Carlo simulation used to determine the efficiencies. Adding in quadrature the three components, we get the final systematic errors summarized in Table 1:

The final results are shown in Table 2 along with a comparison with the previous determinations.

4. Isospin analysis of $`DKK`$ and $`D\pi \pi `$ channels

Final State Interactions (FSI) can dramatically modify the observed decay rates and complicate the comparison of the experimental data with the theoretical predictions. By means of the isospin analysis of the decay channels $`DKK`$ and $`D\pi \pi `$, we can gain some insight on the elastic component of the FSI (pure rotation in isospin space).

Let us consider the $`D\pi \pi `$ transistions: $`D^0\pi ^{}\pi ^+`$, $`D^0\pi ^0\pi ^0`$ and $`D^+\pi ^+\pi ^0`$. The decay amplitudes can be expressed in terms of isospin $`I_f=0(A_0)`$ and $`I_f=2(A_2)`$ amplitudes. The final state with isospin $`I_f=1`$ is forbidden by Bose statistics for an angular momentum zero system of two pions. We denote by $`A^+`$, $`A^{00}`$ and $`A^{+0}`$ the decay amplitudes for the $`D^0\pi ^{}\pi ^+`$, $`D^0\pi ^0\pi ^0`$ and $`D^+\pi ^+\pi ^0`$, respectively. Expressing the decay amplitude in terms of isospin amplitudes, we have :

$`A^+`$ $`=`$ $`+\sqrt{{\displaystyle \frac{2}{3}}}|A_0|\mathrm{exp}(i\delta _0)+\sqrt{{\displaystyle \frac{1}{3}}}|A_2|\mathrm{exp}(i\delta _2)`$ (3)
$`A^{00}`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{3}}}|A_0|\mathrm{exp}(i\delta _0)+\sqrt{{\displaystyle \frac{2}{3}}}|A_2|\mathrm{exp}(i\delta _2)`$ (4)
$`A^{+0}`$ $`=`$ $`+\sqrt{{\displaystyle \frac{3}{2}}}|A_2|\mathrm{exp}(i\delta _2).`$ (5)
Adding the decay amplitudes in quadrature, we find the ratio of the magnitude of isospin amplitudes and their relative phase shift difference in terms of measured branching fractions:222The relationship between the isospin amplitude and the branching fraction is $`\mathrm{\Gamma }^{ij}=\frac{1}{8\pi }\frac{p^{}}{M_D^2}A^{ij}^2`$, where $`p^{}`$ is the center of mass 3-momentum of each final particle .

$`{\displaystyle \frac{A_2}{A_0}}^2`$ $`=`$ $`{\displaystyle \frac{\frac{2}{3}A^{0+}^2}{A^+^2+A^{00}^2\frac{2}{3}A^{+0}^2}}`$ (6)
$`\mathrm{cos}(\delta _2\delta _0)`$ $`=`$ $`{\displaystyle \frac{3A^+^26A^{00}^2+2A^{+0}^2}{4\sqrt{2A^{+0}^2}\sqrt{\frac{3}{2}(A^{00}^2+A^+^2)A^{+0}^2}}}.`$ (7)
The decay rate $`\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)`$ has been determined from our measurement of the branching ratio $`\frac{\mathrm{\Gamma }(\pi ^{}\pi ^+)}{\mathrm{\Gamma }(K^{}\pi ^+)}`$, whereas the other $`D\pi \pi `$ decay rates and lifetimes have been taken from the Particle Data Group compilation .

The results are shown in Table 3. In contrast to $`K\pi \pi `$ decays, where the transitions are dominated by the $`I=1/2`$ amplitude ($`\mathrm{\Delta }I=1/2`$ rule), the $`A_2`$ amplitude in $`D\pi \pi `$ is comparable to the $`A_0`$ amplitude. Furthermore, there is a large phase shift difference between the isospin amplitudes. According to Watson’s theorem , this phase shift cannnot arise from the weak processes alone and thus constitutes direct evidence for FSI .

In the same way, we can consider the two-body $`DKK`$ transitions: $`D^0K^+K^{}`$, $`D^0K^0\text{}^0`$ and $`D^+K^+\text{}^0`$. The decay amplitudes $`A^{ij}`$ for the $`DK^iK^j`$ decay modes can be expressed in terms of $`A_0`$ and $`A_1`$ isospin amplitudes :

$`A^+`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(|A_1|\mathrm{exp}(i\delta _1)+|A_0|\mathrm{exp}(i\delta _0)\right)`$ (8)
$`A^{00}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(|A_1|\mathrm{exp}(i\delta _1)|A_0|\mathrm{exp}(i\delta _0)\right)`$ (9)
$`A^{+0}`$ $`=`$ $`\sqrt{2}|A_1|\mathrm{exp}(i\delta _1).`$ (10)
Using the previous decompositions, we can express the ratio of the magnitudes of the isospin amplitudes and their phase shift difference in terms of the measured branching fractions:

$`{\displaystyle \frac{A_1}{A_0}}^2`$ $`=`$ $`{\displaystyle \frac{A^{0+}^2}{2A^+^2+2A^{00}^2A^{+0}^2}}`$ (11)
$`\mathrm{cos}(\delta _1\delta _0)`$ $`=`$ $`{\displaystyle \frac{A^+^2A^{00}^2}{\sqrt{A^{+0}^2}\sqrt{2A^{00}^2+2A^+^2A^{+0}^2}}}.`$ (12)
The $`\mathrm{\Gamma }(D^0K^{}K^+)`$ decay rate has been determined from our measurement of the branching ratio $`\frac{\mathrm{\Gamma }(K^{}K^+)}{\mathrm{\Gamma }(K^{}\pi ^+)}`$, the $`\mathrm{\Gamma }(D^+K^+\text{}^0)`$ from a previous measurement of FOCUS  and the remaining decay rate from the Particle Data Group compilation .

The results are shown in Table 3. Analogously to the $`D\pi \pi `$ case, the two $`DKK`$ isospin amplitudes are of the same order of magnitude, although the isospin phase shift difference is smaller.

The isospin analysis of the $`DKK`$ and $`D\pi \pi `$ decay channels is summarized in Table 3 (the quoted errors are obtained adding in quadrature the statistical and systematic errors) along with a comparison to previous determinations by CLEO :

These results show that strong interactions, acting on the final particles, play a very important role in $`DKK`$ and $`D\pi \pi `$ decays, modifying the measured $`\frac{\mathrm{\Gamma }(K^{}K^+)}{\mathrm{\Gamma }(\pi ^{}\pi ^+)}`$ ratio.

Another way to see the elastic FSI effect on the $`\frac{\mathrm{\Gamma }(K^{}K^+)}{\mathrm{\Gamma }(\pi ^{}\pi ^+)}`$ branching ratio is to compute the ratio of the sums over the $`D^0`$ isospin rotated decay modes : $`\frac{\mathrm{\Gamma }(K^{}K^+)+\mathrm{\Gamma }(K^0\overline{K}^0)}{\mathrm{\Gamma }(\pi ^{}\pi ^+)+\mathrm{\Gamma }(\pi ^0\pi ^0)}`$. As opposed to the $`\frac{\mathrm{\Gamma }(K^{}K^+)}{\mathrm{\Gamma }(\pi ^{}\pi ^+)}`$ ratio, this ratio is not affected by elastic FSI.

Using these measurements for $`\frac{\mathrm{\Gamma }(K^{}K^+)}{\mathrm{\Gamma }(K^{}\pi ^+)}`$ and $`\frac{\mathrm{\Gamma }(\pi ^{}\pi ^+)}{\mathrm{\Gamma }(K^{}\pi ^+)}`$ and the PDG  values for the other modes, we compute:

$`{\displaystyle \frac{\mathrm{\Gamma }(K^{}K^+)+\mathrm{\Gamma }(K^0\overline{K}^0)}{\mathrm{\Gamma }(\pi ^{}\pi ^+)+\mathrm{\Gamma }(\pi ^0\pi ^0)}}=2.06\pm 0.24.`$ (13)
This ratio is lower than the $`\frac{\mathrm{\Gamma }(K^{}K^+)}{\mathrm{\Gamma }(\pi ^{}\pi ^+)}`$ branching ratio, but still above the expected value of $`1.4`$. Therefore, the elastic FSI cannot account for all the discrepancy between theory and experiments. An inelastic FSI that also allows the transition $`KK\pi \pi `$ seems to be the most reasonable explanation .

4. Conclusions

We have measured the following branching ratios: $`\frac{\mathrm{\Gamma }(D^0K^{}K^+)}{\mathrm{\Gamma }(D^0K^{}\pi ^+)}`$, $`\frac{\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)}{\mathrm{\Gamma }(D^0K^{}\pi ^+)}`$ and $`\frac{\mathrm{\Gamma }(D^0K^{}K^+)}{\mathrm{\Gamma }(D^0\pi ^{}\pi ^+)}`$. A comparison with previous determinations has been shown in Table 2. Our results improve significantly the accuracy of these measurements.

An isospin analysis of the decay channels $`DKK`$ and $`D\pi \pi `$ shows that final state interactions play an important role in these hadronic decay modes.

We wish to acknowledge the assistance of the staffs of Fermi National Accelerator Laboratory, the INFN of Italy, and the physics departments of the collaborating institutions. This research was supported in part by the U. S. National Science Foundation, the U. S. Department of Energy, the Italian Istituto Nazionale di Fisica Nucleare and Ministero della Istruzione Università e Ricerca, the Brazilian Conselho Nacional de Desenvolvimento Científico e Tecnológico, CONACyT-México, and the Korea Research Foundation of the Korean Ministry of Education.


text:
# Effective macroion charge and stability of highly asymmetric electrolytes at various salt conditions

## I Introduction

Electrostatic interactions in colloidal suspensions and their relation to their phase stability have always been attracting much attention because of their occurence in many biological and technological circumstances. It has been known for a long time that presence of the oppositely charged polyions seriously deteriorates the stability of colloidal dispersions. In polyelectrolyte-induced aggregation, a correspondence between the location of the fast aggregation regime and the colloid isoelectric point has been documented repeatedly Gregory ; Miclavic ; Grant ; Borkovec . Furthermore, the electrophoretic mobility of colloidal particles, which reflects the apparent particle charge, has been shown to be tightly correlated with the system’s stability. The attractive interactions due to surface charge heterogeneities have been thought to promote colloidal aggregation in colloid-polyelectrolyte systems (”charge patch flocculation”) Gregory ; Miclavic ; Borkovec . It however only recently became clear that apart from specific binding between different solutes and van der Waals forces, the electrostatic correlations themselves can lead to variety of exciting phenomena in systems contaning oppositely charged polyions such as giant charge inversion or short-range attraction between like-charged polyions (meaning either colloids or polyelectrolytes) Guldbrand ; Gronbech ; RB ; Hribar ; LowenAl ; Kjellander ; Bratko ; PRL ; Shklovskii ; Netz ; Grosberg ; Messina ; Grosberg1 ; Levin ; VlachyRev ; HansenRev . Moreover, computer simulations and theoretical studies have revealed that the correlation induced attraction can cause colloidal aggregation even when just small multivalent ions are present Hribar ; PRL ; Shklovskii ; Messina ; JCP2 . It was shown that the universality of the underlying phenomena is based on the enormous energetic contribution of the Coulombic interactions to the solution’s free energy, which dominates over the entropy and makes the system unsensitive to the internal structure of the counterions Kjellander ; Shklovskii ; Jure ; Marie ; Marie2 ; Carlsson . A phase diagram featuring a phase separation region around the macroion isoelectric point with a consequent re-stabilization was predicted for a purely Coulombic system containing polyions and a sufficient amount of multivalent counterions Grosberg ; Toan .

In the present work, we present a computer simulation study of charge inversion and phase instabilities in an asymmetric electrolyte treated by a multivalent salt. This model is meant as a generic representation of a solution containing two types of oppositely charged polyions. We hope to recover the basic properties of such a system originating from the electrostatics and translational degrees of freedom of the solutes. Our simulation setup resembles a common experimental situation when an initially stable suspension is destabilized by adding oppositely charged multivalent ions (or polyelectrolyte) and aimed to outline the range of activity of the “patch charge attraction” Gregory ; Miclavic ; Borkovec  as well as assess the theoretical predictions Grosberg ; Toan ; Harnau  for phase behaviour of asymmetric Coulombic systems at different electrostatic coupling parameters. This work also complements the recent computer simulation studies of these phenomena JCP2 ; Per1 ; Marie2 ; Carlsson  by the investigation into salt effects and the role of the effective macroion charge.

The paper is organized as follows. In Section II.1, we give the description of the model and parameter settings for the numerical experiment. Section III gives a detailed account of the results for three series of simulations with different charge asymmetries. Section IV presents discussion of the mechanisms of observed phase behaviour and its relation to effective interactions between the colloids. The conclusions are given in the final section.

## II Model and method

### II.1 Model

The systems under consideration are asymmetric electrolytes described within the framework of the primitive model. The reference electrolyte contains two types of spherical charged particles: (i) macroions of diameter $`\sigma _M=40`$ Å and charge $`Z_M=60`$ and (ii) small ions of diameter $`\sigma _I=4`$ Å and charge $`Z_I=+1`$, representing the counterions, whereas the solvent enters the model only by its relative dielectric permittivity $`\epsilon _r`$. The added salt consists of small ions of diameter 4Å: monovalent anions $`Z_a=1`$ and cations of different valencies. We classified the simulated systems by the added cation charge: the series A contains cations with $`Z_c=+1`$, the series B $`Z_c=+3`$, and the series C $`Z_c=+5`$. The two last chosen coupling strengths belong to the instability region on the generic phase diagram for the system including only macroions with their counterions PRL ; Per1 . The initial 60:1 system was previously studied in Per1 ; JCP1  and in deionized system with trivalent counterions a phase separation was found PRL ; JCP2 . We characterize the amount of added salt by a ratio of the overall added cation charge to the overall macroion charge, $`\beta =Z_c\rho _c/(Z_M\rho _M)`$, where $`\rho _i`$ is the number density of the corresponding species.

The interaction between the particles in our model is pairwise additive, and for pair $`ij`$, where $`i`$ and $`j`$ denote either polyion or counterion, it is given by

$$U_{ij}(r)=\{\begin{array}{cc}\mathrm{},& r<(\sigma _i+\sigma _j)/2\\ \frac{Z_iZ_je^2}{4\pi ϵ_0ϵ_r}\frac{1}{r},& r(\sigma _i+\sigma _j)/2,\end{array}$$
(1)
where $`r`$ is the center-to-center separation between the particles.

To characterize the intensity of the electrostatic correlations between counterions on the surface, we use the counterion-counterion coupling parameter $`\mathrm{\Gamma }=Z_I^2l_b/a_Z`$, where $`l_b=(e^2/4\pi \epsilon _0\epsilon _rkT)`$ is the Bjerrum length, $`e`$ is the elementary charge, $`\epsilon _0`$ the dielectric permittivity of vacuum, and $`a_Z=[Z_I/(\sigma /e)]^{1/2}`$ the average distance between two neighboring counterions at the charged surface characterized by the surface charge density $`\sigma `$. It is known that the correlation-induced attraction appears at $`\mathrm{\Gamma }>\mathrm{\Gamma }^{}2`$ RB ; JCP2 ; Per1 . Our systems A, B, and C correspond to $`\mathrm{\Gamma }=0.8`$, 4.1, and 8.8, respectively.

The systems are considered at a fixed macroion number density $`\rho _M=2.510^7`$Å<sup>-3</sup> corresponding to a macroion volume fraction $`\varphi _M=0.0084`$. A temperature of $`T=298`$ K and a relative dielectric permittivity of $`\epsilon _r=78.4`$ were used. For each of the three series, 11 values of $`\beta `$ ranging from 0.1 to 10 were studied.

### II.2 Method and simulation settings

In this study, we applied multiparticle Monte Carlo (MC) simulation method, which was developed in our earlier works JCP1 ; Springer . The canonical NVT ensemble with periodic boundary conditions and the Ewald summation for handling the electrostatic interactions were employed. The simulated bulk systems were composed of 20 macroions, 1200 monovalent counterions, and various amounts of added salt ions. Typically, $`0.5110^5`$ attempted MC moves per particles were made in the production runs, while for samples at $`\beta 1`$ we performed up to $`110^7`$ moves.

For the simulation of solutions with much salt we combined cluster MC moves for macroions JCP1  with swap moves, which is a modified implementation of the cluster moves. In the swap move, we first select the macroion and determine the particles within a concentric sphere around it. Then a second sphere of the same radius is taken centered at the destination vector of the macroion translational move $`\delta 𝐫_M`$ and a central inversion of the content of the whole volume within the two (possibly overlapping) spheres about the center of symmetry at $`\delta 𝐫_M/2`$ is performed. The acceptance ratio is calculated then in the same way as for the standard cluster move JCP1 , which results for a hard-sphere system and with 100% probability of the particle inclusion in the cluster in rejection of the moves that change the number of particles in the cluster (cf. discussion in JCP1 ). This type of move significantly improves the acceptance ratio for the systems with homogeneously distributed ions in the bulk as the macroions always move into an ion-free volume.

To evaluate the effective macroion charge and the osmotic coefficient, we applied the spherical Wigner-Seitz (WS) cell model. The model was solved by canonical MC simulations for one centrally placed macroion with corresponding amount of small ions, while the other settings were the same as for the bulk simulation. The simulation package MOLSIM v3.1 by Per Linse et al. was used throughout. Further details of the simulation protocol can be found in Refs. JCP1 ; Springer .

## III Results

### III.1 Single-macroion properties. Thermodynamics

We start discussing the numerical results from thermodynamic properties of the chosen electrolytes as function of the salt concentration. In Figs. 1-2 we plotted the reduced electrostatic energy of the solution $`U/Nk_BT`$, and the osmotic coefficient $`\mathrm{\Pi }/\rho k_BT`$, where $`N`$ and $`\rho `$ are the total number and the number density of the ionic species, respectively, $`k_B`$ the Boltzmann constant, $`\mathrm{\Pi }`$ the osmotic pressure. The reduced energy was calculated both from bulk simulation and WS cell model. The Fig. 1 shows that the total potential energy is negative at all the studied conditions as a result of strong attraction between the macroions and counterions. Its magnitude is decreasing monotonously with the salt content, which is a simple consequence of increasing the overall number of particles. All the curves have similar shape with varying curvature. An inflection point is observed at about $`\beta =1`$ on each curve. For the osmotic coefficient shown in Fig. 2 we present the cell model results obtained from the total density of ions at the cell boundary $`\mathrm{\Pi }_{cell}/(\rho k_BT)=\rho (R_{cell})/\rho `$. The osmotic coefficient behaves differently for systems A, and B and C. For the system A, it is increasing roughly proportionally to $`\mathrm{ln}\beta `$ up to $`\beta =1`$, and after that is more slowly approaching unity, the ideal solution value. In the two latter series it is again growing proportional to $`\mathrm{ln}\beta `$ at low salt but with slightly higher coefficient. Both curves for series B and C have a peak at $`\beta =1`$ and then decrease proportionally to $`\mathrm{ln}\beta `$, the decay of series C results being faster. The value of the osmotic coefficient at the peak is about 0.8. The proximity of the $`\beta `$ values to unity in all these cases reflects the high fraction of the free ions in the solution. This fraction can be roughly estimated as $`(N_a+N_c+Z_{eff}+N_M)/N`$, which gives at $`\beta =1`$ for complete adsorption of the added cations 121/181, 121/141, and 121/133, or 0.67, 0.86, and 0.91 for the systems A, B, and C, respectively. The further decrease of the osmotic coefficient is a result of association of salt anions with the multivalent cations, which reduces the number of ”free” species as compared to the ideal solution value.

### III.2 Single-macroion properties: Effective charge

The concept of effective polyion charge attracted much attention in literature due to its importance for the reduction problem. It immediately arises whenever the solvent- or ion-averaged interactions between solutes are concerned. The two major trends (excluding the experimental definitions) for two-component systems (macroion and counterions) are represented by the thermal and structural definitions. The different definitions in a spherically symmetric case regard the distance from the macroion center, within which the charge should be counted. The structural definitions are based on the knowledge of the radius, at which the mixed macroion-ion rdf falls below unity, or other criteria. The simplest of the thermal definitions counts the charge inside the macroion-ion binding radius ($`|U(r_{therm})|=kT`$). The more accurate one considers the chemical potential of the counterion in the macroion field. In the cylindrical case it was developed by Manning and Oosawa Manning ; Oosawa  and adopted for spherical geometry by Belloni Belloni . The radius in this approach is measured at the inflection point of the accumulated running charge curve plotted over inverse distance from the macroion center $`1/r`$ (see also refs.Belloni ; Deserno  for discussion of alternative definitions). Obviously, the radius is different for the ions of different valencies.

In the case of multicomponent electrolyte, however, we have to apply different definitions to different ionic species due to qualitatively dissimilar behaviour of the integrated running charge depending on the salt valency. Here, we use MC simulation to solve the WS cell model JCP1  and then extract the effective charge out of the integrated running charge curves. For series A and series B, C at low salt, we apply the inflection point method Belloni . For the At $`\beta >1`$, where the inflection point was not accessible anymore, we evaluated the effective charge at the first extremum of the integrated running charge curves.

The Fig. 3 demonstrates typical ionic distribution profiles for systems A, B and C. While all the curves decay monotonously in the series A, the systems B and C display diverse behaviour for the different ionic species. The number density of monovalent counterions decays similarly to the system A. The multivalent ion density has a high maximum at the polyion surface, while farther in the bulk it is close to zero. The coion distribution is peaked close to the maximum of the multivalent ion distribution because of their attraction to the macroions with the inverted charge and becomes very pronounced at the high salt concentrations Tanaka .

The integrated running charge curves are presented in Fig. 4. There is a distinct qualitative difference between system A and two other systems B, C. While in the former one the total charge is decaying monotonously and reaches zero only at the cell boundary, in the two latter systems the charge exhibit a sharp drop to zero within a very narrow region $`r<1.2\sigma _M`$. At $`\beta =1`$, the running charge remains close to zero up to the cell boundary. At the higher $`\beta `$, the charge is changing its sign and after that decaying slowly. Close to the cell boundary, one can see the cation charge build-up (note that the charge is divided by the negative $`Z_M`$), which originates from repulsion between the charge inverted macroions and cations. Note that the accumulated charge in series A at high salt concentrations reaches zero already at $`r=3\sigma _M`$, which is much less than the mean distance between the macroions.

In Fig. 5, we present the effective macroion charge as a function of $`\beta `$ for all three series. For all the systems $`Z_M^{eff}/Z_M`$ starts from about 0.3 at low added salt (its value for the 60:1 electrolyte JCP1 ) and then monotonously decreases. The effective charge changes its sign for systems B and C and regains roughly the same magnitude (0.17 for series B and 0.36 for series C) on the negative side. It is interesting to note the qualitative similarity of all the three curves. All of them have two relatively flat regions at very low and very high salt concentrations and a region of fast variation with an inflexion point at about $`\beta =1`$. The curve for series A asymptotically approaches but never crosses the zero line. Very slow dependence of the effective charge on the salt content at $`\beta >1`$ follows from logarithmical increase of the counterion entropy with counterion number density Toan . One should note that the decrease of the effective charge with the salt concentration, as it is seen for $`\beta <1`$, is characteristic for dilute and deionized systems, where the effective charge is far from its saturated value Belloni . At the higher monovalent salt content, $`\beta >1`$, the ionic double layer saturates and $`Z_M^{eff}/Z_M`$ starts growing with salt concentration towards unity Belloni ; Trizac . We could not see this upturn due to the failure of the inflection point criterion in this region. The failure is due to the fact that the underlying division of the ionic population into the ”condensed” and ”free” parts is no longer valid. The inflection point can still be detected formally but it refers now to the far field properties of the double layer (the inflection point is located next to the cell boundary) and does not reflect its short-range decay, which is important for the interaction between the macroions.

### III.3 Collective macroion properties: Structure

The object of our main interest here is the macroionic distribution. We therefore analyse only macroion-macroion partial distributions and extract those data that characterize the macroion density fluctuations. Namely, we consider the long-wave limit value of the macroion-macroion partial structure factor $`S_{MM}(q)`$, which would quantify the reduced osmotic compressibility in a corresponding one-component fluid of macroions. We should note that due to the small size of the simulated sample we can access only relatively short-wave region ($`q\sigma _M>0.01`$) and therefore have only rough estimate of $`S_{MM}(0)`$. One should expect stronger deviations from unity for the true $`S_{MM}(0)`$, i.e. even lower values for the stable solutions and the higher ones for the aggregating systems.

The macroion distribution in the bulk solution can be easiest illustrated by simulation snapshots shown in Fig. 6. The snapshot a) corresponds to a low added salt content $`\beta =0.1`$ and reveals a homogeneous spatial distribution of the macroions. In the snapshot b) corresponding to $`\beta =1`$ one can see macroion clusters. All the macroions are in aggregated state. Finally, snapshot c) demonstrates that at high salt $`\beta =5`$ most macroions are free and the solution is back to the stable state. The cluster re-dissolution is in agreement with theoretical predictions Toan  and recent simulations for polyelectrolyte-induced a aggregation of colloids and globular proteins Marie2 ; Carlsson . We should note that even at high $`\beta `$ values occasional particle pairs are formed. To ensure the equilibrium thermodynamic state of the system we performed additional simulations with initial compact configuration of the macroions. The macroions were allowed to move after the equilibration of the ionic subsystem. Then, indeed we observed stable aggregates at $`\beta 1`$ and re-dissolution of the macroions at the higher salt dose. The parameters of the aggregates were drifting slowly towards their values from the simulation with random initial configuration with an increase of simulation length. The problem with matching the two sets of data for the re-dissolution and the cluster-cluster aggregation is caused by the wide energy barriers originating from the large net charge of the clusters Toan2 . Although the cluster charge is not growing proportionally to the number of macroions due to the partial counterion release, the cluster charge yet reached values of $`2Z_M^{eff}`$ in our systems B and C at $`\beta >1`$.

Figure 7 presents partial macroion-macroion radial distribution functions (rdf) for the three systems at different salt concentrations. In Fig. 7a, one can see the results for system A, which show gradual decay of the macroion ordering on increasing salt. The peak of the rdf appears initially at $`r=\rho ^{1/3}`$ and has a value of 1.15, then moves towards shorter separations, while decreasing at the same time and almost disappears at ($`\beta >1`$). The trend is readily understood from decreasing the effective macroion charge and growing the screening ability of the electrolyte at the higher salt contents. For the remaining two series B and C, the trend is very different. Already at $`\beta =0.5`$, there is a peak near the macroion contact distance (or, more exactly, at $`r=\sigma _M+\sigma _c`$), which then grows high upon salt addition and reaches value of 60 at $`\beta =1`$ for the series B, and 150 for the series C. At $`\beta >1`$ it considerably drops in height but still remains quite high even at $`\beta =10`$. A secondary peak at $`r=2\sigma _M`$ corresponding to the second particle layer in a cluster is seen for some of the curves in series B and C at $`\beta >1`$. We would like to note the reappearance of the remote peak at $`r\sigma _M=5`$ for the series C at $`\beta =10`$. We made also another interesting observation for the system C. At the salt concentrations $`\beta 1`$, the main peak of the macroion-macroion rdf shifts towards $`r=\sigma _M`$ and is seen now at $`r=1.02\sigma _M`$ instead of $`r=\sigma _M+\sigma _c=1.1`$. This effect might be related to strengthening the electrostatic depletion contribution to the macroion interaction LowenAl .

We plotted some macroion-macroion partial structure factors for different salt concentrations in Figure 8. The structure factors for the series A again shows decay of the structure on increasing $`\beta `$. For the two remaining series, the structure factors have initially a peak at $`q\sigma _M=2`$ that corresponds to the inverse mean distance between the macroions and a low value in the long-wave limit. On increasing salt, this peak disappears and two new peaks develop instead at $`q\sigma _M0`$ and $`q\sigma _M7`$. The first peak reflects the average macroionic cluster size that grows with the salt addition until $`\beta 1`$ for series B and C. After that, when even more salt is added, the value of the long-wave peak decreases, which indicates the re-dissolution of the clusters. The broad secondary peak at $`q\sigma _M7`$ corresponds to the macroion-macroion density build-up close to contact distance.

The whole set of values of the long-wave structure factor value $`S(0)`$ for all three series is shown in Fig. 9. The curve for system A shows gradual increase between $`\beta =0`$ and 6 while the curves for systems B and C have well expressed maxima around $`\beta =1`$. For the system B, the peak is slightly shifted towards higher salt concentration due to the competition between the adsorption energy and the cation entropy, which leads to incomplete macroion charge compensation at $`\beta =1`$. For the series C, the $`S(0)`$ value returns to the values lower than 1, as it was at $`\beta 1`$. The overall trend for the system C corresponds to its destabilization at $`\beta 1`$ a consequent re-stabilization at high salt that occurs due to recovery of the repulsion between the macroions (cf. the effective charge plot in Fig. 5) Grosberg ; Toan .

## IV Discussion

As we have seen from the above results, the salt conditions are crucial for the stability of asymmetric electrolytes. On the path along the salt concentration axis, two important parameters of macroion interactions are affected: the Debye screening length and the effective charge of the macroions. The effect of the former is rather smooth and can be detected through a simple damping of the interactions upon salt addition. The sensitivity of the asymmetric electrolyte to the added salt depends also on the relative amounts of the “native” and “added” ions. We note that overall properties of the studied systems mark the charge compensation point $`\beta =1`$ either by an extremum (the osmotic coefficient) or an inflection (potential energy, effective charge), which can be understood as a crossover from the counterion dominated ($`\beta <1`$) to the salt dominated screening $`\beta >1`$.

Another interaction parameter, the effective charge, influences mostly the macroion related properties such as the macroion distribution and the stability of the macroion subsystem. In the counterion dominated regime, the effective charge is changing almost linearly with the salt concentration. This makes the system extremely sensitive to the addition of salt. In contrast, at the higher electrolyte strength the effective charge is increasing logarithmically with the amount of added multivalent salt. We note a close similarity between the observed charge behaviour and the macroion electrophoretic mobility vs added polyelectrolyte dose curves reported in Refs.Grant ; Borkovec . Fig. 10 demonstrates the direct relation between the effective charge and the stability of the electrolyte solution. In the plot one can see the relative effective charge curve together with the characteristic of stability $`S_{MM}(0)^1`$ for the series C. It is important to note the approximate symmetry of the curves in Fig. 10 relative to the isoelectric point, which allows one to formulate the result in a symmetric manner. While one can say that on the right-hand side of the diagram the macroions are overcharged by the excess multivalent counterions, on its left-hand side there is excess of macroions themselves and they are overcharging the multivalent ions. The latter gets a richer physical meaning when we imagine polyelectrolyte molecules instead of the multivalent ions. The effect of the internal counterion structure was discussed recently Jure ; Marie2 . As it was shown previously for the asymmetric electrolytes without added salt Per1 , the parameters of importance are the charge asymmetry between the oppositely charged ionic species and the ion correlation parameter $`\mathrm{\Gamma }`$. Extending our conclusions, we should expect that the strongly charged systems with similar energetic parameters have common features irrespectively of the internal structure of the polyions, whether they are polyelectrolytes or colloids.

While changing the salt valency leads to a smooth variation of the system’s properties such as the internal energy, osmotic coefficient, or effective macroion charge at the salt-overdosed conditions far from the isoelectric point, the interactions and structure change drastically in all the series at $`\beta =1`$. Our results confirm that the instability appears at a certain threshold counterion valency and becomes sharper at the higher counterion charge (or counterion coupling $`\mathrm{\Gamma }`$). Therefore, at the higher $`\mathrm{\Gamma }`$ one can expect on one hand a much stronger macroion aggregation at $`\beta =1`$ and, on another hand, a much better expressed overcharging and re-dissolution effects. Thus, the three systems of our choice correspond to (i) the weak coupling region and counterion entropy dominated behaviour (series A), (ii) border region demonstrating the coexistence of condensed and free counterions at $`\beta =1`$ (series B) and (iii) the high coupling (or “low temperature”) regime with completely condensed multivalent counterions (series C). The alllter case was analyzed in detail recently Toan . At the same time, the instability in real solutions should be much enhanced due to van der Waals forces, which due to their short-range character similarly play a role only in the vicinity of the isoelectric point, i.e. the whole sequence must be shifted towards lower $`\mathrm{\Gamma }`$ values Carlsson .

Finally, we would like to report on the simple salt effect on the effective charge and the counterion correlation attraction. The salt effects in similar systems were discussed previously JCP2 . It was confirmed that the added salt screens the repulsions between the macroions and leads to weakening the long-range macroion ordering. In a 60:2 electrolyte, addition of divalent ions caused an increase of the electrostatic attraction. In the 60:3 electrolyte, which is similar to our system B but has no monovalent counterions, no significant effect was noticed. Basing on the above analysis, we can explain the increase of the macroion aggregation in the 60:2 electrolyte by that the added divalent ions in fact lessen the effective macroion charge, which is not compensated completely at $`\beta =1`$ (similar to system B of the present work at $`\beta =1`$) PRE . The apparent salt effect in the 60:3 electrolyte was small because the considered salt concentrations placed the system in the center of the instability region, where it is not very sensitive to the salt. The salt concentration range in JCP2  corresponds to $`1<\beta <4`$ in the current work (cf. the curve for the system B in Fig. 9) Now, we repeat the comparisons of Ref. JCP2  between the solutions B and C at $`\beta =1`$ and corresponding systems with multivalent ions only, namely 60:3 and 60:5 electrolytes.

In the 60:3 system, the absolute potential energy per macroion is about 2.5% smaller than in series B, while the difference is about 3.3% between the 60:5 electrolyte and the system C. Similar relations hold for the partial electrostatic energies: the macroion-cation and cation-cation energies, respectively, are 3.3% and 5.5% higher in the 60:3 solution than in the B series at $`\beta =1`$. For the 60:5 solution and system C, these diffrences are marginal: the macroion-cation energy is 0.3% greater in the former, while the cation-cation energy is the same within the statistical uncertainty. The osmotic coefficient calculated in the cell model in the 60:5 system is equal to zero, which means that all the multivalent ions are adsorbed and the bulk osmotic coefficient is then just $`\mathrm{\Pi }/(\rho _Mk_BT)0.08`$. The simulation for the 60:3 electrolyte showed about 5 times higher cation concentration in the bulk for the corresponding system B, though both concentrations were very small. At the same time, the presence of the monovalent ions leads to a very modest decrease in correlation between the multivalent ions. The cation-cation correlation peak in system B falls by about 10% with respect to 60:3 electrolyte, which corresponds to a slight decrease of the effective correlation parameter $`\mathrm{\Gamma }`$. In the system C compared to the 60:5 system, the analogous decrease is just about 5%. The macroion-cation rdfs display marginal sensitivity to the added electrolyte. The peak height for the salt-free case is only 5% higher in the 60:3 system, and 2% higher in the 60:5 system than for the corresponding samples B and C.

The relative effective charge $`Z_M^{eff}/Z_M`$ for the both charge asymmetries increases in magnitude with the salt addition. With trivalent ions, it rises from $`0.015Z_M`$ in 60:3 system to $`0.03Z_M`$ in system B. For pentavalent ions, it changes from 0 in 60:5 electrolyte to $`0.01Z_M`$ in system C. Note that the presence of the anion layer around the adsorbed cations does not cancel the effect (see Fig. 3). A significant role here is played by the monovalent counterions, which also accumulate close to the macroion surface. The macroion structure displays considerable sensitivity to the presence of monovalent ions. In the 60:3 system, the first peak of the macroion-macroion rdf increases from 60 to 200 as compared to the system B (Fig.7b). The structure factor value $`S_{MM}(0)`$ rises from 2.1 to 13 (Fig.8b). For the case of pentavalent ions, the rdf peak rises from 150 to 230 (Fig.7c), while the initial structure factor peak from 6 to 13 (Fig.8c). One can see from the Fig.7bc and Fig.8bc that the aggregation is much stronger expressed in both 60:3 and 60:5 electrolytes than in the corresponding systems B and C. This drastic change can be attributed to (i) the decrease of the effective correlation parameter $`\mathrm{\Gamma }`$ for the multivalent ions in presence of monovalent ones and (ii) the change in the macroion effective charge, which has slightly higher magnitude for the systems with monovalent ions present. The accumulation of these small but finite macroion charges opposes the cluster-cluster aggregation and prevents formation of very large clusters.

## V Conclusion

We performed a numerical study of a highly asymmetric electrolyte solution containing strongly charged macroions and different amounts of salt. The addition of a simple monovalent salt gradually lowers the effective macroion charge and reduces the solution stability. Addition of small amounts of multivalent salt produces similar effect. The dose of multivalent salt that exceeds the macroion isoelectric concentration causes macroion charge inversion and electrostatically driven macroion aggregation. At the isoelectric point, we observe formation of large macroion clusters. When the inverted macroion charge becomes large enough, the clusters re-dissolve. Enlarging the counterion valency strenghens the instability at the isoelectric point but stabilises the solution at the overdosed conditions. The obtained stability diagram closely resembles polyelectrolyte-induced colloidal aggregation, while the aggregation mechanism is related solely to electrostatic correlation forces rather than the van der Waals attractions.

We are grateful to Per Linse, Boris Shklovskii, Toan Nguyen, Michal Borkovec, Christian Holm, and Peter Schurtenberger for stimulating discussions. This work was financially supported by the Swiss National Science Foundation.