text:
# Parallelization of Cellular Automata for Surface Reactions

## I Introduction

One of the most interesting features of surface reactions is that in a large number of cases produce pattern formation, structures with some well–defined length scale, sometimes with symmetries and temporal behavior, such as oscillations, traveling waves, spirals, Turing patterns, etc rab2000 ; imb1995 . A usual approach to study this pattern formation is reaction–diffusion (RD) equations wal1997 , which simulate the dynamic behavior of chemical reactions on surfaces. However, these partial differential equations give only approximate solutions and in several cases completely wrong results, because they are based on the local mean field approximation, meaning well–mixed reactants at microscopic level, ignoring all the correlation terms between reactants locally, fluctuations and lateral interactions in the adsorbate win2002 . The RD equations describe the coverage which are macroscopic continuum variables which neglect the discrete structure of matter, and do not describe the actual chemical process underlying pattern formation. In fact from experimental studies win2002  it is known that to model correctly, a modified kinetic has to be assumed different from the prescribed for RD.

Based in some general assumptions of the physics processes involved, an esencial master equation can be derived that completely describes the microscopic dynamic of the system kam1981 . An exact method to solve this master equation is the Monte Carlo (MC) method gil1977 ; jan1995 ; luk1998 ; zve2001 . In a Monte Carlo method a sequence of discrete events (reactions, including diffusion) is generated on a 2D lattice which represents the surface sites. These events are generated in general in a random way, looking for possible enabled reactions between nearest neighbors on the lattice, and doing the reactions with some probability in correspondence with some defined reaction rates.

To compare MC simulations with experimental pattern formation it is necessary to fill the gap between the length scale of the individual particles and the diffusion length. The regime where spatio–temporal patterns usually occurs, from $`\mu `$m to mm scales, is orders of magnitude larger than the nm scale of individual particles. However, some new experiments win1997 ; vol1999  show that the fast kinetics processes are typically accompanied by the appearance of nanostructures. Only microscopic simulations could deal with this two–scales behavior. However this regime implies lattices sizes above $`10^6\times 10^6`$ and very large values for the diffusion rates to produce agreement with the experimental observed scales on pattern formation. This would be a very large and slow simulation, due to the huge number of particles involved and the fast diffusion rates which means that most of the simulation time is spend doing diffusion of particles instead of chemical reactions. Fortunately, we do not need to simulate every time such a macroscopic system or experimental diffusion rates. We only need to find scaling laws for lengths and diffusion coefficients, i.e from nm scale to $`\mu `$m or mm scale, and from microscopic to real time. This is the possibility which we explore in this paper by using large simulations within parallel computers.

Although only MC simulations provides solutions to the exact master equations for the surface reactions, they are not suitable for efficient parallelization due to the random selection of lattice sites used. However, there is another important approach to simulate discrete events on lattices, the Cellular Automata (CA) tof1987 ; wei1997 . This approach is fully parallel in the sense that all the lattice sites can be updated simultaneously. This has the advantages that fewer random numbers are required, only those for the reaction probabilities (CA codes are faster than MC), and the global updating is easier to implement in a parallel code.

Under some well–defined conditions a CA is equivalent to a MC simulation (see ref. kor1998  and section II). Then a CA could be used as a MC simulation and its parallelization will provide with the required scaling laws. In this paper we present for first time an attempt to provide a tool to get these scaling laws. The fact that CA are ideal simulation methods for parallelization is shown in a recent special issue dedicated to CA in Parallel Computing rev2000 . In particular a quite interesting paper by J.R. Weimar wei2000  presents an object oriented parallel CA, also based in ref. kor1998 , and explores the possibility of divide the surface in regions to be simulated by RD or CA according to the level of detail required.

In this paper we will describe in detail the implementation of the parallel version of a CA simulating chemical reactions on surfaces, but the ideas discussed here are also applicable to more general CA simulations. In section II we will describe briefly the Cellular Automata method. In section III we explain the key points to implement the parallelization of the algorithm. In section IV we will present results of the application of the parallel algorithm to the A$`+`$B$`0`$ model, and to a realistic model of oxidation of CO on Pt$`(110)`$. Finally in section V we will draw some conclusion.

## II Cellular Automata

A cellular automaton (CA) is a regular array of cells. Each cell can be in one of a set of possible states. The CA evolves in time in discrete steps by changing the states of all cells simultaneously. The next state which a cell will take is based on the previous state of the cell and the states of some neighboring cells. A CA is defined by providing prescriptions for the lattice, the set of states, the neighbors, and the transition rules tof1987 ; wei1997 . The idea of CA deserves by itself a lot of study, since in general the evolution of a CA cannot be predicted other than by executing it. Additionally the amount of possible CAs is quite large. For instance, in the simplest one–dimensional case, with two states and two neighbors, there are $`256`$ possible CA.

Our CA uses the standard square two–dimensional lattice of size $`L\times L`$. The states of each cell represent occupation with some kind of particle. The Margolus neighborhood tof1987  is used instead of the common von Neumann neighborhood. Both definitions of neighborhood are shown in Fig.1. In order to obey the CA laws in von Neumann neighborhood, it is necessary cho1988  disobey the laws of stoichiometry because one particle could participate in more than one reactive pair. Similar problems arise with the diffusion of particles cho1991 . The use of a Margolus neighborhood overcomes these difficulties mai1991 ; mai1992 ; mai1992b ; mai1993 . Using the value of the chemical rates we set up some probabilistic transition rules to change the states of the cells inside each Margolus block.

The Margolus blocks, Fig.1, are used in the following way. The blocks are periodically repeated to build up a tiled mask over the whole lattice, considering periodic boundary conditions. In this way there are four possible tilings as shown in Fig.2. Only neighbor sites belonging to the same block can react, so all the blocks can be accessed in parallel. Inside each block a Monte Carlo update is done. After a full lattice update the whole procedure starts over again using another of the four possible tilings. The dynamics is not confined to blocks, because the boundary between blocks is changing from one global sweep to the next.

For the four tilings shown in Fig.2, we chose randomly one of the four possible sequences of tilings: ($`1`$,$`2`$,$`3`$,$`4`$), ($`2`$,$`3`$,$`4`$,$`1`$), ($`3`$,$`4`$,$`1`$,$`2`$), and ($`4`$,$`1`$,$`2`$,$`3`$). They show always the same clockwise cyclic sequence, but starting at a different tiling. This produces better boundary diffusion and mixing between blocks.

A schematic description of the steps to implement this CA is given at the following:

* Choose randomly one of the four possible sequences of tilings: ($`1`$,$`2`$,$`3`$,$`4`$), ($`2`$,$`3`$,$`4`$,$`1`$), ($`3`$,$`4`$,$`1`$,$`2`$), ($`4`$,$`1`$,$`2`$,$`3`$).
* Choose consecutive tilings from the sequence of step 1.
* Sweep over all the blocks and in each block make a single Monte Carlo update:

+ Choose randomly one pair of neighbors.
    + Choose a reaction $`i`$ from the set of all the possible reactions with a probability proportional to the reaction rate $`k_i`$.
    + Check if the reaction chosen in step 3.2 is possible on the sites chosen in step 3.1. If it is possible do the reaction, otherwise do nothing.
* Increase the time by $`\mathrm{\Delta }t/4`$, and return to step 2 until the sequence of tilings is completed.
* Return to step 1.

The updating scheme inside each block is the same as in the Dynamic Monte Carlo algorithm called Random Selection Method luk1998 ; zve2001 . Consequently, the time increment is the same as one Monte Carlo Step (MCS) luk1998 :

$$\mathrm{\Delta }t=\frac{1}{L^2_ik_i}$$
(1)
where the sum in the denominator corresponds to the sum of all the possible reaction rates.

There is a large number of possible CA prescriptions which could try to reproduce a MC simulation of chemical reactions on surfaces. In fact, an extensive study made by J. Mai mai1991 ; mai1992 ; mai1992b ; mai1993  shows that in general it is difficult produce good CAs reproducing chemical reactions and diffusion. However, in a recent paper Kortlüke kor1998  studies under which conditions a CA could reproduce adequately a MC simulation of chemical reactions on surfaces. He found that the main requirement is using large diffusion coefficients. It is required some sort of compromise between MC and CA: it is necessary use a regular array of blocks as in CA, and a MC update scheme has to be used inside each block. The CA which we use here is a small modification of that used by Kortlüke originally kor1998 . We use blocks of $`4`$–sites (Margolus blocks) instead of $`2`$–sites (Hantels blocks). The main advantage of using this $`4`$–sites instead of $`2`$–sites is that reduces the level of CA noise (see kor1998 ) inside each block. This CA noise is the difference between the diffusion simulated with CA and the correct diffusion simulated with MC. In fact using larger blocks the noise will be smaller.

Additionally we present here a full realization of the CA parallel computing idea by implementing this in a parallel code as is shown in the next section. This is the main reason for use CA instead MC as a solution for large–time and large–sizes simulations. Otherwise the only use of CA instead MC in a serial simulation does not speed up the simulation substantially, in the best cases only $`10`$$`20\%`$. This means that the implementation of an efficient CA is very important.

Despite that we have to accept that our CA is an approximation to exact MC microscopic simulations, we known when the approximation is good and we can always check that the CA is reproducing the MC results. In this way we can think about the CA as the MC realization of the exact master equation for the chemical processes, getting in this way, the possibility of obtain the scaling laws mentioned at the introduction for these physical systems. In this paper we do not discus the problem of quality of the CA approximation to the MC realization. This was done in kor1998  and we have made for our modified CA a similar study with similar results.

The fact that blocks can be updated without referring to other blocks allows us to do CA as a parallel algorithm, because the whole set of blocks can be updated at the same time. This is the subject of the next section.

## III Parallelization

The aim of the parallelization of any code is to distribute the whole simulation over several computer processors, also called nodes. We define $`speed(P)=1/T_P`$, where $`T_P`$ is the computer time used by $`P`$–processors to complete a simulation. The optimum result is that this multiplies the speed of the simulation with the number of nodes, $`speed(P)=P\times speed(1)`$. However, the usual result of parallelize a code is not the optimum. The key point to achieve a good speed up is that the time spent by each node sending and receiving information from/to the other nodes is small in comparison with the time consumed doing computing in each node.

In order to distribute the work, we use here a geometrical division as usual rev2000  in spatial extended systems, i.e., the full lattice is divide in sublattices of equal area. The time spent sending data is proportional to the length of the borders, $`L`$, and the time doing computing is proportional to the system size $`L^2`$. We divide the lattice in strips as is shown in Fig.3. Considering the global periodic boundary conditions, each node has periodicity in one direction and in the other directions has to share information with only two neighbor nodes. Another possibility is divide the lattice in squares. This choice is less convenient. Each node has to share information with 4 nodes instead 2, and it is only possible to use a perfect square number of nodes $`P=4,9,16,\mathrm{}`$. Using the strips sublattices requires sending a factor $`4/\sqrt{P}`$ less data than using the square sublattices note2 , provided to have $`P16`$.

Note that due to the periodicity in the vertical direction inside each node, it is not necessary to interchange border information between nodes when we pass from the tiling $`12`$ and $`34`$. It is only necessary when we pass from $`23`$ and $`41`$. This reduces in half the amount of communication data. From Fig.3 we see the direction in which the first column of data in each sublattice has to be be sent and received in each node. In the change of tilings $`23`$ the data goes from left to right and in $`41`$ from right to left. Instead of shifting the data of each node horizontally we use the first column for sending to or receiving from other nodes. As a consequence, and to get interaction between sites in the same block, we consider the first column neighbor of the last one. For this Blocking CA, where the reactions are confined inside the blocks, this produces an effective periodicity in the horizontal direction. This extra periodic boundary condition inside each sublattice, makes the final parallel code very similar to a full lattice implementation.

We use the SPMD (single program, multiple data) model to implement the parallel simulation, in which every node executes the same code using different sublattice. Also, we use a main node, called node zero, dedicated additionally to distribute and collect the global information needed for the input/output of the simulation and other global computing required. Every node executes the same code, and when some special code has to be executed for only some nodes, is necessary use a node identifier number $`p`$. From the global periodic boundary conditions and the shape of the sublattices, there is periodicity also in the sequence of nodes: the node left with respect to node $`p=0`$ is node $`p=P1`$ (P is the total number of nodes) and the node right with respect to node $`p=P1`$ is node $`p=0`$.

In the following we present the CA code for a single node.

* $`p=0`$: Chose randomly one of the four possible tiling sequences: ($`1`$,$`2`$,$`3`$,$`4`$), ($`2`$,$`3`$,$`4`$,$`1`$), ($`3`$,$`4`$,$`1`$,$`2`$), ($`4`$,$`1`$,$`2`$,$`3`$). Send that choice to the other nodes.
$`p>0`$: Receive the information of which tiling sequences has been chosen.
* $`p`$: Choose consecutive tilings from the sequence of step 1.
* $`p`$: If the previous tiling was $`1`$ or $`2`$ and the new tiling is $`3`$ or $`4`$ then
send the first column to the left node $`p1`$,
receive the data from right node $`p+1`$ and put it in the first column.
* $`p`$: If the previous tiling was $`3`$ or $`4`$ and the new tiling is $`1`$ or $`2`$ then
send the first column to the right node $`p+1`$
receive the data from left node $`p1`$ and put it in the first column.
* $`p`$: Sweep over all the blocks and in each block make a single Monte Carlo update.
* $`p`$: Increase the time by $`\mathrm{\Delta }t/4`$, and return to step 2 until the tiling sequence is completed.
* $`p`$: Return to step 1.

We have modified step $`1`$ and added two new steps $`N1`$ and $`N2`$, to interchange information between the nodes. The rest is basically the same as the single processor code, but applied to the respective sublattices for each node. In order to avoid undesired correlations between different sublattices, special attention should be given to a good random number generator, producing different sequence of random numbers within each node knu1998 .

For implementation of the interprocess communication and synchronization, the message passing interface MPI library has been selected, because it provides source portability to different kind of computers. This library provide, amongst a set of specialized and complete communication routines, a so called six–basic set of routines for interprocessor communication: initialization, termination, getting the set of nodes, getting its own node number, sending data to other nodes, and receiving data from other nodes.

There are computers with connection between nodes using giga–ethernet or several processors inside the same machine, sharing the same memory. These computers represent optimal environments to test parallel algorithms, i.e. high–end supercomputers like Cray T3E, and middle–end supercomputers like Silicon Origin 2000. However, we will show in the next section that the results of running this parallel CA algorithm in a low–end Beowulf cluster of PCs connected only via fast–ethernet produce already almost the ideal speed up.

## IV Results

The parallel algorithm was tested on our local cluster of PCs, (17 Athlon 1.1Ghz/256Mb, fast–ethernet, Linux 2.4.18, MPICH 1.1.2). From the results we see that the improvement of the performance using $`P`$ processors with respect to a single processor is almost the ideal, $`speed(P)=P\times speed(1)`$. In order to test the speed up of the algorithm, we will use two systems from surface catalysis: the A$`+`$B$`0`$ model and a model for CO oxidation on Pt$`(110)`$ surface kuz1998 , which has produced several important results kor1998b ; kor1999 ; kuz1999 ; kor1999b .

The A$`+`$B$`0`$ model has been studied for a long time ovc1978 ; kuz1988 ; kot1996 ; pri1997 ; mar1999 ; arg2001  In the pioneering analitycal paper by Ovchinnikov and Zeldovich ovc1978  it was shown for first time that the kinetic law of mass action is violated in this model, producing incorrect results when standard chemical kinetics is used. An illustrative case, also used in our paper, is a situation with equal concentrations of both reactants $`C_A=C_B=C`$, where the standard kinetics predicts an asymptotic behavior $`Ct^1`$. This prediction correspond to the mean–field approximation and it is only valid for high dimensional systems ovc1978 , $`𝒟`$4. In the low dimensional systems $`𝒟<`$4 with diffusion controlled processes it has been proved, using renormalization group arguments lee1995 , that the correct asymptotic behavior is: $`Ct^{𝒟/4}`$. A qualitative agreement with this behavior was shown for first time using MC simulations in reference tou1983 .

The asymptotic law needs large simulation time $`t_{max}`$. The diffusion length $`\xi (t)=\sqrt{Dt}`$ defines the pattern formation scale. A simulation until $`t=t_{max}`$ needs a lattice length $`L\xi (t_{max})`$, which correspond to a large simulation time of the order of $`t_{max}L^2t_{max}^2`$. This case provides a good example of large–time and large–size systems with pattern formation to test our parallel algorithm described in the previous section.

In the corresponding lattice model we consider two kinds of particles A and B. The only possible chemical reaction is desorption of AB, which happen when two particles A and B are next to each other, creating two empty sites. This process occurs with a rate constant $`k`$. Additionally the particles are allowed to diffuse with rate $`D`$. This happen when a particle is next to an empty site. We simulate the behavior produced for an initial condition without empty sites and where the same number of A and B particles are initially randomly distributed on the surface: $`N_A=N_B=L^2/2`$.

In Fig.4 we present the temporal behavior and also illustrate the segregation process forming regions with high concentration of A or B, which increase in size with time. Also we show how the global concentration, $`C=N_A/L^2=N_B/L^2`$, diminishes with time following the asymptotic power–law $`Ct^{\frac{1}{2}}`$. The system size used is $`L=8192`$ and the parameters are $`k=D=1`$. We present two sets of data in Fig.4. The points correspond to a single processor simulation, and the solid line corresponds to a parallel simulation using $`16`$ processors. Both simulations start with identical random initial distribution of particles. It is noticeable that this initial distribution almost completely determines the following behavior of the system. The snapshots shown as insert in Fig.4 are from the parallel simulation. They are very similar to the ones obtained from the single processor simulation, which uses the same initial conditions but a different sequence of random numbers to simulate the dynamics. Moreover, in order to check quantitatively the agreement between spatial structures between the single and parallel simulations, we present in Fig.5 the radial correlation function. Again the points correspond to the single processor case and the solid line to the parallel case. A correlation length $`r_c`$ could be obtained by fitting these correlation functions with $`\mathrm{exp}[(r/r_c)^2]`$, which we show in Fig.5 by dashed lines. The obtained values of $`r_c`$ are plotted in the insert. This shows that the correlation length for this dynamic also follows a power–law $`r_ct^\alpha `$. By numerical fitting we obtain $`\alpha =0.5122\pm 0.012`$. An analytical asymptotical solution for this correlation functions is given in kot1996 , $`\mathrm{exp}(r^2/4Dt)`$ or $`\mathrm{exp}[c(r/\xi (t))^2]`$. The previous value obtained for $`\alpha `$ means that the diffusion length $`\xi (t)=\sqrt{Dt}`$ define the scale of pattern formation for this reaction model.

The performance or speed up of the parallel algorithm for the A$`+`$B$`0`$ model is shown in Fig.6, using different system sizes $`L=256,512,1024,4096,8192`$ and number of processors $`P=1,2,4,8,16`$. The simulated time for each computation was set to $`t_{max}=500`$. The speed was normalized to the speed of the single processor case. In the insert we can see the behavior of the single processor speed for each system sizes. The advantage of use a large number of processors increases when the system size increases, as we expect from the discussion in the previous section.

The second model used to test the parallel algorithm is a model for CO oxidation on Pt$`(100)`$ and Pt$`(110)`$ surfaces kuz1998 ; kor1998b ; kor1999 ; kuz1999 ; kor1999b . This system shows different types of kinetic oscillations. On Pt$`(100)`$ local, irregular oscillations occur in a wide parameter interval, whereas on Pt$`(110)`$ globally synchronized oscillations exist only in a very narrow parameter interval. Both surfaces exhibit an $`\alpha \beta `$ surface reconstruction, where $`\alpha `$ denotes the hex or $`1\times 2`$ phase on Pt$`(100)`$ or Pt$`(110)`$, respectively. $`\beta `$ denotes the unreconstructed $`1\times 1`$ phase in both cases. Both surfaces have qualitatively quite similar properties with the exception of the dissociative adsorption of O<sub>2</sub>. The ratio of the sticking coefficients of O<sub>2</sub> on the two phases is $`s_\alpha :s_\beta 0.5:1`$ for Pt$`(110)`$ and $`s_\alpha :s_\beta 10^2:1`$ for Pt$`(100)`$ imb1995 . From the experiments imb1995  it is known that kinetic oscillations are closely connected with the $`\alpha \beta `$ reconstruction of the Pt surfaces.

In the model kuz1998 ; kor1998b ; kor1999 ; kuz1999 ; kor1999b , CO is able to absorb onto a free surface site with rate constant $`y`$ and to desorb from the surface with rate constant $`k`$, independent of the surface phase to which the site belongs. O<sub>2</sub> adsorbs dissociatively onto two nearest neighbor sites with rate constant $`(1y)s_\chi `$ with $`\chi =\alpha `$,$`\beta `$. In addition CO is able to diffuse via hopping onto a vacant nearest neighbor site with rate constant $`D`$. The CO$`+`$O reaction occurs with rate constant $`R`$, when CO and O are in nearest neighbor sites desorbing the reaction product CO<sub>2</sub>. The $`\alpha \beta `$ surface phase transition is modeled as a linear front propagation induced by the presence of CO in the border between phases with rate constant $`V`$. Consider two nearest neighbor surface sites in the state $`\alpha \beta `$. The transition $`\alpha \beta \alpha \alpha `$ ($`\alpha \beta \beta \beta `$) occurs if none (at least one) of these two sites is occupied by CO. Summarizing the above transition definitions written in the more usual form of reaction equations gives:

$`\text{CO(g)}+S^\chi \text{CO(a)},`$
$`\text{O}\text{2}\text{(g)}+2S^\alpha 2\text{O(a)},`$
$`\text{O}\text{2}\text{(g)}+2S^\beta 2\text{O(a)},`$
$`\text{CO(a)}+S^\chi S^\chi +\text{CO(a)},`$
$`\text{CO(a)}+\text{O(a)}\text{CO}\text{2}\text{(g)}+2S^\chi ,`$
$`S^\alpha S^\beta ,`$
where $`S`$ stands for a free adsorption site, $`\chi `$ stands for either $`\alpha `$ or $`\beta `$ and (a) or (g) for a particle adsorbed on the surface or in the gas phase, respectively. For additional details see ref. kuz1998 .

Amongst several successful results of this model, we can mention that it was one of the first microscopic models for CO oxidation on Pt including surface reconstruction, which is nowadays widely accepted as the key element in order to get oscillatory behavior. This model reproduce correctly oscillatory regimes for both surfaces Pt$`(100)`$ and Pt$`(110)`$, by changing only one parameter $`s_\alpha `$. The diffusion of CO is consider explicitly and could be applied to the fast diffusion regime without modification. With this model an alternative mechanism for global synchronization of oscillation has been suggested kor1999 , different from the traditional gas–phase coupling. This new mechanism is stochastic resonance, obtained by including a spontaneous nucleation of one surface phase in the other $`\alpha \beta `$ at very low rates. One unique result reproducing experimental observations is the transition into chaotic behavior via the Feigenbaum route or period doubling kor1998b . It is also for this model that the compatibility of both microscopic simulations MC and CA has been study extensively kor1998 .

In Fig.7 we show snapshots from a simulation on Pt$`(110)`$ and a system size $`L=8192`$ using the parameter values, $`D=250`$, $`V=1`$, $`y=0.494`$, $`k=0.1`$, $`R=D`$. In the left part we plot the chemical species: CO particles are dark–grey, O particles are light–grey, and empty sites are black. The right part shows the structure of the surface: $`\alpha `$ phase sites are black, and $`\beta `$ phase sites are white. The pattern formation in this regime shows a spatio–temporal behavior, where a spiral dynamics is the dominant phenomena. It is interesting see the different structures at different spatial scales. For this purpose we include in Fig.7 sections from the upper–left corner with sizes $`4096\times 4096`$, $`1024\times 1024`$, $`256\times 256`$, from top to bottom. This sequence shows that the spiral dynamics occurs on a slowly varying island structure of sizes of the order of $`\sqrt{D/V}`$. The fact that we can see both mesoscopic and microscopic pattern formation is a quite interesting feature of the model, which has some experimental evidence win1997 ; vol1999  and has been studied theoretically in hil2002  by including lateral interactions between adsorbed particles. The model used here is simpler, because does not need that consideration in order to obtain nanostructures.

In Fig.8 we analyze the speed up of the parallel algorithm of this realistic model. We use system sizes, $`L=256,512,1024,4096,8192`$ and number of nodes $`P=1,2,4,8,16`$. The simulated time for each computation was also set to $`t_{max}=500`$. Here we can see that the speed up of the parallel algorithm is good even for small system sizes. This is because the amount of computing in each node for this model is larger than for the A$`+`$B$`0`$ model, while the amount of communication data is the same in both models.

## V Conclusions

In this paper we present a tool to obtain scaling laws connecting experimental system sizes and diffusion coefficients to standard values in microscopic MC simulations. By using a CA equivalent to MC simulation we provide an efficient parallelization algorithm. We have explained in detail how to implement the parallelization. The speed up of the algorithm is almost ideal and it is much better for larger system sizes and more complex models. A full description and analysis of the scaling laws for the second model used here is in preparation sal2002 .

###### Acknowledgements.

We thank J.J. Lukkien and S. Nedea for stimulating discussions. This work was supported by the Nederlanse Organisatie voor Wetenschapperlÿk Onderzoek (NWO), and the EC Excellence Center of Advanced Material Research and Technology (contract N 1CA1-CT-2080-7007). We thank the National Research School Combination Catalysis (NRSCC) for computational facilities.


format: text

text:
# Universally valid reformulation of the Heisenberg uncertainty principle on noise and disturbance in measurement

## I Three formulations of uncertainty principle

The uncertainty principle has been known as one of the most fundamental principles in quantum mechanics. However, there is still ambiguity in formulation; we have at least three different formulations of the “uncertainty principle”.

The Robertson uncertainty relation is generally formulated as the relation

$$\sigma (A,\psi )\sigma (B,\psi )\frac{|\psi |[A,B]|\psi |}{2}$$
(1)
for any observables $`A`$, $`B`$ and any state $`\psi `$, where the standard deviation $`\sigma (X,\psi )`$ of an observable $`X`$ in state $`\psi `$ is defined by $`\sigma (X,\psi )^2=\psi |X^2|\psi \psi |X|\psi ^2`$. This relation was proven mathematically from fundamental postulates of quantum mechanics. Nevertheless, this relation describes the limitation on preparing microscopic objects but has no direct relevance to the limitation of accuracy of measuring devices.

It is a common understanding that the uncertainty principle implies or is implied by a limitation on measuring a system without disturbing it as a position measurement typically disturbs the momentum. However, the limitation has eluded a correct quantitative expression on the trade-off between noise and disturbance. By the $`\gamma `$-ray thought experiment, Heisenberg argued that the product of the noise in a position measurement and the momentum disturbance caused by that measurement should be no less than $`\mathrm{}/2`$. This relation is generally formulated as follows: For any apparatus $`𝐀`$ to measure an observable $`A`$, the relation

$$ϵ(A,\psi ,𝐀)\eta (B,\psi ,𝐀)\frac{|\psi |[A,B]|\psi |}{2}$$
(2)
holds for any input state $`\psi `$ and any observable $`B`$, where $`ϵ(A,\psi ,𝐀)`$ stands for the noise of the $`A`$ measurement in state $`\psi `$ using apparatus $`𝐀`$ and $`\eta (B,\psi ,𝐀)`$ stands for the disturbance of $`B`$ in state $`\psi `$ caused by apparatus $`𝐀`$. We refer to the above relation as the Heisenberg noise-disturbance uncertainty relation. We shall investigate the validity of this relation to solve such questions as: When does this relation hold and when does not? What relation can be considered a universally valid generalization of this relation? How can we experimentally observe the violation of this relation?

Closely related to the above relation, the Heisenberg uncertainty relation for joint measurements is generally formulated as follows: For any apparatus $`𝐀`$ with two outputs for the joint measurement of $`A`$ and $`B`$, the relation

$$ϵ(A,\psi ,𝐀)ϵ(B,\psi ,𝐀)\frac{|\psi |[A,B]|\psi |}{2}$$
(3)
holds for any input state $`\psi `$, where $`ϵ(X,\psi ,𝐀)`$ stands for the noise of the $`X`$ measurement in state $`\psi `$ using apparatus $`𝐀`$ for $`X=A,B`$. This relation was proven mathematically under the joint unbiasedness condition requiring that the (experimental) mean values of the outcome $`𝐱`$ of the $`A`$ measurement and the outcome $`𝐲`$ of the $`B`$ measurement should coincide with the (theoretical) mean values of observables $`A`$ and $`B`$, respectively, on any input state $`\psi `$. It is a common opinion that currently available measuring devices satisfy this relation.

There is a logical relationship between the noise-disturbance relation, (2), and the joint measurement relation, (3). Suppose that the $`A`$ measurement using the apparatus $`𝐀`$ is followed immediately by a measurement of the observable $`B`$ using a noiseless measuring apparatus $`𝐁`$. Then, combining the two apparatuses, we have a new apparatus $`𝐂`$ to jointly measure $`A`$ and $`B`$ on the input state of apparatus $`𝐀`$. Since $`𝐁`$ carries out a noiseless measurement on its own input state, the noise of $`B`$ measurement in the outcome of apparatus $`𝐂`$ is equal to the $`B`$ disturbance caused by apparatus $`𝐀`$. Thus, we have the relations

$`ϵ(A,\psi ,𝐀)`$ $`=`$ $`ϵ(A,\psi ,𝐂),`$ (5)
$`\eta (B,\psi ,𝐀)`$ $`=`$ $`ϵ(B,\psi ,𝐂).`$ (6)
By applying the uncertainty relation for joint measurements to the apparatus $`𝐂`$, we can conclude that the noise-disturbance relation, (2), holds if apparatus $`𝐂`$ satisfies the joint unbiasedness condition for observables $`A`$ and $`B`$. However, there are two deficiencies of the above approach: (i) Even for noiseless measuring apparatus $`𝐀`$ to measure $`A`$ one cannot ensure that the combined apparatus $`𝐂`$ satisfies the joint unbiasedness condition. (ii) The above argument does not give a universally valid trade-off relation between noise and disturbance. Thus, we can conclude that the validity of the noise-disturbance relation, (2), cannot be reduced to the current formulation of the Heisenberg uncertainty relation for joint measurements, (3).

## II Measurement noise and disturbance

Now, let us analyze noise and disturbance in a most general description of measurement in detail. Let $`𝐀`$ be a measuring apparatus with (macroscopic) output variable $`𝐱`$ to measure an observable $`A`$ of a quantum system $`𝐒`$. Then, apparatus $`𝐀`$ measures observable $`A`$ precisely if and only if the (experimental) probability distribution of $`𝐱`$ coincides with the (theoretical) probability distribution of $`A`$; or rigorously the probability of the event that the output $`𝐱`$ of apparatus $`𝐀`$ is in an interval $`\mathrm{\Delta }`$ satisfies the Born statistical formula (BSF) for observable $`A`$, i.e.,

$$\mathrm{Pr}\{𝐱\mathrm{\Delta }\}=\psi |E^A(\mathrm{\Delta })|\psi $$
(7)
on every input state $`\psi `$, where $`E^A(\mathrm{\Delta })`$ stands for the spectral projection of $`A`$ corresponding to interval $`\mathrm{\Delta }`$. Otherwise, we consider apparatus $`𝐀`$ to measure observable $`A`$ with some noise. In order to evaluate the noise, we need to describe the measuring process. The measuring interaction is supposed to turn on at a time $`t`$, the time of measurement, and to turn off at time $`t+\mathrm{\Delta }t`$ between the object $`𝐒`$ and a quantum subsystem $`𝐏`$, called the probe, of the (large) apparatus $`𝐀`$. Denote by $`U`$ the unitary operator representing the time evolution of the composite system $`𝐒+𝐏`$ for the time interval $`(t,t+\mathrm{\Delta }t)`$. We assume that the object and any part of the apparatus do not interact before $`t`$ nor after $`t+\mathrm{\Delta }t`$. At the time of measurement the object is supposed to be in an arbitrary state $`\psi `$ and the probe is in a fixed state $`\xi `$. Thus, the composite system $`𝐒+𝐏`$ is in the state $`\psi \xi `$ at time $`t`$. Just after the measuring interaction has turned off, the probe is subjected to a precise local measurement of an observable $`M`$ of the probe, called the probe observable, and its output is recorded by the output variable $`𝐱`$. In the Heisenberg picture with the original state $`\psi \xi `$ at time $`t`$, we write $`A^{\mathrm{in}}=AI`$, $`M^{\mathrm{in}}=IM`$, $`A^{\mathrm{out}}=U^{}(AI)U`$, and $`M^{\mathrm{out}}=U^{}(IM)U`$. Since the output of this measurement is obtained by the precise measurement of observable $`M`$ at time $`t+\mathrm{\Delta }t`$, the probability distribution of the output variable $`𝐱`$ is given by

$$\mathrm{Pr}\{𝐱\mathrm{\Delta }\}=E^{M^{\mathrm{out}}}(\mathrm{\Delta }),$$
(8)
where $`\mathrm{}`$ stands for $`\psi \xi |\mathrm{}|\psi \xi `$ throughout this paper.

The noise $`ϵ(A,\psi ,𝐀)`$, or $`ϵ(A)`$ for short, of the $`A`$ measurement on input state $`\psi `$ using apparatus $`𝐀`$ is defined as the root-mean-square deviation of the experimental variable $`M^{\mathrm{out}}`$ from the theoretical variable $`A^{\mathrm{in}}`$, i.e.,

$$ϵ(A,\psi ,𝐀)=(M^{\mathrm{out}}A^{\mathrm{in}})^2^{1/2}.$$
(9)
Then, we can prove that $`ϵ(A)=0`$ on any input state $`\psi `$ if and only if apparatus $`𝐀`$ measures observable $`A`$ precisely. This ensures that “precise apparatuses” and “noiseless apparatuses” are equivalent. Let $`B`$ be an observable of $`𝐒`$.

The disturbance $`\eta (B,\psi ,𝐀)`$, or $`\eta (B)`$ for short, of observable $`B`$ on input state $`\psi `$ caused by apparatus $`𝐀`$ is defined as the root-mean-square of the change in the observable $`B`$ during the measuring interaction, i.e.,

$$\eta (B,\psi ,𝐀)=(B^{\mathrm{out}}B^{\mathrm{in}})^2^{1/2}.$$
(10)
Then, we can prove that $`\eta (B)=0`$ on any input state $`\psi `$ if and only if apparatus $`𝐀`$ does not change the probability distribution of the observable $`B`$, i.e., $`E^{B^{\mathrm{in}}}(\mathrm{\Delta })=E^{B^{\mathrm{out}}}(\mathrm{\Delta })`$ for every interval $`\mathrm{\Delta }`$ on any input state $`\psi `$. Thus, apparatuses not disturbing the observable $`B`$ and apparatuses with zero disturbance, $`\eta (B)0`$, are equivalent notions. It should be also noted that the above definitions of noise and disturbance are consistent with the standard formulation for the Heisenberg uncertainty relation for joint measurements, Eq. (3), with Eq. (I). Thus, the above definitions of noise and disturbance can be considered standard.

## III Universally valid uncertainty principle

Under the above general definitions of noise and disturbance, we can rigorously investigate the validity of the Heisenberg noise-disturbance uncertainty relation, Eq. (2). For this purpose, we introduce the noise operator $`N(A)`$ and the disturbance operator $`D(B)`$ by

$`M^{\mathrm{out}}`$ $`=`$ $`A^{\mathrm{in}}+N(A),`$ (11)
$`B^{\mathrm{out}}`$ $`=`$ $`B^{\mathrm{in}}+D(B).`$ (12)
Since $`M`$ and $`B`$ are observables in different systems, we have $`[M^{\mathrm{out}},B^{\mathrm{out}}]=0`$, and hence we have the following commutation relation for the noise operator and the disturbance operator,

$`[N(A),D(B)]+[N(A),B^{\mathrm{in}}]+[A^{\mathrm{in}},D(B)]`$ (13)
$`=[A^{\mathrm{in}},B^{\mathrm{in}}].`$ (14)
Taking the modulus of means of the both sides and applying the triangular inequality, we have

$`|[N(A),D(B)]|+|[N(A),B^{\mathrm{in}}]+[A^{\mathrm{in}},D(B)]|`$ (15)
$`|\psi |[A,B]|\psi |.`$ (16)
Since the variance is not greater than the mean square, we have

$`ϵ(A,\psi ,𝐀)`$ $``$ $`\sigma (N(A),\psi \sigma ),`$
$`\eta (B,\psi ,𝐀)`$ $``$ $`\sigma (D(B),\psi \sigma ),`$
and hence by the Heisenberg uncertainty principle, Eq. (1), we have

$$ϵ(A)\eta (B)\frac{|[N(A),D(B)]|}{2}.$$
(17)
Thus, we obtain the universally valid noise-disturbance uncertainty relation for the pair $`(A,B)`$,

$`ϵ(A)\eta (B)+{\displaystyle \frac{|[N(A),B^{\mathrm{in}}]+[A^{\mathrm{in}},D(B)]|}{2}}`$ (18)
$`{\displaystyle \frac{|\psi |[A,B]|\psi |}{2}}.`$ (19)
The above relation shows that the Heisenberg noise-disturbance uncertainty relation, Eq. (2), holds if the correlation term $`|[N(A),B^{\mathrm{in}}]+[A^{\mathrm{in}},D(B)]|`$ vanishes. In order to characterize a class of measurements satisfying Eq. (2), we introduce the following definition. The measuring interaction is said to be of independent intervention for the pair $`(A,B)`$ if the noise and the disturbance are independent of the object system; or precisely if there is observables $`N`$ and $`D`$ of probe $`𝐏`$ such that $`N(A)=1N`$ and $`D(B)=1D`$. In this case, we have $`[N(A),B^{\mathrm{in}}]=[A^{\mathrm{in}},D(B)]=0`$. Therefore, we conclude that if the measuring interaction $`U`$ is of independent intervention for the pair $`(A,B)`$, the Heisenberg noise-disturbance uncertainty relation for the pair $`(A,B)`$, Eq. (2), holds for any object state $`\psi `$ and any probe state $`\xi `$. The above conclusion was previously suggested in part by Braginsky and Khalili with a limited justification and now fully justified.

The universally valid uncertainty relation shows that for measurements of dependent intervention, the lower bound of the noise-disturbance product depends on the pre-measurement uncertainties (standard deviations) of $`A`$ and $`B`$. In order to obtain the trade-off among the noise $`ϵ(A)`$, the disturbance $`\eta (B)`$, and the pre-measurement uncertainties $`\sigma (A)`$ and $`\sigma (B)`$, we apply the Heisenberg uncertainty principle, Eq. (1), to all terms in the left-hand-side of the universally valid noise-disturbance uncertainty relation, Eq. (18). Then, we now obtain the generalized noise-disturbance uncertainty relation,

$$ϵ(A)\eta (B)+ϵ(A)\sigma (B)+\sigma (A)\eta (B)\frac{|\psi |[A,B]|\psi |}{2}.$$
(20)
The above relation holds for any apparatus $`𝐀`$ (specified by any probe state $`\xi `$, any unitary operator $`U`$, and any probe observable $`M`$), any observables $`A,B`$, and any input state $`\psi `$, and hence ultimately generalizes the Heisenberg noise-disturbance uncertainty relation, Eq. (2), to arbitrary measurements.

Under the finite energy constraint, i.e., $`\sigma (Q),\sigma (P)<\mathrm{}`$, the above relation excludes the possibility of having both $`ϵ(Q)=0`$ and $`\eta (P)=0`$. However, $`ϵ(Q)=0`$ is possible with $`\sigma (Q)\eta (P)\mathrm{}/2`$ and $`\eta (P)=0`$ is possible with $`ϵ(Q)\sigma (P)\mathrm{}/2`$. In particular, even the case where $`ϵ(Q)=0`$ and $`\eta (P)<\epsilon `$ with arbitrarily small $`\epsilon `$ is possible for some input state with $`\sigma (Q)>\mathrm{}/2\epsilon `$, and also the case where $`\eta (P)=0`$ and $`ϵ(Q)<\epsilon `$ is possible for some input state with $`\sigma (P)>\mathrm{}/2\epsilon `$. Such extreme cases occur in compensation for large uncertainties in the input state, while in the minimum uncertainty state with $`\sigma (Q)=\sigma (P)=(\mathrm{}/2)^{1/2}`$, we have

$$ϵ(Q)\eta (P)+\sqrt{\frac{\mathrm{}}{2}}[ϵ(Q)+\eta (P)]\frac{\mathrm{}}{2}.$$
(21)
Even in this case, it is allowed to have $`ϵ(Q)\eta (P)=0`$ with $`ϵ(Q)=0`$ and $`\eta (P)(\mathrm{}/2)^{1/2}`$ or with $`\eta (P)=0`$ and $`ϵ(Q)(\mathrm{}/2)^{1/2}`$.

## IV Violation of the Heisenberg inequality

Now let us consider the problem as to whether one can implement, under the current experimental technique, a good measuring apparatus with small noise-disturbance product beyond the original Heisenberg lower bound. The controversy on the sensitivity limit of gravitational wave detectors suggested that the Heisenberg noise-disturbance uncertainty relation is not universally valid. In fact, based on the Heisenberg noise-disturbance uncertainty relation, Braginsky, Caves, and others claimed that there is a sensitivity limit, called the standard quantum limit (SQL), for monitoring the free-mass position that leads to a quantum mechanical sensitivity limit on interferometer type gravitational wave detectors. However, Yuen proposed the idea of “contractive state measurements” to break the SQL and Ozawa found an explicit Hamiltonian realization of a contractive state measurement that breaks the SQL. Consequently, the above measuring interaction violates the Heisenberg noise-disturbance uncertainty relation. Direct computaions on the position-measuring-niose and momentum-disturbance was given in  showing the violation of the Heisenberg noise-disturbance uncertainty relation.

In what follows, modifying the above interaction in the context of optical quadrature measurement, I will show that the small noise-disturbance product can be achieved beyond Heisenberg’s lower bound by an apparatus carrying out a precise and nearly non-disturbing quadrature measurement with currently available linear optical devices.

## V Backaction evading quadrature amplifiers

Consider the case where the system $`𝐒`$ and the probe $`𝐏`$ are two optical modes represented by annihilation operators $`a`$ and $`b`$, respectively. The quadrature component field operators $`X_a,Y_a,X_b,Y_b`$ are self-adjoint operators satisfying $`a=X_a+iY_a`$ and $`b=X_b+iY_b`$, for which we have $`[X_a,Y_a]=i/2`$ and $`[X_b,Y_b]=i/2`$.

A measuring interaction $`U`$ on $`𝐒+𝐏`$ is called a back-action evading (BAE) quadrature amplifier with gain $`G`$ if we have

$`X_a^{\mathrm{out}}`$ $`=`$ $`X_a^{\mathrm{in}},`$ (23)
$`X_b^{\mathrm{out}}`$ $`=`$ $`X_b^{\mathrm{in}}+GX_a^{\mathrm{in}},`$ (24)
$`Y_a^{\mathrm{out}}`$ $`=`$ $`Y_a^{\mathrm{in}}GY_b^{\mathrm{in}},`$ (25)
$`Y_b^{\mathrm{out}}`$ $`=`$ $`Y_b^{\mathrm{in}}.`$ (26)
In order to measure $`X_a`$, the probe observable is chosen best to be $`M=X_b/G`$. Then we have

$$M^{\mathrm{out}}=X_a^{\mathrm{in}}+\frac{1}{G}X_b^{\mathrm{in}}.$$
(27)
The $`X_a`$-noise operator, the $`X_a`$-disturbance operator, and the $`Y_a`$-disturbance operator are given by

$`N(X_a)`$ $`=`$ $`{\displaystyle \frac{1}{G}}X_b^{\mathrm{in}},`$ (29)
$`D(X_a)`$ $`=`$ $`0,`$ (30)
$`D(Y_a)`$ $`=`$ $`GY_b^{\mathrm{in}}.`$ (31)
The condition $`D(X_a)=0`$ is characteristic of BAE amplifiers. If the probe is prepared nearly in the $`X_b`$ eigenstate $`|X_b=0`$, the measurement is a nearly noiseless ($`ϵ(X_a)0`$) and nondisturbing ($`\eta (X_a)=0`$) measurement of $`X_a`$. From the above, BAE amplifiers are of independent intervention for the pair $`(X_a,Y_a)`$. Thus, the Heisenberg noise-disturbance uncertainty relation for the pair $`(X_a,Y_a)`$ holds. If the $`b`$-mode is prepared in the vacuum, $`\xi =|0`$, the noise and disturbance satisfy $`ϵ(X_a)=1/2`$ and $`\eta (Y_a)=1/2`$, so that the minimum noise-disturbance product attains as $`ϵ(X_a)\eta (Y_a)=1/4`$.

## VI Noiseless quadrature transducers

Consider the following input-output relations

$`X_a^{\mathrm{out}}`$ $`=`$ $`X_a^{\mathrm{in}}X_b^{\mathrm{in}},`$ (33)
$`X_b^{\mathrm{out}}`$ $`=`$ $`X_a^{\mathrm{in}},`$ (34)
$`Y_a^{\mathrm{out}}`$ $`=`$ $`Y_b^{\mathrm{in}},`$ (35)
$`Y_b^{\mathrm{out}}`$ $`=`$ $`Y_b^{\mathrm{in}}+Y_a^{\mathrm{in}}.`$ (36)
In this case, the measuring interaction $`U`$ is called the noiseless quadrature transducer. In order to measure $`X_a`$, the probe observable is chosen to be $`M=X_b`$. The $`X_a`$-noise operator, the $`X_a`$-disturbance operator, and the $`Y_a`$-disturbance operator are given by

$`N(X_a)`$ $`=`$ $`0,`$ (38)
$`D(X_a)`$ $`=`$ $`X_b^{\mathrm{in}},`$ (39)
$`D(Y_a)`$ $`=`$ $`(Y_a^{\mathrm{in}}+Y_b^{\mathrm{in}}).`$ (40)
The condition $`N(X_a)=0`$ is characteristic of the noiseless transducer. Hence, the measurement is a noiseless $`X_a`$-measurement regardless of the probe preparation $`\xi `$. From (39), if the probe is prepared nearly in the $`X_b`$ eigenstate $`|X_b=0`$, the measurement is a noiseless ($`ϵ(X_a)=0`$) and nearly nondisturbing ($`\eta (X_a)0`$) measurement of $`X_a`$. Since $`ϵ(X_a)=0`$, we have

$$ϵ(X_a)\eta (Y_a)=0$$
(41)
for any states $`\psi `$ and $`\xi `$, so that Heisenberg’s lower bound for the noise-disturbance product can be overcome by a noiseless and nearly nondisturbing quadrature measurement, if one can implement a noiseless quadrature transducer. The above model also suggests that the linearity of measuring interaction does not ensure the validity of the Heisenberg noise-disturbance uncertainty relation, despite a prevailing claim that linear measurements, measurements closely connected to linear systems, obey the Heisenberg noise-disturbance uncertainty relation.

## VII Experimental realization of the noiseless quadrature transducers

The noiseless quadrature transducer can be implemented, in principle, as follows. Consider two degenerate modes $`a`$, the signal, and $`b`$, the probe, with frequency $`\omega `$ and orthogonal polarization, which undergo successive parametric interactions in the following 5 steps; see Refs.  for similar implementations of BAE amplifiers.

(i) The two polarization modes undergo a mixing interaction using a polarization rotator which rotates the angle of polarization by $`\theta `$. The operation of the polarization rotator is represented by the mixing operator

$$T(\theta )=\mathrm{exp}[\theta (ab^{}a^{}b)].$$
(42)
(ii) The mixture of the signal and the probe fields will propagate along each of the ordinary and orthogonal extraordinary axis of a KTP crystal pumped by a pulsed intense classical field. This interaction is a nondegenerate parametric amplifier described by the two-mode squeeze operator

$$S(r)=\mathrm{exp}[r(aba^{}b^{})],$$
(43)
where $`r`$ corresponds to the squeezing parameter.

(iii) After the amplification step the fields pass through a second polarization rotator with the mixing angle $`2\theta `$ so that the operation is represented by $`T(2\theta )=T(\theta )^2`$.

(iv) The mixture of the signal and the probe fields undergoes the second parametric amplification described by the two-mode squeeze operator $`S(r)`$.

(v) After the second amplification step the field pass through a third polarization rotator with the mixing angle $`\theta `$ so that the operation is represented by $`T(\theta )`$.

Thus this process defines the measuring interaction operator

$$U(r,\theta )=T(\theta )S(r)T(\theta )T(\theta )S(r)T(\theta ).$$
(44)
We shall determine the parameters $`\theta `$ and $`r`$ for $`U(r,\theta )`$ to realize a noiseless quadrature transducer, Eq. (VI). Suppose that $`\theta `$ and $`r`$ satisfies the relation $`\mathrm{sin}2\theta =\mathrm{tanh}r`$. In this case, it is well-known that the unitary operator $`U_{}=T(\theta )S(r)T(\theta )`$ realizes the BAE quadrature amplifier with $`G=2\mathrm{sinh}r`$, i.e.,

$`U_{}^{}X_aU_{}`$ $`=`$ $`X_a,`$ (46)
$`U_{}^{}X_bU_{}`$ $`=`$ $`X_b+2(\mathrm{sinh}r)X_a,`$ (47)
$`U_{}^{}Y_aU_{}`$ $`=`$ $`Y_a2(\mathrm{sinh}r)Y_b,`$ (48)
$`U_{}^{}Y_bU_{}`$ $`=`$ $`Y_b.`$ (49)
Similarly, the unitary operator $`U_+=T(\theta )S(r)T(\theta )`$ realizes the conjugate BAE quadrature amplifier, i.e.,

$`U_+^{}X_aU_+`$ $`=`$ $`X_a2(\mathrm{sinh}r)X_b,`$ (51)
$`U_+^{}X_bU_+`$ $`=`$ $`X_b,`$ (52)
$`U_+^{}Y_aU_+`$ $`=`$ $`Y_a,`$ (53)
$`U_+^{}Y_bU_+`$ $`=`$ $`Y_b+2(\mathrm{sinh}r)Y_a.`$ (54)
Combining the above equations, we have the input-output relations for the unitary operator $`U(r,\theta )=U_{}U_+`$,

$`U(r,\theta )^{}X_aU(r,\theta )`$ $`=`$ $`X_a2(\mathrm{sinh}r)X_b,`$ (56)
$`U(r,\theta )^{}X_bU(r,\theta )`$ $`=`$ $`(14\mathrm{sinh}^2r)X_b+2(\mathrm{sinh}r)X_a,`$ (57)
$`U(r,\theta )^{}Y_aU(r,\theta )`$ $`=`$ $`(14\mathrm{sinh}^2r)Y_a2(\mathrm{sinh}r)Y_b,`$ (58)
$`U(r,\theta )^{}Y_bU(r,\theta )`$ $`=`$ $`Y_b+2(\mathrm{sinh}r)Y_a.`$ (59)
Thus, if $`\mathrm{sinh}r=1/2`$, i.e.,

$`r`$ $`=`$ $`\mathrm{log}{\displaystyle \frac{1+\sqrt{5}}{2}}`$ (61)
$`\theta `$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}^1{\displaystyle \frac{1}{\sqrt{5}}},`$ (62)
the resulting unitary operator $`U=U(r,\theta )`$ realizes the noiseless quadrature transducer, Eq. (VI).

Thus, if $`r`$ and $`\theta `$ are chosen as Eq. (VII), the parametric process (i)–(v) realizes the noiseless quadrature transducer. Therefore, this process followed immediately by the homodyne detection of the $`X_b`$ component implements a nearly nondisturbing and noiseless measurement of the $`X_a`$ quadrature component that disturbs the conjugate observable $`Y_a`$ much less than the quantum limit set by Heisenberg’s lower bound for the noise-disturbance product, $`ϵ(X_a)\eta (Y_a)1/4`$.

## VIII Conclusion

In this paper, I have proposed new relations, (18) and (20), universally valid for the trade-off between the measurement noise and disturbance. These relations demonstrates that the prevailing Heisenberg’s lower bound for the noise-disturbance product is valid for measurements with independent intervention, but can be circumvented by a measurement with dependent intervention. An experimental confirmation of the violation of Heisenberg’s lower bound is proposed for a measurement of optical quadrature with currently available techniques in quantum optics. The new relation will not only bring a new insight on fundamental limitations on measurements set by quantum mechanics but also advance a frontier of precision measurement technology such as gravitational wave detection and quantum information processing.

## ACKNOWLEDGMENTS

Acknowledgments. The author thanks F. De Martini, S. Takeuchi, and I. Ojima for useful discussions. This work was supported by the programme “R&D on Quantum Communication Technology” of the MPHPT of Japan, by the CREST project of the JST, and by the Grant-in-Aid for Scientific Research of the JSPS.


text:
# CONSTRAINING THE GALACTIC BAR PARAMETERS WITH RED CLUMP GIANTS

## 1 INTRODUCTION

Many observational groups carried/are carrying out microlensing observations (e.g., MACHO: Alcock et al. 1993; OGLE: Udalski et al. 1993; DUO: Alard & Guibert 1997; MOA: Bond et al. 2001; EROS: Aubourg et al. 1993). Over one thousand microlensing events in the local group have been identified (e.g., Alcock et al. 2000; Woźniak et al. 2001; Derue et al. 2001; Alard & Guibert 1997; Bond et al. 2001). In the coming years, $``$ 1000 microlensing events are expected to be discovered by the OGLE III111 http://www.astrouw.edu.pl/~ogle/ogle3/ews/ews.html  and other collaborations every year, many of them in real-time.

Data collected from microlensing experiments have very diverse applications (for reviews see Paczyński 1996; Gould 1996). One of the most important applications is studying the Galactic structure. There is strong evidence that the centre of Galaxy hosts a bar (e.g., de Vaucouleurs 1964; Blitz & Spergel 1991; Stanek et al. 1994, 1997; Kiraga & Paczyński 1994; Häfner et al. 2000 and references therein). However, the parameters of the bar are not well-determined, including its mass, size, and the motion of stars within it.

Data from microlensing experiments provide several ways of probing the structure of the inner Galaxy. For example, the optical depth is roughly proportional to the total mass of the bar while the event time scale distribution probes the mass function and kinematics of the bar. Another important method that was first explored by Stanek et al. (1997) uses red clump giants (RCGs). These bulge stars occupy a distinct region in the colour-magnitude diagram (see, e.g., Stanek et al. 2000 and references therein). They have very small intrinsic widths in their luminosity function, about 0.2 mag for RCGs in the Galactic bulge (see §2 in Stanek et al. 1997; Paczyński & Stanek 1998). The observed luminosity function has a well-defined peak and an apparent width (see Fig. 5 in Stanek et al. 1997). The apparent width depends on both the intrinsic width and the spread caused by the radial depth of the bar. Stanek et al. (1997) used this dependence to constrain the Galactic bar axial ratios and orientation.

In this paper, we use the RCGs to study the kinematics of stars in the Galactic bar. Stars on the front side and far side may have different (tangential and radial) streaming motions. We explore two effects to constrain such motions. The first effect is based on the realisation that red clump stars at the bright slope of the peak of the luminosity function must be (on average) closer to us, i.e. on the near (front) side, while the stars on the faint slope must be (on average) on the far (back) side. Therefore, one can select bright and faint samples of RCGs and examine their differences in the radial velocity and proper motion. The second effect we explore is based on the well-known fact that lensed bulge stars lie preferentially on the far side of the Galactic bar due to the enhanced lensing probability by the stars on the near side. This leads to several observable effects. For example, the lensed RCGs should be fainter compared with all observed red clump stars in the field (Stanek 1995; see also Zhao 1999a, 1999b, 2000 for other effects such as extinction and reddening for microlensing events toward the Large Magellanic Cloud). We will examine systematic differences in the proper motions and radial velocities of the lensed RCGs relative to all observed RCGs in the Galactic bar; a similar bias in the radial velocity for lensing events toward the Large Magellanic Cloud has been discussed by Zhao (1999b).

As this paper is the first feasibility study, we do not attempt to construct a detailed (or even self-consistent) model of the Galactic bulge. Our aim is to provide an order-of-magnitude estimate of the kinematic biases in the hope of motivating observations to be carried out, which may in turn provide strong constraints on the the model of the inner Galaxy and incentives for their refinement. The structure of the paper is as follows. In §2, we outline our model, and in §3 we present our main results. We discuss observational issues to detect these effects in §4. For definiteness, we shall restrict ourselves to studying the red clump stars in Baade’s window ($`l=1^{}`$, $`b=3.9^{}`$).

## 2 MODEL

To explore both effects that we discussed above for RCGs, we need to adopt a bulge model, which is rather uncertain; there are a number of published analytical and numerical models of various complexities (see Han & Gould 1995 for a large selection of analytical models; see Zhao et al. 1996, Häfner et al. 2000, Fux 2001 and references therein for numerical models). For illustrative purposes, we will take two simple analytical models in the literature for the Galactic bulge.

The first toy model we adopt is identical to the one used by Kiraga & Paczyński (1994). The bulge mass density distribution is approximated as axis-symmetric and is based on the study of Kent (1992). Note that there is no bar in this simple model. The bulge velocity distribution is assumed to be isotropic. Each of the three velocity components follows a Gaussian distribution with a dispersion of $`100\mathrm{km}\mathrm{s}^1`$. In addition, the bulge stars on the near side rotate with a velocity of $`+100\mathrm{km}\mathrm{s}^1`$ (in the same direction as the solar rotation) while those on the far side rotate in the opposite sense (see Kiraga & Paczyński for details). In this model, the streaming motion is purely tangential. The second model we consider is somewhat more realistic. It is based on the E2 bar model of Dwek et al. (1995). For the following discussions, it is convenient to establish a coordinate system centred on the bar with $`x,y`$ and $`z`$ corresponding to the major, intermediate and minor axes. The minor axis coincides with the vertical direction while the major axis has a viewing angle of $`\theta _{\mathrm{bar}}=24^{}`$ from the line of sight (for an illustration of the bar geometry, see Figs. 6 and 11 in Stanek et al. 1997). The axial ratios are taken to be $`1:0.42:0.28`$. The luminosity-weighted velocity dispersions can be derived using the virial theorem (Belokurov & Evans 2001) as $`\sigma _x=100\mathrm{km}\mathrm{s}^1`$, $`\sigma _y=80\mathrm{km}\mathrm{s}^1`$ and $`\sigma _z=68\mathrm{km}\mathrm{s}^1`$, respectively. In addition, the stars on the near side have a streaming motion of $`+50\mathrm{km}\mathrm{s}^1`$ parallel to the major axis, while the stars on the far sides move in the opposite sense. Note that the (rotational) sense of the streaming motion is similar to that adopted by Kiraga & Paczyński (1994), but a factor of two smaller in value.

The intrinsic luminosity function of red clump stars is assumed to follow a Gaussian distribution with a width $`\sigma _00.2\mathrm{mag}`$ (Stanek et al. 1997); there is presumably no difference in the luminosity function for RCGs on the near and far sides. We assume that RCGs follow the same bulge mass distribution over a distance range from $`D_{\mathrm{min}}<D_\mathrm{s}<D_{\mathrm{max}}`$; we take $`D_{\mathrm{min}}=0.5R_0`$ and $`D_{\mathrm{max}}=1.5R_0`$, where $`R_0=8\mathrm{kpc}`$ is the distance to the Galactic centre. Other choices of $`D_{\mathrm{min}}`$ and $`D_{\mathrm{max}}`$ only change our results slightly. Note that the effect of extinction and reddening can be minimised by using the extinction-insensitive Wesenheit index, $`W_{VI}`$ (written as $`V_{VI}`$ in Paczyński et al. 1994).

RCGs are usually selected in the colour-magnitude diagram by requiring their colours and magnitudes to be in a certain range (e.g., Fig. 1 in Paczyński & Stanek 1998). Many red giants occupy the same region in the colour-magnitude diagram. As the luminosities of red giants span a broad range they can be at either on the far side or on the near side at any given (apparent) magnitude. Therefore, they dilute the kinematic differences between a bright RCG sample and a faint one as their apparent brightness can not be used as an approximate measure of distance. We model the number counts of red giants as a linear function of the apparent magnitude (see eq. 1 in Paczyński & Stanek 1998). According to the model, within one magnitude of the RCG peak, red giants contribute about 60% of the total number. Therefore, red giants substantially dilute the offset signals if they cannot be separated from the red clump stars; we return to this issue in the discussion.

For the offset in kinematics between lensed and all observed RCGs, we need to calculate the microlensing optical depth and event rate. We consider both bulge self-lensing and lensing by the disc. Therefore, we need to specify the disc mass distribution and kinematics. As the offsets in the radial velocity and proper motion depend only on the ratio of the lensing frequency for stars on the near side and those on the far side, they do not depend on the mass function as the mass cancels out in the ratio.

For the disc mass distribution, we adopt the familiar double exponential form (Bahcall 1986)

$$\rho (r,z)=\rho _{}\mathrm{exp}\left(\frac{rR_0}{3.5\mathrm{kpc}}\frac{|z|}{0.325\mathrm{kpc}}\right),$$
(1)
where $`(r,\varphi ,z)`$ are the usual cylindrical coordinates (in units of kpc), and $`\rho _{}`$ is the local disc density, for which we take a value of $`0.06M_{}\mathrm{pc}^3`$. The velocity dispersions are taken to be $`\sigma _r=30\mathrm{km}\mathrm{s}^1,\sigma _\varphi =20\mathrm{km}\mathrm{s}^1`$, and $`\sigma _z=20\mathrm{km}\mathrm{s}^1`$, which are typical values in all disc models (e.g., Belokurov & Evans 2002). The disc lenses have a systematic rotation of $`v_\varphi =+220\mathrm{km}\mathrm{s}^1`$ at all distances. In addition, we also account for the solar motion relative to the local standard of rest (see e.g., Binney & Merrifield 1998).

Obviously, the optical depth and event rate need to be weighted properly by the source distance (see eqs. 5 and 11 in Kiraga & Paczyński 1994). For the lensing of RCGs, we assume that they are bright enough to be seen with equal probability at all distances (corresponding to $`\beta =0`$ in the formalism of Kiraga & Paczyński 1994).

## 3 RESULTS

Monte Carlo simulations are performed to obtain our results; this is the easiest way to incorporate different source populations and selection strategies. In the following, we first discuss the offset in the radial velocity, $`V_\mathrm{r}`$, and in the longitudinal proper motion, $`\mu _\mathrm{l}`$, between the faint and bright samples of RCGs, and then discuss the offsets between the lensed and all observed RCGs.

### 3.1 Offsets in $`V_\mathrm{r}`$ and $`\mu _\mathrm{l}`$ Between Faint and Bright Sub-samples of RCGs

As we discussed previously, the observed luminosity function of RCGs has a well defined peak in their number counts (see Fig. 5 in Stanek et al. 1997). The RCGs on the brighter side of the peak are (on average) on the front side of the bar and vice versa. We experimented with different ways of selecting ‘bright’ and ‘faint’ samples in order to maximise their differences in the radial velocity and proper motion. We found that the optimum ‘bright’ sample is formed by stars that are brighter than the peak magnitude by 0.2 to 0.4 magnitude. Similarly, the optimum ‘faint’ sample is formed by stars fainter than the peak magnitude by 0.2 to 0.4 magnitude. Notice that each sample also includes a background of red giant stars.

The shift in the radial velocity between the faint and bright samples depends on the model. In the Kiraga & Paczyński model, there is no radial streaming of stars, hence there is no shift in the radial velocity as the stars on the near and far sides have the same random kinematics. This is different for the Dwek et al. model, where we have streaming motions of $`\pm 50\mathrm{km}\mathrm{s}^1`$ mostly along the radial direction as the bar is (by assumption) only $`24^{}`$ from the line of sight.

The left panel in Fig. 1 shows the distributions of radial velocity for the bright (dashed) and faint (solid) samples, respectively. Both distributions are roughly Gaussian with a dispersion of $`110\mathrm{km}\mathrm{s}^1`$. The faint sample is shifted toward more negative radial velocities, with an average shift of $`33\mathrm{km}\mathrm{s}^1`$. This shift can be understood qualitatively. The stars on the near side have an average radial velocity $`(10+50\mathrm{cos}\theta _{\mathrm{bar}})\mathrm{km}\mathrm{s}^1`$, where $`10\mathrm{km}\mathrm{s}^1`$ arises due to the solar motion ($`+10\mathrm{km}\mathrm{s}^1`$) relative to the local standard of rest. Analogously, stars in the faint sample are mostly on the far side, and they have an average radial velocity $`(1050\mathrm{cos}\theta _{\mathrm{bar}})\mathrm{km}\mathrm{s}^1`$. As the faint sample includes more stars on the far side, hence their radial velocities are more negative. The average shift between the bright and faint samples is about one third of the velocity dispersion along the line of sight ($`100\mathrm{km}\mathrm{s}^1`$). We return to the feasibility of observing this shift in the discussion.

As the streaming motion we considered is in the plane of the Galaxy, we do not expect any shift in the latitudinal proper motion, so we will only consider the shift in the longitudinal direction, $`\mu _\mathrm{l}`$. The shift in $`\mu _\mathrm{l}`$ between the bright and faint samples again depends on the model. In the second model, the streaming motion is mostly radial, so the proper motion shift is relatively small. However, for the Kiraga & Paczyński model, there is a tangential streaming motion of $`+100\mathrm{km}\mathrm{s}^1`$ for stars on the near side and $`100\mathrm{km}\mathrm{s}^1`$ on the far side, and so the offset becomes more apparent. The left panel in Fig. 2 shows the distributions of proper motion for the bright (dashed) and faint (solid) samples, respectively. The proper motions for the faint sample is shifted toward a more negative value, by an average amount of $`1.6\mathrm{mas}\mathrm{yr}^1`$; the dispersion for the faint sample is also about 10% smaller. To understand these trends, let us recall that for a star with velocity $`v_l`$ (in units of $`\mathrm{km}\mathrm{s}^1`$) and at distance $`D_\mathrm{s}`$ (in units of kpc) has a proper motion of $`\mu _\mathrm{l}0.2v_l/D_\mathrm{s}\mathrm{mas}\mathrm{yr}^1`$. In the heliocentric frame, stars on the far side move with a streaming motion of $`v_l=300\mathrm{km}\mathrm{s}^1`$ while stars on the near side move with $`v_l=100\mathrm{km}\mathrm{s}^1`$. As the faint sample includes more stars from the far side, their proper motions are therefore more negative relative to the bright sample. The average shift is roughly proportional to the value of the tangential streaming motion; it is about one half of the dispersion in the proper motion of approximately $`3\mathrm{mas}\mathrm{yr}^1`$ for all observed stars. The dispersion of the faint sample is smaller because $`\mu _\mathrm{l}D_\mathrm{s}^1`$ and faint RCGs are on average at slightly larger distances.

### 3.2 Offsets in $`V_\mathrm{r}`$ and $`\mu _\mathrm{l}`$ Between Lensed and All Observed RCGs

The shift in the radial velocity between the lensed and all observed RCGs depends on the models. In the first model where there is no radial streaming of stars, there is no shift in the radial velocity as the stars on the near and far sides have the same random kinematics. This is different for the second model, where we do have radial streaming motions.

The right panel in Fig. 1 shows the distributions of radial velocity for the lensed (solid) and all observed (dashed) RCGs, respectively. The curve predicted for all observed stars has a negative mean value around $`10\mathrm{km}\mathrm{s}^1`$ partly due to the solar motion with respect to the local standard of rest. However, in this model, the stars on the far side move with $`50\mathrm{cos}\theta _{\mathrm{bar}}\mathrm{km}\mathrm{s}^1`$ along the line of sight and those on the near side move in the opposite direction. As the optical depth and lensing event rate are higher for stars on the far side, radial velocities of the lensed sources are shifted toward more negative values. The difference in the average radial velocity is about $`10\mathrm{km}\mathrm{s}^1`$. This is smaller than the offset between the bright and faint samples that we discussed in §3.1. We found that this shift is roughly proportional to the value of the radial streaming motion. In any case, the average shift is fairly small compared with the velocity dispersion along the line of sight ($`100\mathrm{km}\mathrm{s}^1`$).

The shift in the proper motion between the lensed and all observed RCGs again depends on the model. In the second model, the streaming motion is mostly radial, so the shift in proper motion is relatively small. However, for the first model, where there is a tangential streaming motion of $`+100\mathrm{km}\mathrm{s}^1`$ for stars on the near side and $`100\mathrm{km}\mathrm{s}^1`$ on the far side, the offset becomes noticeable. The right panel in Fig. 2 shows the distributions of proper motion for the lensed stars (solid) and all observed stars (dashed), respectively. The curve for the overall population has an average proper motion of $`5.8\mathrm{mas}\mathrm{yr}^1`$ due to the solar motion around the Galactic centre. The curve for the lensed population is shifted toward more negative proper motions. This is because stars on the far side move opposite to the solar rotation, and as they have larger optical depths and are more frequently lensed, their proper motions are shifted toward more negative longitudinal values. The average shift in the proper motion is about $`1.2\mathrm{mas}\mathrm{yr}^1`$, which is smaller than its dispersion, about $`3\mathrm{mas}\mathrm{yr}^1`$. As expected, this offset is roughly proportional to the value of the tangential streaming motion.

## 4 DISCUSSION

In this paper we have studied how to use red clump giants (RCGs) to constrain the kinematics of stars in the Galactic bar. The basic idea is that stars on the near and far sides may have different streaming motions, and if we can select two samples with one preferentially on the near side and the other on the far side, then they should show differences in the radial and tangential streaming motions. We have examined two ways that we can select such samples. The first method selects a bright sample and a faint one with the peak magnitude of the number counts of RCGs as the approximate dividing line. The RCGs in the faint sample is (on average) more on the far side of the bar and vice versa. The second method relies on the fact that stars on the far side are preferentially microlensed due to the enhanced lensing probability by the stars on the near side and in the disc. Hence the lensed stars and all observed stars should behave differently in kinematics.

We illustrate the kinematical differences in the context of two simple models of the Galactic bulge. In the first (axis-symmetric) model, the streaming motion is purely tangential, while in the second (more realistic) bar model the streaming motion is more or less radial. In reality, the streaming motions in the Galactic bar may have both a radial component and a tangential one, so our numbers should be taken as rough order-of-magnitude estimates.

We found at the difference in the radial velocity between the bright and faint samples can reach $`\delta V_\mathrm{r}33\mathrm{km}\mathrm{s}^1`$ for a radial streaming motion of $`50\mathrm{km}\mathrm{s}^1`$. The difference is substantially reduced by the population of red giants that occupy the same part of the colour-magnitude diagram as the RCGs, as they have a large span in luminosity, so a brighter magnitude does not signal that a star is closer to us on average. If there is an effective way of differentiating red giants and red clump stars, then the differences in the radial velocity and proper motion can be much larger. This can be, for example, done by examining colours as RCGs occupy a slightly bluer part of the colour-magnitude diagram. Spectroscopic observations may also ultimately provide features that distinguish red giants and RCGs. For example, if the contribution of red giants can be reduced by one half, then the shift in the proper motions can be as high as $`2.3\mathrm{mas}\mathrm{yr}^1`$, while that for the radial velocity can reach $`50\mathrm{km}\mathrm{s}^1`$. In comparison, we find that the shift between the lensed stars and overall population of stars is more modest, about $`10\mathrm{km}\mathrm{s}^1`$ for a radial streaming motion of $`50\mathrm{km}\mathrm{s}^1`$.

The velocity dispersion along the line of sight is about $`\sigma 100\mathrm{km}\mathrm{s}^1`$, so in order to see this shift at the $`2\sigma `$ level, from Poisson statistics, we only need to obtain the radial velocity for about $`(2\times \sigma /\delta V_\mathrm{r})^2100`$ stars even for $`\delta V_\mathrm{r}=20\mathrm{km}\mathrm{s}^1`$. The required number is moderate, but is within reach of multi-object spectrographs available on many large telescopes. For example, the FORS1 instrument on VLT has 19 slits that can be efficiently used to derive the radial velocities of lensed stars and field stars simultaneously. Other instruments such as 2dF 222http://www.aao.gov.au/2df/ which can take up to 400 spectra simultaneously should also be explored. We note that for RCGs, blending is not a severe problem as they are so bright, so their radial velocities can be measured without much difficulty. In any case, we only require an accuracy of $`1020\mathrm{km}\mathrm{s}^1`$ per star, since the dispersion in the radial velocity for the bulge stars is of the order of $`100\mathrm{km}\mathrm{s}^1`$.

For the shift in the proper motions, we find that a value of $`1.6\mathrm{mas}\mathrm{yr}^1`$ between the bright and faint samples of RCGs for a tangential streaming motion $`\overline{V_l}=100\mathrm{km}\mathrm{s}^1`$. The shift scales roughly linearly with $`\overline{V_l}`$. The difference is about $`1\mathrm{mas}\mathrm{yr}^1`$ between the lensed and overall population of bright stars. These shifts are not much smaller than the dispersion in the proper motion, approximately $`3\mathrm{mas}\mathrm{yr}^1`$. The shifts can be detected with HST or even ground-based instruments. To do this, only a modest number of RCGs need to be monitored to measure the relative shifts in proper motions. The relative proper motion may be feasible to detect even from the ground. Soszyński et al. (2002) demonstrated that proper motions as small as $`4\mathrm{mas}\mathrm{yr}^1`$ can be measured with OGLE II data. This accuracy is achieved through a combination of the new difference image analysis software (Eyer & Woźniak 2001) and a large number of data points for any given star, which reduces the astrometric error through Poisson statistics. The camera system of OGLE III yields much better images than that of OGLE II as the latter operates in the drift scan mode and produced somewhat worse point spread functions. Therefore it seems possible (although still to be demonstrated) that proper motions of many red clump stars, including lensed ones, can be obtained using data from microlensing surveys (e.g., OGLE III), requiring no extra observing resources.

In this paper we have only explored the shifts in the radial velocity and proper motion in Baade’s window. Clearly there must be a spatial dependence of these shifts. The OGLE III collaboration is currently monitoring many fields around the Galactic centre, covering from $`7^{}`$ to $`6^{}`$ in latitude and $`12^{}`$ to $`12^{}`$ in longitude 333For detailed field coordinates, see http://www.astrouw.edu.pl/~ogle/ogle3/ews/gb\_ews.jpg. Red clump stars and microlensing events will become available over large areas, and so it will be very interesting to detect the spatial variation of the shifts in the radial velocity and proper motion. Such observations will likely provide strong constraints on the Galactic bar parameters, and stimulate further theoretical efforts to build a better and self-consistent model of the inner Galaxy. Finally, we point out that if future astrometric missions, SIM and GAIA444http://sim.jpl.nasa.gov; http://astro.estec.esa.nl/GAIA/, perform as planned, they will directly measure the radial depth of the galactic bar by determining parallax’s for a number of stars. This will make our suggestion obsolete, but not sooner than in a decade.

We thank Martin Smith, David Spergel, and Hongsheng Zhao for helpful discussions. S.M. acknowledges travel support by Princeton University. This project was supported by the NSF grant AST-1206213, and the NASA grant NAG5-12212 and funds for proposal #09518 provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.