text:
# Multipartite bound entangled states that violate Bell’s inequality

## Abstract

We study the relation between distillability of multipartite states and violation of Bell’s inequality. We prove that there exist multipartite bound entangled states (i.e. non-separable, non-distillable states) that violate a multipartite Bell inequality. This implies that (i) violation of Bell’s inequality is not a sufficient condition for distillability and (ii) some bound entangled states cannot be described by a local hidden variable model.

Since the celebrated work by Bell , it is evident that quantum mechanics is not compatible with local realist theories. In fact, Gisin  and Gisin and Peres  showed that all bipartite entangled pure states violate the Bell-CHSH inequality . This rules out the existence of a local hidden variable (LHV) model which is capable of describing all statistical correlations predicted by quantum mechanics for such states. These results were readily generalized by Popescu and Rohrlich  to multipartite entangled pure states of arbitrary dimension.

However, for quantum systems in mixed states, the situation is much more complicated and we still lack a complete classification of mixed states into “local” and “non–local” ones. Although structural knowledge about entanglement has increased rapidly in the last decade, many open questions remain to be answered, in particular concerning the relation of certain entanglement properties to the existence of LHV models.

While it is obvious that separable mixed states  can be described by a LHV model, Werner showed in his pioneering paper  that also a certain class of entangled mixed states —now called Werner states— do not violate any Bell–type inequality. This was done by explicitly constructing a LHV model which can simulate the results of any single (i.e. nonsequential) measurement performed at each side. However, some years later Popescu  realized that also sequential measurements can be considered, and showed that most of the Werner states exhibit violation of local realism if sequences of measurements are taken into account. This so called “hidden nonlocality” is revealed by a sequence of two measurements, where the first measurement is used to select a certain subensemble of pairs —those pairs which produce a specific outcome— while the second measurement tests the Bell observable on the subensemble. If the subensemble does not satisfy Bell’s inequality, then one concludes that the initial ensemble violates local realism. In fact, applying a similar reasoning to collective tests of particles , it was shown that all inseparable states of two qubits  violate local realism. More generally, all distillable  states violate local realism. It is an open problem whether violation of Bell’s inequality —which seems to be a rather strong requirement— already implies distillability .

While in systems of two qubits all entangled states are distillable  (and thus violate local realism), this turned out to be false for higher dimensional systems. The Horodecki  discovered states which, although being non–separable and thus entangled, are not distillable. Those states are called bound entangled. The role of bound entangled states is not entirely clear yet. Although they can not be useful for any quantum information processing task directly (since they are non–distillable), it was nevertheless shown that they allow to perform certain processes (e.g. quasidistillation)  which cannot be performed using local operations and classical communication alone. In particular, it is not known whether bound entangled states violate local realism or not. This paper is aimed to shine some light on these questions.

We will, however, tackle this problem not in its bipartite setting but rather in the multipartite setting. To this aim, we consider $`N`$ spatially separated parties. A generalization of the Bell-CHSH inequality to multipartite systems is due to Mermin , and was further developed e.g. in Ref.  (see also Refs.  for some recent developments). Also the notion of separability and distillability can be readily generalize to multipartite systems (see e.g. Refs. ). We have that a $`N`$–partite state $`\rho `$ is called (fully) separable iff it can be written as a convex combination of (unnormalized) product states, i.e.

$$\rho =\underset{i}{}|a_i_{\mathrm{party1}}a_i||b_i_{\mathrm{party2}}b_i|\mathrm{}|n_i_{\mathrm{party}N}n_i|.$$
(1)
If $`\rho `$ is not fully separable, it is entangled. Note that fully separable states are those which can be prepared locally. Due to the fact that many different kinds of multipartite pure state entanglement exist, there are various kinds of distillability . We will however not distinguish between those possible kinds of entanglement but rather say that a state $`\rho `$ is distillable, iff —by means of local operations and classical communication— out of an arbitrary number of identical copies of a state $`\rho `$ some entangled pure state can be created. If no entangled pure state whatsoever can be created, $`\rho `$ is non–distillable. As shown in Ref.  (see also ), there exist also multipartite bound entangled states, i.e. states which are not (fully) separable and hence entangled, but which are non–distillable. Here we investigate multipartite bound entangled states and consider the relation of their distillability to the existence of a LHV model. We show that there exist multipartite bound entangled states which violate Bell’s inequality. This means on the one hand that violation of Bell’s inequality does not imply distillability, and on the other hand that some bound entangled states violate local realism.

In the remainder of this letter, we consider specific $`N`$–qubit mixed states $`\rho _N`$ acting on a Hilbert space $`=(\mathrm{I}\mathrm{C}^2)^N`$. We consider $`N`$ parties, $`A_1,\mathrm{},A_N`$, at different locations, each of them possessing several qubits. The parties possess $`M`$ identical copies of $`\rho _N`$, where $`M`$ can be arbitrary large. Thus the state of all qubits is described by the density operator $`\rho _N^M`$. This ensures that the parties can use distillation protocols  in order to obtain entangled pure states between some of them. The states $`\rho _N`$ we consider are given by

$`\rho _N={\displaystyle \frac{1}{N+1}}\left(|\mathrm{\Psi }\mathrm{\Psi }|+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{N}{}}}(P_k+\overline{P}_k)\right).`$ (2)
We have that $`|\mathrm{\Psi }`$ is a $`N`$-party Greenberger-Horne-Zeilinger (GHZ) state ,

$$|\mathrm{\Psi }=\frac{1}{\sqrt{2}}(|0^N+e^{i\alpha _N}|1^N),$$
(3)
with $`\alpha _N`$ being an arbitrary phase. We denoted by $`P_k`$ a projector on the state $`|\varphi _k`$, where $`|\varphi _k`$ is a product state which is $`|1`$ for party $`A_k`$ and $`|0`$ for all other parties, i.e. $`|\varphi _k=|0_{A_1}|0_{A_2}\mathrm{}|1_{A_k}\mathrm{}|0_{A_{N1}}|0_{A_N}.`$ The projector $`\overline{P}_k`$ is obtained from $`P_k`$ by replacing all zeros by ones and vice versa.

For the states $`\rho _N`$ (2) we show

The states $`\rho _N`$ are bound entangled, i.e. non-separable and non-distillable if the number of parties $`N4`$.

The states $`\rho _N`$ violate the Mermin-Klyshko inequality if the number of parties $`N8`$ and can thus not be described by a LHV model.

We start out by showing $`(𝐢)`$, i.e. $`\rho _N`$ is bound entangled. One readily verifies that $`\rho _N^{T_{A_k}}`$ is a positive operator $`k`$, where $`T_{A_k}`$ denotes partial transposition with respect to party $`A_k`$ . This already implies that $`\rho _N`$ is non–distillable . To see this, assume on the opposite that one can distill some bi- or multipartite entangled pure state. As shown below (Lemma 1), one can always create by means of local operations from any multipartite entangled pure state a maximally entangled bipartite pure state shared between two of the parties, say $`A_i`$ and $`A_j`$. Thus the resulting state has non-positive partial transposition with respect to parties $`A_i`$ and $`A_j`$, while the initial state $`\rho _N`$ has positive partial transposition. Due to the fact that by means of local operations and classical communication, one cannot change the positivity of the partial transposition , we have the desired contradiction, hence $`\rho _N`$ is not distillable.

On the other hand, $`\rho _N`$ can easily seen to be entangled for $`N4`$, e.g. by observing that $`\rho _N^{T_{A_kA_l}}`$ is not a positive operator for $`kl`$ and $`N4`$. This already implies that $`\rho _N`$ is not fully separable, as positivity of all possible partial transpositions is a necessary condition for (full) separability in multipartite systems . It remains to show the announced

Lemma 1: Given a single copy of an $`N`$–party entangled pure state $`|\mathrm{\Phi }`$ of arbitrary dimension, one can always create with non–zero probability of success by means of local operations a maximally entangled bipartite pure state shared among some of the parties.
Proof: One may write $`|\mathrm{\Phi }`$ in its Schmidt decomposition with respect to any party $`A_k`$. Since $`|\mathrm{\Phi }`$ is entangled, there exist at least one party, say $`A_1`$, where one obtains a minimal number of two non-zero Schmidt coefficients (otherwise $`|\mathrm{\Phi }`$ would be a $`N`$-party product state). That is,

$$|\mathrm{\Phi }=\underset{k=0}{\overset{d}{}}\lambda _k|k_{A_1}|\phi _k_{A_2\mathrm{}A_N},$$
(4)
where $`k|k^{}=\phi _k|\phi _k^{}=\delta _{kk^{}}`$ and $`d2`$. For simplicity, let us assume that $`d=2`$ and $`\lambda _0=\lambda _1=1/\sqrt{2}`$ (this can always be accomplished by a filtering measurement in $`A_1`$, e.g. using $`O_{A_1}=\lambda _1|00|+\lambda _0|11|`$). Now, either (a) both $`|\phi _0`$ and $`|\phi _1`$ are product states, or (b) at least one of them, say $`|\phi _0`$, is entangled.

In case of (a), $`|\phi _0`$ and $`|\phi _1`$ have to be locally orthogonal in at least one location, say at $`A_2`$ (this is due to the fact that $`\phi _0|\phi _1=0`$), i.e. $`|\phi _0=|0_{A_2}|\chi _3_{A_3}\mathrm{}|\chi _N_{A_N}`$ and $`|\phi _1=|1_{A_2}|\stackrel{~}{\chi }_3_{A_3}\mathrm{}|\stackrel{~}{\chi }_N_{A_N}`$. There may exist $`l`$ locations $`A_k`$, $`lN2`$ for which $`|\chi _k=|\stackrel{~}{\chi }_k`$. Each of the other parties $`A_k`$ can apply a local filtering measurement of the form $`O_{A_k}=|0\chi _k^{}|+|1\stackrel{~}{\chi }_k^{}|`$, where $`\{|\chi _k^{},|\stackrel{~}{\chi }_k^{}\}`$ is the biorthonormal basis to $`\{|\chi _k,|\stackrel{~}{\chi }_k\}`$, i.e. $`\stackrel{~}{\chi }_k|\chi _k^{}=\chi _k|\stackrel{~}{\chi }_k^{}=0`$ and $`\stackrel{~}{\chi }_k|\stackrel{~}{\chi }_k^{}=\chi _k|\chi _k^{}=1`$. One readily observes that this leads to the creation of a $`Nl`$ party GHZ state (3), from which —by means of local measurements in the basis $`\{|+,|\}`$, where $`|\pm =1/\sqrt{2}(|0\pm |1)`$ at the remaining locations— a maximally entangled bipartite pure state shared between any two out of the remaining $`Nl`$ parties can be created deterministically.

In case of (b), one measures in $`A_1`$ the projector $`P_{A_1}=|00|`$ and is left with an entangled state of $`N1`$ (or less) particles. This situation is similar to the one we started with, however the number of entangled systems decreased. One proceeds in the same vain with the remaining systems until (a) applies, which happens in the worst case if only two entangled particles are left. Finally, one obtains at least a maximally entangled bipartite state shared among two of the parties which concludes the proof of the lemma.

Note that from Lemma 1 follows the non–existence of a LHV model for all multipartite entangled pure states of arbitrary dimension which describes properly also sequences of measurements (see Ref.  for the stronger result of inconsistency with local realism even for single measurements per site). Following the reasoning of Popescu  (see also ) (adopted to the multipartite case), the violation of Bell’s inequality of a subensembles of states (obtained e.g. by local filtering measurements) ensures that also the original ensemble violates local realism. Since, according to Lemma 1, from any multipartite entangled pure state a maximally entangled bipartite pure state can be created —which clearly maximally violates Bell’s inequality— the claim follows.

We now turn to (ii) and show that $`\rho _N`$ violates the Mermin-Klyshko inequality for $`N8`$. Let $`a_j,a_j^{}`$ be two vectors on the unit sphere which indicate two possible measurement directions for party $`A_j`$. The corresponding observables are given by $`O_j=\sigma _{a_j},O_j^{}=\sigma _{a_j^{}}`$. Up to a normalization factor, any $`k`$–qubit Bell inequality involving two observables per qubit can be written as

$$|_k|1,$$
(5)
where $`_k_k(a_1,a_2,\mathrm{},a_k,a_1^{},a_2^{},\mathrm{},a_k^{})`$ is the corresponding Bell operator. We consider the Mermin-Kyshko inequalities  for $`N`$ qubits, whose corresponding Bell operator is defined recursively as

$$_k=\frac{1}{2}_{k1}(\sigma _{a_k}+\sigma _{a_k^{}})+\frac{1}{2}_{}^{}{}_{k1}{}^{}(\sigma _{a_k}\sigma _{a_k^{}}),$$
(6)
where $`_{}^{}{}_{k}{}^{}`$ is obtained from $`_k`$ by exchanging all the $`a_k`$ and $`a_k^{}`$.

We choose the same measurement directions in all $`N`$ locations, $`\sigma _{a_j}=\sigma _x`$ and $`\sigma _{a_j^{}}=\sigma _yj`$, where $`\sigma _x,\sigma _y`$ are Pauli matrices. One readily verifies that in this case $`_N`$ can be written as

$$_N=2^{\frac{N1}{2}}(e^{i\beta _N}|1^N0^N|+e^{i\beta _N}|0^N1^N|),$$
(7)
with $`\beta _N\pi /4(N1)`$. Using that tr$`(_NP_k)`$ = tr$`(_N\overline{P}_k)=0k`$ and tr$`(_N|\mathrm{\Psi }\mathrm{\Psi }|)=2^{\frac{N1}{2}}`$ when fixing the phase $`\alpha _N=\beta _N`$ in $`|\mathrm{\Psi }`$ (3), we have

$$\mathrm{tr}(_N\rho _N)=\frac{1}{N+1}2^{\frac{N1}{2}},$$
(8)
which fulfills $`\mathrm{tr}(_N\rho _N)>1`$ iff $`N8`$. Thus the states $`\rho _N`$ (2) with the choice $`\alpha _N=\pi /4(N1)`$ violate the Mermin-Klyshko inequality for $`N8`$ as announced.

To conclude, we have shown that certain multipartite bound entangled states violate a multipartite Bell inequality. This implies that (i) violation of Bell’s inequality is not a sufficient condition for distillability and (ii) there does not exist a local hidden variable model for certain bound entangled states. Note that the states $`\rho _N`$ for sufficiently large $`N`$ violate the Mermin-Klyshko inequality directly, and no sequence of measurements, eventually performed on a tensor product of the states, is required to rule out the existence of a LHV model as in the case of “hidden nonlocality”. There remain a number of open problems concerning the relation of inseparability to the existance of a LHV model . In particular, it is not known whether all bipartite bound entangled states, those with positive partial transposition as well as the conjectured ones with non–positive partial transposition , can be described by a local hidden variable model or not.

I thank G. Vidal for an essential contribution in the proof of Lemma 1, and J. I. Cirac and G. Giedke for interesting discussions. This work was supported by the Austrian SF under the SFB “control and measurement of coherent quantum systems” (Project 11), the European Community under the TMR network ERB–FMRX–CT96–0087 and project EQUIP (contract IST-1999-11053), the ESF, and the Institute for Quantum Information GmbH.


format: text

text:
# A guided Monte Carlo method for optimization problems

\[

## Abstract

We introduce a new Monte Carlo method by incorporating a guided distribution function to the conventional Monte Carlo method. In this way, the efficiency of Monte Carlo methods is drastically improved. To further speed up the algorithm, we include two more ingredients into the algorithm. First, we freeze the sub-patterns that have high probability of appearance during the search for optimal solution, resulting in a reduction of the phase space of the problem, a concept inspired by the renormalization group equation method in statistical physics. Second, we perform the simulation at a temperature which is within the optimal temperature range of the optimization search in our algorithm. We use this algorithm to search for the optimal path of the traveling salesman problem and the ground state energy of the Anderson-Edwards spin glass model and demonstrate that its performance is comparable with more elaborate and heuristic methods.

\]

Optimization problems arise in areas of science, engineering and other fields. Many of them are NP hard problems, in the sense that the number of computing steps required to solve the problem increases faster than any power of the size of the system. The traveling salesman problem (TSP), the spin glass problem and the Lennard-Jones microcluster problem are some of the examples belonging to this class. Over the years, people have developed heuristic methods that can allow one to possibly obtain optimal solutions in some optimization problems, e.g. the TSP. However, these methods are usually specially designed for the problem that one is seeking to solve. In many cases, these heuristic methods will either not be applicable or will take a much longer time to locate the optimal solution of the problem. In many practical situations, it is more desirable to obtain a near-optimal solution quickly than to find the true optimal solution of the problem. Furthermore, it is also preferrable to have an algorithm that is simple to use and general enough to treat various optimization problems.

Since its introduction about half a century ago, the Monte Carlo (MC) method  has been widely used to treat optimization problems. Despite its simplicity and versatility, it is known that Monte Carlo method cannot usually give one an acceptable solution since it can easliy get trapped in local minima of the problem. In the past decades, people have developed new algorithms based on the original Monte Carlo method. Simulated Annealing(SA)  is one of the algorithms that originates from the Monte Carlo method. Because of its simplicity and versatility, SA has been employed to study varoius statistical systems with success. Despite its improvement over the original Monte Carlo method, SA also suffers from a similar problem, that the algorithm can easliy get trapped in local minima. Progress has been made to improve SA and Monte Carlo methods further by developling more efficient MC sampling algorithms. More recent developments include algorithms such as the multicanonical method  and simuulated tempering . All the above mentioned algorithms are meant to give a more efficient MC sampling over the original method.

In this paper, we introduce a new algorithm by incorporating the Monte Carlo method with a guided distribution function. The advantage of this new algorithm is that it is very simple to use and is generally suitable for most optimization problems, whether they involve continuous or discrete parameter spaces. The essentials of this algorithm is from both biology and theoretical physics. In the following, we will give the main ingredients of this algorithm and apply it to two problems, namely, the traveling salesman problem and the spin glass problem (SGP). The task of the TSP is to find a route through a given set of $`N`$ cities. The spin glass model we are testing here is the three-dimensional Edwards-Anderson Ising spin glass model  in which Ising spins ($`S_i=\pm 1`$) are placed on a three dimensional lattice and there are nearest neighbor interaction only. The corresponding Hamiltonian for this model is

$$H=\underset{<i,j>}{}J_{ij}\sigma _i\sigma _j,$$

where in our simulation, the $`J_{ij}`$ takes on values $`\pm J`$ randomly and $`\sigma _i,\sigma _j`$ can take values +1 or -1. The task here is to find the spin values which minimizes $`H`$. As it turns out, its performance in these two tests is comparable to the more elaborate methods and it can obtain near-optimal solutions rather quickly.

A new substance in this algorithm is the addition of a guided distribution function to the original Monte Carlo method. This introduction of the guided distribution function is motivated from biology. Similar to the subject of evolutionary programming , this guided distribution function is based on the cooperative effort of individuals. To begin with, we have a population of $`M`$ randomly generated solutions. We let each of them evolve independently and according to the conventional Monte Carlo rules, i.e. we accept a new solution over its old solution when the new solution is better fit than the old solution. If a new solution is less fit than the old solution, we pick the new solution over the old solution with respect to a Boltzmann factor. At this point, it is equivalent to $`M`$ individual runs of the Monte Carlo method. In each of the $`M`$ individual runs, we keep record of its best fit solution while the individual solution is evolving.

After a certain number of Monte Carlo steps, a distribution function of the $`M`$ best fit Monte Carlo solutions will be recorded. Let us use the TSP as an example. After we perform a preset number of Monte Carlo steps, each of the $`M`$ independent runs has its own best fit solution, or path. In each of these paths, we record the links between any two cities. There are a total of $`N`$ links in each path and a total of $`MN`$ links of the $`M`$ best fit solutions. One has now a distribution function of the probability of appearance among all possible links that connects any two cities. This completes our first layer of Monte Carlo simulation. In our discussion, a layer of simulation always means a set of $`M`$ individual Monte Carlo runs for a preset Monte Carlo steps plus the evaluation of the distribution function that we just mentioned. The distribution function that we obtain here will be used as a guided distribution function in the next layer of Monte Carlo simulation of $`M`$ individual runs.

In the next layer of Monte Carlo simulation, a set of $`M`$ individual Monte Carlo runs is again performed. An important point here is that we also start with $`M`$ randomly generated solutions at the beginning of the simulation in this layer. There is however one more criterion here—that we pick the links which appear less frequent in the guided distribution function and try to change them into links with a higher probability of appearance. In practice, we pick one of the two cities that are connected with a less probable link and connect it to a city where this new link appears more often in the guided distribution function. The new solution will be compared with the old one to decide whether we keep or disgard it according to the conventional Monte Carlo method. The idea behind this is simple. The more (less) often that a certain link appears among the $`M`$ best fit solutions, the more (less) likely that it would (not) appear in the optimal solution. One would then let a given solution evolve into solutions with links of high probability of appearance while avoiding links of low probability of appearance. The Monte Carlo method still allows enough fluctuation to search for optimal, or near-optimal solutions. After a preset number of Monte Carlo steps is performed, a new guided distribution function will again be obtained for the next layer simulation. In principle, the more Monte Carlo steps one performs in each layer, and the more layers of simulations one has, the better will be the result.

Two more ingredients can be added into this new algorithm. The first is the elimination of irrelevant degrees of freedom. Again taking the TSP as an example. There are a total of $`\frac{N(N1)}{2}`$ possible links among the $`N`$ cities. Some of them have a higher probability of appearing in the best fit solutions while others have very low probability. We can elminate some of the links by freezing a link between two cities that exceeds a certain probability of appearance in the guided distribution function. In practice, after the guided distribution function in a layer is obtained, one freezes links which have a high probability say, 0.8 of appearance before one performs the next layer simulation. In this way, one reduces the available degrees of freedom in the system in the next layer, thus speeds up the search process. This idea is inspired by the renormalization group equation method in statistical physics and is similar in spirit to a recent work  which has proposed to treat optimization problems using the idea of renormalization in statistical physics. By freezing high probability links, one can obtain near-optimal solutions within a much faster CPU time. One can unfreeze the frozen links during the simulation to avoid missing the optimal solution.

The last ingredient that we include in this algorithm is from an observation in a recent work . It was  observed that the average first passage time to reach the global optimum of an optimization problem in Monte Carlo simulation was a U-shape temperature dependent curve. At the optimal temperature, the avereage first passage time to find the global optimal will be the shortest. This means that for the optimization problem one has, it is best to perform simulation at the optimal temperature so that one can locate the optimal, or near-optimal solution in the fastest CPU time possible. The reason for the appearance of such an optimal temperature can be understood as follows. At a given temperature, the system tends to stay in states near the equilibrium value determined by that particular temperature rather than going to the global optimum. As the system relaxes toward this equilibrium value, it would fluctuate around this equilibrium value. The Gaussian like fluctuation would then carry the system from its equilibrium value to the global optimum as long as there are no insurmountable barriers between them. When the simulation is performed at high temperatures, there are too many paths available and it would take the system too long to locate the path that would reach the global optimum. As the temperature lowers, the number of available paths are reduced. When the temperature further decreases, the potential barriers and traps make the system more difficult to move out of the local optima. There is therefore an optimal temperature at which the two factors balance each other and the system approaches the global optimum with the shortest time.

(1) The Traveling Salesman Problem (TSP)

Since we have already used the TSP as an example to explain the implementation of our algorithm above, we do not need to repeat the essential steps here but instead only state the parameters we used and give some comments at the end. In the simulation, we used $`M=50`$ in all the TSP cases that we tested. We performed 150 Monte Carlo steps in each layer at a temperature that is near the optimal temperature in each of the cases which we give in Table I. We have tested five cases from the TSP databank : d198, lin318, pcb442, rat783 and fl1577. These cases are selected from the set of sample problem instances as listed in . We have performed 100 trial runs in each case and computed the average value of the best solutions obtained in each of the 100 trial runs. We also kept record of the best solution among the 100 trial runs and list them in Table I. We used 6 layers for the first three cases and 8 layers for the last two cases respectively. To make our discussion simple, we only consider 2-opt moves in our algorithm. Other moves can be considered but is not our main emphasis here.

Table I. Tests on 5 cases from TSPLIB. $`T`$ is the temperature used for the Monte Carlo simulation, $`RQ`$(=% above the optimal solution) is the relative quality of the average length of the best paths obtained for the 100 trial runs with respect to the optimal solutions given in TSPLIB, $`Best`$ is the shortest path length obtained in the 100 trial runs using the present algorithm and $`Op`$ is the optimal solution given in TSPLIB.

| Case | $`T`$ | $`RQ`$ | $`Best`$ | $`Op`$ |
| --- | --- | --- | --- | --- |
| d198 | 10.0 | 0.0155% | 15780 | 15780 |
| lin318 | 37.5 | 0.30% | 42029 | 42029 |
| pcb442 | 27.5 | 0.67% | 50798 | 50778 |
| rat783 | 4.0 | 0.68% | 8819 | 8806 |
| fl1577 | 2.75 | 0.88% | 22273 | 22249 |

As can be seen in Table I, the algorithm can locate the optimal paths in the two cases with fewest cities, i.e. d198 and lin318 while it can obtain near-optimal solutions for the other three cases within the number of Monte Carlo steps we preset. In general, the more the individual Monte Carlo runs, i.e. the larger the $`M`$, and the more the Monte Carlo steps and layers one uses, the better will be the solution obtained. For example, in the case of d198, 23 out of 100 trial runs found the optimal solution when we perform 150 Monte Carlo steps in each layer. When we increased the Monte Carlo steps in each layer to 200, the average value of the best runs would now be only 0.00906% above the optimal solution with 41 out of 100 trial runs found the optimal solution. As a comparison, we have performed a simulation on the d198 case using the conventional Monte Carlo method, for 100 trial runs with 45,000 Monte Carlo steps per run. The number of Monte Carlo steps here is equivalent to the number of $`M`$ (50 here) times the number of layers (6 in the case of d198) times the 150 Monte Carlo steps in each layer. The average value for the best solution of the 100 trial runs is 1.07% above the optimal value and the best solution is 15843. A further increase of Monte Carlo steps to $`10^5`$ improves the average value slightly to 0.90% above the optimal value. This shows the superiority of the present algorithm over the conventional Monte Carlo method.

It is perhaps worth mentioning that locating the optimal or near-optimal temperature should help one perform simulations in the case of random TSP in a two-dimensional unit square or something similar. As an example, we have obtained an average optimal temperature around 0.02 over 10 configurations of the random TSP with 100 cities, when 100 Monte Carlo steps are performed in each layer of a 6 layer simulation. The optimal temperature range is however affected by the number of Monte Carlo steps and the moves(we use 2opt moves here). A rescaling of the optimal temperature range with respect to the system size, i.e. the number of cities, is also needed.

(2) The Spin Glass Problem (SGP)

Spin glasses have been a subject in statistical physics that is under intensive study. The simplest spin glass model is the Edwards-Anderson model and its Hamiltonian has been given in the above. There are exact methods available  for such spin glass models in two-dimensions but are not effective in three-dimensions. Estimates of the ground state energy are available for the three dimensional case  in the literature for lattice sizes up to $`L=12`$ and could be used for comparison. Our task here is to use the present algorithm to estimate the ground state energy of the three-dimensional Anderson-Edwards model within a reasonable amount of CPU time. Employing our algorithm to this model is straightforward. To begin with, we have a set of $`M`$ randomly generated spin configuration for an $`L\times L\times L`$ lattice with randomly assigned $`\pm J`$ couplings. We again let each of the configurations evolve for a preset Monte Carlo steps and record the $`M`$ best solutions from the $`M`$ individual runs. These $`M`$ best solutions are then used to construct the guided distribution function. Similar to the TSP, we freeze a sub-spin configuration pattern on this three-dimensional lattice if it exceeds a certain preset probability. In this way, we reduce the spin degrees of freedom after each layer of simulation. We again randomly generate $`M`$ spin configurations at the beginning of each new layer. There is one thing one needs to keep in mind in constructing the guided distribution function of the spin glass problem. When all the spins on the three-dimensional lattice changes sign, the energy of the system does not change. Therefore, one needs to fix one of the spins to take on the same value for all the $`M`$ best solutions when one constructs the guided distribution function. For example, we can set the spin value of the (0,0,0) site of the $`L\times L\times L`$ lattice to be 1. If a particular best solution has a value -1 for this site, then all the spins have to flip sign in this best solution before one records its spin configuration.

In our simulation, we have used $`M=15`$. We have carried out simulations for $`L`$ ranging from 4 to 12 and the result is given in Table II. Since our purpose here is to obtain the ground state energy of the model, we will instead perform several runs for each spin configurations and pick the lowest energy as its ground state. The number of layers and the number of Monte Carlo steps in each layer that we used will also depend on $`L`$. For example, for $`L=6`$, we used 2 layers with 150 Monte Carlo steps in each layer and 5 trial runs for each randomly assigned coupling configuration we generated. In the case of $`L=8`$ sites, we used 6 layers with 200 Monte Carlo steps in each layer and 10 trial runs for each randomly assigned coupling configuration before we take the best solution as the ground state energy for that particular spin configuration. The number of Monte Carlo steps and layers for different $`L`$ are given in Table II.

In the simulation of the spin glass model, we have observed that unlike the case of the TSP, the temperature range in the Boltzmann factor that is most efficient to locate the ground state or near ground state energy is between $`T=J`$ and $`T=1.1J`$ for all the $`L`$ that we tested and we have used $`T=1.05J`$ for all our simulation of the spin glass model. That the optimal temperature range is not sensitive to the lattice size here is probably related to the fact that one only considers the nearest neighbor interaction in the model.

Table II. Tests on $`L\times L\times L`$ SG instances. $`m`$ is the number of randomly assigned $`\pm J`$ ($`J=1`$ here) cases, $`MC`$ is the number of Monte Carlo steps in each layer and $`E`$ is the average ground state energy of the $`M`$ spin configurations. We have set $`J=1`$ in our simulation. We also include here the results from Ref for comparison.

| $`L`$ | $`m`$ | $`MC`$ | $`E`$ | Ref | Ref |
| --- | --- | --- | --- | --- | --- |
| 4 | 200000 | 50 | -1.73750(13) | -1.73749(8) | -1.7370(9) |
| 6 | 50000 | 150 | -1.77128(12) | -1.77130(12) | -1.7723(7) |
| 8 | 4000 | 200 | -1.77910(28) | -1.77991(22) | -1.7802(5) |
| 10 | 400 | 400 | -1.7831(8) | -1.78339(27) | -1.7840(4) |
| 12 | 40 | 600 | -1.7847(14) | -1.78407(121) | -1.7851(4) |

As shown in Table II, the ground state energy that we obtained by using our algorithm is comparable with the present estimate of the ground state energy of the model by using other methods. If more layers, Monte Carlo steps and trial runs are used, the result will be improved.

In this paper, we have introduced a new algorithm by incorporating a guided distribution function into the conventional Monte Carlo method. This new algorithm is both biologically and physically motivated. By evolving a population of $`M`$ solutions and record their best individual solutions in each layer, we are relying on the cooperative effort of the population, similar to many evolutionary algorithms such as genetic algorithms. Reducing the computational phase space by eliminating irrelevant degrees of freedom of the problem as we continue our layer by layer simulation is similar to the renormalization group equation method in statistical physics. Together with the use of an optimal temperature in Monte Carlo simulation for optimization problems, which results in a speed up of optimum search processes, it proves to be a powerful new algorithm. It takes about 24 seconds and 9 seconds of a Pentium II 450 MHz CPU time to perform a trial run of the d198 of the TSP and $`L=8`$ of the SGP. Tests of our algorithm in the TSP and SGP are promising. Incorporating heuristic methods into this algorithm will certainly improve the efficiency and power of its optimum search. This algorithm is very general and can be applied to many optimization problems including discrete optimization problems we tested here and optimization problems with continuous parameters. The present algorithm has been applied to the Lennard-Jones microcluster problem and can reproduce all the known ground state energy for the system sizes that the algorithm has so far tested . We believe that this algorithm should be a simple and general enough algorithm that can obtain optimal or near-optimal solutions within a reasonable CPU time for optimization problems.


text:
# Breakdown of Self-Organized Criticality

## Abstract

We introduce two sandpile models which show the same behavior of real sandpiles, that is, an almost self-organized critical behavior for small systems and the dominance of large avalanches as the system size increases. The systems become fully self-organized critical as the system parameters are changed, showing that these systems can make a bridge between the well known theoretical and numerical results and what is observed in real experiments. A simple mechanism determines the boundary where self-organized criticality can or cannot be found, which is the existence of local chaos.

The concept of Self-Organized Criticality (SOC) was introduced by Bak, Tang and Wiesenfeld (BTW) in 1987 to denote a phenomenon in which out of equilibrium, multidimensional systems, drive themselves to a critical state characterized by a power-law distribution of event sizes. Until then, the studies of fractal structures were related to equilibrium systems where this kind of distribution appears only at special parameter values where a phase transition takes place. In that pioneering work, the concept of SOC was illustrated by a model for sandpiles and since then an enormous amount of research on SOC, both theoretically and experimentally, has been done. Among other phenomena in which SOC has been connected with are earthquakes and evolution.

The existence of SOC in an experiment with a quasi-one-dimensional pile of rice was demonstrated by Frette et. al. They found that the occurrence of SOC depends on the shape of the rice. Only with sufficient elongated grains avalanches with a power-law distribution occurred. For more symmetric grains a stretched exponential distribution was seen. Christensen et. al introduced a model for the elongated rice pile experiment in which the local critical slope varies randomly between 1 and 2. They found that their model, known as the Oslo rice pile model, reproduced well the experimental results on the quasi-one-dimensional rice pile. In a recent publication, we introduced a fully deterministic one-dimensional SOC system, which presents the same qualitative and quantitative behavior of the Oslo system. In other words, they belong to the same universality class.

When one goes to sandpiles with geometry of two-dimensions a different picture emerges. That is, the models predict the presence of power-law distributions and the experiments do not display them. The most well known sandpile experiments can be classified in two types: (a) local dropping of sand in the center of the pile and (b) uniform driving, more specifically the rotating drum experiments. In class (a) it was found that small systems present scaling almost consistent with SOC, but in large systems another regime with big avalanches belonging to a different distribution appears. In class (b) one sees small avalanches that seems to display a power-law distribution of limited size. These small avalanches are interwoven with big system-wide avalanches which belong to a different distribution. There are studies of avalanches in natural settings, which also show similar distribution as the ones observed in the sandpile experiments.

Although SOC has not been observed in sandpile experiments is well known that power law distributions does exist in Nature, one of the most well known cases being Gutenberg-Richter law for earthquakes in agreement with SOC. It is clear that there is a missing piece in this puzzle. That is, would there be simple models that would display the observations in real sandpiles, and still present SOC in other parameter regions? To our knowledge it still missing in the literature such models and it is the aim of this letter to introduce them. Our models makes a bridge between the SOC and non SOC behavior, and the boundary that characterizes the separation between the two regimes is the absence or presence of local chaos. However, we will see that local chaos is a necessary but not sufficient condition for SOC to be destroyed.

Our models are inspired by the model introduced in Ref., and is also governed by coupled lattice map (that is, systems characterized by discrete time and continuous values for the space variables). Here we increase the dimensionality and change the drive and the relaxation rules. Having in mind the sandpiles experiments we first introduce a model for the local dropping of sand. We assume a two-dimensional square lattice of linear size $`L`$ and to each site $`i,j`$ in the lattice there is associated to it a variable $`x_{i,j}`$ with $`x[0,+\mathrm{})`$, which is to represent the local slope of the pile. The dynamics of the model is described by the following algorithm:

(1) Start the system by assigning random initial values for the variables $`x_{i,j}`$, so the they are below a chosen, fixed, threshold $`x_{th}`$.

(2) Choose a nearly central site of the lattice and update it slope according to $`x_{i,j}=x_{th}`$.

(3) Check the slope in each element. If an element $`i,j`$ has $`x_{i,j}x_{th}`$, update $`x_{i,j}`$ according to $`x_{i,j}^{}=\varphi (x_{i,j}x_{th})`$, where $`\varphi `$ is a given nonlinear function that has two parameters $`a`$ and $`d`$. Increase the slope in all its nearest neighboring element according to $`x_{nn}^{}=x_{nn}+\mathrm{\Delta }x/4`$, where $`\mathrm{\Delta }x=x_{i,j}x_{i,j}^{}`$ and $`nn`$ denotes nearest neighbors.

(4) If $`x_{i,j}^{}<x_{th}`$ for all the elements, go to step (2) (the event, or avalanche, has finished). Otherwise, go to step (3) (the event is still evolving).

Without losing generality, we can take $`x_{th}=1`$. In our simulation in step (2) we have chosen the site with $`i=j=L/2`$. The nonlinear function we use is

$$\varphi (x)=\{\begin{array}{cc}1dax,\mathrm{if}x<(1d)/a,\hfill & \\ 0,\mathrm{otherwise}.\hfill & \end{array}$$
(1)
The parameter $`d`$ would be associated with the minimum drop in energy after an event involving one single element and $`a`$ would be the parameter associated with the amount of friction between the grains. That is, the smaller the $`a`$, the larger the friction and the smaller the change in the slope of the pile. We have tested several other functions and found that the quantitative and qualitative results we show here are robust. The important ingredient being the shape of $`\varphi (x)`$ in the vicinity of $`x=0`$. In contrast with the one-dimensional case, it is not required here that $`\varphi (x)`$ be periodic in order to find the presence SOC.

We have chosen to evolve the system using parallel dynamics with open boundary conditions. It is beyond the scope of the present letter to study the several possible variations of our models. Further results on these models will be presented in a future publication. The distribution of time duration of the avalanches will also be presented in the future, but our preliminary results show that they are qualitatively similar to the ones for the size distribution.

We display an example of our simulations in Fig. 1, where we show the distribution of events $`P(s)`$ involving $`s`$ update steps, that is, the size of the avalanche. The events that involve all the elements of the system have been excluded from our analysis. However, in this model of local dropping of sand we have observed that nothing very distinct will occur if they are also included in the statistics. In (a) we show $`P(s)`$ for $`L=64`$, $`d=0.1`$ and vary $`a`$ and in (b), we use $`a=1.5`$ and vary $`d`$, keeping $`L=64`$. We notice the existence of two regimes. For small $`d`$, power-law distributions, that is SOC, appears only when $`a1`$, whereas if $`d`$ is large, we see SOC even with $`a>1`$. As in real sandpiles, the almost SOC regime is characterized by a region with an apparent scaling for small events, and the big events belong to a different distribution. To illustrate this, we show in Fig. 1(c) simulations with $`a=1.5`$ and $`d=0.1`$ and varying $`L`$. In small systems $`P(s)`$ can be fit to a scaling form of the type

$$P(s,L)=L^\beta G(s/L^\nu ),$$
(2)
as shown in Fig. 1(d), where we have used $`\beta =6`$ and $`\nu =3`$. The function $`G`$ is not well fit by a power-law, since one can clearly see in the figure that it is curved. We have found that it is consistent with a stretched exponential, as in real sandpiles. The observations of Himalayan avalanches also display the kind of distribution shown in Fig. 1(c) for large systems.

The second model we introduce here is for the rotating drum experiment, which we call uniform drive, since the slope of the pile increases uniformly for all the grains. The algorithm is similar to the one described above, with the exception of step (2), which is now replaced by

(2) Find the element in the lattice that has the largest $`x`$, denoted here by $`x_{max}`$. Then update all the lattice elements according to $`x_{i,j}=x_{i,j}+x_{th}x_{max}`$.

We show examples of $`P(s)`$ for this model in Fig. 2. There, in (a) we fix $`d`$ and vary $`a`$, and in (b) we fix $`a`$ and vary $`d`$. In both cases we have used $`L=64`$ and the events that involve all the elements of the system have been excluded. Distinctly from the model of local dropping, it seems here that there is a power-law distribution for any parameter value. However, the behavior is not exactly SOC. We have found that SOC is only seen if $`a1`$ or $`d`$ is large enough, as in the case of the local dropping. When $`a>1`$ and $`d`$ is smaller than a given value, we see a SOC-like behavior only for small values of $`L`$. As $`L`$ grows, there is a transition to a different behavior, in which the larger the system, the smaller the power-law region, as Fig. 2(c) shows.

System-wide avalanches have been reported in the rotating drum experiment that belong to a different distribution than the one of the small avalanches. This is exactly what we see in this model for $`a>1`$ and small $`d`$. In Fig. 3(a) we show all the events of the system including the ones that involve all the elements (for $`a=3`$, $`d=0.1`$ and varying $`L`$). We see that the scaling region does not get bigger as the system size increases, and a peak related to the events involving all the elements is seen. We have found that the intermediate size events ($`10s100`$) can be fit by a scaling of the type $`P(s,L)=L^\beta G(s/L^\nu )`$ with $`\beta =1`$ and $`\nu =0`$, as shown in Fig. 3(b). As in the rotating drum experiment the function $`G`$ in this case is closer to a power-law than in the case of the experiment with local dropping of sand.

For given $`a`$ and $`d`$ the slope of the power-law distribution seems to be the same for both models, that is, local dropping of sand and uniform drive, but the slope varies with $`a`$, as shown in Fig. 1 and 2. Consequently, the universality class of these models vary with the parameters. In Fig. 4(a) we show the scaling given by Eq. (2) with $`\beta =3.55`$ and $`\nu =2.70`$ for the local dropping (ld) and $`\beta =3.55`$ $`\nu =2.85`$ for the uniform drive (ud). In the limit $`a0`$ or $`d1`$ the slope of the power-law distribution is the same as the one in the BTW model and the conservative OFC model, that is, $`P(s)s^{1.245}`$.

We next investigate what would happen if the relaxation function is just a random number generator. In other words, instead of using $`\varphi `$ in step 3 of the above algorithm we now use $`x_{i,j}^{}=\rho `$, where $`\rho `$ is a random number uniformly distributed in the interval $`[0,1d]`$. We have found SOC for any $`d(0,1]`$ with the same exponents as the BTW model. This is displayed in Fig. 4(b) where we show the size distribution for $`d=0.01`$ in the cases of local dropping (ld) and uniform drive (ld). Therefore nonlinearities in $`\varphi `$ and consequently non ergodicity are necessary for the SOC behavior to be destroyed in these models. In that figure, we also show the case in which $`x_{i,j}^{}=0`$, which corresponds to $`d=1`$. In this limit we recover the OFC model for the case of uniform drive.

The reason why $`a=1`$ determines a special boundary, in which SOC may or may not be present, is due to the fact that it marks the boundary in which local chaos exists. By “local” here we mean that is is at the grain level, no matter what happens at the system level. Chaos is defined as the exponential divergence of trajectories that start with almost the same initial conditions. In fact, if we consider two copies of a given system, copy 1 and copy 2, with all the elements having the same $`x`$ except that in copy 1 the element $`x_{i,j}^1=x^{}`$, whereas in copy 2 $`x_{i,j}^2=x^{}+\mathrm{\Delta }`$. It is not difficult to find that after one iteration by $`\varphi `$, the separation of the two elements instead of $`\mathrm{\Delta }`$ will be $`a\mathrm{\Delta }`$. So, if $`a>1`$ the separation increases (that is, we have local chaos) whereas if $`a<1`$ the separation decreases. We have seen from the above results that, since SOC can exist even with $`a>1`$, if $`d`$ is large enough or $`L`$ is smaller than a given value, this means that local chaos is a necessary but not sufficient condition for SOC to be broken. Consequently, the true boundary between SOC and non SOC depends on three parameters, that is, $`a`$, $`d`$ and $`L`$. Further studies are necessary to understand this boundary more clearly.

In summary, we have introduced a model for sandpiles, or other systems that present avalanche like behavior, that reproduce very well the observation in real sandpiles. We have found that a simple mechanism, i.e local chaos, can explain the breakdown of SOC in those systems. Based on our results, we believe that with an appropriate choice of grains SOC will be seen even in real sandpiles. Grains with large friction would be the best candidates for this.


text:
# Five Years of Magellanic Clouds Research: A Newsletter Editors’ Perspective

## 1 Introduction

In late 1995 the Magellanic Clouds Newsletter (MCNews) was founded by You-Hua Chu’s Magellanic Clouds Working Group at the University of Illinois in Urbana-Champaign (UIUC). Two years earlier a joint “Graduiertenkolleg” (graduate school) for Magellanic Clouds research had been created by the Universities of Bonn and Bochum. Several members of the Graduiertenkolleg went to UIUC as exchange visitors, and the University of Bonn became the European mirror site for MCNews. The Graduiertenkolleg is now reaching the end of its funding period, and we use this opportunity to analyze the topical and demographical evolution of Magellanic Clouds research worldwide over the past five years as reflected in MCNews.

## 2 The Magellanic Clouds Newsletter

MCNews covers all areas of Magellanic Clouds research and publishes abstracts of submitted and accepted refereed papers, PhD theses, conference proceedings, and job and conference announcements. Its editors are Eva Grebel (MPIA) and You-Hua Chu (UIUC). Since September 1997 MCNews appears monthly and is currently sent out electronically in format to $`440`$ subscribers in 31 countries. Fifty five issues have appeared as of May 2001, comprising a total of 532 abstracts of refereed papers (an estimated 80% of the refereed publications in this area), so we estimate that we are reaching $`80`$ % of the researchers active in this field.

## 3 Demographics of Magellanic Clouds research

Assuming that MCNews subscribers are approximately representative of Magellanic Clouds researchers, then Europe (42%) and North America (37%) have the highest concentration of Magellanic Clouds researchers. For the top eight countries the following fractional distribution of subscribers results: USA: 32%, Germany: 16%, Australia: 8%, France: 7%, UK: 6%, Italy: 5%, Chile: 4%, Canada: 3%. Between 8% (Asia) and 24% (South America) are women (i.e., 18% on average).

## 4 Scientific productivity

The number of refereed papers corresponds roughly to the number of Magellanic Clouds researchers per country. The largest total number of publications is found in the US, followed by Germany, France, and the UK. We use the affiliation of the first author to assign a paper to a specific country without considering the rest of the author team. Some of the annual variations are directly correlated with special campaigns or instruments. For example, the pronounced increase in Japanese refereed papers on the Magellanic Clouds is based largely on data from the X-ray satellite ASCA. The activity peaks in Poland are almost exclusively due to the OGLE microlensing experiment.

Dividing the total number of papers per country by the number of subscribers per country gives an approximate measure of activity per researcher in the said country. Countries with few subscribers but very active research groups stand out in this normalization. Poland takes a pronounced lead due to its very productive OGLE group and small number of subscribers. Other countries where researchers have a per capita rate of refereed papers that clearly exceeds 1 are Brazil (especially star clusters), Israel (particularly theory), Japan (X-ray binaries), and the Netherlands. Countries with a large total number of refereed papers (e.g., USA, Germany) are close to one paper per subscriber, illustrating that not all of their subscribers actively work on the Magellanic Clouds. Further, not every person who submits an abstract to MCNews is a subscriber, and the number of papers per subscriber ranges from 0 to $`>20`$.

## 5 Main areas of Magellanic Clouds research

Stellar populations – individual stars such as massive stars, AGB stars, variables, binaries; field populations, clusters and associations – account for the largest number of refereed publications on the Magellanic Clouds. The fractions of the main research areas are as follows: Stellar populations: 45%, interstellar medium: 27%, microlensing: 12%, theory: 7%, distance: 5%, dynamics/kinematics: 3%. Many studies belong to several of these research areas, each of which was credited in these cases. Interestingly, theoretical studies account only for a small fraction of the total. Theoretical research concentrates mainly on evolutionary models, both for individual stars as well as for stellar synthesis and chemical evolution. Many Magellanic Clouds studies are based partially or wholly on optical data (fractional distribution: Optical: 53%, X-ray: 16%, UV: 12%, radio/sub-mm: 11%, infrared: 8%).

## 6 Milestones of Magellanic Clouds Research

The past five years have seen significant progress in many areas of Magellanic Clouds research, with large-scale surveys playing an important role. To highlight just a few, we recall the impact of the Australian H i synthesis maps of the Clouds, which revealed a complex, fractal ISM full of shells, holes, and fragments. The HIPASS multi-beam survey detected the leading arm of the Magellanic stream. The NANTEN CO survey showed the distribution of molecular clouds at unprecedentedly high resolution and correlated their location with the age of star-forming regions. A growing body of high-resolution infrared data (e.g., NICMOS, VLT) is resolving pillars, Bok globules, and pre-main sequence stars in the Clouds. High-resolution studies of the stellar IMF in clusters (mostly with HST) show evidence for Salpeter-like slopes, but also clear indications of mass segregation. Infrared and optical imaging surveys (e.g., DENIS, 2MASS, the Magellanic Clouds Photometric Survey) unveiled the structure of the Magellanic Clouds as traced by their various stellar components and led to comprehensive point source catalogs of special stellar types such as AGB stars and carbon stars. Searches for carbon stars at the periphery of the Clouds, and kinematic studies resulted in extended rotation curves and the identification of kinematic subcomponents. Many studies concentrated on red clump stars as distance and evolutionary indicators. The microlensing surveys (esp. EROS, MACHO, OGLE) yielded a wealth of information on variable stars in the Clouds and contributed to the still unsolved question of the distance to the Clouds. Eclipsing binaries are emerging as one way to solve this question independent of stellar evolutionary assumptions. X-ray surveys (e.g., ROSAT, ASCA) are helping to complete the census of supernova remnants, X-ray binaries, as well as revealing the large-scale distribution of hot gas. Deep HST observations (WFPC2) made possible accurate relative age determinations for the oldest globular clusters in the Magellanic Clouds, and the derivation of detailed star formation histories for their field populations. OGLE and other imaging surveys led to improved star cluster catalogs and cluster age-dating. Abundances of field stars and clusters were derived photometrically (Strömgren, Washington) and through increased spectroscopic samples (Ca ii triplet), revealing a complex age-metallicity relation. High-resolution stellar abundance measurements are also becoming more available; an area that will likely expand with the advent of the new large southern telescopes. The IUE archive, GHRS, and STIS advanced our knowledge of the UV and wind properties of massive stars. GHRS, STIS, ORFEUS, and FUSE increased our knowledge of interstellar abundances in the Clouds and led to the detection of H2. These studies are complemented by radio and sub-millimeter measurements in the Clouds and in the Magellanic bridge. Surprisingly, stellar and gaseous abundances in the bridge turn out to be lower than those of the young field populations in both the LMC and SMC.

We look forward to the exciting results that the coming years and the new large telescopes and satellite missions will bring.


text:
# Security of quantum key distribution using 𝑑-level systems

## Abstract

We consider two quantum cryptographic schemes relying on encoding the key into qudits, i.e. quantum states in a $`d`$-dimensional Hilbert space. The first cryptosystem uses two mutually unbiased bases (thereby extending the BB84 scheme), while the second exploits all the $`d+1`$ available such bases (extending the six-state protocol for qubits). We derive the information gained by a potential eavesdropper applying a cloning-based individual attack, along with an upper bound on the error rate that ensures unconditional security against coherent attacks.

PACS numbers: 03.67.Dd, 03.67.Hk, 89.70.+c

Quantum key distribution is probably one of the most promising concepts in quantum information theory, and has been extensively studied both theoretically and experimentally since its discovery by Bennett and Brassard in 1984. This cryptographic method allows two remote parties to share a secret key by use of a quantum channel supplemented with a public authenticated classical channel (see e.g.  for a review). The impossibility for an eavesdropper to tap the quantum channel without disturbing the communicated quantum data – in a way that can, in principle, be detected using the classical channel – ensures the security of the key distribution. Most of the research effort to date has focused on quantum cryptosystems based on two-dimensional quantum variables (qubits) carried e.g. by the polarization state of individual photons. In particular, the optimal individual attack is now known both for the BB84 protocol (using two mutually unbiased bases) and for the six-state protocol (using all three maximally unbiased bases). Strong bounds have also been derived in the more general case of coherent attacks, which are useful to assess the security of quantum cryptography (see e.g. ). For higher-dimensional systems, however, very few results have been obtained on the resistance to eavesdropping of qudit-based schemes (i.e., schemes based on encoding the key on $`d`$-level systems). The only schemes that have been considered use either two bases for a ququat (4-level system)  or four bases for a qutrit, but their security was only investigated against simple non-optimal attacks.

In this Letter, we investigate more general quantum cryptosystems where the encoding is made into qudits with arbitrary $`d`$, extending on an earlier study by some of us that only considered simple individual attacks. A first protocol we study consists in using two mutually unbiased bases, just as in the original BB84 scheme. The sender Alice sends a basis state in one of these two bases chosen at random, while the receiver Bob makes a measurement in one of these two bases, again at random. The basis used by each party is subsequently disclosed on the public channel, so that Alice and Bob obtain correlated $`d`$-ary random variables if they used the same bases (and if there was no disturbance on the channel), which happens with probability 1/2. The use of mutually unbiased (or complementary) bases implies that if Alice and Bob use different bases, Bob’s measurement yields a random number that is uncorrelated with Alice’s state. The raw secret key is then made out of the correlated data (discarding the uncorrelated data is known as the sifting procedure). This procedures ensures that any attempt by an eavesdropper Eve (oblivious of the chosen basis) to gain information on Alice’s state induces errors in the transmission, which can then be detected by the legitimate parties. The second qudit-based protocol that we study makes use of all the $`d+1`$ mutually unbiased bases that are available in a $`d`$-dimensional Hilbert space, much in the same way as the six-state protocol for qubits. Here, Alice and Bob choose their basis at random among the $`d+1`$ possible bases. This method clearly has a lower yield than the first one since the sifting procedure only keeps one transmission out of $`d+1`$ (instead of 1/2). However, as we shall see, this second protocol is more secure against individual attacks in the sense that a slightly higher error rate is acceptable.

In what follows, we will analyze the security of these two cryptographic protocols against individual attacks (where the qudits are monitored separately) as well as coherent attacks (where several qudits are monitored jointly). For the individual case, we consider a fairly general class of eavesdropping attacks that are based on quantum cloning machines. It is known for qubits that such a cloning-based attack simply is the optimal eavesdropping strategy, that is, the best Eve can do is to clone (imperfectly) Alice’s qubit and keep a copy while sending the original to Bob. An appropriate measurement of the clone (and the ancilla system) after disclosure of the basis enables Eve to gain the maximum possible information on Alice’s key bit. Extending this cloning-based individual attack to higher dimensions results in a lower bound on the information accessible to Eve for a given error rate. Hence, this yields an upper bound on the error rate, which is a necessary condition for security against individual attacks. (Higher error rates do not permit to establish a secret key using one-way communication.) We conjecture that applying the optimal cloner is actually the best strategy for Eve in any dimension, so that this bound is actually tight. For the case of coherent attacks, we consider a situation where Eve interacts with a qudit sequence of arbitrary (but finite) length, and then uses the basis information to extract key information. In particular, we make use of an information-theoretic uncertainty principle to derive a lower bound on Bob’s information, or, equivalently, an upper bound on the error rate. It is a sufficient condition for the protocol to be guaranteed to generate a nonzero net key rate for all attacks.

Let us consider first an individual eavesdropping based on the use of a quantum cloning machine for qudits. The particular cloning machine that is best using depends on whether the protocol uses two bases or $`d+1`$ bases. We focus on the case of two bases first, for which we need to use a cloner that copies equally well two mutually unbiased bases, e.g. the computational basis $`\{|k\}`$, with $`k=0,\mathrm{},d1`$, and its dual under a Fourier transform

$$\overline{|l}=\frac{1}{\sqrt{d}}\underset{k=0}{\overset{d1}{}}e^{2\pi i(kl/d)}|k$$
(1)
with $`l=0,\mathrm{},d1`$. We use a general class of cloning transformations as defined in . If Alice sends the state $`|\psi `$, the transformation reads

$$|\psi _A\underset{m,n=0}{\overset{d1}{}}a_{m,n}U_{m,n}|\psi _B|B_{m,n}_{E,E^{}}$$
(2)
where $`A`$, $`B`$, $`E`$, and $`E^{}`$ stand for Alice’s qudit, Bob’s clone, Eve’s clone and cloning machine, respectively. Here, the amplitudes $`a_{m,n}`$ (with $`_{m,n=0}^{d1}|a_{m,n}|^2=1`$) characterize the cloner, while the states $`|B_{m,n}_{EE^{}}`$ are $`d`$-dimensional Bell states, that is, a set of $`d^2`$ orthonormal maximally-entangled states of two qudits,

$$|B_{m,n}_{EE^{}}=\frac{1}{\sqrt{d}}\underset{k=0}{\overset{N1}{}}\mathrm{e}^{2\pi i(kn/d)}|k_E|k+m_E^{}$$
(3)
with $`m,n=0,\mathrm{},d1`$. Note that the kets must be taken modulo $`d`$ here. The operators $`U_{m,n}`$ defined as

$$U_{m,n}=\underset{k=0}{\overset{d1}{}}\mathrm{e}^{2\pi i(kn/d)}|k+mk|$$
(4)
form a group of qudit error operators, generalizing the Pauli matrices for qubits: $`m`$ labels the shift errors (extending the bit flip $`\sigma _x`$) while $`n`$ labels the phase errors (extending the phase flip $`\sigma _z`$). Tracing the output joint state given by Eq. (2) over $`EE^{}`$ implies that Alice’s state $`|\psi _A`$ is transformed, at Bob’s station, into the mixture

$$\rho _B=\underset{m,n=0}{\overset{d1}{}}|a_{m,n}|^2U_{m,n}|\psi \psi |U_{m,n}^{}$$
(5)
Thus, the state undergoes a $`U_{m,n}`$ error with probability $`|a_{m,n}|^2`$. Note that $`U_{0,0}=𝟙`$, implying that the state is left unchanged with probability $`|a_{0,0}|^2`$. If Alice sends any state $`|k`$ in the computational basis, the phase errors $`(n0)`$ clearly do not play any role in the above mixture since $`U_{m,n}|k=e^{2\pi i(kn/d)}|k+m`$, so Bob’s fidelity can be expressed as

$$F=k|\rho _B|k=\underset{n=0}{\overset{d1}{}}|a_{0,n}|^2$$
(6)
In the complementary basis, we have $`U_{m,n}\overline{|l}=e^{2\pi i(l+n)m/d}\overline{|l+n}`$, so the shift errors ($`m0`$) do not play any role and Bob’s fidelity becomes

$$\overline{F}=\overline{l|}\rho _B\overline{|l}=\underset{m=0}{\overset{d1}{}}|a_{m,0}|^2$$
(7)
For the cloner to copy equally well the states of both bases, we choose a $`d\times d`$ amplitude matrix of the form

$$a=\left(\begin{array}{cccc}v& x& \mathrm{}& x\\ x& y& \mathrm{}& y\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ x& y& \mathrm{}& y\end{array}\right)$$
(8)
with $`x`$, $`y`$, and $`v`$ being real variables satisfying the normalization condition $`v^2+2(d1)x^2+(d1)^2y^2=1`$. Thus, Bob’s fidelity is $`F=v^2+(d1)x^2`$ in both bases, and the corresponding mutual information between Alice and Bob (if the latter measures his clone in the good basis) is given by

$$I_{AB}=\mathrm{log}_2d+F\mathrm{log}_2F+(1F)\mathrm{log}_2\left(\frac{1F}{d1}\right)$$
(9)
since the $`d1`$ possible errors are equiprobable.

Now, the clone kept by Eve can be shown to be in a state given by an expression similar to Eq. (5) but with the amplitudes $`a_{m,n}`$ replaced by

$$b_{m,n}=\frac{1}{d}\underset{m^{},n^{}=0}{\overset{d1}{}}e^{2\pi i(nm^{}mn^{})/d}a_{m^{},n^{}}$$
(10)
that is, the Fourier transform of the $`a_{m,n}`$’s . This corresponds to a matrix similar to Eq. (8) but with

$`x`$ $``$ $`x^{}=[v+(d2)x+(1d)y]/d`$ (11)
$`y`$ $``$ $`y^{}=(v2x+y)/d`$ (12)
$`v`$ $``$ $`v^{}=[v+2(d1)x+(d1)^2y]/d`$ (13)
resulting in a cloning fidelity for Eve given by $`F_E=v^2+(d1)x^2`$. Maximizing Eve’s fidelity $`F_E`$ for a given value of Bob’s fidelity $`F`$ (using the normalization relation) yields the optimal cloner:

$$x=\sqrt{\frac{F(1F)}{d1}},y=\frac{1F}{d1},v=F$$
(14)
The corresponding optimal fidelity for Eve is

$`F_E={\displaystyle \frac{F}{d}}+{\displaystyle \frac{(d1)(1F)}{d}}`$ (16)
$`+{\displaystyle \frac{2}{d}}\sqrt{(d1)F(1F)}`$
Let us see how Eve can maximize her information on Alice’s state. If Alice sends the state $`|k`$, then it is clear from Eq. (2) that Eve can obtain Bob’s error $`m`$ simply by performing a partial Bell measurement (measuring only the $`m`$ index) on $`EE^{}`$. Then, it appears from Eqs. (8) and (14) that, in order to infer Alice’s state, Eve must distinguish between $`d`$ non-orthogonal states (corresponding to all possible values of $`k`$) with a same scalar product $`(dF1)/(d1)`$ for all pairs of states, regardless the measured value of $`m`$. Consequently, Eve’s information $`I_{AE}`$ is simply given by the same expression as Eq. (9) but replacing $`F`$ by $`F_E`$. As a result, Bob’s and Eve’s information curves intersect exactly where the fidelities coincide, that is, at

$$F=F_E=\frac{1}{2}\left(1+\frac{1}{\sqrt{d}}\right)$$
(17)
We now use a theorem due to Csiszár and Körner, which provides a lower bound on the secret key rate, that is, the rate $`R`$ at which Alice and Bob can generate secret key bits via privacy amplification: if Alice, Bob and Eve share many independent realizations of a probability distribution $`p(a,b,e)`$, then there exists a protocol that generates a number of key bits per realization satisfying

$$R\mathrm{max}(I_{AB}I_{AE},I_{AB}I_{BE})$$
(18)
It is therefore sufficient that $`I_{AB}>I_{AE}`$ in order to establish a secret key with a non-zero rate. If we restrict ourselves to one-way communication on the classical channel, this actually is also a necessary condition. Consequently, the quantum cryptographic protocol above ceases to generate secret key bits precisely at the point where Eve’s information attains Bob’s information. In Table I, we have computed the disturbance $`D=1F`$ (or error rate) at which $`I_{AB}=I_{AE}`$ (or $`F=F_E`$), that is, above which Alice and Bob cannot distill a secret key any more by use of one-way privacy amplification protocols. Strictly speaking, since we only conjectured here that the cloning-based attack is optimal for all $`d`$, $`D_2^{\mathrm{ind}}`$ is actually an upper bound on $`D`$ that must necessarily be satisfied to generate secret key bits with one-way protocols. (A tighter upper bound on $`D`$ might exist if the optimal individual attack was not cloning, but we conjecture this is not the case). Interestingly, we note that $`D`$ increases with the dimension $`d`$, suggesting that a cryptosystem based on qudits is more secure for large $`d`$.

Now, we consider the second protocol where all $`d+1`$ bases are used. The cloner that must be used then is an asymmetric universal cloner, characterized by an amplitude matrix of the same form as (8) but with $`x=y`$, the normalization relation becoming $`v^2+(d^21)x^2=1`$. Here, Bob’s fidelity is given by $`F=v^2+(d1)x^2=1d(d1)x^2`$, so that the cloner is characterized by

$$x^2=\frac{1F}{d(d1)},v^2=\frac{(d+1)F1}{d}$$
(19)
As before, Bob’s information is given by Eq. (9). Eve’s clone is characterized by a matrix of the same form, with

$`x`$ $``$ $`x^{}=(vx)/d`$ (20)
$`v`$ $``$ $`v^{}=[v+(d^21)x]/d`$ (21)
so the corresponding fidelity is $`F_E=v^2+(d1)x^2=1d(d1)x^2`$. For deriving Eve’s information, we need first to rewrite the cloning transformation as

$`|k_A{\displaystyle \frac{vx}{\sqrt{d}}}|k_B{\displaystyle \underset{l=0}{\overset{d1}{}}}|l_E|l_E^{}`$ (23)
$`+x\sqrt{d}{\displaystyle \underset{m=0}{\overset{d1}{}}}|k+m_B|k_E|k+m_E^{}`$
After the basis is disclosed, Eve’s strategy is first to measure both $`E`$ and $`E^{}`$, the difference (modulo $`d`$) of the outcomes simply giving Bob’s error $`m`$. Making use of $`vx=x^{}d`$ and expressing $`x`$ and $`x^{}`$ as functions of $`F`$ and $`F_E`$, it is easy to check that the best Eve can do then is to use the state of her clone $`E`$ as an estimate of Alice’s state $`|k`$. If Bob makes no error ($`m=0`$), which happens with probability $`F`$, then it yields the correct value of $`k`$ with probability $`(F+F_E1)/F`$, while it yields any other of the $`d1`$ possibilities $`lk`$ with probability $`(1F_E)/[(d1)F]`$. In contrast, if Bob makes an error ($`m0`$), then Eve obtains the right $`k`$ with probability one. Consequently, the average mutual information between Alice and Eve conditionally on Bob’s error $`m`$ can be written as

$`I_{AE}=\mathrm{log}_2d+(F+F_E1)\mathrm{log}_2\left({\displaystyle \frac{F+F_E1}{F}}\right)`$ (25)
$`+(1F_E)\mathrm{log}_2\left({\displaystyle \frac{1F_E}{(d1)F}}\right)`$
One can check that, for a given $`F`$, $`I_{AE}`$ is slightly lower here than for the 2-bases protocol, which is consistent with the stronger requirement that we put on the cloner. Therefore, the fidelity $`F`$ at which $`I_{AB}=I_{AE}`$ is slightly lower, and the corresponding disturbance $`D=1F`$ is slightly higher. In Table I, we have shown the corresponding upper bound $`D_{d+1}^{\mathrm{ind}}`$ for several values of $`d`$, illustrating that there is a slight advantage in using all $`d+1`$ bases, as for the 6-state protocol for qubits.

Our last result concerns the most general eavesdropping strategy which consists in applying a coherent attack on a sequence of qudits of arbitrary (but finite) size $`n`$. Actually, our reasoning is simpler to state with a single qudit, but it remains valid for qudit sequences. We use an uncertainty principle due to Hall that puts a limit on the sum of Bob’s and Eve’s information: if $`\widehat{B}`$ and $`\widehat{E}`$ are Bob’s and Eve’s observables applied on the qudit sent by Alice, then

$$I_{AB}+I_{AE}2\mathrm{log}_2\left(d\underset{i,j}{\mathrm{max}}|b_i|e_j|\right)$$
(26)
where $`|b_i`$ and $`|e_j`$ are the eigenstates of $`\widehat{B}`$ and $`\widehat{E}`$, respectively. Since Eve has no way of guessing the basis used by Alice, her optimal information is the same for the correct and incorrect bases. Thus, one can bound $`I_{AE}`$ by assuming that Eve measures an observable $`\widehat{E}`$ complementary to $`\widehat{B}`$ (i.e. $`|b_i|e_j|=d^{1/2},i,j`$) :

$$I_{AB}+I_{AE}\mathrm{log}_2(d)$$
(27)
Using the discussion following Eq. (18), we conclude that $`I_{AB}>\mathrm{log}_2(d)/2`$ is a sufficient condition to warrant security against coherent attacks if the key is made out of a large number of independent realizations of $`n`$-qudit sequences (i.e., if the key is much longer than $`n`$). Using Eq. (9), this translates into a lower bound on $`F`$, or, equivalently, an upper bound on $`D`$

$$F\mathrm{log}_2\left(\frac{1}{F}\right)+(1F)\mathrm{log}_2\left(\frac{d1}{1F}\right)<\frac{\mathrm{log}_2d}{2}$$
(28)
which guarantees that one can distill secret key bits. This bound for coherent attacks, shown in Table I, exactly coincides with the well-known 11% bound on the error rate for $`d=2`$ due to Mayers.

In summary, we have extended standard quantum cryptography to protocols where the key is carried by quantum states in a space of arbitrary dimension $`d`$. We have used a general model of quantum cloning in order to calculate the information accessible to an eavesdropper monitoring the qudits individually. This provides an upper bound on the error rate above which the legitimate parties cannot distill a secret key by use of one-way privacy amplification protocols (since $`I_{AB}<I_{AE}`$). We conjectured that this bound is tight (i.e., applying the optimal cloner is the best strategy for Eve to gain the maximum information). Our analysis also suggested that the 2-bases protocol should be preferred to a $`(d+1)`$-bases one since its maximum acceptable error rate is only slightly lower, while the corresponding secret key rate is much larger. Finally, we have derived a very simple security proof of quantum cryptography with qudits that exploits an intuitive information inequality constraining the simultaneous measurement of conjugate observables. This results in an upper bound on the acceptable error rate that is more restrictive than the previous one, but guarantees that a non-zero secret key rate can always be produced (even with coherent attacks on sequences of finite length). In the region between these two bounds, it is unknown whether the security is guaranteed or not. It should be stressed that all the bounds on $`D`$ discussed above tend to 1/2 for $`d\mathrm{}`$, reflecting the advantage of using higher-dimensional spaces. However, practical limitations might be more severe in realistic qudit-based cryptosystems, in particular the influence of the detector’s quantum efficiency and dark count rate. This is discussed in a related paper.

Note: After completion of this work, it was proven in an independent paper that the optimal individual attack for qutrits ($`d=3`$) when using all four mutually unbiased bases exactly coincides with our results based on the optimal cloning machine, as conjectured here.


text:
# The Effect of a Binary Source Companion on the Astrometric Microlensing Behavior

## 1 Introduction

If microlensing occurs in a binary-source system, both source components are gravitationally magnified and the resulting light curve deviates from the standard one of a single source event (Griest & Hu 1992; Han & Gould 1997). However, the photometric binary-source effect becomes important only for the rare cases of the lens trajectory’s passage close to both source stars and in most cases only one source component is significantly magnified. Then, the less magnified source component (the companion) can be treated simply as a blending source, and the resulting light curve may then be approximated by that of a blended single source event (Dominik 1998).

Although microlensing events have till now been observed only photometrically (Alcock et al. 1993; Aubourg et al. 1993, Udalski et al. 1993; Alard & Guibert 1997; Abe et al. 1997), they could also be observed astrometrically by using high-precision instruments that will become available in the near future, e.g. the Space Interferometry Mission (SIM, Unwin et al. 1998), and the interferometers to be mounted on the Keck (Colavita et al. 1998) and VLT (Mariotti et al. 1998). If an event is astrometrically observed, one can measure the displacements in the source star image centeroid with respect to its unlensed position. Once the trajectory of the centroid shifts (the astrometric curve) is measured, the physical parameters (mass and location) of the lens can be better constrained (Miyamoto & Yoshii 1995; Høg, Novikov & Polnarev 1995; Walker 1995; Boden, Shao & Van Buren 1998; Han & Chang 1999).

The binary source companion also affects the astrometric lensing behavior. However, whenever its astrometric effect was considered (Han & Kim 1999; Dalal & Griest 2001), it was always treated as a simple blending source due to the belief that the astrometric binary-source effect would be similar to the photometric one. In this paper, we show that contrary to this belief the effect of the source companion on the astrometric lensing behaviors of a significant fraction of binary-source events cannot be sufficiently well described by the blending approximation, even when the corresponding light curves are well approximated by the blending effect alone.

The paper is organized as follows. In § 2, we describe the basics of binary source and blending effects on the light and astrometric curves of lensing events. In § 3, we the present astrometric curves and the corresponding light curves of example binary-source events to illustrate the relative difficulty in describing the astrometric lensing behaviors by the blending approximation alone. In § 4, we statistically estimate the range of the binary source separation and the companion light fraction where the blending approximation breaks down for the description of the photometric and astrometric lensing behaviors. We conclude with § 5.

## 2 Binary source and Blending Effects

### 2.1 Binary Source Effect

When a single source is lensed, its image is split into two. The two images appear at the positions

$$𝜽_{\mathrm{I}\pm }=\frac{\theta _\mathrm{E}}{2}\left(u\pm \sqrt{u^2+4}\right)\widehat{𝒖},$$
(1)
and have magnifications

$$A_\pm =\frac{1}{2}\left(\frac{u^2+2}{u\sqrt{u^2+4}}\pm 1\right),$$
(2)
where $`\theta _\mathrm{E}`$ is the angular Einstein ring radius, $`𝒖`$ is the lens-source separation vector normalized by $`\theta _\mathrm{E}`$, and $`\widehat{𝒖}`$ represents its unit vector. The angular Einstein ring radius is related to the physical parameters of the lens by

$$\theta _\mathrm{E}=\sqrt{\frac{4Gm}{c^2}}\left(\frac{1}{D_{\mathrm{ol}}}\frac{1}{D_{\mathrm{os}}}\right)^{1/2},$$
(3)
where $`m`$ is the lens mass and $`D_{\mathrm{ol}}`$ and $`D_{\mathrm{os}}`$ represent the distances to the lens and source star, respectively. The lens-source separation vector is related to the lensing parameters by

$$𝒖=\left(\frac{tt_0}{t_\mathrm{E}}\right)\widehat{𝒙}+\beta \widehat{𝒚},$$
(4)
where $`t_\mathrm{E}`$ represents the time required for the source to transit $`\theta _\mathrm{E}`$ (Einstein time scale), $`\beta `$ is the closest lens-source separation in units of $`\theta _\mathrm{E}`$ (impact parameter), $`t_0`$ is the time at that moment, and the unit vectors $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ are parallel with and normal to the direction of the lens-source motion. The two images formed by the lens cannot be resolved, but one can measure the total magnification and the shift of the light centroid of the source star with respect to its unlensed position, which are represented respectively by

$$A=A_++A_{}=\frac{u^2+2}{u\sqrt{u^2+4}},$$
(5)
and

$$𝜹=\frac{A_+𝜽_{\mathrm{I}+}+A_{}𝜽_\mathrm{I}}{A}𝒖=\frac{𝒖}{u^2+2}\theta _\mathrm{E}.$$
(6)
The resulting light curve is symmetric with respect to $`t_0`$ and the trajectory of the centroid shift traces an ellipse with a semi-major axis $`a=(\beta ^2+2)^{1/2}\theta _\mathrm{E}/2`$ and an axis ratio $`q=\beta (\beta ^2+2)^{1/2}`$ (Walker 1995; Jeong, Han & Park 1999).

If lensing occurs in a binary-source system, on the other hand, both source components are gravitationally magnified. For the binary-source event, the lensing behavior involved with each source can be treated as an independent single source event. Then the observed light and astrometric curves are represented by

$$F_{\mathrm{BS}}=A_1F_1+A_2F_2,$$
(7)
and

$$𝜹_{\mathrm{BS}}=\frac{A_1F_1(𝒖_1+𝜹_1)+A_2F_2(𝒖_2+𝜹_2)}{A_1F_1+A_2F_2}\frac{F_1𝒖_1+F_2𝒖_2}{F_1+F_2},$$
(8)
where $`𝒖_i`$ are the separation vectors between the lens and the individual source components, $`A_i`$ and $`𝜹_i`$ are the magnifications and the centroid shifts of the individual single source events with baseline source fluxes $`F_i`$, and the subscripts $`i=1`$ and 2 are used to denote parameters related to the individual single source events ($`i=1`$ for the event with the higher magnification). We note that the reference position of the centroid shift measurements for the binary-source event is not the location of the primary, i.e. $`𝒖_1`$, but the center of light between the unlensed source components, i.e. the second term in eq. (8). In Figure 1, we present the lens system geometry of an example binary-source event.

### 2.2 Blending Effect

If a single source event is affected by the blended light from an unresolved nearby star, the observed light curve is represented by

$$F_{\mathrm{blend}}=AF_0+F_\mathrm{b},$$
(9)
where $`F_0`$ is the unblended baseline source flux and $`F_\mathrm{b}`$ is the amount of blended flux. Since $`F_\mathrm{b}`$ is constant throughout the event, the observed flux is uniformly increased by $`F_\mathrm{b}`$ and the shape of the light curve is still symmetric.

The centroid shift of a blended event is represented by

$$𝜹_{\mathrm{blend}}=\frac{AF_0(𝒖+𝜹)+F_b𝒖_\mathrm{b}}{AF_0+F_\mathrm{b}}\frac{F_0𝒖+F_\mathrm{b}𝒖_\mathrm{b}}{F_0+F_\mathrm{b}},$$
(10)
where $`𝒖_\mathrm{b}`$ is the position vector of the blending source with respect to the lens (Han & Kim 1999; Dalal & Griest 2001). Blending has two effects on the observed astrometric curve. Firstly, it causes the source star image centroid to be additionally shifted towards the blending source during the event. Secondly, due to blending, the reference position of the centroid shift measurements is no longer the position of the lensed source but the center of light between the lensed star and the blending source. As a result, depending on the position and the light fraction of the blending source, $`f=F_\mathrm{b}/(F_0+F_\mathrm{b})`$, the resulting astrometric curves can be seriously distorted from the elliptical one of the unblended event (Han & Kim 1999; Dalal & Griest 2001).

## 3 Blending Approximation

In many cases, the light curves of binary-source events resemble those of single source events. This is because the lens, in general, approaches closely only one of the source components. In this case, the magnification of the companion is small, i.e. $`A_21.0`$, and thus the the resulting light curve can be approximated as

$$A_{\mathrm{BS}}A_1F_1+F_2.$$
(11)
This implies that the observed light curve is well described by that of a single source event where the companion simply acts as a blending source.

Can the blending approximation be equally applied to astrometric curves of most binary-source events? To this question, we find the answer is ‘no’. We demonstrate this in Figure 2, where we present the astrometric (middle panels) and light curves (lower panels) of binary-source events with three different values of source separation, and we compare these curves to those obtained by the blending approximation. One finds that even when the light curve is well approximated by the blending effect alone, the deviation of the astrometric curve from the approximation is considerable.

To see the patterns of both the photometric and astrometric deviations from the blending approximation for more general cases of binary-source events, we construct contour maps of magnification and centroid shift excesses for binary-source systems with various values of the binary separation, $`\mathrm{}`$, and the companion light fraction, $`f`$. The magnification and centroid shift excesses are defined respectively by

$$ϵ=\frac{A_{\mathrm{BS}}A_{\mathrm{blend}}}{A_{\mathrm{blend}}},$$
(12)
and

$$\mathrm{\Delta }𝜹=𝜹_{\mathrm{BS}}𝜹_{\mathrm{blend}},$$
(13)
where $`A_{\mathrm{BS}}`$ and $`𝜹_{\mathrm{BS}}`$ are the exact magnification and centroid shift of the binary-source event, and $`A_{\mathrm{blend}}`$ and $`𝜹_{\mathrm{blend}}`$ are those obtained by the blending approximation. The constructed maps are presented in Figure 3 for the excess magnification, and in Figure 4 for the excess centroid shift, respectively. From a comparison of the two maps, one finds that while significant photometric deviations occur only in a small region around the companion, comparable astrometric deviations occur in a substantially larger area, even for a companion with a large separation and a small light fraction.

The relative difficulty in describing the astrometric curves of binary-source events by just the blending approximation can be understood in the following way. If a source companion is located at a large distance from the lens (i.e. $`u_2\sqrt{2}`$), its contribution to the magnification and the centroid shift are represented by

$$A_21+\frac{2}{u_2^4},$$
(14)
and

$$\delta _2\frac{\theta _\mathrm{E}}{u_2}.$$
(15)
Then, as the separation becomes larger, the photometric contribution falls off rapidly ($`u_2^4`$), while the astrometric contribution decays much more slowly ($`u_2`$) (Miralda-Escudé 1996; Paczyński 1998). As a result, even at the location where the magnification of the companion is negligible (i.e. $`A_21`$), the amount of the centroid shift can be considerable. In this range of $`u_2`$, the centroid shift becomes

$$𝜹_{\mathrm{BS}}\frac{A_1F_1(𝒖_1+𝜹_1)+F_2(𝒖_2+𝜹_2)}{A_1F_1+F_2}\frac{F_1𝒖_1+F_2𝒖_2}{F_1+F_2},$$
(16)
which differs from $`𝜹_{\mathrm{blend}}`$ in eq. (10).

## 4 Validity of the Blending Approximation

In the previous section, we illustrated the relative difficulties in describing the astrometric lensing behavior of binary-source events by the blending approximation, as compared to the photometric behavior. Then, a natural question to ask is: over what range of values of the source separation and the companion light fraction does the blending approximation break down for the description of the photometric and astrometric behaviors of binary-source events? We answer this question as follows.

We proceed by statistically estimating the probabilities of binary-source events whose light and astrometric curves can be approximated by the blending effect as functions of $`\mathrm{}`$ and $`f`$; $`P_{\mathrm{ph}}(\mathrm{},f)`$ for the photometric and $`P_{\mathrm{ast}}(\mathrm{},f)`$ for the astrometric measurements. To determine these probabilities, we first simulate a large number of microlensing events occurred in binary-source systems, characterized by $`\mathrm{}`$ and $`f`$. We then determine the probabilities by estimating the fraction of events whose light and astrometric curves have deviations (from those obtained by the blending approximation) less than preselected threshold values. The lens trajectories of the tested events are selected so that they have random orientations with respect to the binary axis, and so that their impact parameters with respect to the primary source are uniformly distributed in the range of $`0\beta _11`$. We assume that the blending approximation is valid if the light and astrometric curves, measured during the time period $`10t_\mathrm{E}t10t_\mathrm{E}`$, have deviations less than the threshold values of $`ϵ_{\mathrm{th}}=10\%`$ and $`\mathrm{\Delta }\delta _{\mathrm{th}}=0.1\theta _\mathrm{E}`$, respectively.

In Figure 5, we present the resulting probabilities as contour maps in the parameter space of $`\mathrm{}`$ and $`f`$: upper panel for $`P_{\mathrm{ast}}(\mathrm{},f)`$ and lower panel for $`P_{\mathrm{ast}}(\mathrm{},f)`$. From the figure, one finds that for a given $`\mathrm{}`$ and $`f`$, $`P_{\mathrm{ast}}`$ is significantly lower than the corresponding $`P_{\mathrm{ph}}`$. For example, if microlensing occurs in a binary-source system with $`\mathrm{}=6.0\theta _\mathrm{E}`$, we estimate that the probabilities of events whose astrometric curves can be treated by the blending approximation will be only $`P_{\mathrm{ast}}0\%`$, 20%, 50%, and 70% for companions with light fractions of $`f=0.8`$, 0.6, 0.4, and 0.2, respectively; while the corresponding photometric probabilities will be $`P_{\mathrm{ph}}>80\%`$ for all four values of $`f`$. One also finds that the region of low $`P_{\mathrm{ast}}`$ occupies a significantly larger area in the parameter space than the low $`P_{\mathrm{ph}}`$ region, implying that the astrometric blending approximation will not be valid for a significant fraction of binary-source events.

## 5 Conclusion

We have shown that although the blending approximation describes well the photometric lensing behaviors of most binary-source events, the approximation will not be able to adequately describe the astrometric lensing behaviors for a significant fraction of binary-source events. This is because the astrometric effect of the companion endures to relatively large distances where the corresponding photometric effect is negligible. Therefore, it will be important to take the full effect of the binary companion into consideration in determining the lens parameters from the observed source motion.

We thank André Fletcher (Korea Astronomy Observatory) for his help with the preparation of the paper. This work was supported by a grant (1999-2-1-13-001-5) from the Korea Science & Engineering Foundation (KOSEF).