format: text

text:
# SZ and X-ray combined analysis of a distant galaxy cluster, RX J2228+2037

## 1 Introduction

The study of galaxy clusters is one of the clues in our understanding of the universe. Located at the nodes of large-scale structures, they are the object of many different cosmological studies. While it is important to study their statistical properties, it is also important to probe their internal physics in the context of structure formation in the universe. As the main baryonic component of galaxy clusters, the intracluster gas is a rich source of information. It represents between 10 and 30% of the cluster total mass. Largely observed in X-ray through its Bremsstrahlung emission, it can also be detected from submillimeter to centimeter wavelengths. In the last case, the signal is due to the Sunyaev-Zel’dovich (SZ) effect (Sunyaev & Zel’dovich (1972)). The intracluster gas is highly ionized because of its high thermal temperature. The thermal electrons accelerate the low energy photons of the cosmic microwave background (CMB) by inverse Compton scattering. The transfer of energy from the thermal electrons to the photons produces a shift in the CMB spectrum towards higher energies. This characteristic distortion of the CMB spectrum is a well-known signature of a galaxy cluster (Birkinshaw (1999)). Since the SZ flux is independent of the redshift, the study of the SZ effect toward distant clusters is a powerful tool to probe massive clusters at high redshift and to constrain the cosmological density parameter.

The single dish measurements of the SZ effect allows us to determine the absolute intensity of the effect from centimeter to submillimeter wavelengths. When performed on large radio telescopes, they also are the only way to get high resolution SZ observations. Recent SZ studies at millimeter wavelengths using such observing methods and coupled with the use of high performance bolometric detectors have led to significant results. Two independent groups, Pointecouteau et al. (pointecouteau99 (1999, 2001)) with the IRAM-30m$`/`$DIABOLO instrument and Komatsu et al. (komatsu99 (1999, 2001)) with the NRO-45m$`/`$NOBA instrument, have shown surprising results in the direction of the same target, the most luminous X-ray cluster known to date, RX J1347-1145. Discrepancies between the SZ and the X-ray (ROSAT) signal distributions have been pointed out, as well as evidence for an internal structuring of the intracluster gas. The existence of an internal structure is confirmed by the Chandra observation of RX J1347-1145 (Allen et al. (2001)).

In this paper, we present the result of a combined SZ and X-ray analysis in the direction of RX J2228+2037, a distant galaxy cluster ($`z=0.421`$). It was detected in the ROSAT All Sky Survey (RASS). An optical identification of the X-ray source was given by Bade et al. (bade98 (1998)). Voges et al. (voges99 (1999)) identified this source as a galaxy cluster. More recently it has been included in the NORAS galaxy cluster survey (Böhringer et al. (2000)) and in the ROSAT Brightest Cluster sample (Ebeling et al. (2000)). Both catalogs are based on the RASS observations in the \[0.1-2.4\] keV energy band. The characteristics reported for this target in these two papers are in agreement. The results of the NORAS catalog are: a count rate of $`0.2`$ cts s-1, an estimated flux for a temperature of 5 keV of $`F_X=4.2\times 10^{12}`$ ergs s-1 cm-2 and an intrinsic luminosity of $`L_X=20.7\times 10^{44}h_{70}^2`$ ergs s-1.

We have observed RX J2228+2037 through its SZ effect using the Nobeyama Radio Observatory (NRO) 45 m radiotelescope at 21 GHz and combined those data with the X-ray image obtained from the ROSAT/HRI instrument in the \[0.1-2.4\] keV energy band. In the following, we present an extensive study of this target through the combined analysis of the two sets of data. In a first step, we probe the spatial structure of the cluster with the X-ray data (see Sec. 3). Then, we use the direct dependency of the SZ flux with respect to the gas temperature to extract this key parameter and furthermore we deduce the gas mass and the total mass of RX J2228+2037 (see Sec. 4).

Throughout this paper, we used $`H_0=70h_{70}`$ km s-1 Mpc-1, $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. In this cosmology an angular scale of 1 arcmin corresponds to a physical size of 328 kpc at the redshift of RX J2228+2037.

## 2 Observations

### 2.1 The X-ray data

Apart from the RASS data, pointed observations have been conducted in the direction of RX J2228+2037. A total of 30 ks exposure time has been spent with ROSAT/HRI on this target on the 8th and the 12th of December 1996. The pointing coordinates were $`\alpha _{2000}`$=22h28m36s and $`\delta _{2000}`$=+203712′′.

We have processed the data using the EXSAS software (Zimmermann et al. (1993)). In order to limit the background events, we checked each PHA individually. We only selected PHA displaying a significant signal for inclusion in the data processing. Finally, PHA from 1 to 10 were selected. The events were binned into 5 arcsec pixels to generate a raw image of the cluster. The cluster map is presented in Fig 1a.

### 2.2 The SZ data

The SZ effect observations were conducted at the Nobeyama Radio Observatory in December 1999 using the 45 m radio telescope facilities at 21 GHz. The 21 GHz (1.43 cm) detector is a double channel equiped with low-noise HEMT receivers coupled to 2 GHz bandwidth back-ends. Basically the two channels allow detection of a circularly polarized signal. However, we have used them as independent detectors for the same signal. This redundancy allows us to monitor the self consistency of the signal and to avoid unexpected drift and systematics in the signal measurement. It also increases the reliability of our measurements. The beam size is given by the diffraction limit of the telescope at this wavelength: $`\theta _{FWHM}80`$ arcsec, which corresponds to a physical angular scale of 437 kpc at the cluster redshift. During the whole observation, the system temperature was approximately 135 K.

Observing at the 45 m telescope at 21 GHz, we were looking for the extended SZ signal, with the goal of mapping the distribution of the SZ decrement at large angular scales. Our observations covered an area of about $`6\times 6`$ arcmin2. We directly mapped the target in the RA–DEC coordinates referential by using a raster/scanning technique: a single observation is composed by a X raster (the scanning is done in the right ascension direction) and a Y raster (the scanning is done in the declination direction). Their combination is expected to cancel out the scanning effects. Each single map (either a X or a Y raster) had a basic scanlength of 360 arcsec for a total of 9 lines. Two lines are separated by 40 arsec in order to achieve a correct spatial sampling with respect to the beam size. We have monitored the atmospheric emission and done pointing verifications in the direction of point sources lying near the cluster location. Calibration measurements have been done throughout the observing run in the direction of NGC7027 (Ott et al. (1994)). The main steps of the data processing pipeline include a baseline subtraction (a first order polynomial is fitted to the edges of each line of the X and the Y rasters). The signal is cleaned of glitches due to cosmic high energy particles and to the electronic spikes in the readout systems. After the calibration of the signal, the X and the Y rasters are combined using a Fourier transform filtering method to remove the scanning noise (Emerson & Graeve (1988)). Each raster is reprojected on a RA-DEC grid to avoid a difference in the central position between observing sequences. Finally, the SZ map is computed from the sum of the whole set of rasters.

The resulting map is presented in Fig 1b. The sensitivity level we have reached is 0.5 mJy beam-1 for an integration time of 34 h.

## 3  Spatial analysis

### 3.1 X-ray analysis

To characterize the intrinsic emission of RX J2228+2037, we first cleaned the field of view of contaminating sources. The sources were extracted using a maximum likelihood method with a threshold of 8, which roughly corresponds to a 3.5$`\sigma `$ detection level. A total of 27 sources were identified above this level and then removed. The hydrogen column density was fixed to the galactic value $`4.68\times 10^{20}`$ cm-2 (Dickey & Lockman (1990)). The count rate for the cluster at the virial radius is $`0.053\pm 0.003`$ cts s-1.

The X-ray map seems to present an extension in the right ascension direction. Moreover, the shape of the cluster seems to be slightly elliptical and tilted. For this reason, we have chosen to fit different models of the intracluster gas distribution to the data. We checked four different models based on a modified King profile, or so-called $`\beta `$-model (Cavaliere & Fusco-Femiano (1978)) and refered to them as A, B, C and D. Our models also assumed that the intracluster medium is isothermal.

Models A and B are 2D elliptical models and therefore they have been fitted to the X-ray map. Model A is an elliptical $`\beta `$-model. The centroid of this model has been fixed to the position of the X-ray maximum. Model B is also an elliptical $`\beta `$-model, but in this second case, the position of the centroid is left as a free parameter.We apply the following expression of the X-ray surface brightness in radial coordinates $`(\theta ,\varphi )`$:

$$\begin{array}{cc}S_X(\theta ,\varphi )\hfill & =S_X(0)(1+(\mathrm{cos}^2(\varphi \mathrm{\Phi })\hfill \\ & +ϵ^2\mathrm{sin}^2(\varphi \mathrm{\Phi }))\left(\frac{\theta }{\theta _c}\right)^2)^{3\beta +\frac{1}{2}}+b_X\hfill \end{array}$$
(1)
$`S_X(0)`$ is the central X-ray surface brightness and $`b_X`$ the background value. $`\theta _c`$ is the core radius along the major axis of the cluster, and $`ϵ`$ is its ellipticity, defined as the ratio of the major axis core radius over the minor axis core radius. $`\mathrm{\Phi }`$ is the orientation angle, the angle between the cluster major axis and the right ascension axis. In model A, $`S_X(0)`$, $`\theta _c`$, $`\beta `$, $`ϵ`$, $`\mathrm{\Phi }`$ and $`b_X`$ are fitted simultaneously. In model B, the centroid coordinates $`\alpha _{cen}`$ and $`\delta _{cen}`$ are additional free-fit parameters. The 27 contaminating sources have been masked to avoid any bias in the fitting procedure.

For those two models, we strictly followed the 2D fitting procedure of an elliptical $`\beta `$-model described by Neumann & Böhringer (neumann97 (1997)): we smoothed the image of the cluster with a Gaussian filter with $`\sigma =5`$ arcsec . We reproduced this Gaussian filtering on the 2D elliptical $`\beta `$-model. We calculated the corresponding uncertainty taking into account the Gauss filter. The uncertainties on the best fit parameters were obtained via a Monte Carlo method: a random Poisson noise is added to the data. This “modified” data set is fitted the same way as the real data. We reproduced this procedure 100 times for different realizations of the Poisson noise. We took the value of the standard deviation of the 100 fitted values as the uncertainty of each parameter.

Even though this method does not allow us to give the exact value of the $`\chi ^2`$, it allows us to minimize it and it has already shown good results for CL 0016+16, A2218 or 3C295 (see Neumann & Böhringer neumann97 (1997, 1999); Neumann neumann99b (1999)).

Models C and D are “Classical” 1D $`\beta `$-models. They have been fitted to the radial X-ray surface brightness profile using a maximum likelihood method. Model C is centered on the location of the X-ray maximum, and model D has been fitted to the profile folded from the cluster centroid position ($`\alpha _{cen}`$, $`\delta _{cen}`$) determined by model B. For a 1D $`\beta `$-model, the X-ray surface brightness is:

$$S_X(\theta )=S_X(0)\left(1+\left(\frac{\theta }{\theta _c}\right)^2\right)^{3\beta +\frac{1}{2}}+b_X$$
(2)
The radial profiles have been extracted from the X-ray map within a 12 arcmin circular area from the position chosen as the cluster center. Within this area, the X-ray background can be considered as a constant (Ota (2000)). The events of the 27 contaminating sources have been masked, so that the resulting radial profile is free from any contamination (see Fig 2a). The two models include 4 free parameters: $`S_X(0)`$, $`\theta _c`$, $`\beta `$ and $`b_X`$.

Our four models have been tested against the X-ray data. In each case the PSF effects have been taken into account. The results are summarized in Table 1. The best fit parameters are given with their $`1\sigma `$ errors bars. For models A and C, the location of the cluster center, which corresponds to the position of the maximum of the X-ray emission, is: $`\alpha _{cen}=`$22h28m33.3s, $`\delta _{cen}=`$203712.2′′. In the case of models B and D, this position is determined by the best fit position for the centroid of the elliptical $`\beta `$-model in model B: $`\alpha _{cen}=`$22h28m32.7s, $`\delta _{cen}=`$203711.2′′. The core radii for the different models (A, B, C, D) as well as the $`\beta `$ parameters agree within the 2$`\sigma `$ limit. In case of the 1D models, model C presents the best quality fit. In the following, we have used model C as the reference model.

### 3.2 SZ analysis

Before combining the SZ data with the X-ray data, in order to extract the intracluster gas temperature, we decided to check whether our SZ data were consistent with the results obtained from the spatial analysis of the X-ray data.

If the cluster gas distribution follows a modified King profile, then the SZ surface brightness can be expressed as:

$$S_{SZ}(\theta )=S_{SZ}(0)\left(1+\left(\frac{\theta }{\theta _c}\right)^2\right)^{\frac{3}{2}\beta +\frac{1}{2}}+b_{SZ}$$
(3)
where $`S_{SZ}(0)`$ is the central SZ surface brightness and $`b_{SZ}`$ the background emission level.

Unfortunately, the quality of the SZ data is not good enough to allow us to test an elliptical model, so that we limited our analysis to a symmetrical $`\beta `$-model. Using this distribution to model the SZ signal, we have reproduced in it each step of the SZ data processing (convolution with the beam pattern, baseline subtraction, etc) in order to get the most accurate model. Using a maximum likelihood method, we left the SZ central flux, the core radius and the $`\beta `$ parameter as free parameters and fitted the model over the SZ map. The best configuration parameters are obtained for a reduced chi-squared value of 1.1 (96 d.o.f): $`S_{SZ}(0)=3.9\pm 1.5`$ mJy beam-1, $`\theta _c=20.3\pm 3.1`$ arcsec and $`\beta =0.88\pm 0.19`$. The result for $`\theta _c`$ and $`\beta `$ are compatible with the best fit parameter values obtained from the X-ray data analysis within a 1$`\sigma `$ and a 2$`\sigma `$ limit respectively. The large errors (especially in the case of $`\beta `$) are due to the low signal to noise ratio and to the lack of resolution of our SZ data at 21 GHz. This made the precise determination of the core radius and $`\beta `$ very difficult and explained the strong degeneracy between those parameters.

## 4 Physical properties of RX J2228+2037

Knowing that the angular diameter distance of RX J2228+2037 is $`d_A=1127`$ Mpc, the value of the angular core radius, $`\theta _c=18.9\pm 2.3`$ arcsec, corresponds to a physical core radius of $`r_c=103\pm 12h_{70}^1`$ kpc. From the best fit parameters of model C, we have extracted the central electron density, $`n_e(0)=0.014h_{70}^{1/2}`$ cm-3.

### 4.1 Temperature determination

ROSAT/HRI has no spectral resolution so that the X-ray data do not allow us to determine the electron temperature. However, using the high luminosity value in the \[0.1-2.4\] keV energy band reported for RX J2228+2037 in the NORAS catalog $`L_X=20.7\times 10^{44}h_{70}^2`$ ergs s-1(Böhringer et al. (2000)), we can expect a quite high temperature for the intracluster medium. In fact, we can use this value of the luminosity as a lower limit for the bolometric luminosity of the cluster to put a lower limit on the gas temperature. Using the $`L_XT_e`$ relation, we obtained $`k_BT_e>7.7`$ keV (Arnaud & Evrard (1999)) or $`k_BT_e>6.6`$ keV (Fairley et al. (2000)).

Knowing that the SZ signal is directly proportional to the pressure and furthermore to the gas mass integrated along the line of sight, we have used the spatial information of our SZ map to constrain the gas temperature of RX J2228+2037. The central SZ surface brightness can be expressed as follow :

$$S_{SZ}(0)=yf_{SZ}(\nu ,T_e)$$
(4)
where $`f_{SZ}(\nu ,T_e)`$ represents the SZ effect spectral shape. The dependency of this spectral shape with respect to the gas temperature should not be neglected, especially in our case where we expect a high value for the temperature. This dependency is not linear and $`f_{SZ}(\nu ,T_e)`$ has to be computed numerically (Pointecouteau et al. (1998)). $`y`$ is the Comptonization parameter:

$$y=\frac{k_BT_e}{m_ec^2}\sigma _Tn_e(0)_{R_v}^{+R_v}f_g(l)𝑑l$$
(5)
where $`k_B`$, $`m_e`$, $`c`$ and $`\sigma _T`$ are respectively the Boltzmann constant, the electron mass, the speed of light and the Thomson cross-section. $`n_e(0)`$ is the electron density at the center of the cluster. $`R_v=r_{178}=2.9`$ Mpc is the virial radius of the cluster (Bryan & Norman (1998)). $`f_g(l)`$ represents the gas distribution profile along the line of sight.

To determine the electron temperature from the SZ data, we have fixed the $`\beta `$-model parameters to the best fit values of model C, $`\theta _c=18.9`$ arcsec and $`\beta =0.54`$ (see Sec. 3.1). $`n_e(0)`$ has also been fixed so that the temperature was the only remaining free parameter in the model. We have tested this $`\beta `$-model against the SZ map for various temperature values. Taking into account the temperature effect on the SZ spectral shape, we have derived from a maximum likelihood analysis a value of $`k_BT_e=10.4\pm 1.8h_{70}^{1/2}`$ keV ($`\chi ^2=1.02`$). This high temperature value is in agreement with the lower limit value suggested by the \[0.1-2.4\] keV X-ray luminosity. The corresponding central Comptonization parameter is $`y(0)=2.4\times 10^4`$.

From this determination and the X-ray count rate estimates (see Sec. 3.1), we deduced the X-ray flux and luminosity of the cluster in the \[0.1-2.4\] energy band: $`F_X=2.7\times 10^{12}`$ ergs s-1 cm-2 and $`L_X[0.12.4]=16.5\times 10^{44}h_{70}^2`$ ergs s-1. The X-ray bolometric luminosity has been computed from the best fit model: $`L_X^{bol}=58.3\times 10^{44}h_{70}^2`$ ergs s-1. This allows us to check our determination of the temperature against the $`L_XT_e`$ relation expectations: $`k_BT_e=11.9`$ keV (Arnaud & Evrard (1999)) and $`k_BT_e=9.9`$ keV (Fairley et al. (2000)). Our temperature determination fits perfectly on the $`L_XT_e`$ relation. This determination confirms the fact that RX J2228+2037 is indeed a hot cluster.

### 4.2 Point sources

Where the determination of the spatial parameters of the cluster ($`r_c`$ and $`\beta `$) is provided by the X-ray data, the electron temperature is constrained by the SZ data. Within the hypothesis of the hydrostatic equilibrium, the SZ flux and the electron temperature are linearly linked, so that any contamination in the SZ flux measurement will bias the temperature estimation. Such contamination could be inferred by the presence of point sources in the field of view. To prevent such biases, we have investigated for the presence of point sources in our field of view. We have found three radio sources in the NVSS survey (Condon et al. (1998)), that could contaminate the SZ mapping of the cluster. We only know their flux measurements at 1.4 GHz, therefore we are not able to determine their spectra in order to extrapolate their flux at 21 GHz without a few hypothesis. Assuming that the radio emission of the sources is due to synchrotron emission, we can adopt a power law as a model for their spectra. To date many different radio sources have been observed with different spectral shapes. Nevertheless, most of the spectral indexes of models extracted from radio source analysis fit into the interval: $`1<\alpha _r<0.5`$. Using those two limits we have extrapolated the contamination sources fluxes from 1.4 GHz to 21 GHz using a power law spectrum: $`F(\nu )\nu ^{\alpha _r}`$. In order to estimate whether the presence of radio sources affects the determination of $`T_e`$, we have added the point source contribution to the SZ model and performed a new temperature estimation. In a first step, each point source flux was fixed to the value extrapolated from the 1.4 GHz measurements with the two spectral index values of $`0.5`$ and $`1`$. In a second step, we considered their fluxes as free parameters as well as the electron temperature. The characteristics of the radio sources and the results on the temperature determination are summarized in table 2.

In each configuration, the temperature determination is not really affected. The three temperature values are consistent with the value obtained when ignoring the point sources. The $`\chi ^2`$ value is bad in the case of $`\alpha =0.5`$. On the contrary, the fit for a spectral index of $`1.0`$ is satisfying and also agrees very well with the fluxes value we determined from the SZ data (third case) within a 2$`\sigma `$ limit. This can be explained by the low fluxes expected or fitted at 21GHz and by the relatively high noise level on the SZ map. Moreover, in this case the point sources produce a localized effect which is washed out in the temperature determination over the entire map. It appears that the influence of the point sources present in the cluster neighborhood can be neglected in the temperature extraction. Therefore the value of $`k_BT_e=10.4\pm 1.8h_{70}^{1/2}`$ keV can be taken as the isothermal temperature for RX J2228+2037.

### 4.3 Mass estimation

The SZ surface brightness is directly proportion to the pressure integrated along the line of sight. A map of the SZ effect is a direct projected image of the cluster gas mass distribution if the temperature distribution is isothermal. Knowing the gas temperature, the projected gas mass is derived from the SZ flux through a linear transformation : $`M_g(\theta )=C(z,T_e)\times F_{SZ}(\theta )`$. Where $`C`$ is a constant depending on the cluster redshift and on the gas temperature. From our 21 GHz SZ data, we are able to measure the gas mass up to a radius of 0.5 Mpc: $`M_{gas}^{SZ}(r<0.5\mathrm{Mpc})=(2.7\pm 0.6)\times 10^{13}\mathrm{h}_{70}^2\mathrm{M}_{}`$.

Assuming the isothermality of the cluster, we can derive the gas mass and the total mass from our best fit model C to the X-ray data (see Sec. 3.1) :

$$\{\begin{array}{c}M_{gas}=\mu m_p𝑑Sn_e(l)𝑑l\hfill \\ M_{tot}(r)=\frac{k_B}{\mu m_pG}rT_e\left(\frac{d\text{ln}\rho }{d\text{ln}r}+\frac{d\text{ln}T_e}{d\text{ln}r}\right)\hfill \end{array}$$
(6)
where $`\mu m_p`$ is the average weight per massive particle and $`G`$ the gravitational constant. The measured gas mass from the SZ map is in good agreement with the mass predicted by the best fit model C: $`M_{gas}(r<0.5\mathrm{Mpc})=2.9\times 10^{13}\mathrm{h}_{70}^{5/2}\mathrm{M}_{}`$.

Finally, we can estimate the gas mass and the total mass at the virial radius by extrapolating the profile corresponding to the best fit parameters out to the virial radius. Using our best fit model C which parameters are $`\theta _c=18.9`$ and $`\beta =0.54`$ (Sec. 3.1), we obtained the following extrapolated measurements: $`M_{gas}(r<R_v)=(4.0\pm 0.2)\times 10^{14}h_{70}^{5/2}\mathrm{M}_{}`$ and $`M_{tot}(r<R_v)=(1.8\pm 0.4)\times 10^{14}h_{70}^1\mathrm{M}_{}`$.

The errors on the mass determination take into account the uncertainties derived from the best fit parameters. However, our results using $`r_c`$ and $`\beta `$ for the 1D and 2D profiles do not account for a possible oblate or prolate form in the direction of the line-of-sight. Taking these effects into account would certainly increase the error bars. However, such analysis is beyond the scope of this paper.

The corresponding gas fraction is $`f_{gas}=0.22\pm 0.06h_{70}^{3/2}`$. This result is larger than the value derived by Grego et al. (grego01 (2001)) from the SZ measurements of 18 clusters, $`f_g=0.12\pm 0.01`$, but both values agree within the 2$`\sigma `$ limit.

## 5 Conclusion

In this paper we have presented a combined analysis of the SZ effect and X-ray mapping observations. Through this analysis, we have carried out an extensive study of the galaxy cluster RX J2228+2037.

This work led us to the characterization of the cluster physical properties, among which the gas temperature is a key parameter. Our estimation of this parameter, $`k_BT_e=10.4\pm 1.8h_{70}^{1/2}`$ keV, has only been possible through the SZ and X-ray data combination which shows the power of combining this sort of data.

This target will be soon observed by the new generation X-ray satellites, XMM-Newton and Chandra. The high quality spectral data that will be provided by those instruments should confirm and help to improve our results.

Our work also emphasizes how powerful the combination of SZ and X-ray data is. It is a preview of what we can expect for many sources when the full sky coverage of the Planck Surveyor satellite from submillimeter to centimeter wavelengths will be available. Many sources in this catalog are expected to have faint X-ray counterparts. The quality of the X-ray spectral data is expected to be poor in such cases so that the combination with SZ data will be the only path to the extraction of reliable information about the physics of galaxy clusters.

###### Acknowledgements.

The authors are grateful to F. J. Castander, the referee, whose comments helped to improve the quality of this paper. We want to thank the NRO staff members for their help during the observations and T. Kitayama for providing support for the data processing. EP acknowledges the support of the Japanese Society for Promotion of Science.


text:
# Exchange Monte Carlo for Molecular Simulations with Monoelectronic Hamiltonians

## Abstract

We introduce a general Monte Carlo scheme for achieving atomistic simulations with monoelectronic Hamiltonians including the thermalization of both nuclear and electronic degrees of freedom. The kinetic Monte Carlo algorithm is used to obtain the exact occupation numbers of the electronic levels at canonical equilibrium, and comparison is made with Fermi-Dirac statistics in infinite and finite systems. The effects of a nonzero electronic temperature on the thermodynamic properties of liquid silver and sodium clusters are presented.

In the recent years, the physics of materials and complex systems has undergone awesome development in the field of molecular simulation molsim . Significant achievements include ab initio, Car-Parrinello molecular dynamics (MD) cpmd , linear-scaling linescale , or progresses in ergodic techniques such as exchange Monte Carlo (EMC) ptmc , also known as parallel tempering, or the Wang-Landau algorithm wanglandau . Account of electronic structure is often made either implicitely with empirical potentials or more explicitely through one-electron methods, either in the framework of density-functional theory (DFT) or using tight-binding (TB) approximations. In the latter cases, the electronic ground state is described via integer occupation numbers. This is usually adequate for nonmetal systems. However the electronic states of metals can cross the Fermi surface, therefore integer occupations are not appropriate for continuous dynamics. Fractional occupation numbers must also be introduced for insulators or semiconductors, provided that the temperature is high enough for the lowest excited states to be populated.

Building upon a seminal paper by Mermin mermin  who extended the Hohenberg-Kohn theorem of DFT to nonzero electronic temperatures, several authors proposed to combine DFT and fractional occupation numbers weinert ; alavi ; springborg . Most of these works were devoted to metals, and they tried to find suitable forms of the occupation laws for numerical stability purposes and a better sampling in the Brillouin zone. In particular, it was proposed to replace the Fermi-Dirac (FD) function with other expressions springborg .

Even when keeping a FD distribution, the best choice for electronic temperature was shown not to be necessarily related with the nuclear vibrational temperature mehl . Additionally, and strictly speaking, the FD function holds for an infinite system, but should only be considered as an approximation when treating a small molecular system such as a metal cluster.

The goal of this Letter is to show how the true canonical equilibrium of both ionic and electronic degrees of freedom can be simulated through MC methods, for small or large sizes. A specific interest in choosing MC methods over MD simulations is their greater flexibility and wider range of application. They are very convenient for discrete systems, and they can be adapted using statistical biases to accelerate convergence. They also offer a straightforward way to sample grand-canonical ensembles, which is more difficult with MD molsim . As seen below, the Monte Carlo method is well suited to the problem of electronic thermalization, especially for finite systems.

We consider the general class of materials modelled by monoelectronic Hamiltonians. At any given nuclear configuration $`𝐑`$, the total energy $`E`$ depends on the occupation numbers $`\{n_i\}`$. In the Kohn-Sham formalism, these numbers appear explicitely in the expression of the density, hence in the energy. In TB models, the band contribution is the weighted sum of the one-electron energies $`\{\epsilon _i\}`$: $`E(𝐑)=_in_i\epsilon _i(𝐑)`$. The first, most simple MC algorithm consists in treating the nuclear $`(𝐑)`$ and electronic $`(𝐍)`$ degrees of freedom on the same footing, by performing random moves on the generalized coordinates $`𝐐=(𝐑,𝐍)`$. Here $`𝐍=\{n_i\}`$ is the set of instantaneous integer occupation states, which evolve during the simulation accordingly with the level statistics. In the Metropolis scheme, a change from $`𝐐_{\mathrm{old}}`$ to $`𝐐_{\mathrm{new}}`$ due to a change in either $`𝐑`$ or $`𝐍`$ is accepted with probability $`\mathrm{acc}(𝐐_{\mathrm{old}}𝐐_{\mathrm{new}})=\mathrm{min}[1,\mathrm{exp}(\beta \mathrm{\Delta }E)]`$ where $`\mathrm{\Delta }E=E(𝐐_{\mathrm{new}})E(𝐐_{\mathrm{old}})`$. The MC moves involving $`𝐍`$ must keep constant the total number $`𝒩`$ of occupied states. The only moves of this kind that we consider are single exchanges between occupied and unoccupied, neighboring levels. For a given ionic geometry $`𝐑`$, this algorithm clearly converges to the electronic canonical distribution in the ergodic limit.

The local moves involving $`𝐑`$ and $`𝐍`$ are now implemented in the framework of generalized ensembles. In our case, we use exchange Monte Carlo by performing simultaneous simulations at various temperatures. With some probability $`p`$, exchange moves between $`𝐐_i=(𝐑_i,𝐍_i)`$ and $`𝐐_j=(𝐑_j,𝐍_j)`$ at the inverse temperatures $`\beta _i`$ and $`\beta _j`$, respectively, are attempted. They are accepted with the probability

$$\mathrm{acc}(𝐐_i𝐐_j)=\mathrm{min}[1,\mathrm{exp}(\mathrm{\Delta }\beta \mathrm{\Delta }E)].$$
(1)
In this equation, $`\mathrm{\Delta }\beta =\beta _i\beta _j`$ and $`\mathrm{\Delta }E=E(𝐐_i)E(𝐐_j)`$. With probability $`1p`$, the usual local moves are attempted for all trajectories.

We test this algorithm on a model liquid silver, described by a TB Hamiltonian tbag . In Fig. 1, we have plotted the average radial density $`g(r)`$ at $`T=5000`$ K, from three different methods, for a system with 108 Ag atoms at constant density $`\rho =10.510^3`$ kg.m-3. The first method only considers nuclear moves, the electrons being frozen at $`T_e=0`$. The second one is the previous MC method with adjacent trajectories at 2500 K and 7500 K, and $`p=0.1`$.

Finally, we perform a single trajectory MC calculation using only nuclear moves, but fractional occupation numbers given by the Fermi-Dirac distribution at $`T_e=T=5000`$ K. For each atomic configuration, the chemical potential which normalizes $`𝒩`$ is found by a Newton-Raphson minimization of the error function $`\chi ^2(\mu )=[_kn_k^{\mathrm{FD}}(\mu )𝒩]^2`$ with $`n_k^{\mathrm{FD}}(\mu )=[1+\mathrm{exp}(\beta (\epsilon _k\mu ))]^1`$, starting from the distribution $`𝐍`$ of the previous MC step. In this case, an entropic correction ($`TS`$) to the energy is included with the usual form given by Mermin mermin :

$$S=k_\mathrm{B}\underset{k}{}[n_k\mathrm{ln}n_k+(1n_k)\mathrm{ln}(1n_k)].$$
(2)
All simulations consist of $`10^7`$ nuclear MC cycles. In the second method, $`𝒩`$ electronic moves are attempted for each vibrational move.

The good agreement between the two calculations with $`T_e0`$ shows that the present MC algorithm is able to equilibrate both the electronic and nuclear degrees of freedom, for a quite large system whose electronic statistical distribution can be safely represented by a Fermi-Dirac distribution at the same temperature.

We now turn to finite atomic metal clusters. Finding the electronic average occupation numbers for any given nuclear geometry is a combinatorial task, which cannot be solved exactly in the canonical ensemble, except for very few $`(20)`$ levels. In the bulk limit, the grand-canonical ensemble is relevant and gives the FD distribution. For finite, intermediate sizes, the previous MC algorithm can solve this problem numerically. However, as is well known for discrete spin systems, exchange “flip” moves between occupied and unoccupied neighboring states can be expected to be mostly rejected at low temperatures $`T<\mathrm{\Delta }E/k_\mathrm{B}`$, where $`\mathrm{\Delta }E=E_{\mathrm{LUMO}}E_{\mathrm{HOMO}}`$ is the energy gap between the lowest unoccupied and highest occupied orbitals. The convergence can be greatly accelerated using the kinetic Monte Carlo (KMC) method bkl  of Bortz, Kalos and Lebowitz. We use the KMC method in conjunction with the local exchange moves, starting from the $`T_e=0`$ electronic distribution where only the lowest levels are populated. Tests on small clusters (up to 16 atoms) have shown that about $`10^5`$ KMC steps are necessary to ensure convergence towards the exact statistical

population. In Fig. 2 we plot the average occupation number at $`T=500`$ K versus level energy for the ground state geometry of Na40 described by a TB Hamiltonian tbna . While the general shape is that of a Fermi-Dirac type, the actual FD distribution at electronic temperature $`T_e=500`$ K does not match the result of the KMC calculation, which is best fitted by a FD law at effective temperature $`T_e^{\mathrm{eff}}=335`$ K.

We can combine the KMC algorithm for the electrons with the usual MC moves for the nuclei. To save computational time, the electronic problem is solved periodically, once every $`M`$ steps. Again, the EMC strategy is used to improve global convergence and reduce quasi-ergodicity. However, one must be careful when attempting exchange moves between trajectories at different temperatures, because the energy depends explicitely on the nuclear coordinates, but also on temperature via the average occupation numbers. The same problem would hold for any other temperature-dependent potential, such as the effective potentials with quantum corrections used in liquids theory jph .

More precisely, the acceptance probability of an exchange between configurations $`𝐑_i`$ and $`𝐑_j`$ initially at the respective inverse temperatures $`\beta _i`$ and $`\beta _j`$ is now

$$\mathrm{acc}(𝐑_i𝐑_j)=\mathrm{min}[1,\mathrm{exp}(\beta _j\mathrm{\Delta }E_j+\beta _i\mathrm{\Delta }E_i)],$$
(3)
where $`\mathrm{\Delta }E_k=E(𝐑_j,\beta _k)E(𝐑_i,\beta _k)`$, $`k=i`$ or $`j`$. When the energies are temperature-independent, Eq. (1) is recovered. However, in the present case, they include the entropic correction of Eq. (2) and should be considered as free energies. The canonical equilibrium of both nuclear and electronic coordinates can be simulated using this MC method. Since we are dealing with temperature-dependent energies, useful analysis techniques such as the histograms methods histo  cannot be applied here. In addition, corrective terms to the thermodynamic properties appear due to this dependence. Unfortunately, because the occupation numbers are obtained numerically, their derivatives with respect to $`\beta `$ are hard to get. We can reasonably assume that the corresponding effects are small at low temperature, since in bulk metals the electronic heat capacity is only a small quantity with a weak (linear) dependence upon temperature solid .

The previously described MC method is used to simulate the solidlike-liquidlike phase change in small sodium clusters. Recent theoretical works calspi1  emphasized the significant role of geometric (nuclear) effects on the caloric curves. However, the possible effects of electrons thermalization have not been considered yet. The influence of electronic temperature on shape was discussed by Yannouleas and Landman landman , but the jellium model used by these authors does not provide information about phase changes within an atomistic description. From the experimental point of view, the results by Haberland and coworkers hab  are still far from a complete understanding, especially the complex size-dependence of the thermal properties.

We first compare in Fig. 3 the heat capacity curves of Na8 obtained without consideration of electron temperature, nor use of sophisticated sampling method calspi2 , with the ones obtained with the present Monte Carlo algorithms. Here we can compute the exact average occupation numbers by solving the combinatorial problem for each configuration, or we can employ the numerical KMC procedure, both within the EMC scheme. 35 simultaneous trajectories were propagated with $`10^7`$ cycles each, the occupation numbers being calculated for each configuration. For comparison, the curves obtained at $`T_e=0`$ but with EMC moves are also plotted on Fig. 3. The resulting heat capacities show a very good agreement between the two present simulations with nonzero electronic temperature. Thus the KMC scheme provides a good approach to electronic thermalization. The thermodynamic curves obtained without using exchange MC, or assuming frozen electrons, are very similar. Therefore, in this case, we can conclude that (i) the electronic temperature plays only a small role; and (ii) EMC does not bring such an improvement over conventional Monte Carlo.

This situation becomes somewhat different for larger sizes. Na20, Na40, Na$`{}_{}{}^{+}{}_{59}{}^{}`$ and Na$`{}_{}{}^{+}{}_{93}{}^{}`$ are studied using the same

methods except the one involving the exact calculation of the fractional numbers. Again, we also compare the heat capacities with our previous calculations. The whole results are shown in Fig. 4. The overall behavior resembles much that of Ref. calspi2, , with a major peak in the heat capacity which marks the onset of the solidlike-liquidlike phase change. However, two significant differences can be noted. First, the melting peak appears rather clearly, without any strong premelting feature (shoulder or peak) at low temperature. This must be contrasted with most other theoretical studies calspi1 ; calspi2 ; alonsoblundell  which emphasized multistep melting in sodium clusters described by various models, but is consistent with experiments hab . Second, the melting temperature indicated by the top of the peak is shifted to lower temperatures by about 10–30 K depending on size. Again this brings the present results closer to experiments, as our previous ones with the same tight-binding model were seen to overestimate melting points in these clusters calspi2 . Electronic temperature does not have such a large effect, since the Monte Carlo results with $`T_e=0`$ nearly match the ones with $`T_e=T`$. This is not quite surprising, since the temperature required for electronic transitions is large, even if it decreases as the cluster grows. However, at least two reasons can be invoked to explain the differences between the present results and the previously published thermodynamical curves. Exchange Monte Carlo is known to be a convenient method in reducing quasi-ergodicity and accelerating convergence ptmc2 . Also, by adding the possibility of changing electronic surface, barrier crossing is further favored. Thus we expect that the MC methods developed here are much more efficient than the simple, single trajectories simulations. The new results are then consistent with a lower melting point calguet .

Extensions of the present algorithms to mean-field monoelectronic Hamiltonians other than tight-binding is possible. A physical limitation is the relevance of the calculated excited levels as single-particle states and the neglect of many-body electron interactions. Another more important limitation is the difficult combination with molecular dynamics, due to the non-explicit dependence of the average occupation numbers on the nuclear structure, at least for finite systems. The use of Car-Parrinello dynamics including entropic corrections wentzcovitch  requires some practical approximations, such as the Fermi-Dirac form for the occupation numbers, or the assumption that these numbers do not vary much during a short time scale.

The KMC scheme makes the present method naturally suited for use with Monte Carlo sampling of the nuclear degrees of freedom. Accelerating procedures calspi2 ; mousseau  for the total energycalculations in tight-binding models can also be a valuable improvement for large systems. The methods introduced in this Letter have a wide range of applications, for both finite and infinite systems, metals, insulators or semiconductor materials. They provide a numerically accurate way of calculating the fractional occupation numbers, and we gave evidences that the Fermi-Dirac statistical distribution is not appropriate for small clusters. Combined with advanced Monte Carlo techniques such as parallel tempering or the multicanonical ensemble sampling, these methods enable one to investigate the equilibrium thermodynamics of large complex systems which exhibit various kinds of phase transitions due to structural isomerization or electronic excitations. In this respect, clusters close to the insulator-metal crossover or having magnetic properties offer good candidates for further investigations.


text:
# New measurements of the 𝐷⁰ and 𝐷⁺ lifetimes

## Abstract

A high statistics sample of photoproduced charm particles from the FOCUS (E831) experiment at Fermilab has been used to measure the $`D^0`$ and $`D^+`$ lifetimes. Using about $`\mathrm{210\hspace{0.17em}000}D^0`$ and $`\mathrm{110\hspace{0.17em}000}D^+`$ events we obtained the following values: $`409.6\pm 1.1\text{(statistical)}\pm 1.5`$ (systematic) fs for $`D^0`$ and $`1039.4\pm 4.3\text{(statistical)}\pm 7.0`$ (systematic) fs for $`D^+`$.

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

The study of the charm hadron lifetimes has been fundamental for our understanding of the heavy quark decays. The most important contribution is the spectator diagram which contributes equally to the widths of all hadrons of a given flavour . In the early days of charm physics it was quite a surprise when the experiments measured a large value for $`\tau _{D^+}/\tau _{D^0}`$. It is generally believed that this large ratio ($`2.5`$) is mainly due to the destructive interferences between different quark diagrams that contribute only to $`D^+`$ decays. The increasingly precise measurements of the heavy quark lifetimes have stimulated the development of theoretical models, like the Heavy Quark Theory , which are able to predict successfully the rich pattern of charm hadron lifetimes, that span one order of magnitude from the longest lived ($`D^+`$) to the shortest lived ($`\mathrm{\Omega }_c^0`$).

In this letter we present the most accurate measurement to date of the lifetimes of the $`D^+`$ and $`D^0`$. Although the accuracy reached by the previous experiments is remarkable, for example the $`D^0`$ lifetime is known with an uncertainty of $`1\%`$ , we think that a more precise determination of the $`D^0`$ lifetime would be needed. For example it would allow a more accurate check for the determination of the lifetime difference in the neutral $`D`$-meson system (to evaluate the parameter $`y=\mathrm{\Delta }\mathrm{\Gamma }/2\mathrm{\Gamma }`$ of the $`D^0\overline{D}^0`$ mixing ).

Charmed particles were produced by the interaction of high energy photons, obtained by means of bremsstrahlung of electron and positron beams (with typically $`300`$ GeV endpoint energy), with a beryllium oxide target. The mean energy of the photon beam was approximately $`180`$ GeV. The data were collected at Fermilab during the 1996–97 fixed-target run. More than $`6.3\times 10^9`$ triggers were collected from which more than 1 million charmed particles have been reconstructed.

The particles from the interaction are detected in a large-aperture magnetic spectrometer with excellent vertex measurement, particle identification and calorimetric capabilities. The vertex detector consists of two systems of silicon microvertex detectors. The upstream system consists of 4 planes interleaved with the experimental target  (2 target slab upstream then 2 silicon planes and the replica of this setting), while the downstream system consists of 12 planes of microstrips arranged in three views. These detectors provide high resolution separation of primary (production) and secondary (decay) vertices with an average proper time resolution of $`35\mathrm{fs}`$. The momentum of the charged particles is determined by measuring their deflections in two analysis magnets of opposite polarity with five stations of multiwire proportional chambers. Kaons and pions in the $`D`$-meson final states are well separated up to $`60\mathrm{GeV}/c`$ of momentum using three multicell threshold Čerenkov counters.

The final states are selected using a candidate driven vertex algorithm . A secondary vertex is formed from the reconstructed tracks and the momentum vector of the $`D`$ candidate is used as a seed to intersect the other tracks in the event to find the primary vertex. Once the production and decay vertices are determined, the distance $`\mathrm{}`$ between them and the relative error $`\sigma _{\mathrm{}}`$ are computed. Cuts on the $`\mathrm{}/\sigma _{\mathrm{}}`$ ratio are applied to extract the $`D`$ signals from the prompt background. The primary and secondary vertex are required to have a confidence level greater than $`1\%`$.

The vertices (primary and secondary) have to lie inside a fiducial volume111The reason for this cut is the presence of a trigger counter just upstream of the second microstrip device, therefore we define the fiducial volume as the region between the first slab of the experimental target and this trigger counter. and the primary vertex must be formed with at least two reconstructed tracks in addition to the seed track. The Čerenkov particle identification cuts used in FOCUS are based on likelihood ratios between the various stable particle identification hypotheses. These likelihoods are computed for a given track from the observed firing response (on or off) of all cells within the track’s ($`\beta =1`$) Čerenkov cone for each of our three Čerenkov counters. The product of all firing probabilities for all cells within the three Čerenkov cones produces a $`\chi ^2`$-like variable $`W_i=2\mathrm{ln}(\mathrm{Likelihood})`$ where $`i`$ ranges over the electron, pion, kaon and proton hypotheses (see Ref.  for more details). We require $`\mathrm{\Delta }_KW_\pi W_K>1`$, called kaonicity, for the tracks reconstructed as a kaon. Analogously the tracks reconstructed as pions have a pionicity, $`\mathrm{\Delta }_\pi W_KW_\pi `$, exceeding 1.

In Fig. 1a and 1b we show the invariant mass plots obtained with this set of cuts and with $`\mathrm{}/\sigma _{\mathrm{}}>9`$ for the decay modes $`D^0K^{}\pi ^+`$ and $`D^0K^{}\pi ^+\pi ^+\pi ^{}`$ respectively (throughout this paper the charge conjugate state is implied). Fig. 1c is the invariant mass plot for the decay mode $`D^+K^{}\pi ^+\pi ^+`$ with $`\mathrm{}/\sigma _{\mathrm{}}>14`$. This set of cuts is chosen to optimize the yield and the background underneath the signal ($`S/N`$ ratio).

From a binned maximum likelihood fit we find $`\mathrm{139\hspace{0.17em}433}\pm 520D^0K^{}\pi ^+`$, $`\mathrm{68\hspace{0.17em}274}\pm 360D^0K^{}\pi ^+\pi ^+\pi ^{}`$ and $`\mathrm{109\hspace{0.17em}877}\pm 385D^+K^{}\pi ^+\pi ^+`$ candidates. The plots are fit with two Gaussians with the same mean but different widths to take into account the different resolution in momentum of the tracks passing through one or two magnets (see Ref.  for more details) of our spectrometer plus a $`2^{\text{nd}}`$ order polynomial. The low mass region is excluded in the fit to avoid possible contamination due to other hadronic charm decays involving an additional $`\pi ^0`$.

The lifetime is measured using a binned maximum likelihood fitting technique . A fit is made to the reduced proper time distribution in the signal region. The reduced proper time is defined by $`t^{}=(\mathrm{}N\sigma _{\mathrm{}})/(\beta \gamma c)`$ where $`\mathrm{}`$ is the distance between the primary and the secondary vertex, $`\sigma _{\mathrm{}}`$ is the resolution on $`\mathrm{}`$ and $`N`$ is the minimum “detachment” cut required to extract the signal. Our vertexing algorithm provides very uniform reduced proper time acceptance even at very low reduced proper times. If absorption and acceptance corrections are small enough that they can be neglected, and if $`\sigma _{\mathrm{}}`$ is independent of $`\mathrm{}`$, one can show that the $`t^{}`$ distribution for decaying charmed particles will follow an exponential distribution. These assumptions are very nearly true in FOCUS .

The signal region reduced proper time distributions (indicated with dashed lines in Fig. 1) are formed from events with invariant mass within $`\pm 2\overline{\sigma }`$ of the mean $`D`$ mass; $`\overline{\sigma }`$ stands for weighted sigma because the invariant mass plots were fitted with two Gaussians. The dependence of the lifetime measurement on the choice of the signal and background region is discussed later in the text. The binned maximum likelihood method allows direct use of the proper time distribution of the data above and below the $`D`$ mass peak to represent the background underneath the signal instead of using a background parametrization. We have chosen two sidebands starting $`4\overline{\sigma }`$ above and below the mean $`D`$ mass, each half as wide as the signal region (indicated with dotted lines in Fig. 1). The signal and background reduced proper time distributions are binned in proper time wide bins ($`200`$ fs) spanning about $`10`$ nominal lifetimes.

The observed numbers of events in a reduced proper time bin $`i`$ (centered at $`t_i^{}`$) in the signal and sideband histograms are labeled $`s_i`$ and $`b_i`$ respectively. The predicted number of events $`n_i`$ in a reduced proper time bin is given by:

$$n_i=(N_sB)\frac{f(t_i^{})\mathrm{exp}(t_i^{}/\tau )}{_if(t_i^{})\mathrm{exp}(t_i^{}/\tau )}+B\frac{b_i}{b_i}$$
(1)
where $`N_s`$ is the total number of events in the signal region, $`B`$ is the total number of background events in the signal region and $`f(t_i^{})`$ is a correction function. The fit parameters are $`B`$ and $`\tau `$. The $`f(t_i^{})`$ correction function, derived from a Monte Carlo simulation, corrects the reduced proper time evolution of the signal for the effects of geometric acceptance, reconstruction efficiency, analysis cuts, hadronic absorption and decay of charm secondaries. The use of a multiplicative $`f(t_i^{})`$ correction, rather than an integral over a resolution factor, is justified since our reduced proper time resolution ($`35\mathrm{fs}`$) is much less than the $`D^0`$ or $`D^+`$ lifetime.

A separate $`f(t^{})`$ correction function is used for each of the three decay modes. Our Monte Carlo simulation includes the Pythia  model for photon-gluon fusion and incorporates a complete simulation at the digitization level of all detector and trigger systems and includes all known multiple scattering and particle absorption effects. The Monte Carlo was run with $`15`$ times the statistics of the experiment.

The plots $`a)`$ of Fig. 2, 3 and 4 show the correction function $`f(t^{})`$ for the three decay modes in bins of reduced proper time. The $`f(t^{})`$ function is obtained by dividing the simulated reconstructed charm yield in each bin by the input decay exponential. The fall off in $`f(t^{})`$ for the $`D^+`$ case is due to the exclusion of long lived events with vertices downstream of the fiducial volume1.

A factor $`_{\text{bg}}`$ is included in the likelihood function in order to relate $`B`$ to the number of background events expected from the side band population. The background level is thereby jointly determined from the invariant mass distribution and from the reduced proper time evolution in the side bands. The likelihood function is then given by:

$$=_{\text{signal}}\times _{\text{bg}}$$
(2)
where

$$_{\text{signal}}=\underset{i=1}{\overset{\text{bins}}{}}\frac{n_i^{s_i}\mathrm{exp}(n_i)}{s_i!}$$
(3)
and

$$_{\text{bg}}=\frac{B^{N_{\text{bg}}}}{N_{\text{bg}}!}\mathrm{exp}(B)$$
(4)
with $`N_{\text{bg}}=_ib_i`$ (we assume a linear background because the $`2^{\text{nd}}`$ order term of the polynomial is negligible).

The plots $`b)`$ of Fig. 2, 3 and 4 show the predicted events (histogram) superimposed on the observed events, the background events $`b_i`$ are also superimposed. In plots $`c)`$ of Fig. 2, 3 and 4 a pure exponential function with the fitted lifetime is superimposed on the background subtracted and $`f(t^{})`$ corrected $`t^{}`$ distribution.

The measured lifetimes are: $`408.75\pm 1.42`$ fs for $`D^0K^{}\pi ^+`$, $`411.25\pm 1.95`$ fs for $`D^0K^{}\pi ^+\pi ^+\pi ^{}`$ and $`1039.42\pm 4.28`$ fs for $`D^+K^{}\pi ^+\pi ^+`$.

Our lifetime measurements have been tested by modifying each of the vertex and Čerenkov cuts individually. For example in Fig. 5 one can see the measured lifetimes versus the $`\mathrm{}/\sigma _{\mathrm{}}`$ detachment cut. The measured lifetimes of the three decay modes are stable with respect to $`\mathrm{}/\sigma _{\mathrm{}}`$. Our measured lifetimes show no significant variation with the cuts employed to extract the signal.

We check our $`D^0`$ lifetime evaluation partitioning the total sample into $`D^{}`$tag and *no-tag* according to their origin. The obtained lifetimes are in very good agreement with the reported values.

To further check our lifetime measurements we have used tight cuts in order to extract a signal with virtually no background. Fig. 6 shows the invariant mass plots of the three decay modes for this set of cuts. The lifetime measurements from these samples, $`411.26\pm 3.11`$ fs for $`D^0K^{}\pi ^+`$ $`413.10\pm 4.80`$ fs for $`D^0K^{}\pi ^+\pi ^+\pi ^{}`$ and $`1036.66\pm 8.00`$ fs for $`D^+K^{}\pi ^+\pi ^+`$, are in very good agreement with our previous determination.

Systematic uncertainties for these lifetime measurements can arise from several sources. We performed a detailed study to analyze these sources.

There is an uncertainty due to the absolute time scale which was determined by studying the absolute length and momentum scale in the experiment (see Ref.  for more details). We estimate an uncertainty of $`\pm 0.11\%`$ for this source.

Another source of systematic uncertainty is linked to the detector and reconstruction efficiency. The $`f(t^{})`$ corrects the reduced proper time distribution for these effects, but an uncertainty could originate if there is a mismatch between the Monte Carlo simulation and the data. We have verified that our Monte Carlo accurately reproduces the distributions of several relevant variables, such as the longitudinal and transverse momenta, the multiplicity of the production vertex, the measured decay length and the estimated error on the reconstructed proper time. In order to estimate this uncertainty we split our total sample into independent subsamples depending on $`D`$ momentum, particle versus anti-particle and the different periods in which the data were collected. The splits into $`D`$ momentum and charge conjugation are the natural tests to reveal a possible mismatch between data and Monte Carlo because they probe the response of the detector. The main reason for the period dependence is the insertion of the upstream silicon system (which improved the resolution) in the target region during the 1997 fixed-target run period. A technique, employed in FOCUS and in the predecessor experiment E687, modeled after the S-factor method from the Particle Data Group , was used to try to separate true systematic variations from statistical fluctuations. The lifetime is evaluated for each of the $`8(=2^3)`$ statistically independent subsamples and a scaled variance is calculated; the split sample variance is defined as the difference between the reported statistical variance and the scaled variance if the scaled variance exceeds the statistical variance. This contribution to the systematic error is reported as split sample in Table 1.

The reported lifetimes are obtained with a particular set of fitting conditions. For example the width of the bins or the range of the $`t^{}`$ distribution. This is a particular choice and the lifetime should be independent of it. We investigated if this could be a possible source of uncertainty by varying the width of the bins, the upper limit of the $`t^{}`$ distribution, the location and the width of the sidebands, and the width of the signal region. In addition we studied the effect of using only the low or only the high mass sideband as well as the effect of eliminating the background term in the likelihood (the second term in equation $`2`$). For all these fit variants the sample variance is used as an estimate of this uncertainty because the various measurements are all taken as a priori likely.

A further source of systematic error can be due to uncertainties in the target absorption corrections. Two effects are present: hadronic absorption of decay daughters which would increase the fitted lifetime if neglected and absorption of the $`D`$ in the target which would tend to decrease the fitted lifetime if not taken into account. In our Monte Carlo all known particle absorption effects (values from Particle Data Group ) have been simulated for the decay daughters. For the charm hadron absorption our simulation assumes $`1/2`$ of the cross section for neutrons. We estimate the uncertainty of this contribution by varying the charm cross-section by $`50\%`$ and the daughter particle interaction cross-sections by $`25\%`$ in the Monte Carlo simulation. We verify that these estimates, reported in Table 1, are consistent with a determination of this contribution obtained comparing the lifetimes of decays with the $`D`$ produced in the upstream half of each target with those produced in the downstream half of the same target (see Ref.  for more details on the target setup). Each partition represents a different mixture of hadronic absorptions of decay daughters and $`D`$ mesons.

The acceptance could be another source of uncertainty. We analyzed this effect determining the $`D^0`$ lifetime without the correction function $`f(t^{})`$ and removing the fiducial volume cut from the set of analysis cuts (this makes the correction function almost flat). We obtained results in good agreement with the reported values. This check is not possible for the $`D^+`$ because of the longer lifetime, a geometric acceptance correction is always needed. A study was performed in FOCUS (see Ref.  for more details) comparing the acceptance part of the Monte Carlo correction with the high statistics $`K_S^0\pi ^+\pi ^{}`$ decays. The result of this study showed an excellent agreement between the acceptance observed in the data and the acceptance simulated by the Monte Carlo; however we assess a $`2\%`$ uncertainty due to the finite statistics. This $`2\%`$ uncertainty in the $`f(t^{})`$ correction function gives a $`0.21\%`$ and $`0.52\%`$ uncertainty in the lifetime of $`D^0`$ and $`D^+`$ respectively.

The finite Monte Carlo statistics give a negligible contribution to the systematic uncertainty.

Table 1 shows the contributions of each of these sources to the total systematic uncertainty. For the combined $`D^0`$ lifetime the systematic error (also shown in Table 1) is obtained combining the individual sources of systematic uncertainty from $`D^0K^{}\pi ^+`$ and $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$. We assume the absolute time scale and the absorption correlated. To obtain the final systematic error the uncertainties from the different sources are then added in quadrature.

The final lifetime values are $`409.62\pm 1.15\text{(statistical)}\pm 1.55`$ (systematic) fs for $`D^0`$ (weighted average) and $`1039.42\pm 4.28\text{(statistical)}\pm 6.97`$ (systematic) fs for $`D^+`$. A final check was performed for the $`D^0`$ lifetime. We compute the lifetime using a combined likelihood, that is forming a global likelihood for $`K^{}\pi ^+`$ and $`K^{}\pi ^+\pi ^{}\pi ^+`$. The fit parameters are the two distinct backgrounds and one lifetime. The lifetime from the combined likelihood, $`409.62\pm 1.15`$ (statistical), is identical to the reported value.

This measurement of the $`D^0`$ lifetime value is in very good agreement with the result we obtained in our lifetime difference paper .

In conclusion we have measured the lifetimes of the $`D^0`$ and $`D^+`$ mesons. Our results are reported in Table 2 along with a comparison with the most recent published measurements.

Our results will significantly decrease the errors on the current world average values for the $`D^0`$ and $`D^+`$ lifetimes.

From our measurements of the $`D^0`$ and $`D^+`$ lifetimes we can update the determination of the ratio $`\tau (D^+)/\tau (D^0)`$: $`2.538\pm 0.023`$. This result and the inclusive semileptonic branching ratios , $`D^+eX=(17.2\pm 1.9)\%`$ and $`D^0eX=(6.75\pm 0.29)\%`$, show that the $`D^0`$ and $`D^+`$ semileptonic decay widths are nearly equal:

$$\frac{\mathrm{\Gamma }(D^0eX)}{\mathrm{\Gamma }(D^+eX)}=\frac{B(D^0eX)}{B(D^+eX)}\times \frac{\tau (D^+)}{\tau (D^0)}=1.00\pm 0.12$$
(5)
This implies that differences in the total decay widths between $`D^0`$ and $`D^+`$ must be due to differences in the hadronic decay sector.

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


text:
# ENTANGLED STATES AND LOCAL MEASUREMENTS

## Abstract

We show that the maximum entanglement in a composite system corresponds to the maximum uncertainty and maximum correlation of local measurements.

It has been recognized that entanglement is one of the most symbolic, suggestive, and, at the same time, enigmatic concepts in quantum mechanics. On the one hand this concept touches on the fundamental notions of locality and reality in quantum physics  and on the other hand plays a key role in quantum information theory .

The physical essence of entanglement consists in the existence of quantum correlations between the individual parts of a composite system that have interacted once in the past but are no longer interacting. Formally, these correlations are caused by the combination of the superposition principle in quantum mechanics with the tensor product structure of the space of states .

In the usual treatment, the entanglement is associated with the nonseparability of states of a composite quantum system. At the same time, the nonseparability of states by itself is not a sufficient condition of either maximum or distilled entanglement . We also note that the violation of Bell’s inequalities and other conditions of classical realism, which has attracted a great deal of interest in the context of entanglement (e.g., see  and references therein), demonstrate nothing but nonexistence of hidden variables in quantum mechanics and therefore cannot be considered as a criterion of entanglement. Moreover, there are a certain examples of mixed states , which do not violate Bell’s inequalities but manifest entanglement and therefore can be used for quantum teleportation .

Definitely the nonseparability of entangled states has a certain informational aspect . In particular, there is an assumption that the essence of entanglement consists in the fact that the information resides mostly in the correlations then in properties of individual parts of a composite system .

In fact, the correlations peculiar to entanglement manifest themselves via local measurements. Since Shannon’s definition of an information channel implies a correlation between input and output, it is possible to say that an entangled state represents a virtual information channel between the individual parts of the system. The main objective of this note is to show that the maximum entanglement, corresponding to the maximum channel capacity, is associated with the maximum correlation between the local measurements, which is equivalent to the maximum uncertainty of these measurements.

As an example of some considerable interest, we investigate the case of a composite system, consisting of the two parts with two quantum degrees of freedom each. Such a system is defined in the Hilbert space

$`=𝐂^2𝐂^2,`$ (1)
where each subspace is spanned by the basis

$`|e_{1,2}=|\pm .`$ (2)
A general state defined in (1) has the form

$`|\psi ={\displaystyle \underset{i,j=1}{\overset{2}{}}}a_{ij}|e_i|e_j`$ (3)
$`=a_{11}|+++a_{12}|++a_{21}|++a_{22}|,`$ (4)
where

$`\psi |\psi ={\displaystyle \underset{i,j}{}}|a_{ij}|^2=1.`$ (5)
This state (2) describes the information channel between the two parts of the system under consideration. It is seen that (2) is more general than the conventional Schmidt form

$`|\psi _S={\displaystyle \underset{i}{}}b_i|e_i|e_i`$ (6)
usually discussed in the context of entanglement .

In either individual part of the composite system, the local measurement is represented by the set of infinitesimal generators of the $`SL(2)`$ algebra represented by the Pauli matrices

$`\sigma _1`$ $`=`$ $`|+|+|+|,`$ (7)
$`\sigma _2`$ $`=`$ $`i(|+||+|),`$ (8)
$`\sigma _3`$ $`=`$ $`|++|||.`$ (9)
The physical realization of these operators depends on the concrete definition of the system under consideration.

It is then clear that the correlation between the local measurements at different individuals, described by the correlation function

$`\sigma _{\mathrm{}}^{(1)};\sigma _{\mathrm{}^{}}^{(2)}\sigma _{\mathrm{}}^{(1)}\sigma _{\mathrm{}^{}}^{(2)}\sigma _{\mathrm{}}^{(1)}\sigma _{\mathrm{}^{}}^{(2)},`$ (10)
achieves the maximum at the following condition

$`i,\mathrm{}\sigma _{\mathrm{}}^{(i)}=0.`$ (11)
The averaging in above formulas is performed over the state (2).

In turn, the uncertainty of local measurement described by the variance

$`\left(\delta \sigma _{\mathrm{}}^{(i)}\right)^2=\left(\sigma _{\mathrm{}}^{(i)}\right)^2\sigma _{\mathrm{}}^{(i)}^2`$ (12)
also achieves the maximum under the same condition (7).

We now note that the opposite condition of minimum uncertainty of the $`SL(2)`$ generators corresponds to the so-called atomic coherent states have been discussed in . It is known that the coherent states describe almost classical behavior of the corresponding system. Therefore, all one can expect is that the states with maximum uncertainty should describe the pure quantum behavior of the system.

In fact, the condition (7) together with the normalization condition (3) impose on the matrix

$`A=a_{ij}`$ (13)
of the coefficients in (2) the following constraints

$`|a_{11}|^2+|a_{12}|^2={\displaystyle \frac{1}{2}},|a_{22}|=|a_{11}|,|a_{21}|=|a_{12}|,`$ (14)
$`\mathrm{arg}a_{11}+\mathrm{arg}a_{22}\mathrm{arg}a_{12}\mathrm{arg}a_{21}=\pm \pi +2k\pi ,`$ (15)
where $`k=0,1,\mathrm{}`$. These conditions specify an infinite subset of states (2). In particular, the choice of either $`a_{11}=0`$ or $`a_{12}=0`$ in (10) lead to the two sets of mutually orthogonal EPR states

$`|\psi _{EPR}={\displaystyle \frac{1}{\sqrt{2}}}(|++e^{i\varphi }|+)`$ (16)
and

$`|\phi _{EPR}={\displaystyle \frac{1}{\sqrt{2}}}(|+++e^{i\varphi ^{}}|),`$ (17)
respectively. Here

$`\varphi =\mathrm{arg}a_{21}\mathrm{arg}a_{1,2},`$ (18)
$`\varphi ^{}=\mathrm{arg}a_{22}\mathrm{arg}a_{11}.`$ (19)
Another example is provided by the orthogonal states

$`|\psi `$ $`=`$ $`{\displaystyle \frac{1}{2}}[i(|+++|)+|++|+],`$ (20)
$`|\psi ^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[|+++|+i(|++|+)].`$ (21)
By definition, the states (2) are nonseparable and hence entangled. Let us now show that the condition (7) implies the maximum entanglement of the states (2). For this, we consider a standard way based on the use of reduced density matrices

$`\rho ^{(i)}=Tr_i\rho ,\rho =|\psi \psi |.`$ (22)
It is easily seen that

$`\rho ^{(1)}=AA^+,\rho ^{(2)}=A^+A.`$ (23)
It is a straightforward matter to arrive at the conclusion that the conventional conditions of internal consistency of quantum mechanics  together with the condition of maximum of reduced entropy

$`iS_i=Tr(\rho ^{(i)}\mathrm{ln}\rho ^{(i)})=\mathrm{ln}2`$ (24)
yield, on allowing of (3), the constraints (10) and (11). In other words, they are equivalent to the condition (7).

It is no wonder that the requirements of maximum uncertainty of local measurements and maximum of reduced entropy lead to the same condition (7). In fact, the group of symmetry of the system under consideration is

$`G=SU(2)SU(2)`$ (25)
and $`SU(2)SL(2)`$. Since $`SL(2)`$ represents the complexification of $`SU(2)`$, the Pauli operators (5) are responsible for the dynamical symmetry in the system under consideration. By means of the transformation

$`g=g_1g_2,gG`$ (26)
applied to the matrix (9) we get

$`Ag_1Ag_2^+,A^+g_2Ag_1^+.`$ (27)
The only quantity, which has the physical meaning for the virtual information channel described by the state (2), is represented by the invariant with respect to (13), that is, by the trace of matrices $`AA^+`$ and $`A^+A`$:

$`I=Tr(AA^+)Tr\rho ^{(1)}=Tr(A^+A)Tr\rho ^{(2)}.`$ (28)
Use of (14) allows this to be cast into the form

$`I=Tr(g_1AA^+g_1^+)=Tr(g_2A^+Ag_2^+),`$ (29)
which completes the prove.

In particular, the equality (15) can be used to show that the Pauli matrices (5) commute with the reduced density matrices

$`i,\mathrm{}[\sigma _{\mathrm{}},\rho ^{(i)}]=0,`$ (30)
if the elements of (9) obey the condition (7).

Let us stress that the obtained results for the channel defined by the relations (2) and (7) are compatible with the quantum counterpart of Shannon’s coding theorem developed by Holevo .

Similar results can be established for the general GHZ states as well. In this case, the system is defined in the space

$`_3=𝐂^2𝐂^2𝐂^2`$ (31)
and the group of symmetry has the form

$`G^{}=SU(2)SU(2)SU(2).`$ (32)
A general state defined in (16) has the following form

$`|\psi ={\displaystyle \underset{i,j,k=1}{\overset{2}{}}}x_{ijk}|e_i|e_j|e_k,|e_{1,2}=|\pm ,`$ (33)
and obeys the normalization condition. It is straightforward to show that condition (7) with averaging over the state (18) provides the maximum of reduced entropy as above. In particular, condition (7) defines the conventional GHZ states

$`|\psi _{GHZ}={\displaystyle \frac{1}{\sqrt{2}}}(|+++\pm |).`$ (34)
Besides that, there are infinitely many other maximum entangled states defined by the relations (18) and (7). A nontrivial example is provided by

$`|\psi ={\displaystyle \frac{1}{\sqrt{8}}}[|++++|i(|+++|+)`$ (35)
$`+|+++|++i(|+++|+)].`$ (36)
In general, the maximum entangled states of an arbitrary multi-component system can be determined by the condition (7) in the similar fashion.

To conclude, we have shown that the maximum entanglement corresponds to the maximum correlation between the local measurements and to the maximum uncertainty of these measurements. The quantitative criterion is represented by the condition (7) which enables us to fairly simplify the analysis of entangled states. The novelty is that the condition of maximum entanglement is defined directly in terms of what can be measured for the individual parts of the system.

Our consideration so far have applied to the composite systems with individuals defined in the two-dimensional space. As a matter of fact, the increase of the number of quantum degrees of freedom per individual part leads to a certain change of local measurement. For example, the case of free degrees of freedom per component assumes that the measurement is represented by the eight Hermitian operators associated with the generators of the $`SU(3)`$ algebra. This operators are similar, in a certain sense, to the Stokes operators, describing the polarization of multipole photons . The maximum entangled states, in this case, can be defined in direct analogy to the above result by the requirement of maximum uncertainty of local measurements. The detailed investigation of this case deserves special consideration.


text:
# Effect of a wide binary companion to the lens on the astrometric behavior of gravitational microlensing events

## 1 Introduction

Light curves of microlensing events caused by a single point-mass lens are characterized by non-repeating, single peak, and symmetric curves (Paczyński, 1986), known as Paczyński curves. On the other hand, if events are caused by a binary lens, the resulting light curves deviate from those of single lens events. The deviation caused by the binary nature of the lens becomes prominent when the angular separation between the two lens components (binary separation) is comparable to the angular Einstein ring radius of the lens system, which is related to the physical parameters of the system by

$`\theta _\mathrm{E}`$ $`=`$ $`\left({\displaystyle \frac{2R_{\mathrm{Sch}}}{D}}\right)^{1/2}`$ (1)
$``$ $`0.72\text{mas}\left({\displaystyle \frac{D_\mathrm{S}}{D_\mathrm{L}}}1\right)^{1/2}\left({\displaystyle \frac{D_\mathrm{S}}{8\text{kpc}}}\right)^{1/2}\left({\displaystyle \frac{M}{0.5\text{M}\text{}}}\right)^{1/2}`$
where $`D^1=D_\mathrm{L}^1D_\mathrm{S}^1`$, $`D_\mathrm{L}`$ and $`D_\mathrm{S}`$ are the distances to the lens and the source from the observer respectively, and $`R_{\mathrm{Sch}}`$ is the Schwarzschild radius corresponding to the total mass $`M`$ of the lens system. As the binary separation becomes larger, however, the individual lens components begin to act as if they are two independent single lenses. For these wide binary lens systems, it is likely that the source trajectory will pass close only to one of the lens components and the resulting light curve will look very similar to a standard Paczyński curve. If the Galactic binaries have a similar distribution of binary separations to that of binaries in the solar neighborhood (Duquennoy & Mayor, 1991), the majority of the binaries located along the lines of sight towards the Galactic bulge will have separations substantially larger than their Einstein ring radii whose typical values are several hundred micro- to a few milli-arcseconds. Therefore, for most of Galactic binary lens events, it will be difficult to identify the existence of the companion to the lens from the analysis of the light curves alone except for the rare case of either caustic crossing or the source star’s passage near both lens components (a repeating event; Di Stefano & Scalzo, 1999).

As an alternative method of microlensing observations, it was proposed to measure the lensing-induced displacement of the (brightness-weighted) image centroid position (usually referred to centroid shifts) by high precision astrometry (Miyamoto & Yoshii, 1995; Høg, Novikov, & Polanarev, 1995; Walker, 1995; Miralda-Escudé, 1996; Paczyński, 1998; Boden, Shao, & Van Buren, 1998; Han & Chang, 1999; Dominik & Sahu, 2000). When a source star is microlensed, it becomes split into several images whose number depends on the structure of the lens system, e.g., two images for a single lens and three or five images for a binary lens system. The typical separation between the images is an order of $`\theta _\mathrm{E}`$, and thus it is impossible to resolve the individual images even with the highest angular resolution achieved so far. However, several planned high-precision astrometric satellite experiments, e.g., Space Interferometry Mission (SIM), Double Interferometer for Visual Astrometry (DIVA), and Global Astrometric Interferometer for Astrophysics (GAIA), and optical/NIR interferometers on 10-m class ground-based telescopes, e.g., the Keck interferometer, the Very Large Telescope Interferometer (VLTI), and the Large Binocular Telescope (LBT), will have precisions in position measurements of the order of tens of $`\mu `$-arcsecs or better, and thus will enable one to measure the lensing-induced centroid shifts.

If an event is caused by a point-mass lens, the trajectory of the centroid shifts traces an ellipse (see § 2), known as an astrometric ellipse (Walker, 1995; Jeong, Han, & Park, 1999; Dominik & Sahu, 2000). If an event is caused by a binary lens, on the other hand, the centroid shift trajectory deviates from the elliptical one of a single lens event, which is characterized by distortions, twistings, and jumps (Han, Chun, & Chang, 1999).

In this paper, we investigate how a wide binary companion of the lens affects the observed photometric and astrometric microlensing behaviors. From this investigation, we find that the effect on the trajectory of the image centroid shifts is significantly more important than on the light curve. Hence, we argue that a careful consideration of the effect of a wide binary lens will be required in analyzing the measured centroid shifts even if the corresponding light curve appears to be a standard Paczyński curve.

## 2 Basics of Microlensing

If a source located at $`𝒓_\mathrm{S}`$ (on the projected plane of the sky) is lensed by a $`N`$-point-mass lens system, where the individual components’ masses and locations are $`m_i`$ and $`𝒓_{\mathrm{L},i}`$, the positions of the images $`𝒓`$ are obtained by solving the equation of lens mapping (the lens equation);

$$𝒓_\mathrm{S}=𝒓\theta _\mathrm{E}^2\underset{i=1}{\overset{N}{}}\frac{m_i}{M}\frac{𝒓𝒓_{\mathrm{L},i}}{|𝒓𝒓_{\mathrm{L},i}|^2},$$
(2)
where $`M=_i^Nm_i`$. Since the lensing process conserves the surface brightness (i.e. does not create nor destroy photons), the flux ratio between the lensed image and the unlensed source is simply given by the surface area ratio between the image and the unmagnified source, i.e., the magnification. For a point source limit, the magnifications of the individual images are given by the inverse of the Jacobian determinant of the lens mapping evaluated at each image position $`𝒓_j`$;

$$A_j=\left(\frac{1}{|detJ|}\right)_{𝒓=𝒓_j};detJ=\left|\frac{𝒓_\mathrm{S}}{𝒓}\right|$$
(3)
and the total magnification is given by the sum of the magnifications of the individual images, i.e., $`A_{\mathrm{tot}}=_j^{N_\mathrm{I}}A_j`$, where $`N_\mathrm{I}`$ is the total number of images. The location of the apparent image centroid corresponds to the magnification weighted mean of the individual image positions, and therefore, the centroid shift is given by

$$𝜹=\underset{j}{\overset{N_\mathrm{I}}{}}\frac{A_j}{A_{\mathrm{tot}}}𝒓_j𝒓_\mathrm{S}.$$
(4)
For a single point-mass lens ($`N=1`$), the lens equation (eq. ) is easily solvable and solving the equation yields two image positions ($`N_\mathrm{I}=2`$). The total magnification and the centroid shift of the single lens event are expressed in analytical forms of

$$A=\frac{u^2+2}{u\sqrt{u^2+4}};$$
(5)
$$𝜹=\frac{\theta _\mathrm{E}𝒖}{u^2+2},$$
(6)
where $`𝒖(𝒓_\mathrm{S}𝒓_\mathrm{L})/\theta _\mathrm{E}`$ is the dimensionless lens-source separation vector normalized by $`\theta _\mathrm{E}`$. Since the source and the lens move relative to each other, $`𝒖`$ changes with time. Under the approximation that this relative motion is rectilinear, it is represented by

$$𝒖=\left(\frac{tt_0}{t_\mathrm{E}}\right)\widehat{𝒙}+\beta \widehat{𝒚},$$
(7)
where $`t_\mathrm{E}`$ represents the time required for the source to transit $`\theta _\mathrm{E}`$ (Einstein time scale), $`\beta `$ is the closest lens-source separation in units of $`\theta _\mathrm{E}`$ (impact parameter), $`t_0`$ is the time at that moment, and the unit vectors $`\widehat{𝒙}`$ and $`\widehat{𝒚}`$ are parallel and perpendicular to the direction of the relative lens-source motion, respectively. The trajectory of the centroid shifts is an ellipse with a major axis and an eccentricity of $`a=(\theta _\mathrm{E}/2)\times (\beta ^2+2)^{1/2}`$ and $`e=[(\beta ^2/2)+1]^{1/2}`$ and its semi-major axis lies parallel to $`\widehat{𝒙}`$.

For a multi-lens system ($`N2`$), on the other hand, the lens equation cannot be solved algebraically. Hence, to obtain the image positions, one has to solve the equation numerically. In particular, for the case of a binary lens ($`N=2`$), the lens equation can be expressed as a fifth-order complex polynomial (Witt, 1990) and there exist three or five images depending on the source position with respect to the lens positions.

## 3 Effects of a Wide Binary Lens Companion

The existence of a binary companion to the lens affects both the photometric and astrometric behaviors of a microlensing event. The most important effect of a wide binary companion of the lens is that it makes the effective position of the primary111Here, the primary denotes the binary component which lies closer to the projected path of the source on the sky, and the use of ‘the companion’ is reserved for the other component of the binary. We note, however, that the basic results of the subsequent discussion do not change under the usual definition of the primary as the more massive binary component. move towards it or, equivalently, the effective source position move with respect to the primary (Di Stefano & Mao, 1996). From the lowest order approximation of equation (2), one finds that this effect is well characterized by the use of an effective lens-source separation vector, which is expressed by

$$\stackrel{~}{𝒖}𝒖_1\frac{q}{d}\frac{𝒓_{\mathrm{L},2}𝒓_{\mathrm{L},1}}{|𝒓_{\mathrm{L},2}𝒓_{\mathrm{L},1}|},$$
(8)
where the subscripts ‘1’ and ‘2’ represent the primary and the companion respectively, $`q=m_2/m_1`$ is the mass ratio, and $`𝒖_1=(𝒓_\mathrm{S}𝒓_{\mathrm{L},1})/\theta _{\mathrm{E},1}`$, and $`d=|𝒓_{\mathrm{L},2}𝒓_{\mathrm{L},1}|/\theta _{\mathrm{E},1}`$ are the distances of the source and the companion from the primary lens normalized by the Einstein ring radius of the primary $`\theta _{\mathrm{E},1}`$, respectively. We note that this effect is not generally measurable because it requires an a priori information on the actual positions of the lens components on the sky.

As an additional effect of the binary lens companion, one may expect that the observed light curve and the centroid shift trajectory would be the superposition of those resulting from the two events in which the component lenses behave as two independent lenses. However, this superposition effect of the lens companion is important only for the trajectory of the centroid shifts but not for the light curve of an event. This can be seen from the magnification and the centroid shift caused by the companion when the source is located close to the primary (i.e., $`𝒓_\mathrm{S}𝒓_{\mathrm{L},1}`$), which are, respectively, represented by

$$A_21+2\frac{q^2}{d^4}(\mathrm{for}d2),$$
(9)
$$\delta _2\theta _{\mathrm{E},1}\frac{q}{d}(\mathrm{for}d\sqrt{2}).$$
(10)
From equations (9), one finds that as the binary separation becomes larger, the photometric superposition effect falls off rapidly ($`d^4`$): faster than the decrease of $`|\stackrel{~}{𝒖}𝒖_1|`$ ($`d^1`$; see eq. ). On the other hand, the astrometric effect endures up to a very large separation at which the photometric effect is negligible (the same order effect as the positional shift; see eq. ).

To illustrate the larger astrometric superposition effect of the lens companion than its photometric effect, we present the light curves (Fig. 1) and the trajectories of the image centroid shifts (Fig. 2) of events caused by a wide binary lens with the mass ratio $`q=0.3`$ and the separation $`d=5.7`$. One finds that the light curves of these events are well approximated by that of a single lens event after accounting for the effective positional shift of the primary. By contrast, the corresponding centroid shift trajectories deviate significantly from an elliptical one expected for a single lens event.

## 4 Detection of Wide Binary Lens Companions

In order to assess how well one can distinguish light curves of wide binary lensing events from the ones of single point-mass events, we define an magnification excess by

$$ϵ=\frac{A\stackrel{~}{A}_0}{\stackrel{~}{A}_0},$$
(11)
where $`A`$ is the exact magnification of the wide binary lens event and $`\stackrel{~}{A}_0`$ is that of the single lens event with a mass equal to the primary and located at the effective position defined by equation (8),

$$\stackrel{~}{A}_0=\frac{\stackrel{~}{u}^2+2}{\stackrel{~}{u}\sqrt{\stackrel{~}{u}^2+4}}.$$
(12)
Figure 3 shows the magnification excess map of the region close to the primary lens for wide binary lens systems with various values of $`d`$ and $`q`$. From the figure, one finds that the region of significant deviations shrinks rapidly as $`d`$ increases. Noticeable deviations (e.g., $`ϵ2`$%) occur only if the source passes very close to the primary, implying that the presence of the companion to the lens can be detected effectively only if the source crosses the caustic associated with the primary or passes the region of the cusp-influence magnification. However, both the size of the caustic and the strength of the cusp converge to those of a point-mass lens much faster than the effect of the positional shift. As a result, even with the photometric precision as good as 0.01 mag and sufficiently dense sampling during the peak magnification, the binary detection efficiency is expected to be low for most wide binary microlensing events with $`d10`$.

In order to compare the photometric deviation to the one expected from astrometric observations, we also define the deviation of centroid shift by

$$\mathrm{\Delta }𝜹=𝜹\stackrel{~}{𝜹}_0,$$
(13)
where $`𝜹`$ is the centroid shift of the wide binary lensing event while $`\stackrel{~}{𝜹}_0`$ represents the centroid shift caused by the primary alone placed at its effective position,

$$\stackrel{~}{𝜹}_0=\frac{\theta _{\mathrm{E},1}\stackrel{~}{𝒖}}{\stackrel{~}{u}^2+2}.$$
(14)
Figure 4 shows the contour maps of $`\mathrm{\Delta }\delta `$ of the same region and for the same binaries as in Figure 3. When the source is located near the primary, the centroid shift deviation – mostly due to the superposition effect – is an order of $`(q/d)\theta _{\mathrm{E},1}`$ (see eq. ). Therefore, we draw contours at the levels which are scaled by $`(q/d)\theta _{\mathrm{E},1}`$. Figure 5 show the maps in a wider region where the contours are drawn at 6 different levels ranging $`1\%`$ to $`25\%`$ of $`\theta _{\mathrm{E},1}`$. Considering that the size of the angular Einstein ring radius is an order of a few milli- to several hundred micro-arcseconds for events caused by Galactic G, K, and M dwarfs with masses in the range of $`0.1\text{M}\text{}`$$`1\text{M}\text{}`$ (see eq. ), the deviations for the events caused by a binary lens even with $`q/d0.01`$ are expected to be still several micro-arcseconds, which is comparable to the expected astrometric precision of the new generation interferometers.222The astrometric precision of SIM, for instance, will be as good as $`5\mu \text{as}`$ for stars brighter than $`V=20`$ (Unwin, Boden, & Shao, 1997). Moreover, large deviations are found throughout the map, rather than concentrating near the primary as in the map of magnification excess. Notably, quite complex patterns of deviation are present over the region of the circle around the primary with the radius of about $`\theta _{\mathrm{E},1}`$ even for the case of $`d10`$. This implies that the astrometric deviation will typically last for $`t_\mathrm{E}`$. From astrometric microlensing observations of events, therefore, one can detect the presence of binary companions to the lens for a broader range of separations with an increased efficiency than from conventional photometric observations, possibly even with less dense sampling.

As a quantitative comparison between the detection efficiencies of wide binary companions to the lens expected from the photometric and the astrometric lensing observations, we perform a Monte-Carlo simulation. In the simulation, we produce a large number of wide binary lens events and fit their light curves and image centroid shift trajectories to those of single point-mass lens events (eqs.  and ). We assume the lens binarity is detected if the resulting $`\chi ^2`$ of the fit is larger than a statistically acceptable value (see below). Tables 1 and 2 list the photometric and astrometric detection rate determined by the simulation, each for the same wide binaries as Figures 3 and 4 under two different sampling strategies and three different measurement errors. The resulting detection rate in each entry is based on 5,000 trial simulations. We note that, since we are more interested in the possibility of the inference to the presence of the lens companion, rather than in the determination of the absolute detection efficiency, we exclude the repeating events and the caustic-crossing events for which the binarity of the lens can be immediately inferred. Trial events are generated with random trajectories for each fixed configuration of the binary (i.e, at fixed $`d`$ and $`q`$), but restricted only to high-magnification,333The maximum magnification due to the primary is greater than two, i.e. $`\mathrm{min}\{|𝒓_\mathrm{S}𝒓_{\mathrm{L},1}|^2)\}<[(4/\sqrt{3})2]\theta _{\mathrm{E},1}^2`$. non-repeating,444The magnification due to the companion is not bigger than two, i.e. $`\mathrm{min}\{|𝒓_\mathrm{S}𝒓_{\mathrm{L},2}(𝒓_{\mathrm{L},1}𝒓_{\mathrm{L},2})/d^2|^2\}[(4/\sqrt{3})2]\theta _{\mathrm{E},2}^2`$. and non-caustic-crossing events.555The number of data points when the source is inside the caustic is, at most, one. To simulate light curves, we fix the unlensed source flux at an arbitrary value and include also random amounts of blended light of less than 20% of the unlensed source flux. As for the observation, we examine two different types of sampling. In the first type, we assume 200 photometric and astrometric measurements during $`t_\mathrm{E}t_{\mathrm{obs}}t_03t_\mathrm{E}`$ (‘Broad’ sampling). In the second type, the same number of measurements are performed during a shorter period of time around the photometric peak of $`0.5t_\mathrm{E}t_{\mathrm{obs}}t_01.5t_\mathrm{E}`$ (‘Dense’ sampling). For simplicity, the sampling rate is taken to be uniform and the uncertainty associated with each measurement is the same throughout the measurements (i.e., the deviations distribute as Gaussian with a single variance). Finally, the best-fitting single lens model and $`\chi ^2`$ of a simulated event are determined through the $`\chi ^2`$ minimization by successive linearized fitting. For light curves, we fit five parameters: $`t_0`$, $`t_\mathrm{E}`$, $`\beta `$, and the unlensed source flux and the blended light while, for the image motion, we fit seven parameters: $`t_0`$, $`t_\mathrm{E}`$, $`\beta `$, $`\theta _\mathrm{E}`$, and the direction of the relative lens-source motion and the lens position on the sky (which is characterized by two parameters). The detection is defined by

$$\sqrt{2\chi ^2}\sqrt{(2\mathrm{dof}1)}3,$$
(15)
where ‘$`\mathrm{dof}`$’ denotes the degrees of freedom of the fitting. The criterion was chosen to be the 3$`\sigma `$ deviation of $`\chi ^2`$ from the expectation value if the underlying model (i.e., single point-mass lens event) were correct. However, due to the positive skewness of the $`\chi ^2`$-distribution, the criterion actually corresponds to the rejection of the single lens model with the confidence level of approximately 99.8%.

As expected, the result shows a higher binary detection rate for the 10-$`\mu `$as astrometry ($`\theta _\mathrm{E}0.5\text{mas}`$) than for the 1% photometry. For example, the astrometric detection efficiency of a equal-mass binary with a separation of $`d10`$ stays more than 50% even for events with $`\theta _\mathrm{E}0.5\text{mas}`$, while the photometric detection rate falls below 10% (2% photometry) to 30% (0.5% photometry). Furthermore, even in the regime where the photometric deviations are basically indistinguishable from random errors, the astrometric deviations can be unambiguously detected for a substantial fraction (5$``$10%) of events under a reasonable accuracy and a proper sampling strategy.

Another interesting result one finds from Tables 1 and 2 is that ‘Broad’ sampling strategy has an advantage in detecting wide binary lens companions over ‘Dense’ sampling strategy if the binary effect is strong (i.e., smaller $`d`$ and/or larger $`q`$) and/or the measurement precision is high. On the other hand, ‘Dense’ sampling strategy becomes advantageous over ‘Broad’ sampling strategy as the binary signal becomes weaker (i.e., larger $`d`$ and/or smaller $`q`$) and/or the measurement precision is poor. It is also interesting to see that the relative advantage of ‘Dense’ sampling strategy for weak binaries and/or poor measurements is more obvious for the photometry than the astrometry. This is more or less expected by the examination of the deviation maps in Figures 3 and 4, i.e., the region of substantial deviation is much more wide-spread in the astrometric deviation map than in the photometric one, and thus the photometric deviations will be likely to be shorter-lived than the astrometric deviations.

## 5 Discussion

The relatively significant effect of the wide binary lens companion on the astrometric lensing behavior compared to the effect on the photometric behavior, on one side, calls for careful considerations of the lens binarity when one attempts to analyze and interpret the future astrometric lensing data. In the current lensing experiments based on photometric observations, undetected wide binary lens companions do not pose as a serious problem because the resulting light curve is, in most cases, well approximated by that of a single lens event and more importantly the resulting time scale of the event (the only physically interesting lensing parameter for a single lens event) is hardly affected by the existence of the wide binary lens companion. However, as demonstrated in the previous sections, one generally cannot ignore the effect of the wide binary lens companion on the observed centroid shift trajectory even if the corresponding light curve has little signature of the lens binarity.

The large astrometric effect of the wide binary lens companion, on the other side, makes the future astrometric lensing experiment a useful method to detect binary lens companions over a broader range of separations. Especially, in searches for sub-stellar-mass and/or sub-luminous companions, astrometric lensing can provide a much more effective channel to detect them than the lensing light curve analysis. In addition, the study of the binary frequency can benefit from the astrometric lensing observations since it is sensitive to the existence of the companion to the lens over a wider range of separations and consequently suffers less severe bias towards resonant separation than the conventional lensing technique.

One possible complication involving the interpretation of binary-affected centroid shift trajectories is the parallax effect. The parallax effect becomes important due to the long-range duration of the centroid shifts up to a large lens-source separation. Therefore, proper interpretation of the observed centroid motion of a wide binary lens event will require careful consideration of the parallax effect. We defer systematic studies of the parallax effect for a future work.

## 6 Conclusion

By investigating the effect of wide binary lens companions on the photometric and astrometric behaviors of microlensing events, we find that the signature of the companion to the lens is significantly more important in the centroid shifts than in the light curve. On one side, this implies that, in analyzing the centroid shifts to be measured from the prospective astrometric lensing experiments, one should consider the effect of wide binary lens companions with more caution, which is generally not required in the analysis of light curves obtained from the current-type photometric lensing experiments. On the other side, this signifies the importance of the future astrometric lensing observations in providing a useful channel to search for binary companions to the lens over a broader range of separations.

This work was supported by a grant R01-1999-00023 from the Korea Science & Engineering Foundation (KOSEF). Work by J. An is supported by the Presidential Fellowship from the Graduate School of the Ohio State University.