text:
# Kondo effect in the presence of itinerant-electron ferromagnetism studied with the numerical renormalization group method

(Submitted: April 16, 2003)

## Abstract

The Kondo effect in quantum dots (QDs) - artificial magnetic impurities - attached to ferromagnetic leads is studied with the numerical renormalization group method. It is shown that the QD level is spin-split due to presence of ferromagnetic electrodes, leading to a suppression of the Kondo effect. We find that the Kondo effect can be restored by compensating this splitting with a magnetic field. Although the resulting Kondo resonance then has an unusual spin asymmetry with a reduced Kondo temperature, the ground state is still a locally-screened state, describable by Fermi liquid theory and a generalized Friedel sum rule, and transport at zero temperature is spin-independent.

Introduction. – The prediction glazman  and experimental observation of the Kondo effect in artificial magnetic impurities – semiconducting quantum dots (QDs) goldhaber ; kondo-odd  – renewed interest in the Kondo effect and opened new opportunities of research. The successful observation of the Kondo effect in molecular QDs like carbon nanotubes cobden ; bachtold  and single molecules molecular  attached to metallic electrodes opened the possibility to study the influence of many-body correlations in the leads (superconductivity superconducting  or ferromagnetism) on the Kondo effect. Recently the question arose whether the Kondo effect in a QD attached to ferromagnetic leads can occur or not. Several authors have predicted sergueev ; zhang ; bulka ; lopez  that the Kondo effect should occur. However, it was shown recently martinek1  that the QD level will be spin-split due to the presence of ferromagnetic electrodes leading to a suppression of the Kondo effect, and that the Kondo resonance can be restored only by applying an external magnetic field. The analyses of Refs. sergueev ; zhang ; bulka ; lopez ; martinek1  were all based on approximate methods.

In this Letter we resolve the controversy by adapting the numerical renormalization-group (NRG) technique costi1 ; hewson-book , one of the most accurate methods available to study strongly-correlated systems in the Kondo regime, to the case of a QD coupled to ferromagnetic leads with parallel magnetization directions. We find that in general the Kondo resonance is split, similar to the usual magnetic-field-induced splitting hofstetter2 ; costi . However, we find that by appropriately tuning an external magnetic field, this splitting can be fully compensatedand and the Kondo effect can be restored RKKY\_vs\_K  (confirming Ref. martinek1 ). We point out that precisely at this field the occupancy of the local level is the same for spin up and down, $`n_{}=n_{}`$, a fact that follows from the Friedel sum rule langreth . Moreover, we show that the Kondo effect then has unusual properties such as a strong spin polarization of the Kondo resonance and, just as for ferromagnetic materials, for the density of states (DOS). Nevertheless, despite of the spin assymetries in the DOS of the QD and the leads, the symmetry in the occupancy $`n_{}=n_{}`$ implies that the system’s ground state can be tuned to have a fully compensated local spin, in which case the QD conductance is found to be the *same* for each spin channel, $`G_{}=G_{}`$.

The Model. – For ferromagnetic leads electron-electron interactions in the leads give rise to magnetic order and spin-dependent DOS $`\rho _r(\omega )\rho _r(\omega )`$, $`r=L,R`$. Magnetic order of typical band ferromagnets like Fe, Co, and Ni is mainly related to electron correlation effects in the relatively narrow $`3d`$ sub-bands, which only weakly hybridize with $`4s`$ and $`4p`$ bands nolting . We can assume that due to a strong spatial confinement of $`d`$ electron orbitals, the contribution of electrons from $`d`$ sub-bands to transport across the tunnel barrier can be neglected tsymbal . In such a situation the system can be modeled by noninteracting interaction  $`s`$ electrons, which are spin polarized due to the exchange interaction with uncompensated magnetic moments of the completely localized $`d`$ electrons. In mean-field approximation one can model this exchange interaction as an effective molecular field, which removes spin degeneracy in the system of noninteracting conducting electrons, leading to a spin-dependent DOS. The Anderson model (AM) for a QD with a single energy level $`ϵ_\mathrm{d}`$, which is coupled to ferromagnetic leads, is given by

$`\stackrel{~}{H}`$ $`=`$ $`{\displaystyle \underset{rk\sigma }{}}ϵ_{rk\sigma }c_{rk\sigma }^{}c_{rk\sigma }+ϵ_\mathrm{d}{\displaystyle \underset{\sigma }{}}\widehat{n}_\sigma +U\widehat{n}_{}\widehat{n}_{}`$ (1)
$`+{\displaystyle \underset{rk\sigma }{}}(V_{r,k}d_\sigma ^{}c_{r,k\sigma }+V_{r,k}^{}c_{r,k\sigma }^{}d_\sigma )g\mu _BBS_z.`$
Here $`c_{rk\sigma }`$ and $`d_\sigma `$ ($`\widehat{n}_\sigma =d_\sigma ^{}d_\sigma `$) are the Fermi operators for electrons with momentum $`k`$ and spin $`\sigma `$ in the leads ($`r=L/R`$), and in the QD, $`V_{rk}`$ is the tunneling amplitude, $`S_z=(\widehat{n}_{}\widehat{n}_{})/2`$, and the last term is the Zeeman energy of the dot. In general, all information about spin asymmetry in the leads can be modeled by the spin-dependent hybridization function $`\mathrm{\Gamma }_{r\sigma }(\omega )=\pi _k\delta (\omega ϵ_{k\sigma })V_{r,k}^2=\pi \rho _{r\sigma }(\omega )V_r^2`$, where $`V_{r,k}=\mathrm{const}V_r`$ and $`\rho _{r\sigma }(\omega )`$ is the spin-dependent DOS.

In order to understand the Kondo physics of a QD attached to two identical ferromagnetic electrodes with parallel configurations, it suffices to study, instead of the above general model \[Eq. (1)\], a simpler one, which captures the same essential physics, namely the fact that $`\mathrm{\Gamma }_r(\omega )\mathrm{\Gamma }_r(\omega )`$ will generate an effective local magnetic field, which lifts the degeneracy of the local level (even for $`B=0`$). A simple (but not unique) way of modeling this effect is to take the DOS in the leads to be constant and spin-independent, $`\rho _{r\sigma }(\omega )\rho `$, the bandwidths to be equal $`D_{}=D_{}`$, and lump all spin-dependence into the spin-dependent hybridization function, $`\mathrm{\Gamma }_{r\sigma }(\omega )`$, which we take to be $`\omega `$-independent, $`\mathrm{\Gamma }_{r\sigma }(\omega )\mathrm{\Gamma }_{r\sigma }`$. By means of a unitary transformation glazman  the AM \[Eq. (1)\] can be mapped onto a model in which the correlated QD level couples only to one electron reservoir described by Fermi operators $`\alpha _{k\sigma }`$ with strength $`\mathrm{\Gamma }_\sigma =_r\mathrm{\Gamma }_{r\sigma }`$,

$`H`$ $`=`$ $`{\displaystyle \underset{k\sigma }{}}ϵ_k\alpha _{k\sigma }^{}\alpha _{k\sigma }+ϵ_\mathrm{d}{\displaystyle \underset{\sigma }{}}\widehat{n}_\sigma +U\widehat{n}_{}\widehat{n}_{}`$
$`+{\displaystyle \underset{k\sigma }{}}\sqrt{{\displaystyle \frac{\mathrm{\Gamma }_\sigma }{\pi \rho }}}(d_\sigma ^{}\alpha _{k\sigma }+\alpha _{k\sigma }^{}d_\sigma )g\mu _BBS_z.`$
Finally, we parametrize the spin-dependence of $`\mathrm{\Gamma }_\sigma `$ in terms of a spin-polarization parameter $`P(\mathrm{\Gamma }_{}\mathrm{\Gamma }_{})/\mathrm{\Gamma }`$, by writing $`\mathrm{\Gamma }_{()}\frac{1}{2}\mathrm{\Gamma }(1\pm P)`$, where $`\mathrm{\Gamma }\mathrm{\Gamma }_{}+\mathrm{\Gamma }_{}`$.

In the model of Eq. (Kondo effect in the presence of itinerant-electron ferromagnetism  studied with the numerical renormalization group method), we allowed for $`\mathrm{\Gamma }_{}\mathrm{\Gamma }_{}`$ but not for $`D_{}D_{}`$, as would be appropriate for real ferromagnets, whose spin-up and down bands always have a Stoner splitting $`\mathrm{\Delta }DD_{}D_{}`$, with typical values $`\mathrm{\Delta }D/D_{}20\%`$ (for Ni, Co, Fe). However, no essential physics is thereby lost, since the consequence of taking $`D_{}D_{}`$ is the same as that of taking $`\mathrm{\Gamma }_{}\mathrm{\Gamma }_{}`$, namely to generate an effective local magnetic field Stoner .

The occurrence of the Kondo effect requires spin fluctuations in the dot as well as zero-energy spin-flip excitations in the leads. Indeed, a Stoner ferromagnet without full spin polarization $`1<P<1`$ provides zero-energy Stoner excitations yosida , even in the presence of an external magnetic field.

Method. – The model \[Eq. (Kondo effect in the presence of itinerant-electron ferromagnetism  studied with the numerical renormalization group method)\] can be treated by Wilson’s NRG method. This method, with recent improvements related to high-energy features and finite magnetic field hofstetter2 ; costi , is a well-established method to study the Kondo impurity (QD) physics. It allows one to calculate the level occupation $`n_\sigma \widehat{n}_\sigma `$ (a static property), the QD spin spectral function, $`\mathrm{Im}\chi _s^z(\omega )=\left\{i\mathrm{\Theta }(t)[S_z(t),S_z(0)]\right\}`$, where $``$ denotes the Fourier transform, and the spin-resolved single-particle spectral density $`A_\sigma (\omega ,T,B,P)=\frac{1}{\pi }\mathrm{Im}𝒢_{d,\sigma }^R(\omega )`$ for arbitrary temperature $`T`$, magnetic field $`B`$ and polarization $`P`$ \[where $`G_\sigma ^R(\omega )`$ denotes a retarded Green function\]. From this we can find the spin-resolved conductance $`G_\sigma =\frac{e^2}{\mathrm{}}\frac{2\lambda }{(\lambda +1)^2}\mathrm{\Gamma }_\sigma _{\mathrm{}}^{\mathrm{}}𝑑\omega A_\sigma (\omega )(\frac{f(\omega )}{\omega })`$ with $`f(\omega )`$ denoting the Fermi function and $`\lambda =\mathrm{\Gamma }_{L\sigma }/\mathrm{\Gamma }_{R\sigma }`$ describing the asymmetry in the couplings. We choose $`\lambda =1`$ below.

Generation and restoration of spin splitting. – In this Letter, we focus exclusively on the properties of the system at $`T=0`$ in the local moment - Kondo regime, where the total occupancy of the local level, $`n=_\sigma n_\sigma 1`$. The occurence of charge fluctuations of course broadens and shifts the position of the QD levels (for both spin up and down), and hence changes their occupation. For $`P0`$, the charge fluctuations and hence level shifts and level occupations become *spin-dependent*, causing the QD level to split martinek1  and the dot magnetization $`n_{}n_{}`$ to be finite (Fig. 1). As a result, the Kondo resonance, too, is spin-split hofstetter3 ; exception  and weakened (Fig. 2), similarly to the effect of an applied magnetic field hofstetter2 ; costi . This means that Kondo correlations are reduced or even suppressed in the presence of ferromagnetic leads. However, for any fixed $`P`$, it is possible to compensate the splitting of the Kondo resonance \[Fig. 2(c) and 3(c)\] by fine-tuning the magnetic field to an appropriate value, $`B_{\mathrm{comp}}(P)`$, defined as the field which maximizes the height of the Kondo resonance. This field is found to depend linearly linear  on $`P`$ \[Fig. 1(c)\] (as predicted in martinek1  for $`U\mathrm{}`$). Remarkably, we also find (throughout the local moment regime) that at $`B_{\mathrm{comp}}`$ the local occupancies satisfy $`n_{}=n_{}`$ \[Fig. 1(b)\]. The fact that this occurs simultaneously with the disappearance of the Kondo resonance splitting suggests that the local spin is fully screened at $`B_{\mathrm{comp}}`$.

Spectral functions. – We henceforth fix the magnetic field at $`B=B_{\mathrm{comp}}(P)`$. To learn more about the properties of the corresponding ground state, we computed the spin spectral function $`\mathrm{Im}\{\chi _s^z(\omega )\}`$ for several values of $`P`$ at $`T=0`$ (Fig. 3). Its behavior is characteristic for the formation of a local Kondo singlet: as a function of decreasing frequency, the spin spectral function shows a maximum at a frequency $`\omega _{\mathrm{max}}`$ which we associate with the Kondo temperature \[i.e. $`k_BT_\mathrm{K}\mathrm{}\omega _{\mathrm{max}}`$ at $`B=B_{\mathrm{comp}}(P)`$\], and then decreases linearly with $`\omega `$, indicating the formation of the Fermi liquid state hewson-book . By determining $`T_\mathrm{K}(P)`$ (from $`\omega _{\mathrm{max}}`$) for different $`P`$-values, we find that $`T_\mathrm{K}`$ decreases with increasing $`P`$ \[Fig. 3(b)\]. For metals like Ni, Co, and Fe, where $`P=0.24`$, $`0.35`$, and $`0.40`$ respectively, the decrease of $`T_\mathrm{K}`$ is rather weak, so the Kondo effect should still be experimentally accessible. Remarkably, both $`\mathrm{Im}\{\chi _s^z(\omega )\}`$ and the spectral function $`A_\sigma (\omega )`$ collapse rather well onto a universal curve if plotted in appropriate units \[Figs. 3(a) and 3(c)\]. This indicates that an applied magnetic field $`B_{\mathrm{comp}}`$ restores the universal behavior characteristic for the isotropic Kondo effect, in spite of the presence of spin-dependent coupling to the leads. Fig. 4(a) shows that the amplitude of the Kondo resonance is strongly spin dependent for ferromagnetic leads, which is unusual and unique. The nature of this asymmetry is related to the asymmetry of the DOS in the leads, and its value is exactly proportional to $`1/\mathrm{\Gamma }_\sigma `$. As a result the total conductance $`G_\sigma `$ is not spin dependent \[Fig. 4(b)\]. This indicates the robustness of the Kondo effect in this system: if the external magnetic field has been tuned appropriately, it is able to compensate the presence of a spin asymmetry in the leads by creating a proper spin asymmetry in the dot spectral density, thereby conserving a fully compensated local spin and achieving perfect transmission.

Friedel sum rule. – Further insights can be gained from the Friedel sum rule, an exact $`T=0`$ relation langreth  that holds for arbitrary values of $`P`$ and $`B`$ friedel . The interacting Green’s function can be expressed as hewson-book  $`𝒢_{d,\sigma }^R(\omega )=[\omega ϵ_{\mathrm{d}\sigma }+i\mathrm{\Gamma }_\sigma \mathrm{\Sigma }_\sigma (\omega )]^1`$, with spin-dependent $`ϵ_{\mathrm{d}\sigma }`$ and $`\mathrm{\Gamma }_\sigma `$; the former due to Zeeman splitting ($`ϵ_{\mathrm{d}\sigma }=ϵ_\mathrm{d}1/2\sigma g\mu _BB`$) and latter due to the ferromagnetic leads. Here $`\mathrm{\Sigma }_\sigma (\omega )`$ denotes the spin-dependent self energy. Now, the Friedel sum rule langreth  implies that at $`T=0`$, the occupancy $`n_\sigma `$ and spectral functions can be written as

$`n_\sigma =\varphi _\sigma /\pi ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\pi }}\mathrm{tan}^1\left({\displaystyle \frac{ϵ_{\mathrm{d}\sigma }ϵ_\mathrm{F}+\mathrm{\Sigma }_\sigma ^\mathrm{R}(ϵ_\mathrm{F})}{\mathrm{\Gamma }_\sigma }}\right),`$ (3)
$`A_\sigma (ϵ_\mathrm{F})={\displaystyle \frac{\mathrm{sin}^2(\pi n_\sigma )}{\pi \mathrm{\Gamma }_\sigma }},`$ (4)
where $`\mathrm{\Sigma }_\sigma ^\mathrm{R}(\omega )\mathrm{Re}\mathrm{\Sigma }_\sigma (\omega )`$, and $`\varphi _\sigma (\omega )`$ is the phase of $`𝒢_{d,\sigma }^R(\omega )`$. Since $`\mathrm{\Sigma }_{}^\mathrm{R}(ϵ_\mathrm{F})\mathrm{\Sigma }_{}^\mathrm{R}(ϵ_\mathrm{F})`$, an equal spin occupation, $`n_{}=n_{}`$, is possible only for $`(ϵ_\mathrm{d}ϵ_\mathrm{F}+\mathrm{\Sigma }_{}^\mathrm{R}(ϵ_\mathrm{F}))/\mathrm{\Gamma }_{}=(ϵ_\mathrm{d}ϵ_\mathrm{F}+\mathrm{\Sigma }_{}^\mathrm{R}(ϵ_\mathrm{F}))/\mathrm{\Gamma }_{}`$, which can be obtained only for an appropriate external magnetic field $`B=B_{\mathrm{comp}}`$. For the latter, in the local moment regime ($`n1`$) we have $`n_{}=n_{}0.5`$, so that $`\varphi _{}=\varphi _{}\pi /2`$, which implies that the peaks of $`A_{}`$ and $`A_{}`$ are aligned. Thus, the Friedel sum rule clarifies why the magnetic field $`B_{\mathrm{comp}}`$ at which the splitting of the Kondo resonance disappears, coincides with that for which $`n_{}=n_{}`$.

For $`B=B_{\mathrm{comp}}`$, the spin-dependent amplitude $`A_\sigma (ϵ_\mathrm{F})`$ of Eq. (4) and the conductance $`G_\sigma \mathrm{\Gamma }_\sigma A_\sigma (ϵ_\mathrm{F})`$ agree well with the abovementioned NRG results \[Fig. 4(a), 4(b)\].

In conclusion, applying the NRG technique, we showed that the Kondo effect in a QD attached to ferromagnetic leads is in general suppressed, because the latter induce a spin splitting of the QD level, which leads to an asymmetry in the occupancy $`n_{}n_{}`$. Remarkably, the Kondo effect may nevertheless be restored by applying an external magnetic field $`B_{\mathrm{comp}}`$, tuned such that the splitting of the Kondo resonance is compensated and the condition $`n_{}=n_{}`$ is fulfilled. Although the Kondo resonance is strongly spin polarized, it then features a locally-screened state, a spin-independent conductance, and a Kondo temperature which decreases with increasing spin asymmetry.

Acknowledgements. – We thank R. Bulla, T. Costi, L. Glazman, W. Hofstetter, H. Imamura, B. Jones, S. Maekawa, A. Rosch, M. Vojta, and Y. Utsumi for discussions. This work was supported by the DFG under the CFN, ’Spintronics’ RT Network of the EC RTN2-2001-00440, Project PBZ/KBN/044/P03/2001, OTKA T034243 and the Emmy-Noether program.

Note added. – After submission of our paper, a preprint \[M. S. Choi et al. cond-mat/0305107\] studying a similar problem using the NRG technique appeared, with conclusions consistent with our’s.


format: text

text:
# Comments on the 2001 run of the EXPLORER/NAUTILUS gravitational wave experiment

## 1 Comments on the experimental results

At present, the 2001 scientific run performed by our gravitational wave (GW) detectors EXPLORER and NAUTILUS  constitutes the most sensitive experiment for the detection of GW bursts ever carried out. When we analyzed the data collected, we became convinced of the importance of communicating that unprecedented sensitivities were reached by our detectors and that the use of powerful tools in analysis of the data (energy consistency of the events and sidereal time analysis) led to a coincidence excess centered at sidereal hour 4, which was defined in the paper “an interesting indication”.

Criticisms to this analysis have been recently published . A first criticism regards the mere existence of a coincidence excess. It is stated that the potential significance of the peak centered at sidereal hour 4 must be “diluted” because we did not declare before the analysis (i.e. “a priori”) that we were searching for sources giving excess in that sidereal hour (i.e. galactic sources). In other words, the significance of the peak was recognized “a posteriori”. As a consequence, the probability should be decreased of roughly a factor 24 because a priori the peak could have been found in any of the 24 sidereal hours.

To be quantitative about what we called an “interesting indication”, let us report (see Table 1) the Poisson probabilities of the observed number of coincidences with respect to the corresponding number of accidentals in intervals of increasing duration centered at sidereal time 4.

These probabilities are, at most, at the level of a few percent. Instead, going up to the entire 24 hour period, the overall number of coincidences (31) with respect to the accidentals (25) gives a Poisson probability of 0.14.

We are well aware that extraordinary claims need extraordinary evidence. By simply looking at these Poisson probabilities, we see that the coincidence excess we reported was certainly not strong enough to claim a detection, and in fact no such claim was made in our paper. However, two remarkable facts need to be underlined and give sense to the words “interesting indication”:

* the peak is centered at a physically significant sidereal hour, corresponding to the most favorable orientation of the detectors with respect to sources in the Milky Way. The Galaxy is certainly the privileged place of the sources attainable by present GW detectors and we think that the experiment described in  should be considered as based on the “a priori” hypothesis of signals originating in the Galaxy. This was clearly indicated in our previous paper  (pag. 248), “No extragalactic GW signals should be detected with the present detectors. Therefore we shall focus our attention on possible sources located in the Galaxy.” Recent work by Paturel and Baryshev , quantitatively indicated the signature expected from galactic sources in bar detectors of different sensitivity, and in particular the presence of a peak centered at sidereal hour 4 for EXPLORER and NAUTILUS. So we don’t think that the reported Poisson probabilities should be further “diluted” by a factor taking into account all the possible peak positions in sidereal time.
* there is a strong energy correlation for the coincidence events during the sidereal hours interval 3 to 5 (see figure 8 of reference ). For those events the output of the two detectors have a correlation coefficient of 0.96 and the slope of the linear regression line for NAUTILUS energy on EXPLORER energy is 1.18, which, within the accuracy of the detector calibration, is in agreement with the hypothesis of having equal signals on the two parallel bar detectors. On the contrary, the events at the other sidereal hours exhibit a correlation coefficient of -0.19 (a very poor correlation, as expected for random values).

Another criticism is that, since a bar detector has a rather broad antenna pattern, one should expect a correspondingly broad peak in the sidereal time distribution of the coincidences, instead of the relatively narrow peak reported in . This argument is very weak, because, given a source distribution, the geometrical antenna pattern of a single detector is not representative of the distribution of the coincidences between two detectors. Each detector has its own detection efficiency, which depends on the detector noise, threshold, and on the signal level. The coincidence distribution depends on the product of the efficiencies of the two detectors, and may give a relatively narrow peak if, for instance, the signals are near the thresholds. The fact that the width of the peak depends on the signal level and detector threshold was already shown for a single detector in reference , where the case of EXPLORER was explicitly considered.

For these reasons, we disagree with the objections to our paper, even if we agree that no claim of detection can be made with our data.

Last comment: we were told that some statements in our paper about the possible contribution of real GW signals to our data led the reader to think that we were claiming a discovery. We admit that, in spite of our conservative attitude, in a couple of phrases the possible presence of gw signals in our data was considered. On the other hand it should be recognized by any experimentalist that it is difficult to take data with detectors of unprecedented sensitivity, containing the reported indication, without at least being open to this possibility.

## 2 The future

In order to calculate the probability, or the degree of belief, of having observed real GW signals (and not simply the Poisson probability that a background fluctuation produces $`n_c`$ or more counts), one should introduce astrophysical models for the sources. We deliberately left out this step in our paper, for two reasons: because we wanted to keep separate the observations from the possible astrophysical interpretations and, moreover, because we thought that more data were needed to evaluate the (unexpected) possibility that many gw bursts at the level of $`h210^{18}`$ are bathing the Earth.

The exercise presented by Astone  at the workshop GWDAW 2002 outlines how we can proceed in the future: updating the degree of beliefs adding new pieces of information as more data become available. That paper is an important contribution to the debate (now starting) on how to evaluate and compare the results reported by different gravitational wave experimental collaborations.

In conclusion, let us say that we are proud to have crossed the benchmark of “nihil obstat” upper limit on the strength of the wave: we find that the amplitudes and rate of our candidate events are not only compatible with existing experimental upper limits  but also they are permitted by theoretical upper limits on gravitational wave strengths based on cherished beliefs about the astrophysical structure of the Galaxy and about the physical laws governing gravitational radiation . This was not the case of the events observed by Weber with its first generation bars at room temperatures, which amplitudes were at least two orders of magnitude larger and which can be explained only by invoking unconventional hypotheses, like strong beaming by sources near the galactic center, or today being a very special time in the evolution of the Galaxy .

We remark that the indication we reported is not destined to remain an intriguing observation (and source of trouble …) for long: the existence of bursts bathing the Earth at the level of $`h210^{18}`$ and at a rate of many per year would be such an unexpected gift from nature that it will easily be confirmed or denied soon with interferometers and bars operating at their expected sensitivity.

## References


text:
# Recent Microlensing Results from the MACHO Project

## 1. The MACHO Microlensing Experiment

The Milky Way is one of the two most massive galaxies in the Local Group and so is likely to represent a typical product of cosmological evolution of matter trapped in a relatively deep potential well. The detailed formation processes of our Galaxy put constraints on cosmology, star formation history, environmental influence on galactic evolution etc. The key to understanding the formation process of the Galaxy is to understand its structure and stellar populations. There are still many questions related to the structure of the Galaxy that have not been satisfactory answered. “What is the Galactic dark halo made off?” is probably the most momentous question but the others related to the inner structure of our Galaxy or the origin of different components also have profound consequences for our understanding of the galaxy formation.

Microlensing is one of the preferred techniques to study the Milky Way because of its sensitivity to both luminous and dark matter. Microlensing results from the bending of light by massive objects. The signal depends sensitively on the alignment between the observer, the source of light (at a distance $`D_s`$), and a lens located between them (at a distance $`D_l`$). Strong signal is observed if at some time the projected distance between the source and the lens is of order of the Einstein radius, $`r_E`$, defined as:

$$r_E=\sqrt{\frac{4Gm_{\mathrm{lens}}}{c^2}\frac{(D_sD_l)D_l}{D_s}},$$
(1)
and the event produces significant magnification on a time scale given by

$$t_E=\frac{r_E}{v_{}},$$
(2)
where $`m_{\mathrm{lens}}`$ is the mass of the lens, and $`v_{}`$ a relative transverse velocity. The observable quantity, Einstein radius crossing time, $`t_E`$, is a degenerate combination of $`m_{\mathrm{lens}}`$, $`v_{}`$, $`D_l`$ and $`D_s`$. Fortunately, microlensing phenomenology goes much beyond the simple point lens approximation. The effects of ’parallax’ (Gould 1992; Alcock et al. 1995), binary caustic crossing (Mao & Paczyński 1991; Afonso et al. 2000), or finite source (Gould 1994; Alcock et al. 1997b) allow one to break degeneracies present in the simplest cases and provide useful constraints on stellar physics. Below we will use the effect of parallax, from which one may infer the projected transverse velocity:

$$\stackrel{~}{v}=v_{}\frac{D_s}{D_sD_l}.$$
(3)
Substitution of equations (2) and (3) into (1) results in:

$$M=\frac{c^2}{4G}t_E^2\stackrel{~}{v}^2\frac{D_sD_l}{D_sD_l}\frac{c^2}{4G}t_E^2\stackrel{~}{v}\mu _{\mathrm{rel}},$$
(4)
where in the second equity we simply recognized the fact that the relative proper motion of the lens and the source was equal to $`\mu _{\mathrm{rel}}=v_{}/D_l`$. When we have $`\stackrel{~}{v}`$ from the parallax fit to the microlensing light curve, then the only uncertainty in the mass comes from the lack of knowledge of the distance to the lens (distance to the source is oftentimes approximately known). With the additional measurement of the relative proper motion, one solves for the mass with no ambiguity.

The amount of mass between the source and observer and its distribution along the line of sight is described by the optical depth, which is the probability that microlensing happens on a given star (typically of order $`10^710^6`$). The actual count of events (microlensing rate) is proportional to the optical depth and to the number of sources. Therefore, microlensing surveys target concentrations of stars to maximize the number of detected events. Microlensing observations toward the LMC and SMC test the existence and composition of the lensing population in the Galactic dark halo. The microlensing optical depth toward the bulge is a function of the matter content of the inner Galaxy and is sensitive to the characteristics of the Galactic model.

Less than two decades ago, microlensing was still at a stage of a purely theoretical concept (Paczyński 1986, 1991; Griest 1991). It became such a powerful astrophysical tool thanks to the experiments like the MACHO Survey, which had a major impact on the development and application of the microlensing technique to Galactic studies. The MACHO experiment collected images of the Galactic bulge and Magellanic Clouds from 1992 through 1999. All observations were taken with the 1.3 m Great Melbourne Telescope with a dual-color wide-field camera. The MACHO camera consisted of two sets of four 2k x 2k CCDs that collected blue ($`B_M`$) and red ($`R_M`$) images simultaneously using a dichroic beam splitter. A single observed field covered an area of 43’ by 43’. We observed 82 fields in the LMC, 21 in the SMC and 94 in the Galactic bulge region, and collected about 100000 images. The fields with high priority were observed on most photometric nights in the season. The LMC observations were carried out all year long, whereas the bulge season lasted from March till October. The resulting light curves often contain several hundred points with individual fields sampled on average every few days. Details of the MACHO telescope system are given by Hart et al. (1996) and of the camera system by Stubbs et al. (1993) and Marshall et al. (1994). Details of the MACHO imaging, data reduction and photometric calibration are described by Alcock et al. (1999).

## 2. Where are the source stars of the LMC events?

We (Alcock et al. 2000b) used two sets of microlensing candidates: conservative set A of 13 events and inclusive set B of 17 events to argue that Massive Compact Halo Objects (MACHOs) contribute of order of 20% to the mass of the dark halo of the Milky Way. Using maximum likelihood technique, we showed that the most likely mass of a typical lens is $`0.5M_{}`$. Upper limits on the halo contribution of MACHOs consistent with our results were also obtained by the EROS collaboration (Lasserre et al. 2000) based on 3 events in the LMC or 4 events in both Magellanic Clouds. Neither of the groups had enough statistical power to reveal the nature or the location of the lenses along the line of sight. Figure 1 summarizes the MACHO results for a standard model of the Galactic dark halo.

In all cases when we have some additional constraints on the location of the sources, they seem to belong to the known stellar systems (Afonso et al. 2000; Alcock et al. 2001a; Alcock et al. 2001c), which makes some to believe that all of the events must be of this nature. This probability argument is however faulty to some extent. Most effects that give the indication of the lens position occur preferentially when the lens is close to either observer (parallax) or the source (xallarap). Therefore, they provide little information about the location or population membership of the remaining events.

Being unable to conclusively decide on the location of the lenses from available microlensing information, we will now attempt to obtain some indirect bounds. The idea we want to entertain here is the following: find the location of the sources first and then infer the location of the lenses. A useful demonstration of how this works would be to consider two situations with somewhat oversimplified conclusions: 1) sources are in the LMC $``$ lenses are in the Milky Way, 2) sources are behind the LMC $``$ lenses are in the LMC. The standard case 2) was considered by Zhao, Graff, & Guhathakurta (2000), who placed the source population $`7`$ kpc behind the main body of the LMC and assumed uniform LMC extinction of $`E(BV)=0.13`$. We will not analyze all the possible locations of the sources and the implications for the lens population here, but we note that this issue is discussed in detail by Alcock et al. (2001b) and Nelson et al. (2003).

The sky appearance at the MACHO observing site at Mount Stromlo is very similar to the one portrayed in “Starry Night”, the oil canvas painted by Vincent van Gogh in Saint–Rémy in 1889111Color reproduction available at http://www.vangoghgallery.com/painting/p\_0612.htm.. With median seeing above 2 arc-seconds, the individual stars in the LMC are typically not resolved. The objects observed are blends of several stars of different brightness (Alcock et al. 2001b,d). Therefore, we use Hubble Space Telescope (HST) images to construct true color-magnitude diagrams (CMDs) and Difference Image Analysis (DIA) to establish magnitudes and colors of our source stars. The HST and MACHO frames are put on the same astrometric system, and then the centroid of the source star is recovered based on the motion of the centroid of light during the event. The HST star closest to the recovered centroid is assumed to be the source star. This procedure is illustrated in Figure 2. The DIA analysis applied to the entire set A uniquely determines the positions of 12 out of 13 microlensing events. In one case the recovered centroid of source’s light falls exactly in the middle between the two HST stars. Fortunately, they have very similar colors and apparent magnitudes and so we choose one of them at random.

Having obtained the unblended characteristics of the source stars of microlensing events we now have to model the source star populations MACHO experiment is sensitive to. We assume that our source stars can be drawn from two population: the bar and disk of the LMC and a background population behind the LMC. We assume that the background population does not differ intrinsically from the one in the LMC. The difference in the CMD appearance comes from additional reddening (here taken to be the mean reddening of stars in the LMC) and the shift in the distance modulus. We start with the HST CMD and convolve it with the MACHO efficiency for detecting stars, which goes to almost zero at $`V=22`$, and cuts off a substantial fraction of the HST CMD. In this way we reproduce the unblended LMC population observed by the MACHO experiment. The route to mock the background population is more complicated. Each star in the HST CMD must be first shifted by $`\mathrm{\Delta }(VI)=E(VI)=1.38E(BV)=0.18`$ mag and $`\mathrm{\Delta }V=A_V+\mathrm{\Delta }_{\mathrm{BKG}}=0.43+\mathrm{\Delta }_{\mathrm{BKG}}`$ mag, where $`\mathrm{\Delta }_{\mathrm{BKG}}`$ is the distance modulus shift of the background population with respect to the LMC. Only then we modify the HST CMD by “filtering” the number of stars through our efficiency curve as in the first case. Based on these pure LMC and background populations, we construct a series of CMDs with a fraction $`f_{BKG}`$ stars from the background and $`(1f_{BKG})`$ from the main body of the LMC. The two extreme cases of $`f_{BKG}=0`$ and $`f_{BKG}=1`$ are shown in Figure 3. Spherical symbols with central dots are the sources of the 13 microlensing events, superposed on the underlying population.

We perform a two-dimensional Kolmogorov-Smirnov (KS) test to tell, which models are most consistent with the data. In one dimensional case, a KS test of two samples with number of points $`N_1`$ and $`N_2`$ returns a statistic $`D`$, defined to be the maximum distance between the cumulative probability functions at any ordinate. Associated with $`D`$ is a corresponding probability $`P(D)`$ that if two random samples of size $`N_1`$ and $`N_2`$ are drawn from the same parent distribution a worse value of $`D`$ will result. This is equivalent to saying that we can exclude the hypothesis that the two samples are drawn from the same distribution at a confidence level of $`1.0P(D)`$. If $`N_2N_1`$, then this is also equivalent to excluding at a $`1.0P(D)`$ confidence level the hypothesis that sample 1 is drawn from sample 2. The concept of a cumulative distribution is not defined in more than one dimension. However, it has been shown that a good substitute in two dimensions is the integrated probability in each of four right-angled quadrants surrounding a given point (Peacock 1983; Fasano & Franceschini 1987). The integrated probability of each quadrant for each distribution is the fraction of data from the distribution which lies in that quadrant. The two dimensional KS statistic $`D`$ is now taken to be the maximum difference (ranging both over all data points and all quadrants) between the integrated probability of distributions $`N_1`$ and $`N_2`$. The statistic $`D`$ and the corresponding $`P(D)`$ are subject to the same interpretation as in the one dimensional case (see Press et al. 1992 for details).

We show the resulting $`P(D)`$ versus $`f_{BKG}`$ in Figure 4, where we used the distance modulus shift $`\mathrm{\Delta }_{\mathrm{BKG}}=0.3`$. More models are considered by Nelson et al. (2003). Because the creation of the efficiency convolved CMD is a weighted random draw from the HST CMD, the model population created in each simulation differs slightly. This in turn leads to small differences in the KS statistics. The error bars indicate the scatter about the mean value for 20 simulations. We find that the 2-D KS-test probability is highest at a fraction of source stars behind the LMC $`f_{BKG}0.00.2`$. Most of the features of considered CMDs are almost vertical. As a result, the significant preference for low $`f_{BKG}`$ arises mostly from the reddening term (see discussion below). In a strict statistical sense, the Kolmogorov-Smirnov test can only be used to reject or confirm the hypothesis that two datasets have a common parent distribution. KS test results do not allow one to draw conclusions about the relative probability, e.g., a KS test probability of 90% does not mean that it is three times as likely that two distributions are identical than if they had returned a KS test probability of 30%. Therefore, Figure 4 should be interpreted with caution. Based on the models analyzed by Nelson et al. (2003), we conclude that the general shape of the KS probability versus $`f_{BKG}`$ is quite insensitive to the distance modulus displacement.

We rule out a model in which the source stars all belong to some background population at more than a 95% confidence level (e.g., the KS test probability for the distance modulus displacement of $`\mathrm{\Delta }_{\mathrm{BKG}}=0.3`$ for $`f_{BKG}>0.9`$ is $`0.025`$). The less extreme models of the LMC spheroid or stellar shroud self-lensing are not excluded. The allowed region of the $`f_{BKG}`$ plot ($`P>10`$%) is consistent with the expected location of the source stars in both the Milky Way lensing and LMC disk+bar self-lensing geometries. However, detailed modeling of the LMC disk+bar self-lensing suggests that it contributes at most 20% of the observed optical depth (Gyuk, Dalal, & Griest 2000). Therefore, the results of the KS test taken together with external constraints, could suggest that the lens population comes mainly from the Milky Way dark halo. This statement, however, would be too strong. Our results are not completely conclusive due to the intrinsic problem of the method. Its serious limitation is the uncertainty in the model of the LMC reddening. The significance of our test will decrease substantially if the average reddening of the LMC is significantly lower or the reddening is very patchy.

## 3. Direct identification of a microlens

The HST Wide Field Planetary Camera (PC) observations described in the previous section included imaging of the LMC-5 event (Alcock et al. 2001c). Image of this event in HST V, R, and I filters was taken on May 13, 1999, 6.3 years after the peak on February 5, 1993. The system, composed of two objects displaced by 0.134”, was resolved with PC pixels of 0.046” (Figure 5). One of these objects was shown to be a main sequence source stars in the LMC with $`V=21.02\pm 0.06`$, $`(VR)=0.30\pm 0.09`$, and $`(VI)=0.68\pm 0.10`$. Its neighbor, marked with an arrow, is a very red, faint object with $`V=22.67\pm 0.10`$, $`(VR)=1.60\pm 0.12`$, and $`(VI)=3.18\pm 0.11`$. A chance superposition of such a red object on the LMC source star is of order of $`10^4`$. Therefore, it is tempting to assume that the red object is the lens. The light curve of event LMC-5 showed very large maximum flux amplification $`A_{\mathrm{max}}`$, which corresponds to very small impact parameter. As a result, we expect that at the peak, the angular distance between the lens and the source was negligible compared to the measured separation. Indeed, for large amplification events the angular distance between the lens and the source is at the peak equal to $`4.03610^4`$$`\times A_{\mathrm{max}}^1\sqrt{1/x1}\sqrt{m_{\mathrm{lens}}/M_{}}`$, where $`x=D_l/D_s`$. For LMC-5, the term $`A_{\mathrm{max}}^1\sqrt{1/x1}\sqrt{m_{\mathrm{lens}}/M_{}}`$ is very likely to be $`<1`$, and, as a result, the lens-source separation at the peak was of order of or smaller than $`10^3`$ of the current separation. Consequently, the current angular separation is an excellent measure of the relative proper motion between the source and the lens, $`\mu _{\mathrm{rel}}=0.0214\pm 0.0007`$ arcsec/yr with the direction, $`\mathrm{\Theta }_{\mathrm{rel}}=92^{}`$, given by the line connecting the light centroids.

Interestingly, LMC-5 has a significant detection of the microlensing parallax effect. We obtained the projected transverse velocity $`\stackrel{~}{v}=18`$ km/s, the direction of lens proper motion, $`\mathrm{\Theta }_{\mathrm{par}}=97^{}`$, and duration $`t_E=22.5`$ days. The agreement between $`\mathrm{\Theta }_{\mathrm{par}}`$ and $`\mathrm{\Theta }_{\mathrm{rel}}`$ suggests that the red object is indeed the lens for the microlensing event LMC-5. Equation (4) can therefore be used to estimate the mass of the lens $`m_{\mathrm{lens}}=0.036_{0.004}^{+0.009}M_{}`$, which argues for a sub-stellar object at a very high significance. The errors, however, are not Gaussian and the $`3\sigma `$ limit on the mass based on the parallax fit is $`m_{\mathrm{lens}}0.097M_{}`$. In natural units, the relative parallax of the event is $`\pi _{\mathrm{rel}}=\mu _{\mathrm{rel}}/\stackrel{~}{v}`$. For a nearby lens and a source in the LMC, $`\pi _{\mathrm{rel}}\pi _{\mathrm{lens}}`$. This allows us to derive the distance to the lens: $`d_{\mathrm{lens}}=200_{30}^{+40}`$ pc. Assuming no attenuation of light due to Galactic extinction, we obtain the absolute magnitude of $`M_V^{200\mathrm{p}\mathrm{c}}=22.675.0\mathrm{log}(200\mathrm{p}\mathrm{c}/10\mathrm{p}\mathrm{c})=16.16`$.

We took the spectrum of the LMC-5 microlensing system with the European Southern Observatory Very Large Telescope on February 2, 2001 . The separation of the lens-source system at this time was about 0.2”, which was unresolved and resulted in obtaining a composite spectrum (Figure 6). This spectrum becomes more and more dominated by the lens flux as we move redwards of 7000 Å. The presence of KI, NaI, the absence of CsI, RbI, and the TiO band at 7100 Å coupled with the absence of the VO band at 7450 Å lead to conclusion that the lens is of spectral type M4-5V. Our HST colors are consistent with this spectral classification (Bessell 1991). According to Cassisi et al. (2000) this spectral type implies $`m_{\mathrm{lens}}0.0950.013M_{}`$, which is in $`3\sigma `$ disagreement with the parallax result.

If the lensing star is an M-dwarf, then it should follow the empirical relation

$$M_V=2.89+3.37(VI)_0$$
(5)
derived by Reid (1991). Assuming zero reddening, we obtain $`M_V=13.61\pm 0.55`$ and $`d_{\mathrm{lens}}=650\pm 190`$ pc, again in disagreement with parallax result.

Can the reddening be the reason of this discrepancy? In general:

$$M_V=V5\mathrm{log}\frac{d}{10\mathrm{p}\mathrm{c}}A_V$$
(6)
where $`d`$ is the distance and $`A_V`$ is the visual extinction, which can be written as $`A_V=R_{VI}E(VI)`$. Color excess is defined as $`E(VI)(VI)(VI)_0`$. Combining equations (5) and (6) for our gravitational lens yields:

$$d_{\mathrm{lens}}(\mathrm{pc})=10^{0.2(5.0+V2.893.37(VI)+(3.372.50)E(VI))}=10^{2.815+0.174E(VI)},$$
(7)
where we assumed a standard value for the selective extinction coefficient $`R_{VI}=2.50`$. Equation (7) shows that a non-zero extinction would only increase the disagreement between the parallax fit and the spectroscopically motivated estimate of the lens distance (Figure 7). If the lens is not at the very unusual line of sight, then $`E(VI)0.1`$ and correction due to reddening is small, well within the errors of our estimate based on the assumption of zero extinction.

A heliocentric radial velocity of the lens computed using potassium and sodium lines is equal to $`v_{\mathrm{rad}}=49\pm 10`$ km/s. Assuming a distance to the lens and the proper motion inferred from the HST images we solve for the space velocity. At 200 pc, the velocity of the lens relative to the local standard of rest is $`U=27`$ km/s, $`V=43`$ km/s, and $`W=4`$ km/s, toward the Galactic center, along rotation, and up out of the plane, respectively. Therefore, the lens has disk-like kinematics. Motivated by the surprisingly faint absolute magnitude implied by the solution based on the parallax fit, Alves & Cook (2001) considered a possibility that the LMC-5 lens is a halo subdwarf. However, since metal-poor stars with disk-like kinematics are very rare, there is a strong indication that the lens belongs to the Galactic disk.

The inconsistencies between solutions based on the HST-based proper motion coupled with the MACHO parallax fit on one side and HST photometry coupled with the VLT spectrum on the other might have arisen from difficulty of measuring parallax effect for such a short event. This conflict will be resolved with observations we took on HST Advanced Camera for Surveys. Thanks to a very small pixel size of 0.027”, two epochs of data will allow us to verify the proper motion of the lens and measure the parallactic motion. Despite the described inconsistencies we believe that the red object detected in our images is the lens for the LMC-5 event and that it is a normal low-mass disk star.

## 4. Microlensing optical depth toward the Galactic bulge

After almost a decade of surveys, the constraints from the Milky Way microlensing observations imposed on the dark matter were so tight (Binney & Evans 2001) that they challenged the standard cuspy dark halo profiles (Navarro, Frenk, & White 1997; Moore et al. 1999). On the other hand, these cuspy profiles are very likely a generic feature of the cold dark matter cosmology. Below we briefly describe the measurements of the optical depth toward the Galactic bulge that led to this situation and proceeded the study described here. In brief, the studies discussed below yield the microlensing optical depth toward the Galactic bulge in the range $`\tau _{\mathrm{bulge}}=`$3$``$4 $`\times 10^6`$, 2$``$3 times higher than implied by other observational constraints and theoretical models. Udalski et al. (1994) found 9 events in the first two year of the Optical Gravitational Lensing Experiment (OGLE) data. They set the lower limit on the optical depth to the Galactic bulge at $`\tau _{\mathrm{bulge}}=(3.3\pm 1.2)\times 10^6`$. The uncertainties in this study were related to the detection efficiency analysis as well as small number statistics. Alcock et al. (1997a) described a set of 45 bulge events. They computed the so called sampling efficiencies, which are a good approximation only for bright events. Therefore, unbiased analysis could have been done only for 13 clump giants, which resulted in large uncertainties of the optical depth ($`\tau _{\mathrm{bulge}}=3.9_{1.2}^{+1.8}\times 10^6`$). Udalski et al. (2000) presented a catalog of over 200 microlensing events from the last 3 seasons of the OGLE-II bulge observations. Woźniak et al. (2001) released another catalog of about 500 events selected using image subtraction technique. Unfortunately, no efficiency analysis has been done for these two samples of OGLE events so the information that can be extracted from them is very limited. Alcock et al. (2000a) performed Difference Image Analysis of three seasons of bulge data in 8 frequently sampled MACHO fields and found 99 events. They determined $`\tau _{\mathrm{bulge}}=3.23_{0.50}^{+0.52}\times 10^6`$ at $`(l,b)=(2.^{}68,3.^{}35)`$. This was a major development in bulge microlensing. The DIA technique resulted in a substantial improvement in photometry, so this analysis was less vulnerable to uncertainties in the parameter determination. However, there are two potential problems with this analysis. First, the detection efficiency estimate suffers from the fact that a deep HST luminosity function was available for only one field. Second, the lensed sources are not guaranteed to be in the Galactic bulge. The measured optical depth is, therefore, converted to the optical depth toward the bulge using a fudge factor. Popowski (2002) showed that the reasonable modification of this factor alone can lower the estimate of the optical depth by almost 20%. To cure this uncomfortable situation we would like to select the event that have sources in the bulge. Events with clump giants as sources are excellent candidates, since clump giants are among the brightest and most numerous stars in the bulge. Being bright, clump sources have an additional nice feature – they are almost unaffected by blending.

Blending is a major problem in any analysis of the microlensing data involving point spread function photometry. The bulge fields are crowded, so that the objects observed at a certain atmospheric seeing are blends of several stars. At the same time, typically only one star is lensed. This complicates a determination of the events’ parameters and the analysis of the detection efficiency of microlensing events. If the sources are bright one can avoid these problems. First, a determination of parameters of the actual microlensing events becomes straightforward. Second, it is sufficient to estimate detection efficiency based on the sampling of the light curve alone. This eliminates the need of obtaining deep luminosity functions across the bulge fields. Therefore, in this discussion we concentrate on the events where the lensed stars are clump giants (Popowski et al. 2001a,b). We are using the results of the analysis performed on five seasons of data (1993–1997) in 77 fields.

The events with clump giants as sources have been selected from the sample of all events. The procedure that leads to a selection of microlensing events of general type consists of several steps. First, all the recognized objects in all fields are tested for any form of variability. Second, a microlensing light curve is fitted to all stars showing any variation and the objects that meet very loose selection criteria (cuts) enter the next phase. Here, this selection returns almost 43000 candidates. These candidates undergo more scrutiny and are subject to more stringent cuts, most of which test for a signal-to-noise of the different parts of the light curve. Here, this last procedure narrows a list of candidate events to $`280`$. The question, which of those sources are clump giants, is investigated through the analysis of the global properties of the color-magnitude diagram in the Galactic bulge. Using the accurately measured extinction towards Baade’s Window (Stanek 1996 with zero point corrected according to Gould, Popowski, & Terndrup 1998 and Alcock et al. 1998) allows us to locate bulge clump giants on the dereddened color-magnitude diagram. Such diagram can be then used to predict the positions of clump giants on the color-apparent magnitude diagram for fields with different extinction. Based on Baade’s Window data we conclude that unreddened clump giants are present in the color range $`(VR)_0(0.45,0.65)`$, where they concentrate along a line $`V_0=14.35+2.0(VR)_0`$. We assume that the actual clump giants scatter in $`V_0`$-mag around these central values, but by not more than 0.6 mag toward both fainter and brighter $`V_0`$. This defines the parallelogram-shaped box in the upper left corner of Figure 8. With the assumption that the clump populations in the whole bulge have the same properties as the one in the Baade’s Window, the parallelogram described above can be shifted by the reddening vector to mark the expected locations of clump giants in different fields. The solid lines are the boundaries of the region where one could find the clump giants in fields with different extinctions. There are a few more $`V`$-mag and $`(VR)`$-color cuts that determine the final shape of the clump region. Several assumptions that went into creating this region should be carefully reviewed. For example, the assumed spread in $`V`$ magnitudes can be either bigger or smaller or asymmetric around the central value, clump giants in different fields may have different characteristics etc. In retrospect, we conclude that the effect of the uncertainty in the selective extinction coefficient $`R_{VR}`$ is probably the hardest to deal with222All the selection-related problems are treated in depth in the analysis of 8 years of microlensing data.. The clump region from Figure 8 contains 52 unique clump events from the five bulge seasons discussed here. Their location in Galactic coordinates is shown in Figure 9, superposed on the outlines of the analyzed MACHO fields. The size of each circle is proportional to the duration of the event.

Based on this sample of events we estimate the bulge optical depth from the following formula:

$$\tau =\frac{\pi }{2NT}\underset{\mathrm{all}\mathrm{events}}{}\frac{t_E}{ϵ(t_E)},$$
(8)
where $`N`$ is the number of observed stars (here about 2.1 million clump giants), $`T`$ is the total exposure (here about 2000 days) and $`ϵ(t_E)`$ is an efficiency for detecting an event with a given $`t_E`$. The sampling efficiencies were obtained with the pipeline that has been previously applied to the LMC data (for a description see Alcock et al. 2001d). In brief, artificial light curves with different parameters have been added to 1% of all clump giants in the 77 considered fields and the analysis used to select real events was applied to this set. For a given duration of the artificial event, the efficiency was computed as a number of recovered events divided by a number of input events. The efficiencies used in this analysis are global efficiencies averaged over clump giants in all 77 fields. The optical depth is reported at the central position that is an average of positions of 1% of the analyzed clump giants. We obtain:

$$\tau _{\mathrm{bulge}}=(2.0\pm 0.4)\times 10^6\mathrm{at}(l,b)=(3.^{}9,3.^{}8).$$
(9)
with the error computed according to the formula given by Han & Gould (1995). In Figure 10 we plot the spatial distribution of the optical depth. The variation of the optical depth is dominated by the Poisson noise. The gradient of the optical depth is stronger in $`b`$ than in $`l`$. About 40% of the optical depth is in the long events with $`t_E>50`$ days. Such a population of long events is not easily explained by the standard models of Galactic structure and kinematics. Field 104 at $`(l,b)=(3.^{}1,3.^{}0)`$ has the largest number of clump events (10) and the optical depth of $`(1.4\pm 0.5)\times 10^5`$.

There is a high concentration of long-duration events in this field (5 out of 10 events longer than 50 days are in 104, including the longest 2). We investigate how statistically significant is this concentration. The analysis of event durations uncorrected for efficiencies provides a lower limit on the difference between the frequently-sampled field 104 and all the remaining clump giant fields. We use the Wilcoxon’s number-of-element-inversions statistic (see Popowski et al. 2001a for a description). The Wilcoxon’s statistic is equal to 320, whereas the expected number is 210 with an error of about 43. Therefore the events in 104 differ (are longer) by $`2.55\sigma `$ from the other fields. That is, the probability that events in 104 and other fields originate from the same parent population is of order of 0.011. Popowski’s (2002) analysis of the DIA events from Alcock et al. (2000a) confirms this duration difference. We conclude that field 104 is anomalous both in terms of the optical depth and duration distribution. Both effects can be explained simultaneously by the concentration of mass along this line of sight.

The currently available Galactic models are smooth on the scales of a single MACHO field. Therefore, they cannot account for such a localized anomalous behavior as in field 104, and it seems justified to remove this field as an outlier. With this modification, the optical depth decreases to $`\tau _{\mathrm{bulge}}=(1.4\pm 0.3)\times 10^6`$, which is $`2/3`$ of its original value, and becomes fully consistent with infrared-based models of the Galactic bar (Bissantz & Gerhard 2002, Evans & Belokurov 2002). This lower optical depth is supported by the very recent determination from the EROS group by Afonso et al. (2003), who report $`\tau _{\mathrm{bulge}}=(0.94\pm 0.29)\times 10^6`$ based on 16 clump events in fields centered on $`(l,b)=(2.^{}5,4.^{}0)`$. Surprisingly, another recent determination of the optical depth by the MOA collaboration (Sumi et al. 2002), based on 28 events recovered with DIA technique, gives the high optical depth of $`\tau _{\mathrm{bulge}}=3.36_{0.81}^{+1.11}\times 10^6`$. This value compares very well with the DIA result of Alcock et al. (2000a). In conclusion, there is no longer a general discrepancy between the microlensing observations and Galactic models. Instead, there is a puzzling gap between the optical depth measurements based on clump events and the samples from the Difference Image Analysis technique.

## 5. Stellar remnants as gravitational microlenses

Equations (1) and (2) show that durations of microlensing events could in principle be a good probe of lens masses. However, event’s duration depends also on the kinematics of populations involved in the lensing process. As nicely showed by Gould (2000; Figure 1), the duration distribution is very wide even for a population of lenses with a single mass. Therefore, despite the fact that $`20`$ % of events may be due to remnants in the form of white dwarfs, neutron stars and black holes, one cannot use just durations to identify this population. Events showing parallax effect (equation 4) are in high demand as they provide very stringent constraints on the lens masses.

We present the analysis of the longest timescale microlensing events with $`t_E>70`$ days (Bennett et al. 2002). We used the MACHO experiment data described in §1 as well as the follow up data from the MACHO/GMAN observations on the CTIO 0.9 m telescope (Becker 2000) and from the MPS Collaboration on the Mount Stromlo 1.9 m telescope (Rhie et al. 1999). We performed parallax fit to over 300 events detected by MACHO toward the Galactic bulge. We considered parallax signal to be reliable when the difference in the goodness of fit between the standard microlensing curve and the parallax fit was $`\mathrm{\Delta }\chi ^2200`$. This selection criterion returned 6 events: 104-C, 96-BLG-5, 96-BLG-12, 98-BLG-6, 99-BLG-1, and 99-BLG-8. The light curves and microlensing fits for the two microlensing events with likely black hole lenses, 96-BLG-5 and 98-BLG-6, are presented in Figure 11.

The CMD analysis of the events indicates that MACHO-104-C and MACHO-96-BLG-12 source stars are bulge clump giants. Source star for MACHO-99-BLG-8 is likely a bulge giant. MACHO-96-BLG-5 is likely a blended main sequence star. The CMD locations of MACHO-98-BLG-6 and MACHO-99-BLG-1 source stars suggest that they can be either bulge subgiant stars or red clump giants in the Sgr dwarf galaxy. However, their radial velocities of $`v_{\mathrm{rad}}=65\pm 2`$ km/s and $`v_{\mathrm{rad}}=64\pm 2`$ km/s, respectively are not consistent with the radial velocity of Sgr dwarf $`v_{\mathrm{rad}}^{\mathrm{SGR}}=140\pm 10`$ km/s (Ibata et al. 1997). Therefore, we concluded that they must have bulge sources.

The 6 parallax events have $`\stackrel{~}{v}`$ in the range $`3080`$ km/s. If we were to assume that the bulge and disk velocity dispersions are negligible in comparison with the Galactic rotation of about 200 km/s, then, for the disk lenses, equation (3) would imply that $`D_l=D_s\stackrel{~}{v}/(200\mathrm{km}/\mathrm{s}+\stackrel{~}{v})`$. For the sources in the bulge at $`D_s8`$ kpc and our measured $`\stackrel{~}{v}`$ this would result in the distances of the lenses between about 1 kpc and 2.3 kpc. In reality, the bulge and disk velocity dispersions modify this result to some extent, but this example shows how we can obtain constraints on the lens distance (and consequently mass) from kinematics of Galactic populations. Here we assume that the stellar remnants have the same kinematics and density distributions as the observed stellar populations. To construct likelihood functions of the lens distance given its $`\stackrel{~}{v}`$, we assume for the density profiles a standard double exponential disk and barred bulge of Han & Gould (1996). Our model disk has an isotropic velocity dispersion of 30 km/s, and a flat rotation curve of 200 km/s. The bulge is non-rotating and has an isotropic velocity dispersion of 80 km/s. The resulting likelihood curves are shown in Figure 12 as dotted lines for the two parallax events with the highest predicted lens masses. The solid line is a graphic representation of equation (4) under the assumption that the sources are in the bulge at $`D_s=8`$ kpc. Following the Bayesian method and assuming a uniform prior, we may interpret our likelihood function as a probability distribution of a certain lens mass, which allows us to determine the uncertainty in mass estimation.

We may obtain an additional constraint on the lens distances and masses if we assume that they are main sequence stars. This is achieved by finding the upper limit on the lens brightness. To this end, we use the blending fraction derived from the parallax fit. For our two best candidates: the main sequence lens is extremely unlikely for MACHO-96-BLG-5; for MACHO-98-BLG-6, the main sequence lens is not excluded and its most likely distance is about 6.1 kpc. These additional constraints are depicted in Figure 12 by short-dashed lines.

In Figure 13 we present the mass estimates for 6 bulge events with strong detections of the parallax signal. The region of likely black holes, i.e., the lenses with masses larger than $`2M_{}`$ (Akmal, Pandharipande, & Ravenhall 1998) is shaded. The 95% lower limits on the masses of two candidates that fall into this region are: 1.64 $`M_{}`$ for MACHO-96-BLG-5 and 0.94 $`M_{}`$ for MACHO-98-BLG-6. Therefore, in principle both of them could be neutron stars and for MACHO-98-BLG-6 even the main sequence star is a possibility. But we note that the measured neutron stars masses are close to the Chandrasekhar mass, $`M_{NS}=1.35\pm 0.04M_{}`$ (Thorsett & Chakrabarty 1999). Such masses are excluded for MACHO-96-BLG-5 at more than 95% confidence and for MACHO-98-BLG-6 at more than 90% confidence, which makes them excellent black hole candidates. We obtain similar conclusions when we generalize our likelihood function to include the lens mass function prior as suggested by Agol et al. (2002), but the results are very sensitive to the assumed mass function.

Together with MACHO-99-BLG-22/OGLE-1999-BUL-32 (Mao et al. 2002), the lenses of these events are the first candidates for black holes detected without observations of matter bound to the black hole. These 2 (3) events come from the considered sample of around 300 events, thus constituting around 1% of all events. However, they have very long durations, and contribute substantially more to the total amount of mass along the line of sight. If all three events are causes by black hole lenses, then the stellar mass fraction of black holes may be as high as 10%. This in turn would suggest that most of the black holes are not observed as X-ray binary systems.

## 6. Summary

We described four separate investigation into the nature of microlenses and Galactic structure based on the microlensing events detected by the MACHO collaboration. Briefly:

1. We compared the color-magnitude location of the source stars for the LMC events with composite CMDs representing different lensing scenarios. If extinction is not too patchy, then it is unlikely that most of the microlensing events in the LMC have source stars behind the LMC.
2. We performed a detailed analysis of event LMC-5 from 1993, which in 1999 was resolved into 2 stars in the HST images. We identified one of them as a source and another as a lens. The lens is a low-mass member of the Galactic disk. Its mass $`m_{\mathrm{lens}}`$ is somewhat uncertain, but likely in the range between $`0.04`$ and $`0.13M_{}`$.
3. We described the problem of the microlensing optical depth toward the Galactic bulge. From clump giant events, we obtained $`\tau _{\mathrm{bulge}}=(2.0\pm 0.4)10^6`$ when we averaged over 77 MACHO fields and $`\tau _{\mathrm{bulge}}=(1.4\pm 0.3)10^6`$ when a single highly anomalous field was excluded. Therefore, we showed that the results from clump giants are no longer in disagreement with infrared-based models of the Galactic bar. The situation is, however, still complicated due to the fact the Difference Image Analyses result in substantially higher optical depths. This difference is not understood.
4. We found several events that show clear signs of parallax effect in their light curves, and we selected the six strongest candidates. The average mass of the lenses of these events is 2.7 solar masses. The two longest events may be due to black holes with $`M_{\mathrm{bh}}`$ of order of 6 solar masses. It is very likely that the long-duration events toward the Galactic bulge are caused by massive remnants.

### Acknowledgments.

This work was performed under the auspices of the US Department of Energy, National Nuclear Security Administration, by the University of California, Lawrence Livermore National Laboratory, under contract W-7405-ENG-48. PP thanks the conference organizers for invitation to this exciting meeting.

## References

Agol, E., Kamionkowski, M., Koopmans, L.V.E., & Blandford, R.D. 2002, ApJ, 576, L131

Afonso, C., et al. 2000, ApJ, 532, 340

Afonso, C., et al. 2003, A&A, accepted (astro-ph/0303100)

Akmal, A., Pandharipande, V.R., & Ravenhall, D.G. 1998, Phys. Rev. C, 58, 1804

Alcock, C., et al. 1995, ApJ, 454, L125

Alcock, C., et al. 1997a, ApJ, 479, 119

Alcock, C., et al. 1997b, ApJ, 491, 436

Alcock, C., et al. 1998, ApJ, 494, 396

Alcock, C., et al. 1999, PASP, 111, 1539

Alcock, C., et al. 2000a, ApJ, 541, 734

Alcock, C., et al. 2000b, ApJ, 542, 281

Alcock, C., et al. 2001a, ApJ, 552, 259

Alcock, C., et al. 2001b, ApJ, 552, 582

Alcock, C., et al. 2001c, Nature, 414, 617

Alcock, C., et al. 2001d, ApJS, 136, 439

Alves, D.R. & Cook, K.H. 2001, AAS, 199.1613

Becker, A.C. 2000, Ph.D. thesis, Univ. Washington

Bennett, D.P., et al. 2002, ApJ, 579, 639

Bessell, M.S. 1991, AJ, 101, 662

Binney, J.J, & Evans, N.W. 2001, MNRAS, 327, L27

Bissantz, N., & Gerhard, O. 2002, MNRAS, 330, 591

Cassisi, S., Castellani, V., Ciarcelluti, P., Piotto, G. & Zoccali, M. 2000, MNRAS, 315, 679

Evans, N. W., & Belokurov, V. 2002, ApJ, 567, L119

Fasano, G., & Franceschini, A. 1987, MNRAS, 225, 155

Gould, A. 1992, ApJ, 392, 442

Gould, A. 1994, ApJ, 421, L71

Gould, A., Popowski, P., & Terndrup, D.T. 1998, ApJ, 492, 778

Gould, A. 2000, ApJ, 535, 928

Griest, K. 1991, ApJ, 366, 412

Gyuk, G., Dalal, N., Griest, K. 2000, ApJ, 535, 90

Han, C, & Gould, A. 1995, ApJ, 449, 521

Han, C, & Gould, A. 1996, ApJ, 467, 540

Hart, J., et al. 1996, PASP, 108, 220

Ibata, R.A., Wyse, R.F.G., Gilmore, G., Irwin, M.J., & Suntzeff, N.B. 1997, AJ, 113, 634

Lasserre, T., et al. 2000, A&A, 355, L39

Mao, S., & Paczyński, B. 1991, ApJ, 374, L37

Mao, S., et al. 2002, MNRAS, 329, 349

Marshall, S.L., et al. 1994, in IAU Symp. 161, Astronomy From Wide Field Imaging, ed. H.T. MacGillivray, et al., (Dordrecht: Kluwer)

Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J., Tozzi, P. 1999, ApJ, 524, L19

Navarro, J.F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493

Nelson, C.A., et al. 2003, AJ, submitted

Paczyński, B. 1986, ApJ, 304, 1

Paczyński, B. 1991, ApJ, 371, L63

Peacock, J.A. 1983, MNRAS, 202, 615

Popowski, P., et al. 2001a, in ASP Conference Series 239, Microlensing 2000: A New Era of Microlensing Astrophysics, eds. J.W. Menzies and P.D. Sackett (San Francisco: Astronomical Society of the Pacific), 244 (astro-ph/0005466)

Popowski, P. et al. 2001b, in ASP Conference Series 245, Astrophysical Ages and Times Scales, eds. T. von Hippel, C. Simpson, & N. Manset (San Francisco: Astronomical Society of the Pacific), 358 (astro-ph/0202502)

Popowski, P. 2002, MNRAS, submited (astro-ph/0205044)

Press, W.H., Teukolsky, S.A., Vetterling, W.T., & Flannery, B.P. 1992, Numerical Recipes in C, Second Edition

Reid, N. 1991, AJ, 102, 1428

Rhie, S.H., Becker, A.C., Bennett, D.P., Fragile, P.C., Johnson, B.R., King, L.J., Peterson, B.A., & Quinn, J. 1999, ApJ, 522, 1037

Stanek, K.Z. 1996, ApJ, 460, L37

Stubbs, C.W., et al. 1993, in Proceedings of the SPIE, Charge Coupled Devices and Solid State Optical Sensors III, ed. M. Blouke, 1900.

Sumi, et al. 2002, ApJ, accepted (astro-ph/0207604)

Thorsett, S.E. & Chakrabarty, D. 1999, ApJ, 512, 288

Udalski, A., et al. 1994, Acta Astron., 44, 165

Udalski, A., et al. 2000, Acta Astron., 50, 1

Woźniak, P., et al. 2001, Acta Astron., 51, 175

Zhao, H., Graff, D.S., & Guhathakurta, P. 2000, ApJ, 532, L37


text:
# Evidence of a higher-order singularity in dense short-ranged attractive colloids

## Abstract

We study a model in which particles interact through a hard-core repulsion complemented by a short-ranged attractive potential, of the kind found in colloidal suspensions. Combining theoretical and numerical work we locate the line of higher-order glass transition singularities and its end-point — named $`A_4`$ — on the fluid-glass line. Close to the $`A_4`$ point, we detect logarithmic decay of density correlations and sub linear power-law increase of the mean square displacement, for time intervals up to four order of magnitudes. We establish the presence of the $`A_4`$ singularity by studying how the range of the potential affects the time-window where anomalous dynamics is observed.

Recently, a great deal of interest has grown around dynamical phenomena arising in dense colloidal suspensions when particles interact through a hard-core repulsion followed by a rather short-ranged attractive potential. Some of these unusual phenomena, initially predicted theoretically, have been observed experimentally and in numerical simulations. The novelty, as compared to the well-studied hard-sphere-colloids case, comes from the possibility of generating structurally arrested states with an additional physical mechanism, driven by the range of the attractive interaction. Indeed, in the case of hard-sphere, structural arrest is controlled by the well-known excluded volume cage effect, where particle motion is hindered by the presence of neighboring ones. This mechanism generates a localization length of the order of 10 % of the particle diameter $`\sigma `$. In the additional mechanism, the motion of the particles is restricted by the adhesiveness, i.e. the formation of non-permanent bonds due to the attractive part of the inter-particle potential; particles are confined by bonds and the corresponding localization length is fixed by the attraction range $`\mathrm{\Delta }`$. When $`\mathrm{\Delta }\sigma `$, the interplay of the two localization mechanisms creates liquid states for packing fractions higher than those possible for a pure hard-sphere system. This means that the attraction stabilizes the liquid and leads to a reentry phenomenon, where this new liquid state can arrest into a glass by cooling as well as heating. The experimental and numerical confirmation of this prediction suggests that the mode-coupling theory (MCT) , which was the basis of the cited theoretical work , can contribute to explain the slow dynamics in short-ranged attractive colloids.

These systems are characterized by three control parameters — the packing fraction $`\varphi `$, the ratio $`\theta `$ of the thermal energy $`k_BT`$ to the typical well depth $`u_o`$ and the range $`\mathrm{\Delta }`$ of the attractive potential. Within MCT, the phase diagram of this three-dimensional control parameter space (shown in Fig.1) is organized around a critical point $`(\varphi ^{},\theta ^{},\mathrm{\Delta }^{})`$, referred to as a type $`A_4`$ higher-order glass transition singularity in MCT classification. $`A_4`$ is the end-point of a line of higher-order singularities (of type $`A_3`$). From a physical point of view, $`A_4`$ is characterized for being the only singularity point accessible from the liquid phase. For $`\mathrm{\Delta }>\mathrm{\Delta }^{}`$ no singular points are predicted by the theory, while for $`\mathrm{\Delta }<\mathrm{\Delta }^{}`$ the $`A_3`$ singularity points are buried in the glass phase and their presence can be only observed indirectly. Near the $`A_4`$ singularity, MCT predicts a structural relaxation dynamics, which is utterly different from that known for conventional glass-forming liquids. It is ruled by logarithmic variations of correlators as function of time and subdiffusive increase of the mean squared displacement (MSD). Theory makes precise predictions for the time-interval where this unusual dynamics is expected, as well as for its variation with changes of the $`\varphi `$, $`\theta `$ and $`\mathrm{\Delta }`$ parameters. Indications of logarithmic dynamics have been reported, but neither a systematic study of the dependence on the control parameters has been published nor the $`A_4`$ point has been located numerically. It is in fact difficult to discriminate between straight lines caused by some inflection points in the usual logarithmic representation of correlators and the logarithmic decay laws generated by the higher-order singularity. In this Letter, we establish the existence of the location of the $`A_4`$ point for a well defined model, amenable of simultaneous numerical and theoretical treatment. Reported data — with control parameters explicitly chosen close to the $`A_4`$ point — exhibit the mentioned laws over time intervals up to four orders of magnitude. More importantly, we show that the decay patterns vary with changes of the control parameters and wave-vectors as expected, properly testing the theoretical predictions.

We simulate a $`50\%50\%`$ binary mixture of $`N=700`$ hard spheres of mass $`m`$, with diameters $`\sigma _{_{AA}}`$ and $`\sigma _{_{BB}}`$ and ratio $`\sigma _{_{AA}}/\sigma _{_{BB}}=1.2`$. The hard core between particles of different type $`\sigma _{_{AB}}=0.5(\sigma _{_{AA}}+\sigma _{_{BB}})`$. The hard core potential is complemented by an attractive square well potential of depth $`u_o`$, independent on the particle type. The small asymmetry in the diameters is sufficient to prevent crystallization at high values of $`\varphi `$. We focus on two specific systems, which we label $`S_1`$ and $`S_2`$, differing in the width of the square well $`\mathrm{\Delta }_{ij}`$. $`S_1`$ and $`S_2`$ have $`\mathrm{\Delta }_{ij}=0.031\sigma _{ij}`$, and $`\mathrm{\Delta }_{ij}=0.043\sigma _{ij}`$ respectively ($`i=A,B`$). The system $`S_1`$ has been extensively studied in previous simulations. As discussed in the following, $`S_2`$ is chosen to coincide with the critical amplitude parameter within the accuracy of our calculations. Averages over five independent realizations have been performed to reduce noise. Equilibration has been carefully checked. We note that our results refer to Newtonian dynamics, but are relevant also to Brownian dynamics in the structural relaxation time window .

To estimate the predicted location of the $`A_4`$ point for the binary mixture square well model (and consequently select the parameters to be used in the simulations) we proceed as follows. First, we calculate the glass line $`\varphi _c(T)`$ for the mixture $`S_1`$ extrapolating diffusivity $`D`$ data according to $`D(\varphi \varphi _c(T))^{\gamma (T)}`$, for eight different $`T`$. Second, we calculate the MCT ideal glass transition line. The partial structure factors, the only input needed, are calculated with the Percus-Yevick approximation, numerically solving the Ornstein-Zernike equation. Third, we map the theoretical glass line on the simulation line, following the procedure first used by Sperl. Indeed MCT cannot reproduce the numerical values for $`\varphi _c`$. As shown in Fig.1, the linear transformation $`\varphi 1.897\varphi 0.3922`$ and $`T0.5882T0.225`$, allows to superimpose the MCT result with the simulation data, in the studied region of $`T`$ and $`\varphi `$. By calculating the MCT ideal glass transition line for several values of $`\mathrm{\Delta }`$, the theoretical location of the $`A_4`$ point for the binary mixture square-well model is obtained. Assuming that, to a first approximation, the same linear transformation holds also for the critical value of the well width, we find that the $`A_4`$ location maps to $`\mathrm{\Delta }_{ij}^{}=0.043\sigma _{ij}`$, $`\varphi ^{}=0.611`$ and $`\theta ^{}=0.416`$. This allows us to select a state point with parameters close to the critical ones and perform a simulation (in thermodynamic equilibrium) close to $`A_4`$. We choose $`\theta =0.4`$ and $`\varphi =0.6075`$ and compare, for the same $`T`$ and $`\varphi `$, the dynamics for $`\mathrm{\Delta }_{ij}^{}`$ (system $`S_2`$) and for the close-by value $`\mathrm{\Delta }_{ij}=0.031\sigma _{ij}`$ ($`S_1`$).

Fig.2 shows the MSD for the larger species $`\delta r_A^2`$ for both systems. We notice the presence of a subdiffusive regime at intermediate times, i.e. a variation according to $`ln(\delta r_A^2)=f+aln(t/\tau )`$ for about three decades for $`S_1`$, where $`a_1=0.44`$. The logarithmic regime extends by more than a decade if the attraction range is increased to that for system $`S_2`$. Simultaneously the exponent $`a`$ decreases to $`a_2=0.28`$. Fig. 3 shows the density-fluctuation auto-correlation functions $`\mathrm{\Phi }_q(t)`$ for a representative set of wave vectors $`q`$. The decay curves in Fig. 3 do not show the two steps scenario with a plateau characteristic of conventional glass forming liquids. Indeed it is impossible to fit these curves for large $`q`$ with the standard stretched exponential function. Instead, there is a region of clear logarithmic decay at $`q_1^{}\sigma _{_{BB}}=23.5`$ for $`S_1`$ and at $`q_2^{}\sigma _{_{BB}}=16.8`$ for $`S_2`$. The time-intervals over which logarithmic decay is observed are of similar size as those for the MSD. For $`q<q^{}`$, the $`\mathrm{\Phi }_q`$ vs $`lnt`$ curves are concave, and for $`q>q^{}`$ convex.

To show that the above described features are consistent with the ones predicted by MCT, we cite the general asymptotic decay law for the correlation function of a generic variable $`X`$, near the higher-order singularity

$$\mathrm{\Phi }_X(t)=f_Xh_X\left[B^{(1)}ln(t/\tau )+B_X^{(2)}ln^2(t/\tau )\right].$$
(1)
Here $`\tau `$ abbreviates a time-scale which diverges if the state approaches the singularity. The formula is obtained by asymptotic solution of the MCT equations, using the parameter differences $`\varphi ^{}\varphi `$, $`\theta ^{}\theta `$, $`\mathrm{\Delta }^{}\mathrm{\Delta }`$ as small quantities, say, of order $`ϵ`$. The coefficient $`B^{(1)}`$ is of order $`\sqrt{ϵ}`$ while $`B^{(2)}`$ is of order $`ϵ`$. The amplitude $`h_X`$ is independent of $`ϵ`$. The first term $`f_X`$ is the sum of the non-ergodicity parameter of variable $`X`$ at the singularity and a correction of order $`ϵ`$. Terms of order $`ϵ^{3/2}`$ are neglected. We can interpret $`ln\delta r_A^2`$ also according to Eq. (1). Hence, the straight line for $`S_1`$ with slope $`a_1`$ shown in Fig. 2 is consistent with the assumption that $`S_1`$ is close to the singularity. Changing from $`S_1`$ to $`S_2`$, the slope $`a=h_{_{MSD}}B^{(1)}`$ has to decrease and the range of validity of the leading order description has to expand. This is demonstrated impressively by the data.

To estimate the possibility of describing study the time dependence of $`\mathrm{\Phi }_q(t)`$ according to Eq. 1, we fit the density autocorrelation functions to a quadratic polynomial in $`log(t/\tau )`$ for different $`q`$ values. Fits are reported in Fig. 3. The fitting time window extends from about $`2.5`$ decades for the state point $`S_1`$ to four decades for $`S_2`$. The fitting parameter $`f_q`$, shown in Fig.4, provides an estimate of the non-ergodicity parameter at the $`A_4`$ point. We find that $`f_q`$ does not depend on the state point. This confirms the preceding conclusion that the studied state points are very close to the singularity and that the order $`ϵ`$ correction in $`f_q`$ cannot be detected.

The fit parameters for the coefficient $`H_q^{(1)}h_qB^{(1)}`$ are reported in the middle panel of Fig. 4. The decrease of $`H_q^{(1)}`$ upon changing from $`S_1`$ to $`S_2`$ is in agreement with the prediction that $`B^{(1)}`$ tends to zero with $`ϵ`$ approaching zero. As shown in Eq.1, $`H_q^{(1)}`$ factorizes into a control-parameter independent factor $`h_q`$, depending on $`q`$, and a control-parameter dependent factor $`B^{(1)}`$, independent of $`q`$. This implies that the $`q`$-dependence of $`H_q^{(1)}`$ should be the same for $`S_1`$ and $`S_2`$. As shown in Fig.4, this property is verified by the data. The same property does not apply to $`H_q^{(2)}h_qB_q^{(2)}`$ because of the $`q`$-dependence of $`B_q^{(2)}`$. Moreover, the fitting results confirm that $`H_q^{(2)}`$ is smaller than $`H_q^{(1)}`$, as expected being the first of order $`ϵ`$ and the second of order $`\sqrt{ϵ}`$. The wave-vector value where $`B_q^{(2)}=0`$ allows to identify the characteristic length scale associated to the pure logarithmic decay (see Fig.3). Such length scales are much shorter than the typical first neighbor shell.

To explicitly test the mentioned factorization, Fig. 5 shows a set of correlators rescaled to $`\widehat{\mathrm{\Phi }}_q(t)=(\mathrm{\Phi }_q(t)f_q)/H_q^{(1)}`$. MCT predicts that in leading order $`ϵ`$ all $`\widehat{\mathrm{\Phi }}_q(t)`$ vs $`ln(t)`$ curves collapse on the logarithmic decay law $`ln(t/\tau )`$. This is indeed shown for a large time interval that expands upon approaching the singularity, as it is demonstrated by comparing results for $`S_1`$ and $`S_2`$.

Results reported in this Letter have an intrinsic value associated to the observation, in a particularly simple system — a square well potential —, of a particularly complex dynamics over more than four decades in time. They show that MCT, a theory essentially developed to address the problem of the excluded volume glass transition is able — without any modification — to handle the logarithmic dynamics and to provide an interpretative scheme in term of $`A_4`$ point . The numerical results and the comparison with the theoretical predictions do constitute, in fact, a stringent verification that the logarithmic dynamics in the density autocorrelation functions and the subdiffusive behavior in the MSD can be fully rationalized by MCT. Indeed, data in Figs.3 and 4 are in agreement with MCT predictions (Fig. 12 and 13 of Ref.). The simplicity of the model and the complexity of the dynamics suggest that the short-range attractive colloids have the potentiality to become a benchmark for the development of extended theories of the glass transition.

We acknowledge support from MIUR PRIN and FIRB and INFM PRA-GENFDT. We thank W. Götze for illuminating discussions, G. Foffi and M. Sperl for comments and S. Buldyrev for providing us the molecular dynamics code for square-well binary systems.