text:
# New light on the baryon fraction in galaxy clusters

## 1 Introduction

One of the primary goal of cosmology is the determination of the mean density of the universe. Theoretical arguments such as the inflation predict a flat universe with $`\mathrm{\Omega }_0=1`$. However, most of the observational data do not seem to provide a strong evidence for a critical universe.
Rich clusters of galaxies are the largest virialized systems in the universe, they therefore provide us with useful constraint on $`\mathrm{\Omega }_0`$. The evolution of the abundance of X-ray clusters of a given mass (temperature) is one of the most powerful test of $`\mathrm{\Omega }_0`$ (Oukbir & Blanchard 1992) that has been developed during these last years (Viana & Liddle, 1996, Henry 1997, Eke et al. 1998, Sadat et al. 1998), but results are still contradictory (Blanchard et al. 2000). Clusters are the most massive objects for which both the luminous baryonic mass (dominated by the X-ray emitting intracluster gas) and the total gravitating mass can be determined and therefore an upper limit on the baryon fraction $`f_b`$ can be estimated. As primordial nucleosynthesis can provide reliable constraints on the value of the baryonic density parameter $`\mathrm{\Omega }_b0.030.08h_{50}^{2.}`$ ($`h_{50}`$ being hereafter the Hubble constant in unit of 50 km/s/Mpc; for numerical applications we used $`h_{50}=1.`$), one can derive an upper limit to the matter density parameter: $`\mathrm{\Omega }_0=\mathrm{\Omega }_b/f_b`$. White et al. (1993), estimated the baryon fraction in Coma cluster of $`f_b0.15h_{50}^{3/2}`$ at the Abell radius and several recent analyses of large samples of clusters have confirmed such a high baryon fraction $`f_b1520\%`$ (David et al, 1995, Evrard 1997, Mohr et al. 1999, Ettori & Fabian 1999, EF99 hereafter). This fact suggests that the baryon content observed in rich clusters exceeds the universal value expected in an $`\mathrm{\Omega }_0=1`$ universe. Since numerical simulations have shown that baryons cannot segregate significantly in clusters, this result implies that the mean density of the universe is lower than the critical value.

In this paper, we examine the issue of the baryon fraction in clusters, focussing on its distribution from both observational and numerical simulations sides.

## 2  The baryon fraction: Observations versus simulations.

Let us first examine the main results obtained from theoretical side, mainly from numerical simulations, and compare them to observations. The fate of gas in clusters is expected to be settled during the gravitational collapse of the mixture of baryons and dark matter. In the early phase of the collapse, both gas and dark matter are expected to follow the free fall trajectories until the gas is shock-heated which may lead to a segregation in the denser parts of clusters (corresponding more or less to the virialized region) between the gas and dark matter due to different physics governing the two fluids (Teyssier et al., 1997). In the case where only gravity is the source of heating of the gas, the baryon fraction is expected to follow a scaling law, depending only on the density contrast (Kaiser 1986). Hydro-dynamical simulations help to study quantitatively such a behavior. As a fiducial model we use the results of the simulations developed by different groups (the Santa Barbara project) and nicely summarized by Frenk et al. (1999, F99 hereafter). The theoretical gas fraction distribution for a value of the global gas fraction of 16% (corresponding to the range of what is generally inferred from the observations) is shown in Fig 1, as well as the range in which $`2/3`$ of the results from the simulations lie, taken as our 68% confidence interval. On the observational side, the gas fraction at a given radius is inferred from the gas mass derived from the X-ray gas emissivity (see section 4) and the gravitating mass which can be determined from either the hydrostatic $`\beta `$-model or from the universal dark matter profile (Navarro et al., 1995, NFW) normalized to the temperature-mass relation calibrated from numerical simulations. Isothermal hydrostatic mass estimates are systematically lower than numerical simulations calibrated masses: Roussel et al. (2000, hereafter RSB00) found gas fractions 20% smaller when virial masses are estimated from scaling laws calibrated from numerical simulations by Evrard, Metzler and Navarro, (1996, EMN96 hereafter). Taking into account for the temperature profile, hydrostatic equation would lead to masses nearly twice smaller (Nevalainen et al. 1999). Moreover, Balland and Blanchard (1997) have shown that even when accurate measurements of the temperature profile are available, hydrostatic mass estimates could suffer from a large uncertainty as later confirmed by Hughes (1998). To ensure consistency in the comparison with the simulations we restrict ourselves to the mass-temperature relations from numerical simulations and use the following two normalizations : $`T_X=4.75(M_{200}/10^{15}M_{})^{2/3}\text{keV}`$ (EMN96) and $`T_X=3.81(M_{178}/10^{15}M_{})^{2/3}`$ keV (Bryan and Norman 1998, BN98), $`M_\delta `$ being the mass enclosed in a region with a density equal to $`(1+\delta )`$ times the critical density. These two values for the normalization can be considered to be the extreme values among existing numerical simulations, BN98 normalization leading to virial masses nearly 40% higher.
We use the gas fractions estimated at various $`M_\delta `$ values from a compilation of X-ray data of a large sample of groups and galaxy clusters which are published in RSB00. For comparison, we also use the gas fractions at $`r_{500}`$ and $`r_{200}`$ (corresponding to $`\delta =500`$ and 200 respectively) quoted by Arnaud & Evrard (1999, AE99) who used a similar approach to determine the virial mass. The resulting gas fraction distribution is shown in Fig. 1. We have also included the point at $`r_{500}`$ given by EF99; while their original value was 17% (based on hydrostatic mass estimates), we found that this value translates to $`f_g13.7`$ at a radius of contrast density $`\delta 620`$, when one is using EMN96 mass estimator.
Two striking features can be immediately noticed when we compare the gas fraction distribution derived from observations to the simulations. Firstly, in the central regions where the X-ray data are of better quality, a strong disagreement exists between the data and numerical predictions for primordial fraction of the order of $`16\%`$. Secondly, the gas fraction profile inferred from observations shows a strong trend of increasing from the inner parts to the outskirts, as also found by David et al. (1995), contrary to theoretical predictions which profiles flatten to approach more or less the asymptotic value close to the universal value. Clearly, this discrepancy calls for caution when one is using the baryon (gas) fraction to set upper limit on the mean density of the universe.

## 3 Is the baryon discrepancy robust?

Let us now examine what could be the origin of this discrepancy. A first possibility is that the numerical simulations are suffering from numerical limitations that forbid a correct description of the baryon fraction shape. A second possibility is that the physics involved in the modeling of clusters is not adequate. A third possibility is that the strengh of the conclusions is not warranted by the data. In this section we will examine the first two possibilities and comment in the next section on the third one.
Numerical simulations have been extensively used in the modern cosmology. Several groups in the world have designed hydrodynamic numerical codes based on different techniques to simulate the gas dynamics in the cosmological context. F99 provides a comparison between these different codes allowing to test the ability of numerical simulations to reproduce the characteristics of cosmological formation of a cluster, assuming no cooling and heating processes. This comparison has been performed on the results of simulated X-ray clusters. F99 found that all the simulations agree to better than 10% on the gas mass at the virial radius and that the mean behavior of the gas fraction appears quite similar in all the simulations. Although noticeable dispersion in the absolute value of the gas fraction exists, the general behavior is always similar: the gas fraction raises up in the central part and flattens in the outer part. This suggests that numerical simulations are reliable on this issue given the assumed physics and thus the origin of the above difference cannot be due to numerical problems.

A second possibility, which seems more likely, lies in the physics of the formation of clusters: most of the numerical simulations assume adiabatic gas physics while there are convincing evidences that non gravitational processes are playing a significant role, otherwise the observed $`LT`$ relation could not be recovered. Heating of the IGM before or after the formation process is the often advocated mechanism (Cavaliere et al. 1998). The consequence of heating is to affect the gas distribution inside the cluster and consequently the radial distribution of $`f_b`$. One may think that this effect could provide a natural explanation to the steep increase of $`f_b`$ observed in the data. However, we believe that there are two reasons why this is not the most probable explanation (of course there is enormous freedom in designing a reheating scenario, so excluding totally this possibility is not possible) : i) even when extreme models with heating are considered, the asymptotic slope of the gas fraction is not as steep as observed. Fig. 1, shows the theoretical gas fraction at $`\delta =500`$ and 200 resulting from the simulations which include energy injection by SNe-driven winds (Metzler & Evrard, 1997). As we can see, the effect on the asymptotic slope is not very important. Even in the extreme cases of heating as illustrated by the small dashed curve in Fig. 1 (corresponding to low temperature clusters for which heating is expected to be very efficient), the predicted shape does not show the steep rising of the observed distribution of the gas fraction; ii) in the presence of strong heating large departure from the scaling laws is expected : the gas fraction as well as the asymptotic slope of the gas profile are expected to vary significantly with temperature. Although there is some debate on these questions, several groups have found no clear evidence of strong variations of the gas fraction with temperature (RSB00, Mohr et al. 1999, Vikhlinin et al. 1999, VFJ99 hereafter). Moreover, it has been shown that the emission from the gas follows roughly the scaling law in the outer regions (Neumann & Arnaud, 1999, VFJ99, RSB00) for clusters with $`T110`$ keV. Despite the wide range of possibilities on heating mechanism, and the limited number of cases that have been investigated, it seems to us difficult to significantly change the asymptotic slope and preserve at the same time the scaling of the gas fraction over a rather wide range of contrast densities.
These two arguments provide serious evidence that the apparent rising of the gas distribution in the outer parts of clusters is unlikely to be the result of some heating process.

## 4 Some remarks about gas estimation.

In this section issues concerning the reliability of the gas fraction as estimated from the observations are discussed.

### 4.1 From observations to gas masses

The X-ray gas masses are usually derived by fitting the X-ray surface brightness $`S(r)`$ to a $`\beta `$-model, $`S_o(1+(\theta /\theta _c)^2)^{3\beta +1/2}`$ where $`\theta `$ is the projected angular distance to the center. The mean value of the slope $`\beta `$ is generally found to be around 0.65 smaller than the theoretical asymptotic value derived from the simulations $`<\beta >=0.81`$. However, the X-ray gas emission is rarely traced out to very large radii, so that gas fractions are measured only up to an X-ray limiting radius $`R_{Xlim}`$ for which the signal-to-noise is good enough. Most of the time, it was necessary to extrapolate to the virial radius $`r_{200}`$ or to some other outer radius like $`r_{500}`$. However, it is essential to realize by doing so, that gas masses are estimated at a radius at which the actual emission is not well known, poorly constrained and with a value of the order of the X-ray background or less. Furthermore, the emissivity being dominated by the central region, a parametric fit will be rather insensitive to the outer profile. The extrapolation of a typical $`\beta `$-model to asymptotic regions ($`rr_v`$) leads to $`\rho _{gas}r^{3\beta }r^{1.95}`$ (for $`<\beta >0.65`$) which is significantly less steep than the universal dark matter profile (NFW) found in the simulations $`\rho _{DM}r^{2.7}`$. These two facts may actually explain the trend of steep increasing in the $`f_b`$ profile (Fig. 1). It is therefore vital to examine whether the extrapolation of the $`\beta `$-model up to the virial region does or not introduce a bias in the gas estimations. Recently, a large sample of X-ray ROSAT images in which the data trace the emission up to very large radii have been analyzed by VFJ99. They found that the $`\beta `$-model does not provide an accurate description of the surface brightness over the whole range of radii, the slope in the outer regions (0.3$`r_{200}`$ to $`r_{200}`$) being actually steeper than with the one found when the $`\beta `$-model is applied to the entire cluster. This result is very important as it shows that a dedicated treatment of the outer parts of the cluster leads to a different gas distribution than what is generally assumed and provides direct evidence that the extrapolation to the virial radius of the $`\beta `$-model could be unsafe.

### 4.2 The clumping of the gas

Most of the previous studies of the gas and mass distribution assume that clusters are spherically symmetric and that the ICM is uniform. In reality, most of the clusters exhibit a certain level of density fluctuations at small and large scales due to accretion and merger events which continuously occur inside clusters. This clumping may introduce a further potential bias in gas mass estimates and therefore in the inferred gas fractions. Furthermore it is reasonable to assume that the clumping of the gas is more pronounced in the outer part, as the gas may not have enough time to relax. If the gas is locally inhomogeneous, it is useful to introduce the clumping factor $`C=<\rho _g^2>/<\rho _g>^2`$ where $`\rho _g`$ is the gas density. In presence of clumping, gas mass estimates are overestimated by roughly $`C^{1/2}`$. Moreover, this factor is likely to increase in the outer parts. Studying the effects of density fluctuations on small scale and large scale substructures on the ICM mass measurement of a set of simulated clusters, Mathiesen et al. (1999, MEM99) have quantified the clumping effect and found that the estimation of the gas mass from models assuming a uniform density carries a bias of $`\sqrt{C}1.16`$ at $`\delta =500`$. It is worth noticing that this factor is probably still poorly known, as its measured value might be affected by the numerical simulation limitation. Interestingly enough is to notice that although the various simulations reported in the Santa Barbara project do find a nearly constant baryon profile in the outer regions of the cluster, much larger dispersions are found in the emissivity at large radii (more than a factor of 2, see their figure 19), which may either reflect the numerical limitations or the differences in the clumping factor treatment.

## 5 The baryon fraction revisited

Here we evaluate quantitatively the consequence of the above biases. From what we have seen, robust gas mass estimates require that the emissivity is actually detected up to the radius at which gas masses are given (notice that this may still introduce a kind of Malmquist bias). To our knowledge, only VFJ99 sample satisfies this criteria: their published results allow us to compute the gas mass content enclosed within two different radii, $`R_{1000}`$ and $`R_{2000}`$ (see their equation 2): $`M_g(\delta )=4\pi /3(1+\delta )\rho _bR_\delta ^3`$ (using their $`\rho _b`$). We then computed the mass enclosed within the same physical radius assuming an universal profile (NFW) with a concentration parameter $`c=5`$. As they found a slight temperature dependence we work at their median temperature $`T=6`$ $`\mathrm{keV}`$. We obtain $`f_g(\delta =770)=0.1064`$ and $`f_g(\delta =343)=0.1194`$ with EMN calibration and $`f_g(\delta =974)=0.0842`$ and $`f_g(\delta =437)=0.0934`$ with BN98 normalization. These values are similar to those derived by RSB00, and smaller than those of AE and EF by nearly 2%.
We then correct for the clumping effect of the gas: not much is known about its properties, we therefore restrict ourselves to correct only the outer value at $`\delta =500`$, using the factor given by MEM99. In the inner parts of the cluster we use RSB00 gas fractions, assuming that no correction has to be applied to them. The resulting values are shown as big stars and circles in Fig. 2. By matching the point at $`\delta =500`$, we obtain the following primordial gas fractions : $`f_g=0.0875\pm 0.0075(68\%)`$ (BN98) and $`f_g=0.108\pm 0.0092`$ (EMN96), the uncertainties being due to numerical simulations. Two important features arise: i) The gas fraction now appears to be consistent over a large range of $`\delta `$ with numerical simulations predictions: in particular we can see that the asymptotic flattening is now recovered. Note that, the observations fall off slightly below the predicted profile in the inner parts ($`\delta >10^4`$), this could well be due to energy injection, ii) the inferred gas fraction at $`r_{500}`$ is now lower than what is previously found (Fig. 1). The uncertainty being due to the uncertainty in the calibration of the $`MT`$ relation from numerical simulation.

In order to estimate the total baryon fraction we add the stellar contribution, taking a mean value of the order of $`1\%`$ (RSB00). Our final values are : $`f_b0.10\pm 0.01`$ (with BN98 calibration) and $`f_b0.12\pm 0.01`$ (with EMN96 calibration).
It is worth noticing that the baryon fraction shape inferred in this present work (Fig. 2) corresponds to a steeper gas distribution asymptotic slope than the one derived from the standard $`\beta `$-model and even steeper than VFJ99 slopes. Plugging this gas profile into the isothermal hydrostatic equilibrium would lead to masses which are much higher than previous estimates and thus more consistent with masses calibrated from numerical simulations.

## 6 Implications on $`\mathrm{\Omega }_0`$

The arguments presented above clearly show that previous mean baryon fraction estimations at the virial radius are probably biased due to the extrapolation of data and to the effect of gas clumping. This leads us to re-examine the constraint on $`\mathrm{\Omega }_0`$ set by the baryon fraction. Following the notation of White et al. (1993) :

$$\mathrm{\Omega }_0=\gamma (r_{500})\frac{\mathrm{\Omega }_b}{f_b(r_{500})}$$
(1)
where $`\gamma (r_{500})0.91`$ is the baryon depletion as predicted by the simulations. $`\mathrm{\Omega }_b`$ is set by the primordial nucleosynthesis theory through the baryon-to-photon parameter $`\eta _{10}`$ which is constrained by the measurement of the abundance of light elements. The main uncertainty comes from the deuterium abundance (D/H) measurement. There are two conflicting estimates of D/H, a high value which leads to $`\eta _{10}2`$ (Rugers & Hogan 1996) and a low value corresponding to $`\eta _{10}5.6`$ ( O’Meara et al., 2000). This latter value is consistent with that measured by Wadsley et al. (1999) using the He Lyman–$`\alpha `$ forest and the amount of baryons in order to produce the observed Lyman–$`\alpha `$ systems (Rauch et al., 1997; Weinberg et al., 1997) and are favored by CMB fluctuations properties (Le Dour et al, 2000; Tegmark & Zaldarriaga, 2000). Taking $`\mathrm{\Omega }_bh_{50}^20.08`$ (assuming that all baryons in clusters are seen) and using our estimated primordial baryon fraction obtained with BN98 calibration ($`f_b=0.10\pm 0.01`$), we estimate the density of the universe: $`\mathrm{\Omega }_0=0.8\pm 0.1`$ in agreement with the value derived from the evolution of the X-ray clusters abundance (Blanchard et al. 2000).
If we use the primordial baryon fraction obtained with EMN96 calibration we find $`\mathrm{\Omega }_00.65`$.

## 7 Conclusion

Baryon fractions derived from X-ray clusters of galaxies are generally in the range (15% - 20%) (for $`h_{50}=1.`$). This exceeds the value predicted by Big Bang nucleosynthesis in an $`\mathrm{\Omega }_0=1`$ universe which has led to the conclusion that we are living in a low-density universe with $`\mathrm{\Omega }_00.20.4`$. However, we found that apparent gas fraction profiles, which obey reasonably well to scaling laws, differ in their shape from what is found in numerical simulations: these gas fraction profiles seem to continuously increase from the inner to the outer parts, while it flattens asymptotically, approaching the universal value at the virial radius in the numerical simulations.
Using updated X-ray data on the gas distribution out to large radii from VFJ99, binding mass based on numerical simulations and correcting for the clumping effect, we obtain the following results: i) the overall shape of the baryon fraction is roughly consistent with numerical simulations. ii) This agreement is obtained provided that the global baryon fraction is of the order of 10%, while higher values of the order of 15–20% are clearly ruled out. iii) This revised value of the baryon fraction, combined with current limits on $`\mathrm{\Omega }_b`$ from primordial nucleosynthesis, is consistent with a high density of the universe. Thus we conclude that the apparent discrepancy between the observations and numerical simulations in the baryon fraction distribution inside clusters can be erased once the above systematics are taken into account. This work has also raised the following question: is the baryon fraction a strong argument against $`\mathrm{\Omega }_0=1`$? Several effects and systematics may act in concert to create the apparent gap between the current estimates of the baryon fraction in clusters and the universal value, as it is shown in this work, therefore it seems to us that the baryon fraction is no longer a strong and “clean” argument in favor of a low density universe as is too oftenly claimed. It is therefore important to better constrain and control all the systematics acting in the gas fraction measurements, this will be soon possible with XMM, AXAF and SZ experiments. On the other hand numerical simulations of X-ray clusters are also needed to quantify the clumping effect as well as the role of non-gravitational physics.

###### Acknowledgements.

We would like to thank the anonymous referee for his useful remarks which have significantly improved this work. R. Sadat acknowledges the support by the R. Chrétien International Research Grant.


format: text

text:
# Quantum computing in a Spin ordered phase

## Abstract

We show that quantum computation can be performed in a system at thermal equilibrium if a spontaneous symmetry breaking occurs. The computing process is associated to the time evolution of the statistical average of the qubit coherence operator. This average defines logical states which evolve under the action of impulsive gate perturbations. Non trivial quantum coherence properties appear in the ordered phase characterized by a spontaneous symmetry breaking. Qubits are associated to spins and interaction with the lattice vibrations is approximated by an effective spin-spin interaction. A finite decoherence time is obtained only if spin correlations are taken into account in the framework of the Ising model as corrections to the Mean Field approximation. Decoherence is slowed down as the order parameter saturates.

In the quantum computing process information is stored in qubits i.e. in a quantum state of a two levels system. Coupling with the environment produces an “entanglment” of the qubits with the degrees of freedom of the environment. This “entanglement” is the origin of a progressive loss of information, a process known as decoherence. Hence, systems at thermal equilibrium have not been considered suitable for a computing process. We show that quantum computation is however possible even in thermal equilibrium in presence of a spontaneous symmetry breaking. Logical states of computer will be associated to suitable thermal averages. The 0-logical state is associated to thermal equilibrium and 1-logical state to local excitation determined by an impulsive perturbation. Decoherence is then associated to the local excitation relaxation toward equilibrium. We shall see on a specific model that relaxation vanishes in the Mean Field Approximation (MFA) and is slowed down as magnetization saturates when corrections to MFA are taken into account.

For sake of simplicity we work with one of the most simple system which exhibits symmetry breaking: the Ising model $`H=_{i,j}J_{i,j}S_i^zS_j^z`$. This model is related to the more familiar spin - boson model where spins interact with environment vibrations. If up and down local spin states are degenerated in energy, lattice degrees of freedom can be eliminated by a canonical transformation. Ising model is then obtained as an approximation if lattice modes polarized in the spin quantization direction have smaller frequencies than perpendicular ones. Another physical realization is given by quantum dots interacting with the oscillation modes of the substrate. In this case the two states ( $`|>`$ and $`|>`$ ) can be associated to the presence or absence of an extra electron. In Coulomb blockade condition, i.e. when the two states are degenerated in energy, the coupling of quantum dots with phonons gives rise to an effective interaction between spins equivalent to the Ising model. Another possible physical realization of this model is a system of bosonic neutral atoms localized in the sites of a lattice interacting with repulsive interaction. Possible realization of such a system in an optical lattice has been discussed in .

We define the i-th qubit as the state $`|\theta _i,\varphi _i>=\mathrm{cos}(\frac{\theta _i}{2})|_i>+e^{i\varphi _i}\mathrm{sin}(\frac{\theta _i}{2})|_i>`$ associated to the corresponding site of a regular array ( i=0,..,N-1 ). The local quantum coherence operator is the projection operator on a given qubit state.

$$\pi _i(\theta ,\varphi )=|\theta _i,\varphi _i><\theta _i,\varphi _i|$$
(1)
The projection operator in terms of spin operators reads as

$$\pi _i(\theta ,\varphi )=\frac{1}{2}+n_iS_i$$
(2)
where $`n`$ is the unitary vector associated to the $`\theta `$, $`\varphi `$ direction and $`S`$ is the vector spin operator. The statistical average of the quantum coherence operator is the qubit coherence probability (QCP)

$$\rho _i=\frac{Tr(e^{\beta H}\pi _i)}{Tr(e^{\beta H})}$$
(3)
$`\rho _i`$ is the probability of finding the spin at the i-th site, in the state $`|\theta _i,\varphi _i>`$, when the system is in equilibrium with an external bath. It is immediately seen from (3) that as a consequence of the spin inversion symmetry of the Ising Hamiltonian $`H(\{S_i\})=H(\{S_i\})`$ we have $`\rho _i=\frac{1}{2}`$. All directions $`\theta _i,\varphi _i`$ are then equally probable and coherence, as expected,is completely lost at equilibrium. Thus a “symmetric” thermodynamic state is obviusly useless as a reference state in a computation process. However, if a spontaneous symmetry breaking occurs below a critical temperature ($`T_c`$), we obtain a non trivial QCP. In our case spin inversion symmetry breakdown implies $`m_i=<S_i^z>0`$ ( angular brackets stand for thermodynamic average ). As a consequence the QCP becomes

$$\rho _i=\frac{1}{2}+m_i\mathrm{cos}\theta _i$$
(4)
We see in (4) that coherence is partially preserved at equilibrium even if the dependence on the phase $`\varphi _i`$ is lost. In the states space the direction $`\theta _i=0`$ is the most probable and the average of any even function of $`\theta _i`$ does not vanish. Hence, the thermodynamic state associated to the ordered spin phase is suitable as the reference state of a computing process. It is convenient to consider a reduced QCP as the difference between the ordered and disordered phase QCPs $`\stackrel{~}{\rho }_i=<\stackrel{~}{\pi }_i>=<\pi _i\frac{1}{2}>`$. We find convenient to identify the 0-logical state ( 1-logical state ) of the i-th site with the dependence of $`\stackrel{~}{\rho }_i`$ on $`\mathrm{cos}\theta _i`$ $`(\mathrm{sin}\theta _i)`$. In other words different logical states are associated with different properties of a quantity related to the statistical average of the quantum coherence instead of different quantum states.

Universal quantum computer needs 1-qubit and 2-qubits gates. The first is a rotation of the 0-logic state in a superposition of 0 and 1 logical states. For this reason we need a perturbation that modifies the local QCP introducing a term with depends on $`\mathrm{sin}\theta _i`$. The second modifies the two sites state introducing an entanglement between the two logical states. In the corresponding 2-qubits QCP a $`\mathrm{sin}\theta _i\mathrm{sin}\theta _j`$ will appear as a result of the perturbation acting on a two sites 0 logical initial state.

For sake of simplicity we consider only impulsive perturbations.

$$H_t=H+ϵ\delta (t0^+)H_I$$
(5)
After the action of a gate the QCP is given by

$$\stackrel{~}{\rho }_i(t)=\frac{Tr(e^{\beta H}\stackrel{~}{\pi }_i(t))}{Tr(e^{\beta H})}$$
(6)
Where $`\stackrel{~}{\pi }_i(t)`$ is the reduced quantum coherence operator in the Heisenberg representation.

As 1-bit gate we consider a perturbation perpendicular to the z direction acting on a given site i=0. We choose to work with

$$H_I=S_0^x$$
(7)
If the Ising model is associated with localized particles ( quantum dots, optical lattices, etc. ), the perturbation could be obtained by exchanging particles with a reservoir while in the spin case a magnetic field could be used. The new logical state is characterized as

$$\stackrel{~}{\rho }_0(t)=\stackrel{~}{\rho }_0\mathrm{cos}ϵ+A_0(t)\mathrm{sin}\theta _0\mathrm{sin}ϵ$$
(8)
where $`A_0(t)`$ is given by

$$A_0(t)=\mathrm{sin}\theta _0<S_0^z\mathrm{sin}(2h_0t+\varphi _0)>$$
(9)
and $`h_0=_jJ_{0j}S_j^z`$ is the internal field. The strength of the perturbation determines the amplitude of the rotation from the 0-logical state $`(ϵ=0)`$ and 1-logical state $`(ϵ=\frac{\pi }{2})`$. Note that only the 1-logical state component of (8) will evolve in time and then will be affected by decoherence ( $`A_i(t)0`$ ) while permanent modification are found in the 0-logical state amplitude. A simple oscillating behaviour is obtained in the MFA when internal field fluctuations are neglected, i.e. $`h_i=_jJ_{ij}m_j`$ and $`m_i=\frac{1}{2}\mathrm{tanh}\left(\beta h_i\right)`$. It will be shown that a finite decoherence time arises as a consequence of non vanishing spin correlations $`c_{ij}=<S_iS_j><S_j><S_j>`$ ( neglected in the MFA ).

The 2-qubits gate can be realized by switching on a particles transfer between two sites. A simple result is obtained in the case of sites i,j belonging to two non-interacting adiacent layers of opposite magnetization. The Quantum coherence of the 2-qubits states associated with sites i and j is conveniently characterized by the statistical average of the product of reduced two sites quantum coherence operators $`\stackrel{~}{\rho }_{ij}=<(\pi _i\frac{1}{2})(\pi _j\frac{1}{2})>`$. At equilibrium the two sites 0-logical states is characterized as

$$\stackrel{~}{\rho }_{ij}=<S_i^zS_j^z>\mathrm{cos}\theta _i\mathrm{cos}\theta _j$$
(10)
Note that in the MFA, where spin correlation are neglected, i.e. $`<S_iS_j>=<S_j><S_j>`$ $`(ij)`$ the two sites reduced QCP is equal to the product of one site reduced QCP. The impulsive perturbation hamiltonian associated with sites i,j is then

$$H_I=S_i^+S_j^{}+S_j^+S_i^{}$$
(11)
Where i and j belong to two different non-interacting layers of opposite magnetization. In the time evolution of the reduced two sites QCP, following the impulsive perturbation we distinguish a constant term and a rotated one associated only with $`|_i,_j>`$ and $`|_i,_j>`$ states.

$$\stackrel{~}{\rho }_{ij}(t)=\stackrel{~}{\rho }_{ij}(0)+A_{ij}\mathrm{sin}\left(2ϵ\right)\mathrm{sin}\theta _i\mathrm{sin}\theta _j$$
(12)
where $`A_{ij}`$ is given by

$$A_{ij}=\frac{1}{4}<\left(S_i^zS_j^z\right)\mathrm{sin}\left[\left(h_ih_j\right)t+\varphi _i\varphi _j\right]>$$
(13)
The impulsive local interaction produces the 1-logical state on two sites which oscillate with the interaction strenght. The peculiar superposition of 0 and 1-logical states is a consequence of the quantum entanglement induced by the interaction. In the MFA limit, the internal fields $`h_i`$ and $`h_j`$ become non fluctuating quantities to be determined self consistently. To evaluate decoherence corrections to MFA must be taken into account in the framework of a perturbation expansion in the interaction range . If the interaction range is small with respect to the system size the relaxation toward equilibrium can be seen as a quantum diffusion process wich spreads the perturbation over the whole lattice. We use the technique ref. of to calculate an auxiliary free energy functional where a single site time dependent field appears

$$G^\sigma (\lambda ,t)=\mathrm{ln}\left\{Tr\left[exp\left(\beta H+\underset{k}{}(2i\sigma tg_k+\lambda _k)S_k^z\right)\right]\right\}$$
(14)
Here $`\sigma `$ is a dicotomic quantity ($`\sigma =\pm `$) and $`g_k=J_{0k}`$ in the one site ( i=0 ) gate case and $`g_k=(J_{0k}J_{1k})`$ in the two sites ( i=0, j=1 ) gate. In terms of this dynamical functional we can express both the one (9) and the two sites (13) QCPs variations due to one and two sites gates

$$A_0=\frac{i}{2}\underset{\sigma =\pm }{}\sigma exp\left(G^\sigma (t)G^\sigma (0)\right)\frac{G^\sigma }{\lambda _0}|_{\{\lambda _j\}=0}$$
(15)
$$A_{01}=\frac{i}{8}\underset{\sigma =\pm }{}\sigma exp\left(G^\sigma (t)G^\sigma (0)\right)e^{i\sigma (\varphi _0\varphi _1)}\left(\frac{G^\sigma }{\lambda _0}\frac{G^\sigma }{\lambda _1}\right)|_{\{\lambda _j\}=0}$$
(16)
As a simple generalization of procedure of ref. , we introduce a Legendre transform at fixed time and temperature

$$F^\sigma (m,\beta ,t)=G^\sigma (\lambda ,\beta ,t)\underset{k}{}\lambda _k(\beta ,t)m_k$$
(17)
Where $`m_k=\frac{G^\sigma }{\lambda _0}`$. The dynamical functional of (17) is obtained adding an imaginary time dependent term to the static functional.

$$F^\sigma (m,\beta ,t)=F(m,\beta )+\underset{k}{}2i\sigma tg_k$$
(18)
The time independent part of the (18) $`F(m,\beta )`$ is given in ref. . It is conveniently splitted in two terms  $`F(m,\beta )=F_0(m)+F_1(m,\beta )`$ . The first term is associated to the non interacting system with fixed magnetization and the latter takes into account sites interaction.

$$F_0(m)=\underset{k}{}\left[(m_k+\frac{1}{2})\mathrm{ln}(m_k+\frac{1}{2})+(\frac{1}{2}m_k)\mathrm{ln}(\frac{1}{2}m_k)\right]$$
(19)
$$F_1(m,\beta )=+\beta \underset{lk}{}J_{lk}m_lm_k+\beta ^2\underset{lk}{}J_{lk}^2\left(\frac{1}{4}m_{k}^{}{}_{}{}^{2}\right)\left(\frac{1}{4}m_{l}^{}{}_{}{}^{2}\right)+\mathrm{}$$
(20)
The main point is to study the extremum condition $`\frac{F^\sigma }{m_k}=0`$ in the complex plane.

$$\frac{F(m,\beta )}{m_k}=2i\sigma tg_k$$
(21)
For short time ( $`J_{0k}t1`$ ) we can express the solution of the extremum problem as $`m_k^\sigma (t)=m_k+\delta m_k^\sigma `$ where $`m_k`$ is the solution at t = 0. Expanding to the first order equation (21) we obtain

$$\underset{l}{}\frac{^2F(m,\beta )}{m_km_l}\delta m_l^\sigma =2i\sigma tg_k$$
(22)
From (22) we see that $`\delta m_k^\sigma =it\sigma \delta \overline{m}_k`$ where $`\delta \overline{m}_k\frac{1}{D}`$ and, then, vanishes in the MFA limit. The same expansion gives for the dynamical function of (17)

$$F^\sigma (m,\beta ,t)F(m,\beta )+\underset{k}{}\left[2i\sigma tg_km_kt^2\delta \overline{m}_kg_k\right]$$
(23)
From (22) we can derive the time dependence of $`F(m)`$

$$F^\sigma (m,\beta ,t)F(m,\beta )+\underset{k}{}2i\sigma tg_km_k\mathrm{\Gamma }t^2$$
(24)
The quantity $`2it\sigma _kg_km_k`$ is related to the oscillating behaviour of the MFA while $`\mathrm{\Gamma }=_kg_k\delta \overline{m}_k`$ is real quantity related to the decoherence time. Because of the stability condition on the extremum of $`F(m)`$, $`\mathrm{\Gamma }>0`$. If the number of next neighbors of site 0 ( for one bit gate ), or sites 0 and 1 ( for two bits gate ), which are next neighbors each other do not scale with the dimensionality, then $`\delta \overline{m}_k`$ to first order in $`\frac{1}{D}`$ is equal to $`\delta \overline{m}_k=g_k\left(\frac{1}{4}m_k^2\right)`$ and $`\mathrm{\Gamma }`$ becomes

$$\mathrm{\Gamma }=\underset{k}{}g_k^2\left(\frac{1}{4}m_k^2\right)$$
(25)
As expected $`\mathrm{\Gamma }`$ is proportional to $`\frac{1}{D}`$ and vanishes in the limit of infinite dimensionality ( MFA limit ). It is important to observe that decoherence decreases as magnetization approaches saturation $`(m^2\frac{1}{4})`$.

In conclusion we studied quantum computation in a spin system at thermal equilibrium in the presence of a spin symmetry spontaneous breaking. With respect to the conventional approach based on out of equilibrium systems two main results seem to be of interest.

a) It is possible to define logical states as suitable thermal averages instead of quantum states.

b) Decoherence, due to the conservative dynamics associated to energy operator wich defines statistical equilibrium, can be treated as a perturbation effect with respect to MFA and it is even possible to reduce it at will by increasing the order parameter.

We think that the main results of our work should be valid for other systems which exhibit spontaneous symmetry breaking.


text:
# Study by MOA of extra-solar planets in gravitational microlensing events of high magnification

## 1 Introduction

A remarkable diversity of planetary systems has been uncovered in recent years through studies of extra-solar planets (Mayor & Queloz 1995, Marcy & Butler 1998, Perryman 2000). Many further studies will be required to obtain a full understanding of this diversity of planetary systems and their formation processes. In this paper we report a study of extra-solar planetary systems using a gravitational microlensing technique.

Mao and Paczynski (1991) pointed out that extra-solar planets could be detected using gravitational microlensing111Liebes (1964) may have been the first to consider the detection of planets by microlensing and also the first to consider high magnification microlensing events of the type described in this paper, because the characteristic scale in gravitational microlensing happens to coincide quite closely with the characteristic size of the solar system. The characteristic scale of gravitational microlensing is the radius of the Einstein ring, $`R_\mathrm{E}`$, which is

$$R_\mathrm{E}=\sqrt{\frac{4GM_\mathrm{L}}{c^2}}\times \sqrt{\frac{D_{\mathrm{OL}}D_{\mathrm{LS}}}{D_{\mathrm{OS}}}}.$$
(1)
Here $`M_L`$ denotes the mass of the lens, and $`D_{\mathrm{OL}}`$, $`D_{\mathrm{LS}}`$ and $`D_{\mathrm{OS}}`$ the distances from the observer to lens, lens to source and observer to source respectively. For microlensing events that occur in the galactic bulge, $`M_\mathrm{L}`$ is typically $`0.3M_{}`$, and $`D_{\mathrm{OL}}`$, $`D_{\mathrm{LS}}`$ and $`D_{\mathrm{OS}}`$ are typically about 6 kpc, 2 kpc and 8 kpc respectively. These values imply $`R_\mathrm{E}`$ 1.9 AU. Consequently, if other planetary systems have the same characteristic size as ours, planets orbiting a lens star may perturb the lensing, because light from the source star may pass close to one or more of the planets at some time during a microlensing event. Gould & Loeb (1992) and Bolatto & Falco (1994) showed that Jupiter-like planets could be detected with about 20% efficiency in typical microlensing events if they are monitored at about hourly intervals with a photometric precision of a few per cent throughout each event. The PLANET collaboration has used this technique to set an upper limit on the abundance of Jupiter-like planets of order 30% (Albrow et al. 2001). This result is consistent with results obtained independently by radial velocity and transit measurements (Marcy & Butler 1998, Charbonneau et al. 2000, Henry et al. 2000, Jha et al. 2000). Bennett & Rhie (2000) emphasized the detectability of terrestrial planets by microlensing where they showed that Earth-mass planets could be detected with about 5% efficiency in typical events if they are monitored hourly at a photometric precision of better than 1%.

A modification of the gravitational microlensing plant search strategy was pointed out by Griest & Safizadeh (1998) in which planets, including those less massive than Jupiter, could be detected with high efficiency in microlensing events of high magnification. This occurs when the distance of closest approach of the lens star to the line-of-sight to the source star, $`u_{\mathrm{min}}`$, is much smaller than the Einstein radius, i.e. when $`u_{\mathrm{min}}1`$ where $`u_{\mathrm{min}}`$ is expressed in units of $`R_\mathrm{E}`$. In these events a circular, or near-circular, image generally222The exceptional situation where “caustic crossings” occur is discussed in §5. forms around the lens with radius $`R_\mathrm{E}`$ at the time of maximum amplification, and the maximum amplification $`A_{\mathrm{max}}1/u_{\mathrm{min}}1`$. If a planet is orbiting the lens star at a projected radius of about a few AU at the time of an event, it will perturb the image, although the perturbation will be small because most of the image will lie far from the planet and remain unperturbed. For planets with masses as low as that of Earth, simulations show (see §5 below) that the perturbation may nevertheless be detectable with appropriate observations. The aim of the present study was to exploit this sensitivity.

The first analysis of a high amplification event in terms of planets was reported by Rhie et al. (2000). These authors reported a low-mass planetary signal (few Earth masses to few Uranus masses) at a marginal level of confidence for the event MACHO 98–BLG–35, and excluded Jovian planets over an extensive region surrounding the lens star. They also pointed out that further such results in additional events would provide statistically significant constraints on the abundance of Earth-mass planets, but that more accurate planetary parameters can be obtained in events in which a ”planetary caustic” is crossed or approached, which generally occurs at lower magnifications.

In this paper we report a re-analysis of MACHO 98–BLG–35, where the data in Rhie et al (2000) obtained by MOA have been analysed using an improved technique based on image subtraction and the inclusion of additional data. We also report analyses of two more events of high magnification, MACHO 99–LMC–2 and OGLE 00–BUL–12. We note that the event MACHO 99–LMC–2 coincides with the event OGLE 99–LMC–1 that was found independently by the OGLE group.

## 2 Observational Datasets

Data from the MOA (Abe et al. 1997), OGLE (Udalski et al. 1997), MPS (Bennett et al. 1999) and PLANET (Albrow et al. 2001) microlensing groups333MOA website: www.vuw.ac.nz.scps/moa/
OGLE website: bulge.princeton.edu/~ogle/
MPS website: bustard.phys.nd.edu/MPS/
PLANET website: thales.astro.rug.nl/~planet/ were used in the present analysis. The datasets are summarized in Table 1. Wherever possible, we attempted to obtain actual images and analyse these using an image subtraction procedure to achieve optimum photometric precision in the dense stellar fields in which microlensing is observed. The reduced datasets so obtained are posted at the MOA website.

Observations of MACHO 98–BLG–35 and MACHO 99–LMC–2 by MOA were made with the cameras MOA-cam1 and MOA-cam2 respectively (Yanagisawa et al. 2000). We note here that MOA-cam2 was operating in a slightly non-linear mode during 1999, and that the MOA data on MACHO 99-LMC-2 were corrected for this on a pixel-by-pixel basis using linearly calibrated stars from observations made in 2000.

Datasets on MACHO 98–BLG–35 by the MPS and PLANET groups were included in the present analysis. These datasets were reduced by these groups using the DoPHOT procedure of Schechter et al. (1993) and posted at their websites. In the case of the MPS dataset, it is the same dataset that was used previously by Rhie et al. (2000). In the case of the PLANET dataset, the data used here were obtained by extracting the unpublished, graphical data at the PLANET website, and allowing for a time offset of 0.007 day (A. Gould, private communication, 2000). The accuracy of the graphical extraction was much better than the error bars on the data. This may be seen by comparing the PLANET data shown in Fig. 10 of the present paper (see below) with Fig. 19 of their subsequent publication (Gaudi et al 2002). The need for the time offset, which was presumably caused by a heliocentric correction and/or simple clock errors of a few minutes, became clearly evident from comparisons of the raw PLANET and raw MPS/MOA datasets for the event

## 3 Image Subtraction Analysis

The observations listed in Table 1 by the MOA and OGLE groups were reduced by an image subtraction procedure. Data from each telescope in each passband were treated as separate datasets. For each dataset, the flatfielded images were first geometrically aligned to an astrometric reference chosen from amongst the best seeing images. A reference image for the image subtraction process was then formed by stacking the very best seeing and signal-to-noise images. The convolution kernel which relates the seeing on the reference image to that for a particular image was then determined using our own implementation of the technique of Alard & Lupton (1998) for modeling the kernel along with the modification of Alard (1999) which models the spatial variations of the kernel across the CCD field of view.

Photometry on the subsequent dataset of subtracted images was carried out by first constructing an empirical point-spread-function (PSF) on the reference image using 5–6 bright, isolated stars near the event. A PSF for a particular image was then formed by convolving the reference PSF with the appropriate kernel. This PSF was fitted to the flux profile on the subtracted image at the position of a star to obtain a “delta-flux” measurement, $`\mathrm{\Delta }F`$. The fit was done by re-aligning the centroid of the empirical PSF to that of the $`\mathrm{\Delta }F`$ profile on the subtracted image, and performing a least squares fit to the pixel by pixel cross-plot of the two profiles. Typical fits are shown in Fig. 1.

The scatter in fits such as those shown in Fig. 1 provides an indicator of the quality of the image subtractions. We used the standard deviation from the best fit straight line as a measurement of the subtraction quality associated with a given $`\mathrm{\Delta }F`$ measurement. In a good quality subtraction, the profile of an object on the subtracted image should match that of the PSF on the unsubtracted image, and thus have a low standard deviation. Uncertainties in the $`\mathrm{\Delta }F`$ measurements for any object were determined empirically. The profile fitting technique was applied to the positions of all resolved stars in the field. The frame-to-frame statistical errors for any particular object were then determined from the scatter in the $`\mathrm{\Delta }F`$ values for stars with similar profile standard deviations.

Stars close to the microlensing events were selected as checks. Fig. 2 shows stars with similar statistical profile errors as MACHO 98–BLG–35, and Fig. 3 shows their light curves. Figs. 46 show similar stars from the OGLE database on OGLE 00–BUL–12. It is evident from these that the systematic errors in the image subtraction reduction procedure are significantly less than the statistical errors, and that the statistical errors are realistic.

## 4 Single Lens Fitting

For the datasets that were obtained by image subtraction the fluxes in any passband were fitted to the function

$$\mathrm{\Delta }F(t)=F_\mathrm{b}A(u(t))F_\mathrm{r}$$
(2)
where $`F_\mathrm{b}`$ denotes the baseline flux of the source star, and $`F_\mathrm{r}`$ is the source flux on the reference image used in the image subtraction process, which in general contains some lensed flux444For datasets obtained by DoPHOT a similar equation holds, viz. $`F(t)=F_\mathrm{b}A(u(t))+F_\mathrm{u}`$, where $`F_\mathrm{u}`$ denotes the flux from nearby unresolved stars that are not recognized by the DoPHOT reduction procedure as being separate stars, and which are unlensed. This equation was used to incorporate the data of the MPS group in the analysis of MACHO 98–BLG–35 presented in §5.1..

The well known amplification factor, $`A(u)`$, expressed in terms of the distance $`u`$ of the lens star from the line-of-sight to the source star expressed in units of $`R_\mathrm{E}`$, is given by (Paczynski 1986)

$$A(u)=\frac{u^2+2}{u\sqrt{u^2+4}},$$
(3)
where

$$u(t)=\sqrt{u_{\mathrm{min}}^2+\left(\frac{tt_{\mathrm{max}}}{t_\mathrm{E}}\right)^2}.$$
(4)
Here $`u_{\mathrm{min}}`$ is the minimum value of $`u`$ during an event and $`t_{\mathrm{max}}`$ is the time of maximum amplification, i.e. the time when $`u=u_{\mathrm{min}}`$. The quantity $`t_\mathrm{E}`$ is the time scale that characterizes an event. It is equal to $`R_\mathrm{E}/v_\mathrm{T}`$ where $`v_\mathrm{T}`$ is the transverse velocity of the lens with respect to the line-of-sight to the source. In fitting a single lens curve to photometry taken in one or more passbands for a given event, there are therefore three microlensing parameters and two flux parameters for each passband.

The light curves for the three events are shown in Figs. 79 together with the best single lens fits to the data. For MACHO 98–BLG–35 only the data obtained by the MOA group have been included, since these were the only data to be reduced by image subtraction. The light curve for this event appears to exhibit substantial deviations from the single lens fit near its peak, but not at other times, suggestive of the presence of a planet (or planets) orbiting the lens star. Other possible causes of a deviation near the peak include the source star being a binary or having starspots. The binary interpretation seems unlikely (Rhie et al 2000). Also it would seem that starspots are unlikely to produce perturbations of the same strength as planetary deviations except in extreme situations (Rattenbury et al 2002, submitted to MRAS). For MACHO 99-LMC-2 no data were available to us at the peak, and in this case no apparent deviation from the single lens light curve is seen. For OGLE 00-BUL-12 there appears to be a possible deviation near the peak. The parameters and statistics of the single lens fits are displayed in Table 2.

We note that the maximum magnification determined for OGLE 00-BUL-12 by image subtraction, $`A_{\mathrm{max}}=159`$, is considerably higher than the value of 50 first reported by the OGLE Early Warning System. It should be noted that the first analysis was carried out using DoPHOT and that there is a significant degree of blending present in this event. This would underestimate the fitted value for the peak magnification. On the other hand, the light curve presented here is based on image subtraction photometry and should be unaffected by blending.

## 5 Planetary Modeling

We used the inverse ray shooting technique of Wambsganss (1997) to simulate lensing by planetary systems of arbitrary complexity, by firing photons from the telescope through the lens system to the plane of the source star. For every component of the lens system a deflection $`4GM/bc^2`$ was applied, where $`b`$ denotes impact parameter. Photons that hit the source star were retained; those that missed were discarded. In this way the finite size of the source star was allowed for, but the lens was treated in the thin-lens approximation. All the lens components were treated as if they were in a plane perpendicular to the ray direction. This treatment is expected to be sufficiently accurate because the dimensions of the lens system are much smaller than the distances between the lens and the observer, and the lens and the source.

The general procedure that was followed for planetary modeling was to calculate the $`\chi ^2`$ values of the data for a broad range of planetary models to find the model with the smallest $`\chi ^2`$. The light curve for any particular planetary model was calculated using the inverse ray shooting technique described above. Unless otherwise stated, the computations for all events corresponded to the following assumed values of the lens and source star parameters: $`M_L=0.3M_{}`$, $`D_{\mathrm{OL}}`$=6 kpc, $`D_{\mathrm{LS}}`$=2 kpc and $`R_S=R_{}`$. These values correspond to $`D_{\mathrm{OS}}=8`$ kpc and $`R_\mathrm{E}`$ = 1.9 AU. The coordinate system used for the computations is depicted in Fig. 10.

For the planetary modeling, the ratio of planet-to-lens mass $`ϵ`$ was initially allowed to vary from $`10^7`$ to $`10^3`$ in 33 approximately equally logarithmically spaced steps for each event. Similarly, the projected coordinates of a planet $`x_\mathrm{p}`$ and $`y_\mathrm{p}`$ at the time of peak magnification were allowed to vary from $`2R_\mathrm{E}`$ to $`+2R_\mathrm{E}`$ in 129 equally spaced steps. Thus a total of 549,153 planetary configurations was initially trialed for each event. The trialing was done on a super-cluster computer described by Rattenbury et al. (2002). For each trial, the $`\chi ^2`$ was calculated over the time interval $`t_\mathrm{E}`$ to $`+t_\mathrm{E}`$. This corresponds to the source star traversing the diameter of the Einstein ring. Maps of $`\chi ^2`$ over the lens plane for the 33 planetary mass fractions were examined to determine approximate positions and masses of possible planets. Further minimization, described below, was subsequently carried out to determine the most likely planetary models.

### 5.1 MACHO 98–BLG–35

Planetary modeling of MACHO 98–BLG–35 was carried out using the combined dataset from MOA, MPS and PLANET. A typical $`\chi ^2`$ map is shown in Fig. 11. This corresponds to an Earth-mass planet ($`ϵ=10^5`$). The minima on these maps indicate possible positions of planets at the time of the microlensing event. Maps for heavier (lighter) planets are similar, with the minima displaced further from (closer to) the Einstein ring. The Einstein ring is generally depicted clearly on these maps, either as a region where planets are strongly excluded, or where they may be present. The maps also depict a degeneracy that exists in the microlensing method, namely that planets at projected radii $`a`$ and $`1/a`$ (expressed in units of the Einstein radius) yield similar light curves, and hence are indistinguishable in this method (Griest & Safizadeh 1998).

The $`\chi ^2`$ map shown in Fig. 11 for MACHO 98–BLG–35 shows six possible positions for planets orbiting the lens star, two at position A, two at position B, and two at position C, where the doubling is caused by the degeneracy noted above. We have performed model fitting for each combination of one planet, two planet, and three planet configurations. For each model we carried out a $`\chi ^2`$ minimization allowing all parameters to vary over small ranges near their initial values. The total number of free parameters was 8 for one-planet models, 11 for two-planet models, and 14 for the three-planet models. The minimization was achieved using the Simplex procedure. The results are shown in Table 3.

Among the single planet models, model A gives the best improvement over the single lens with a renormalized $`\mathrm{\Delta }\chi ^2=63.2`$. Model B+C is the most favoured of the two planet models, although the $`\chi ^2`$ value is not a significant improvement over model A. The improvement in $`\chi ^2`$ for the triple planet model A+B+C is also not significant and as such there is no need to introduce a third planet. Models A and B+C appear as the most favourable solutions for the configuration of the planetary system. While it may seem that there is no need to introduce a two planet model, we consider model B+C as an alternative because it is physically distinct from the single planet model A. The ambiguity arises here because the intensive coverage and image subtraction photometry of the event, shown in Fig. 7, did not cover the entire full-width at half-maximum (FWHM) of the peak. Full coverage of the FWHM would reduce the ambiguities (Rattenbury et al. 2002).

We have drawn the light curve for model A in Fig. 12 and for model B+C in Fig. 13 together with the observational data. The data are plotted in terms of the fractional deviation from the single lens model due to the stellar component defined as

$$\delta =\frac{A_{\mathrm{SP}}(t_{\mathrm{max}},t_\mathrm{E},u_{\mathrm{min}},x_1,y_1,ϵ_1,\mathrm{})A_\mathrm{S}(t_{\mathrm{max}},t_\mathrm{E},u_{\mathrm{min}})}{A_\mathrm{S}(t_{\mathrm{max}},t_\mathrm{E},u_{\mathrm{min}})}$$
(5)
Here $`A_\mathrm{S}`$ is the amplification due to a single lens given by Eqns. 2–3. This depends only on the three single lens parameters $`t_{\mathrm{max}}`$, $`t_\mathrm{E}`$, $`u_{\mathrm{min}}`$. The amplification due to the star+planet system, $`A_{\mathrm{SP}}`$, depends upon the single lens parameters together with the positions and mass ratios of the component planets. It should be noted that the three parameters comprising the single lens component are fitted separately for each model. This affects the shape of the profile of the fractional deviation light curve.

In Fig. 14, we depict an additional view of model A, this time in the source plane. This shows the contour of infinite magnification for the event (i.e. the “caustic”). It is seen that the source star does not pass between the star and the planet. Hence the perturbation is negative as can be seen in Fig. 12. It is also seen that the source star does not intersect the caustic at any time.

Rhie et al. (2000) estimated a radius of 1–3$`R_{}`$ for the source star in this event. We have repeated the above $`\chi ^2`$ minimization process for model A using a source radius of 0.004$`R_\mathrm{E}`$($`2R_{}`$). We obtain best fit values for the planetary parameters of $`ϵ=1.3\times 10^5`$ at position $`(0.08,1.17)`$ with $`\mathrm{\Delta }\chi ^2=36.1`$. The effect of increasing the source size serves to partially wash out the fine structure in the microlensing light curve near the peak. The improvement in $`\chi ^2`$ is then not as good as those for smaller source radii. These results indicate a likely source star radius $`2R_{}`$.

Model B may be seen to correspond to the one previously found by Rhie et al. (2000) on the basis of DoPHOT analyses of the MPS and MOA datasets. The inclusion of the PLANET dataset has served to lower the planetary mass-fraction of this model from $`7.0\times 10^5`$ to $`2.8\times 10^5`$. Model B+C may be considered a generalization of the original model of Rhie et al. (2000) to which a second, nearby, very-low-mass planet has been added. The inclusion of a second planet as in model B+C gives a significant decrease in $`\chi ^2`$ over the single planet model B with $`\mathrm{\Delta }\chi ^2=30.0`$. While model A is the simplest, we have included model B+C to illustrate the potential of microlensing to map multi-planet systems..

Thus far the discussion has been in terms of $`\mathrm{\Delta }\chi ^2`$ values because these are both convenient for comparing different models, and because they are less sensitive than absolute $`\chi ^2`$ values to the uncertainties in the measurements. Visual inspection of Figs. 7 and 12 shows that model O is not a good fit to the data, and that model A is significantly better. Assuming that model A is actually correct, a likelihood for model O may be estimated by forcing the $`\chi ^2`$ of model A to be the same as the number of degrees of freedom by renormalizing all the errors by a constant factor.555 A similar procedure was adopted by Albrow et al. (2001). This yields $`\chi ^2=358.0`$ with 298 degrees of freedom for model O which corresponds to a deviation of $`2.5\sigma `$ for which the probability is $`<1`$%. A more satisfactory procedure for interpreting MACHO 98-BLG-35 would be to analyze all the datasets using difference imaging, determine all the errors self-consistently, and redo the planetary modelling.

For the two planetary configuration models that we favour here, model A and model B+C, the renormalized $`\mathrm{\Delta }\chi ^2`$ was around 60 and the number of degrees of freedom were around 290. Gaudi et al (2002) subsequently proposed setting an alternative detection threshold for low-mass planets, namely $`\mathrm{\Delta }\chi ^2>60`$ irrespective of the number of degrees of freedom. Howver, as the chance probability depends on the number of degrees of freedom, we advocate consideration of this factor. Gaudi et al also noted that heavy planets should be detected before light planets unless the planetary mass function is steep. We note that, in the absence of observational information on the planetary mass function, a steep function cannot be excluded. Gaudi et al analysed PLANET data for several microlensing events, including MACHO 98–BLG–35, and reported no evidence for planets. Our representation of their data in Figs. 12 and 13 is consistent with this. However, it is clear from Figs. 12 and 13 that the PLANET data, derived using profile fitting photometry, are not in conflict with the MOA data, derived using the more precise image subtraction photometry. As such the PLANET data do not rule out the planetary models we considered in our study.

Exclusion regions for giant planets orbiting the lens star of MACHO 98–BLG–35 were also computed using the inverse ray shooting technique. The $`\chi ^2`$ maps for one-planet models for mass ratios of ($`ϵ=2.8\times 10^3`$), ($`ϵ=8.5\times 10^4`$) and ($`ϵ=1.3\times 10^4)`$ were computed. These correspond to Jupiter, Saturn, and Uranus mass planets orbiting a one third solar mass star. The contours were found where $`\chi ^2`$ exceeds its value for the single lens model by 90. This value was chosen, under the assumption that any excluded model should have $`20`$ or more consecutive measurements deviating systematically by $`(12)\sigma `$ or more from it. This threshold would appear conservatice when compared with the $`\mathrm{\Delta }\chi ^2`$ values for the planetary models. The higher value was chosen to allow for the fact that the other parameters $`t_{\mathrm{max}}`$, $`t_\mathrm{E}`$, and $`u_{\mathrm{min}}`$ were not allowed to float but were fixed at the single lens model values in this computation. The exclusion regions so obtained are shown in Fig. 15. It is clear that a large region surrounding the lens star of MACHO 98–BLG–35 is devoid of gas-giants like Jupiter or Saturn. Jovian planets with projected radii $``$0.2–15 $`R_\mathrm{E}`$ or $``$0.4–30 AU are excluded. A similar result was found previously by Rhie et al. (2000), but the inclusion of further data and the use of subtraction photometry have served to enlarge the exclusions regions.

### 5.2 MACHO 99–LMC–2

The microlens event MACHO 99-LMC-2 was unusual in that it is the only high magnification event observed to date that occurred towards the Magellanic Clouds. It was found independently by both the MACHO and the OGLE groups. It was included in the observing program of the MOA group because it afforded an opportunity to search for a planet in an external galaxy, i.e. for an extra-galactic planet, under the assumption that the event is an example of the “self-lensing” process of Sahu (1994). In this process, a foreground star in the LMC lenses a background star in the LMC. We note, however, that the question of the most likely location of lenses observed towards the LMC has not yet been settled (see, for example, Alves & Nelson 2000, Alcock et al. 2001).

The typical value of the Einstein radius for self-lensing is the same as for galactic bulge events, i.e. $`2`$ AU. This follows from Eqn. 1 assuming a typical lens mass $`0.3M_{}`$ and values of order 48 kpc, 2 kpc and 50 kpc for $`D_{\mathrm{OL}}`$, $`D_{\mathrm{LS}}`$ and $`D_{\mathrm{OS}}`$. Consequently, high magnification events observed towards the LMC offer similar prospects for planet detection as they do towards the galactic bulge.

The data for MACHO 99-LMC-2 which are displayed in Fig. 8 and Table. 2 show no clear evidence for a deviation from the light curve of a single lens. Consequently, only exclusion regions for giant planets surrounding the lens were computed. This was done with the procedure used above for MACHO 98–BLG–35. The exclusion regions so obtained are displayed in Fig. 16. Jovian planets with projected radii $``$0.4–10 AU are excluded. These results are, to our knowledge, the first limits placed on planetary companions to an extra-galactic star.

### 5.3 OGLE 00-BUL-12

The data on OGLE 00-BUL-12 in Fig. 9 and Table 2 show possible evidence of a deviation from the behaviour of a single lens. Consequently, they were subjected to the same analysis as that accorded to MACHO 98–BLG–35. The initial search with the super-cluster computer yielded several possible planetary models, but none of these were statistically significant improvements over the single lens fit. However, one of them showed interesting behaviour, and this is reproduced below in Figs. 17 and 18 to illustrate the potential of the high magnification technique.

The light curve shown in Fig. 17 shows the characteristic double-spiked behaviour of a caustic crossing. This is illustrated in Fig. 18 which shows the source plane view of the model. Here a source radius of $`1R_{}`$ has been assumed. The combination of high magnification and caustic crossing occurring simultaneously in one event leads to enhanced sensitivity for planet detection. The planetary mass-fraction for the light curve in Fig. 17 ($`ϵ=0.055\times 10^5`$) corresponds to a planet of mass similar to that of Mercury orbiting a one-third solar mass star. The improvement over the single lens fit corresponds to $`\mathrm{\Delta }\chi ^2=13.7`$. We note that further data were obtained during the peak of this event (Sackett 2001). These could be analysed by the image subtraction technique and included in a future analysis of the event, thus reducing the uncertainties in the present analysis.

Given the limited coverage of this event and the low statistical significance, we do not present this as a planetary detection. However, the analysis raises interesting possibilities. The striking feature of this event is its very high peak magnification of about 160. Light curves such as that shown in Fig. 17 for a Mercury mass planet should apply to other events with similarly high magnifications. It can be seen that such planets can produce deviations from the single lens by about 1%. Such deviations should be detectable by a network of 1-m class telescopes providing continuous and complete coverage of the peaks of high magnification events.

Exclusion regions for gas-giant planets orbiting the lens of OGLE 00-BUL-12 were also computed. These are shown in Fig. 19. Jovian planets with projected radii $``$0.2–30 AU are excluded.

## 6 Discussion and Conclusion

The original proposal of Griest & Safizadeh (1998) to study extra-solar planets in gravitational microlensing events of high magnification has received support from the present work. It has been demonstrated that terrestrial-mass planets at orbital radii $`2`$ AU may be detected with 1-m class telescopes. The essential requirement is the relentless observation of the peaks of high magnification events with coordinated telescopes encircling the globe, a relatively straightforward task. The dense stellar fields in which microlensing is necessarily observed require special photometric techniques to take account of the blending of stellar images. Difference imaging analysis appears to be well suited to this task.

If a concerted effort was made by the microlensing community to detect and monitor relentlessly future events of high magnification, it could reasonably be anticipated that a first approximate measurement of the abundance of terrestrial planets could be obtained in a few years. If $`10`$ high magnification events can be detected per year, a rough measurement or upper limit on the abundance of Earth-mass planets orbiting $``$0.3$`M_{}`$ stars at projected radii $``$1.5–2.5 AU could be obtained in a few years (Bond et al 2002). As demonstrated here, one might also obtain new information on planets as light as Mercury, on extra-galactic planets, and on gas-giants. Also, multi-terrestrial-planet systems should be able to be mapped in some cases (Rattenbury et al 2002), rather like the multi-gas-giant systems presently being mapped by the radial velocity community (Butler et al. 1999). Ultimately, a dedicated wide-angle space-borne telescope such as that planned for the proposed GEST mission would provide excellent sensitivity towards all solar-system analogues except Mercury and Pluto (Bennett et al 2002).

It appears likely that none of the three lens objects studied here contains a Jovian planet. This is consistent with previous measurements by Marcy & Butler (1998) and of Albrow et al. (2001) of other systems. The dearth of Jupiter-like planets enhances the sensitivity of microlensing studies of terrestrial planets, merely through the lack of “background” signals from gas-giants.

For all three events further data by microlensing groups exist that could be analysed using image subtraction and incorporated in the planetary modeling. This could be expected to produce results of greater precision, and to diminish ambiguities that presently exist in the planetary modeling.

The assumptions made in the present analysis on the masses, sizes and distances of the lens and source stars may be able to be avoided in the future. When the Next Generation Space Telescope comes into operation, and as the lens and source stars in today’s microlensing events begin to diverge, it should be possible to measure the lens and source stellar parameters in these events by elementary photometry (Rattenbury et al 2002). This being the case, the planetary modeling would be entirely free of undetermined parameters666We have examined a 50-year old image of the region containing MACHO 98–BLG–35, but unfortunately the lens and source stars were not identifiable.

We conclude that coordinated observations using 1-m class telescopes of gravitational microlensing events of high magnification are capable of making valuable contributions to the study of planets orbiting other stars.

## Acknowledgments

We are very grateful to the OGLE collaboration for kindly making their images available to us. In particular we thank Andrzej Udalski and Karol Zebrun for providing images. We are also indebted to Dave Bennett and Sun Rhie for introducing us to the high magnification technique, Joachim Wambsganss for general discussion on microlensing, and Peter Dobcsanyi for introducing us to super-cluster computing. The Marsden Fund of New Zealand and the Ministry of Education, Science, Sports and Culture of Japan are thanked for the financial support that made this work possible.