text:
# The infinite-range quantum random Heisenberg magnet.

## I Introduction

The physics of disordered magnets is a fascinating subject of condensed matter physics. Traditionally, two main ingredients are singled out as crucial to set the physical behavior of these systems: strong interaction and frustration. It is well known that the interplay between them leads to a rich variety of magnetically ordered phases, including conventional commensurate or incommensurate spin-density waves as well as the more exotic spin-glass state. The latter is characterized by an ordered magnetic state with permanent magnetic moments in the microscopic scale, but randomly oriented producing a vanishing net macroscopic magnetization. While the concept of spin itself is purely quantum, it is often maintained that quantum fluctuations are not important for spin-glass physics . However, recent experiments have put the role of quantum fluctuations  at the center of the stage.

One example of real systems with a spin-glass phase at low temperature is the compound Li1-xHoxYF4 which is a dipolar coupled random magnet  and has been recently the focus of beautiful experiments  that introduced quantum fluctuations by means of an external transverse magnetic field. Another notable example is the LiV2O4  compound in which the magnetic V atoms are placed at the vertices of a phyrochlore-like structure which produces strong frustration and a spin-glass state. The V $`d`$-electrons also form a very narrow conduction band that behaves as a strongly renormalized Fermi liquid, with parameters comparable to those of the so called heavy fermion compounds that usually involve only $`f`$electrons. The likely connection between this observation and the spin-glass state remains a challenging open question. Finally, we have the cuprate superconductors, with vast experimental evidence that a glassy phase exists at low temperatures within a narrow range of doping concentrations between the antiferromagnetic and superconducting phases .

Much effort has been dedicated to investigate the spin-glass physics, and many ingenious and insightful theoretical ideas have contributed to our current understanding . Nevertheless, despite this effort, many fundamental questions remain unsolved and the role of quantum fluctuations is only beginning to receive due attention. Among the many long-standing unresolved issues of spin-glass physics we can mention the intriguing behavior of the specific heat. In fact, experiments show that this quantity systematically has a maximum well above the spin-glass freezing temperature $`T_g`$ . It is usually claimed that quantum mechanical effects are not essential for the physical phenomena related with the spin-glass phase . However, we shall later argue that precisely quantum effects might be at the origin of this long-standing puzzle.

Theoretical progress in the description of the glassy phase is usually prevented by technical difficulties. In particular, most of the analytical and numerical approaches rely on the so called replica trick . Unfortunately, this clever technique becomes usually impractical within the glassy phase when replica symmetry breaking occurs. In a recent paper  we introduced the use of the method of numerical exact diagonalization of finite size clusters to investigate models of quantum random magnets at $`T=0`$. This approach has two main advantages: (i) Averages over disorder can be directly performed, avoiding the use of replicas. In fact, the same effort is required to tackle both the disordered and the spin-glass ordered phases. (ii) The dynamical response is directly calculated on the real frequency axis and, unlike other numerical techniques such as quantum Monte Carlo, no analytical continuation from the imaginary axis is necessary. Another advantage is that unlike usual exact diagonalization calculations that give just a few poles in the response functions, in our case we obtain smooth functions due to the average over the disorder. However, the main drawback of this numerical method is that systems with a rather small number of spins are tractable. Nevertheless, this technical obstacle can be overcome by a careful finite size analysis. In fact, we found that most of the relevant physical quantities exhibit a smooth behavior as a function of the system size, and reliable extrapolations to the thermodynamic limit can be obtained . Moreover, as we have complete knowledge of the system for every realization, including the ground state wave function, we can look at its structure to try to gain new insights. This has indeed turned out to be a useful idea, as we obtain an appealing physical picture of the low energy excitations of the spin-glass state, that would have not emerged from classical model numerical calculations.

In this work we consider the SU(2) Heisenberg model for a system of $`S=1/2`$ spins with random infinite-range exchange random interactions. It is defined by the hamiltonian

$$H=\frac{1}{\sqrt{N}}\underset{i,j=1}{\overset{N}{}}J_{ij}𝐒_i𝐒_j,$$
(1)
where $`i,j`$ labels sites of a lattice with $`N`$ spins, $`𝐒`$ denote the SU(2) spin $`1/2`$ operators and the infinite-range exchange constants $`J_{ij}`$ are normally distributed with variance $`J^2`$ that we set to unity. The phase diagram of this model was first outlined by Bray and Moore , who argued that a glassy phase should exist at low temperature for all values of $`S`$. These authors proposed an approach which is formulated on the imaginary time axis and uses the replica trick to obtain a set of self-consistent dynamical mean-field equations. The exact numerical solution of these equations was later obtained with quantum Monte Carlo techniques in the paramagnetic phase , and the solution was found to become unstable towards spin-glass order at a critical temperature $`T_g0.14J`$. The study of SU(M) extensions of this model has been treated in the limit of M$`\mathrm{}`$ , and spin-glass and spin-liquid phases where found at sufficiently low temperatures.

We have recently investigated the dynamical response of the hamiltonian (1) at zero temperature with exact diagonalization techniques . We were able to describe in detail the structure of the ground state and the nature of the elemental excitations within the glassy phase. We found compelling evidence that the dynamical spin response behaves as $`\chi _{loc}^{\prime \prime }(\omega )\mathrm{q}\delta (\omega )+\chi _{reg}(\omega )`$ in the thermodynamic limit. We estimated $`q0.06`$ for the value of the Edwards-Anderson order parameter.

The aim of this paper is to extend the methodology based in exact diagonalization to investigate the physical behavior of the disordered Heisenberg model (1) at finite temperature. Our main goal is to achieve a detailed understanding of the different kind of excitations which occur in the dynamical spin response at different temperatures, within and above the glassy phase. We shall show that the regular part $`\chi _{reg}(\omega )`$ is in fact dominated by spin excitations due to quantum disorder. This contribution to the response function can be qualitatively understood in the framework of an heuristic mean-field theory that we present. Moreover, its functional form is very similar to that of the quantum spin-liquid discussed in Refs.  for the SU(M) generalization of the Heisenberg model in the regime of large M and small S. We also investigate the behavior of the specific heat $`C_v`$ as a function of temperature. We find that this quantity displays a smooth maximum at a temperature $`T_M`$ well above the freezing temperature $`T_g`$, and we link this fact to the presence of strong quantum fluctuations. Interestingly, the unconventional behavior of the specific heat is in accordance with the observation of this effect in real materials that have a spin-glass state at low temperatures and remained unexplained so far.

The paper is organized as follows. In Section II we explain the numerical method and some technical details. Section III contains the study of the dynamical spin response. Results for the specific heat $`C_v`$ are shown and discussed in Section IV. The summary and conclusions are presented in the Section V.

## II Methodology.

The general strategy is to take samples from the random ensemble of systems of size $`N`$ and exactly diagonalize the ensuing hamiltonians (1). The different physical quantities are computed for each realization and then averaged over the number of samples. Finite size effects are analyzed and results are extrapolated to the thermodynamic limit. Typically, systems with up to $`N=17`$ spins are solved at $`T=0`$ and up to $`N=12`$ at finite $`T`$. Averages are performed over several thousands to hundreds of thousands of disorder realizations. A typical run demands up to a week for the larger systems on an 8 node parallel cluster.

The ground state and the dynamical correlation functions at $`T=0`$ are calculated by the Lanczos method . It is convenient to take advantage of the SU(2) symmetry of the model. The selected basis of states belongs to the $`S^z=0`$ representation, where $`S^z`$ is the quantum number corresponding to the $`z`$-component of the total spin. A projector is used to find the ground state within each subspace with total spin $`S`$ . The local spin susceptibility is obtained from

$$\chi _{loc}(\omega )=\frac{1}{M}\underset{m=1}{\overset{M}{}}\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{1}{3}\underset{\alpha }{}\mathrm{\Phi }_0^{(m)}|S_i^\alpha \frac{1}{\omega H^{(m)}}S_i^\alpha |\mathrm{\Phi }_0^{(m)},$$
(2)
where $`M`$ is the number of realizations of disorder, $`|\mathrm{\Phi }_0^{(m)}`$ denotes the ground state for the $`J_{ij}`$ set corresponding to the $`m`$-th realization, and $`\alpha =x,y,z`$ labels the three components of the spin operator. Although we deal with systems having a finite number of poles for each realization, the average over disorder naturally produces smooth response functions without need of introducing an artificial broadening as in usual exact diagonalization methods. In some cases, we found useful to use a logarithmic discretization of the $`\omega `$-axis to obtain accurate results due to the large number of poles occurring at low frequencies.

To study the physical behavior at finite $`T`$, we exactly diagonalized the full hamiltonian matrix $`H^{(m)}`$. Each subspace belonging to the different $`S^z`$ representations is separately diagonalized in order to optimize the use of memory. In this case, the local spin susceptibility is obtained from the spectral function

$`\chi _{loc}^{\prime \prime }(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{1}{Z^{(m)}}}{\displaystyle \underset{k,l}{}}|\mathrm{\Phi }_k^{(m)}|S_i^z|\mathrm{\Phi }_l^{(m)}|^2[\mathrm{exp}(\beta E_k^{(m)})\delta (\omega (E_k^{(m)}E_l^{(m)}))`$ (3)
$``$ $`\mathrm{exp}(\beta E_l^{(m)})\delta (\omega (E_l^{(m)}E_k^{(m)}))],`$ (4)
where $`\chi _{loc}^{\prime \prime }(\omega )=2/\pi Im[\chi _{loc}(\omega )]`$, while $`E_k^{(m)}`$ is the eigenenergy of the eigenstate $`|\mathrm{\Phi }_k^{(m)}`$ corresponding to the $`m`$-th realization of disorder, $`\beta =1/T`$ is the inverse of the temperature and $`Z^{(m)}`$ is the ensuing partition function.

The spin-glass phase is signaled by the divergence of the spin-glass susceptibility $`\chi _{SG}`$, which is related to the local-spin susceptibility by

$$\chi _{SG}=\frac{\chi _{loc}^2}{1J^2\chi _{loc}^2},$$
(5)
where $`\chi _{loc}=Re[\chi _{loc}(\omega =0)]`$. Thus, the condition

$$J\chi _{loc}=1$$
(6)
indicates the instability of the system towards a spin glass state.

In previous papers  it was demonstrated the accuracy of the method by reproducing several known results for the infinite-range Ising model with random exchange interactions and transverse magnetic field $`\mathrm{\Gamma }`$ . In particular, an accurate estimate for the critical value of the transverse field $`\mathrm{\Gamma }_c`$, at which a quantum transition between the spin-glass and the paramagnetic phases takes place, was obtained. At finite $`T`$ we should also test the accuracy of this approach for the current model we are studying. There are not many well known results that can be used for benchmark of this numerical method. To the best of our knowledge, the only quantitative prediction is $`T_g0.14J`$ for the critical temperature at which the spin-glass to paramagnetic phase transition occurs. This estimate was obtained by quantum Monte Carlo numerical solution  of the mean-field equations derived by Bray and Moore .

Results for the behavior of $`\chi _{loc}`$ as a function of $`T`$ are shown in Fig. 1. The curve with circles correspond to extrapolations to the thermodynamic limit obtained from data for systems with up to $`N=12`$ spins. At high temperatures, the local susceptibility obeys a Curie law, as expected. The asymptotic behavior $`\chi _{loc}=\beta /4`$ is plotted in dashed lines for comparison. As the temperature decreases, quantum fluctuations cause the reduction of the effective local magnetic moment and the response at $`\omega =0`$ becomes smaller than that of classical spins. The condition (6) is satisfied at the critical temperature $`T_g0.13`$ where the systems begin to freeze. The diamond over the curve indicates the corresponding result obtained with quantum Monte Carlo  that shows very good agreement and validates our approach at finite $`T`$.

In Fig. 2 we show the behavior of $`\chi _{loc}`$ as a function of $`1/N`$ for a selected set of temperatures. A linear fit was performed to obtain the extrapolated values at $`1/N=0`$. Error bars in Fig. 1 indicate the corresponding standard deviation. Interestingly, an even-odd effect is observed at temperatures below $`T_g`$, which suggests a change in the scaling law as the system enters the ordered phase.

## III The dynamical spin response.

### A $`\chi _{loc}^{\prime \prime }(\omega )`$ at $`T=0`$.

Let us begin with the study of the spin response at $`T=0`$. In this subsection we shall summarize and extend our recent results of Ref. .

The spectral function $`\chi _{loc}^{\prime \prime }(\omega )`$ is plotted in Fig. 3 for systems of different sizes. As we discussed in our previous paper , the most prominent feature is the existence of a $`\delta (\omega )`$ piece plus a low frequency hump mounted on a regular contribution. Here, we shall present a more detailed analysis of the various contributions to $`\chi _{loc}^{\prime \prime }(\omega )`$, including some new excitations that were not previously noticed. We shall present just an heuristic, by no means rigorous, mathematical description of the different contributions to $`\chi _{loc}^{\prime \prime }`$ that will provide us with an appealing physical picture.

To gain insight on the nature of the dynamical response we analyzed individual realizations of disorder. We distinguished four different contributions which we shall describe in detail throughout this section:

$$\chi _{loc}^{\prime \prime }(\omega )=K\delta (\omega )+\chi _{low}^{\prime \prime }(\omega )+\chi _{high}^{\prime \prime }(\omega )+\chi _{inc}^{\prime \prime }(\omega ),$$
(7)
where $`\chi _{low}^{\prime \prime }(\omega )`$ corresponds to the low frequency hump, $`\chi _{high}^{\prime \prime }(\omega )`$ indicates a feature at high frequency $`𝒪(J)`$ carrying a small spectral weight, while $`\chi _{inc}^{\prime \prime }(\omega )`$ denotes the contribution of excitations of incoherent nature which conform a quantum disordered paramagnetic background.

We will discuss the physical nature of all these contributions next. We begin with the latter that is due to incoherent excitations. We found that it can be described with the following simple expression

$$\chi _{inc}^{\prime \prime }(\omega )=C\mathrm{exp}[\frac{\omega ^2}{2J^2S(S+1)}],$$
(8)
with $`S=1/2`$ and $`C`$ a constant. This function is plotted in thick line in Fig. 3 along with our results for $`\chi _{loc}^{\prime \prime }(\omega )`$.

We now show that this contribution can be heuristically described as independent spins in the presence of effective randomly fluctuating magnetic fields $`𝐡`$. We stress that we were not able to derive any proper demonstration of the validity of this procedure and we include it here with the sole hope of motivating the reader and perhaps providing new insights. The fields are the molecular fields at each site due to the action of all other spins,

$$𝐡=\frac{1}{\sqrt{N}}\underset{i,j}{}J_{ij}𝐒_j.$$
(9)
Following standard arguments , we assume that for a given realization of disorder, they point to any spatial direction with equal normal probability

$$P(h)\mathrm{exp}[\frac{h^2}{2\sigma ^2}],$$
(10)
where the variance $`\sigma `$ can be taken to be $`\sigma ^2=𝐡𝐡J^2S(S+1)`$.

We now recall that a spin under a magnetic field has a dynamical response given by

$$\chi _o^{\prime \prime }(\omega ,h)=\frac{1}{3}[\frac{1}{4}\delta (\omega )+\frac{1}{2}\delta (\omega h)],$$
(11)
where the finite frequency transverse part depends only on the magnitude of the field $`h`$.

Finally, we may combine (10) and (11) to obtain the incoherent background part of the dynamical spin response function within this heuristic picture

$$\chi _{inc}^{\prime \prime }(\omega )_{infty}^+\mathrm{}𝑑hP(h)\chi _o^{\prime \prime }(\omega ,h),$$
(12)
that leads to (8) where the constant $`C`$ = 0.175 is set to fit the numerical data. We carefully examined systems of several sizes and found that $`C`$ varies negligible with $`N`$. A point we would like to make is that the averaging (12) involves a one-dimensional integral. The reason is that we are describing the quantum disordered regular contribution, therefore these fields cannot be assumed to have a fixed direction in space.

Another interesting observation is that this contribution to the response function, is very similar to one found in the spin-liquid phase of the SU(M) generalization of the model at large M and in the quantum disordered regime (small $`S`$) .

In order to study the remaining three pieces of $`\chi _{loc}^{\prime \prime }(\omega )`$ in more detail, we substract $`\chi _{inc}^{\prime \prime }(\omega )`$ from this quantity and define

$$\delta \chi _{loc}^{\prime \prime }(\omega )=\chi _{loc}^{\prime \prime }(\omega )\chi _{inc}^{\prime \prime }(\omega ).$$
(13)
This quantity exhibits a strong dependence on the system size $`N`$ and a careful finite size analysis is necessary in order to extract reliable conclusions.

The $`K\delta (\omega )`$ part is a direct consequence of the SU(2) rotational invariance of the Hamiltonian. To demonstrate this we sorted out the contributions to $`\delta \chi _{loc}^{\prime \prime }(\omega )`$ coming from the different $`S`$-sectors. In Fig. 4 we show results for the $`S=0`$, $`S=1`$ and $`S=2`$ sectors. A remarkable observation is that the $`\delta (\omega )`$ part is absent in the $`S=0`$ case, while present in the remaining $`S0`$ sectors. Moreover, the systematic analysis of histograms of $`S`$ (cf. Fig. 5a) for increasing $`N`$, reveals that $`S\sqrt{N}`$ as shown in Fig. 5b. Therefore, although the mean total spin $`S0`$ in the thermodynamic limit (in fact diverges), the local magnetization, i.e. the magnetization per site $`S/N0`$. We have thus traced the origin of this $`\delta (\omega )`$ feature to a “soft mode” connecting the $`2S+1`$-fold degenerate ground state that persists in the $`N\mathrm{}`$ limit.

The low frequency hump $`\chi _{low}^{\prime \prime }(\omega )`$ evolves towards a narrowing feature at $`\omega =0`$ that becomes a $`\delta `$-function in the thermodynamic limit. In fact, the spectral weight of this part contributes together with $`K`$ to the Edwards-Anderson order parameter $`q`$. To elucidate the origin of this hump, we looked carefully to the structure of the ground state wave function for realizations of disorder that yielded spectral weight in the energy range of the hump $`\chi _{low}^{\prime \prime }(\omega )`$. We found that in those realizations the ground state wave function is dominated by some large amplitudes of just a few pairs of states. The states conforming a pair are related by a time-reversal operation . They appear in a symmetric (when $`N/2`$ is even) or antisymmetric (when $`N/2`$ is odd) linear combination. Moreover, we also observed that most of the individual spin states in each of the two states of a given pair are in a configuration compatible with the particular realization of $`J_{ij}`$ interactions. Thus we can think of this large fraction of spins as defining an unfrustrated cluster of size $`N_c`$ ($`𝒪(N)`$). Therefore, several unfrustrated clusters are defined for each realization of disorder. In fact, a number equal to the number of dominating pair of amplitudes occurring in the ground state wave function. An appealing picture to visualize the ground state is that the clusters undergo a quantum mechanical tunneling back and forth between the two corresponding pair of states. As many pairs coexist in any given ground state, then the picture is that many unfrustrated cluster simultaneously undergo these quantum oscillations.

With this physical picture of the ground state in mind it should not be hard to anticipate that the excitations that contribute to $`\chi _{low}^{\prime \prime }(\omega )`$ (the hump) are nothing but wave functions where one (any) of the clusters has its pair of associated amplitudes appearing with the opposite symmetry respect to that in the ground state (a “flipped” cluster). Since many coexistent clusters occur in a single ground state for a given realization, then there are many corresponding excited states contributing to the hump. We can think of these as real collective excitations, as they involve the simultaneous change in the state of a large number of individual spins. We have confirmed this qualitative picture by direct inspection of the amplitudes of the wave functions of the excited states in many realizations of disorder.

An important observation that we made is that the difference in energy between the ground state and the excitations is small (we found it to be $`𝒪(1/N)`$). Thus, these collective excitations become degenerate with the ground state and contribute with a sizable weight to the dynamical spin response at $`\omega =0`$ in the thermodynamic limit. In other words, the hump collapses into a $`\delta `$-function contribution in that limit. A side consequence of this observation is that even within the S=0 sector, which also shows a hump in the dynamical response, the T=0 state is that of a spin-glass with local magnetizations with long range order in time.

To perform a quantitative analysis of the behavior of $`\chi _{low}^{\prime \prime }(\omega )`$ as a function of $`N`$, we used the following parametrized form

$$\chi _{low}^{\prime \prime }(\omega )=A\mathrm{exp}[\frac{\omega ^2}{\mathrm{\Gamma }^2}].$$
(14)
and tracked the systematic evolution of the parameters as function of temperature and system size.

Typical curves are shown in Fig. 6 for two different sizes. The fitting parameters $`A`$ and $`\mathrm{\Gamma }`$ have a strong dependence on the system size. Their behavior as functions of $`1/N`$ is shown in Fig. 7. The width $`\mathrm{\Gamma }`$ obeys a linear extrapolation to $`0`$, while the integral of $`\chi _{low}^{\prime \prime }(\omega )`$ (squares in Fig. 7) remains approximately constant. This provides further quantitative support to the claim that $`\chi _{low}^{\prime \prime }(\omega )`$ evolves to a $`\delta (\omega )`$ contribution in the thermodynamic limit. The integrated spectral weight provides an estimate for the Edwards-Anderson order parameter $`q`$. With this procedure it is found $`q0.04`$, which improves upon our previous estimate $`q0.06`$ . We shall later discuss the origin of this difference.

To complete this subsection, we turn to discuss a last and small contribution to $`\chi _{loc}^{\prime \prime }(\omega )`$, the high frequency part $`\chi _{high}^{\prime \prime }(\omega )`$. This piece of the spectral function is related with the existence of a spin-glass order and it is essentially the high energy counterpart of $`\chi _{low}^{\prime \prime }(\omega )`$. Similar to the former case, $`\chi _{high}^{\prime \prime }(\omega )`$ is due to excitations of the clusters. But the key difference is that now the excitations are incoherent in nature, involving the unbinding of single spins out of the unfrustrated clusters. The spins will then independently revolve around its local field with a fast precession frequency. In other words, the relevant excitations contributing to $`\chi _{high}^{\prime \prime }(\omega )`$ correspond to the flipping of individual spins out of the large unfrustrated clusters.

An heuristic mean field picture can also be presented in this case. To describe these processes let us think of an originally unfrustrated cluster of $`N_c`$ spins as being made of a single spin plus the remaining cluster of $`N_c1`$ correlated spins. The latter would produce an effective magnetic field pointing to some fixed (or very slowly moving) direction which couples to the single spin that will undergo a fast precession around the slow effective field due to the cluster. To estimate the magnitude of the effective field $`h=|𝐡|`$, we assume it to be normally distributed around $`h_0=|_{ij}J_{ij}𝐒_j|/\sqrt{N}`$, which is a quantity of $`𝒪(J)`$. Therefore, similarly as we did before for $`\chi _{inc}^{\prime \prime }(\omega )`$, we can perform the average over the disorder as an average over the distribution of effective fields, to obtain the contribution of $`\chi _{high}^{\prime \prime }(\omega )`$ from

$$\chi _{high}^{\prime \prime }(\omega )_0^+\mathrm{}𝑑hh^2P(h)\chi _o^{\prime \prime }(\omega ,h),$$
(15)
with $`\chi _o^{\prime \prime }(\omega ,h)`$ of the form (11) where we point out that in the present case, the directions of the effective fields are now frozen and thus we include the $`h^2`$ term of the Jacobian of $`d`$h.

The integral leads to the following functional form for this high frequency incoherent contribution,

$$\chi _{high}^{\prime \prime }(\omega )=B\omega ^2\mathrm{exp}[\frac{(\omega h_0)^2}{2\sigma ^2}].$$
(16)
We were not able to determine $`h_0`$ and the variance $`\sigma `$ of the probability distribution analytically. These parameters and the constant $`B`$ can be chosen to fit the numerical data, however, the sizable finite size effects observed and the small spectral intensity did not allow us for a very accurate determination. We found that the choice of $`h_0=1`$ and $`\sigma =0.5`$ is consistent with the behavior of $`\chi _{high}^{\prime \prime }(\omega )`$ as the system grows to the large size limit, as shown in Fig. 8. We note that the spectral weight of this contribution is extremely small, less that 5 percent of the total weight of $`\chi ^{\prime \prime }(\omega )`$, which implies that the agreement remains only qualitative.

A final issue that we would like to discuss briefly concerns the difference in the estimated value for the Edwards Anderson order parameter $`q`$, in comparison with the result we reported in Ref. . This has a simple explanation: the $`T=0`$ sum rule $`𝑑\omega \chi ^{\prime \prime }(\omega )=1/4`$ is used to numerically compute the order parameter from the expression

$$q=\frac{1}{4}_0^{\mathrm{}}𝑑\omega \chi _{reg}^{\prime \prime }(\omega ),$$
(17)
where $`\chi _{reg}^{\prime \prime }(\omega )`$ denotes the regular part of the spectral function. We have already discussed the many different contributions to the $`\chi _{reg}^{\prime \prime }(\omega )`$, including the small last one from $`\chi _{high}^{\prime \prime }(\omega )`$, which went unnoticed in our previous work . The spectral weight carried by $`\chi _{high}^{\prime \prime }(\omega )`$ is only $`0.015`$ and is in fact the origin and amount of the difference. We have now achieved a better and more accurate description of the regular part of the frequency dependent spin-spin response function that through (17) enabled us to obtain a more accurate estimate for the order parameter $`q`$.

### B $`\chi _{loc}^{\prime \prime }(\omega )`$ at finite $`T`$.

We now turn to the behavior of the dynamic response at finite temperature. We recall the definition of the spectral density $`\rho (\omega )`$,

$$\rho (\omega )=\frac{1}{e^{\beta \omega }1}\chi _{loc}^{\prime \prime }(\omega ).$$
(18)
It is more convenient to use $`\rho (\omega )`$ instead of $`\chi _{loc}^{\prime \prime }(\omega )`$. While both quantities coincide at $`T=0`$, the former obeys a simple sum-rule which allows for a clearer interpretation of the transfer of spectral weight. In fact, the sum-rule reads

$$_0^{\mathrm{}}\rho (\omega )𝑑\omega =\frac{1}{4}.$$
(19)
The heuristic derivation of a response function for the incoherent degrees of freedom that we presented before can be easily extended to finite temperatures. We obtain

$$\chi _{inc}^{\prime \prime }(\omega ,T)=C\mathrm{exp}[\frac{\omega ^2}{2JS(S+1)}]tanh(\frac{\beta \omega }{2}).$$
(20)
and correspondingly,

$$\rho _{inc}(\omega ,T)=C\mathrm{exp}[\frac{\omega ^2}{2JS(S+1)}]tanh(\frac{\beta \omega }{2})\frac{1}{e^{\beta \omega }1}.$$
(21)
where the constant $`C`$ is now $`C=C(T)`$. Rather remarkably, the analysis of our results demonstrate that the functional form for $`\rho _{inc}(\omega ,T)`$ remains an excellent fit for the incoherent part of the spectra. Moreover, we found that within the temperature range investigated ($`0`$ to $`2T_g`$) the numerical value of the coefficient $`C`$ remains a constant within our numerical precision. We should mention, however, that at high enough temperature, at least the numerical value of $`C`$ has to change for the fit to comply with the sum-rule. The rather uninteresting high temperature regime has been left out from the scope of the present work.

Similarly as we did before, we use the expressions above to guide our analysis of the numerical data. We define

$$\delta \rho (\omega ,T)=\rho (\omega ,T)\rho _{inc}(\omega ,T).$$
(22)
Typical curves for $`T=0.25`$ are depicted in Fig. 9 for different system sizes. As in the case of $`T=0`$ the response is split in a low frequency feature and a higher frequency one. The low frequency peak corresponds to the $`\delta `$-function contribution discussed at $`T=0`$ that now acquires a finite width ($`𝒪(T)`$) due to the slow diffusive modes that appear at finite temperatures.

A last and smaller contribution is also apparent at higher frequencies of $`𝒪(J)`$. It is nothing but the signature of the high frequency contribution that we discussed at the end of the previous subsection. It corresponds to the fast precession of single unbinded spins around the now slowly varying magnetic fields. While the physical picture is qualitatively clear, we were not able to perform reliable extrapolations to the thermodynamic limit to extract quantitative detailed results upon the evolution of the spectrum as $`T`$ increases. In particular, it was not possible to obtain an accurate estimate of the Edwards-Anderson order parameter as a function of $`T`$. While no technical impediment exists a priori with the method, it turns out that the statistics that we collected were not sufficient. These details will have to wait for longer runs or more powerful computational resources. A final point we would like to consider is the transfer of spectral weight. Although the size of our systems is limited to 12 sites, this seems to be sufficient for a qualitative discussion the evolution and transfers of spectral weight of $`\rho (\omega )`$ as a function of temperature.

In Fig. 10 we show results for $`\delta \rho (\omega ,T)`$ at temperatures in the range of 0.04 to 0.25. One may think of this contribution to the spectra as due to the slow degrees of freedom that were originally frozen at $`T=0`$ and then gradually melt as $`T`$ is increased. Our results indicate that the width of the low frequency peak is linear in $`T`$, as may be expected for diffusive modes . Interestingly, the small higher frequency part is a rather broad feature that does not show any significant shift nor broadening with temperature, but merely looses spectral intensity with increasing $`T`$. At finite $`T`$, the effect of the gradual melting of the clusters of frozen spins leads to a slow (difussive) motion of the directions in which they point. However, as long as the drift of those clusters is slow compared to the fast precession of the unbinded spins, the higher frequency contribution should remain approximately unchanged, except for its spectral intensity which decreases with $`T`$.

This brings us to a last issue worth mentioning, namely, the transfer of spectral weight towards the background when $`T>T_g`$. This feature is better seen in the integrals

$$I(\omega ,T)=_{\mathrm{}}^\omega 𝑑\omega ^{}\rho (\omega ^{},T),$$
(23)
which are shown in Fig. 11 for the same set of temperatures as in Fig. 10.

It is clear in Fig. 11 that as the temperature increases, the distinction between low and high degrees of freedom gradually disappears and that the width of the main feature grows towards the single bare energy scale of the model $`J`$ that must control the physics at high temperatures.

## IV The specific heat.

In real materials with a spin-glass phase below $`T_g`$, the specific heat $`C_v`$ exhibits a smooth maximum well above the critical temperature . The origin of this behavior remains an open question as it cannot be explained from the predictions of classical models. In particular, for the Sherrington Kirpatrick model , which is the classical version of the Heisenberg model, it is found that $`C_v`$ has a maximum exactly at $`T_g^{SK}`$ with a small discontinuity in $`C_v/T`$ . In this section we shall argue that the observed behavior of $`C_v`$ can be understood as a result of quantum effects. We shall discuss the role of quantum fluctuations in the behavior of the specific heat

$$C_v=\frac{E}{T}.$$
(24)
and show that it leads to a shift of the maximum of $`C_v`$ towards temperatures higher than $`T_g`$.

The mean value of the energy per spin $`E`$ is computed from

$`E`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{M}}{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{Z^{(m)}}}{\displaystyle \underset{k}{}}[e^{\beta E_k^{(m)}}E_k^{(m)}]`$ (25)
$`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega \omega \chi _{loc}^{\prime \prime }(\omega ),`$ (26)
and (24) is obtained by numerically differentiating this quantity. The energy $`E`$ as a function of $`T`$ is shown in Fig. 12 for several system sizes. In the spirit of the previous analysis of the different contributions to $`\chi _{loc}^{\prime \prime }(\omega )`$ and using expression (26), we may think of $`E`$ as resulting from two different contributions: $`E=E_{inc}+E_{coh}`$. The energy $`E_{inc}`$ comes from the fraction of degrees of freedom that remain unfrozen down to $`T=0`$ due to quantum fluctuations, and $`E_{coh}`$ is the gain in energy due to the formation of the glass. As it turns out, the former is the largest contribution of the two and is obtain by replacing (20) in (26)

$$E_{inc}=C_{\mathrm{}}^+\mathrm{}𝑑hh\mathrm{exp}[\frac{h^2}{2J^2S(S+1)}]tanh(\frac{\beta h}{2}).$$
(27)
This energy is indicated in thick line in Fig. 12. At low temperatures it underestimates (in absolute value) the total $`E`$ as it lacks the sizeable contribution from the frozen degrees of freedom. On the other hand, at higher temperatures its qualitative behavior is similar to that of $`E(T)`$.

The results for $`C_v(T)=E/T`$ are shown in Fig. 13 for systems of various sizes. The contribution of the incoherent background is also plotted in thick line for comparison. It is remarkable that $`C_v`$ exhibits an smooth cusp with a maximum at a temperature $`T_M0.25`$ which is significantly larger than $`T_g`$. Interestingly, the contribution of the incoherent background exhibits a similar maximum at approximately $`T_M`$. The actual contribution to $`C_v`$ coming from the frozen part is relatively small compared to the quantum disordered part. Thus, the latter is found to dominate the behavior of $`C_v`$ and in fact it predominantly determines the actual position of the maximum.

Therefore, our results strongly suggest that in the present model incoherent mean-field like excitations due to quantum fluctuations set the behavior of $`C_v`$ as a function of $`T`$. In real systems the quantum fluctuations (controlled by, for instance, anisotropy) may play a more or less important role than in the current model but should always be present. One may thus expect that its effect could produce larger or smaller shifts of $`T_M`$ always to higher temperatures respect to $`T_g`$, as is in fact observed experimentally.

To give further support to our claims, we investigate the finite size effects on $`C_v`$ in more detail. We applied the same numerical procedure to the classical Sherrington Kirpatrick model , which is defined by the hamiltonian

$$H_{SK}=\frac{1}{\sqrt{N}}\underset{i,j=1}{\overset{N}{}}J_{ij}S_i^zS_j^z.$$
(28)
where $`S_i^z`$ denote Ising spins now. This model has been extensively studied and many well established results have been reported in the paramagnetic and glassy phases . Results for the behavior of $`C_v`$ as a function of $`T`$ obtained with exact diagonalization are shown in Fig. 14 for systems of various sizes. Although these data does not allow us to capture a subtle issue such as the small cusp in $`C_v`$ at $`T_g^{SK}`$ in the thermodynamic limit, it is clear that the curves become sharper as the size of the system increases. More remarkable is the fact that for these finite systems the maximum $`T_M`$ occurs below $`T_g`$, in marked contrast with the behavior of $`C_v`$ in the random Heisenberg model that is also depicted in the figure for comparison. Note that the finite size effects on $`T_M`$ in the quantum model are negligible.

We therefore conclude that our results show strong evidence that the smooth maximum of $`C_v`$ at $`T_M>T_g`$ is a genuine feature of the quantum model rather than a finite size effect.

## V Summary and discussion

In summary, we have studied the random infinite-range Heisenberg model at $`T=0`$ and at finite temperature with exact diagonalization techniques by carrying out a careful analysis of the finite size effects. We have first shown that our approach is able to reproduce the sole well established quantitative result known to us in this model, namely, the quantum Monte Carlo estimate of the critical temperature $`T_g`$  for spin glass ordering.

At $`T=0`$ we found that the dynamical spin response in the thermodynamic limit behaves as $`\chi _{loc}^{\prime \prime }(\omega )=q\delta (\omega )+\chi _{reg}^{\prime \prime }(\omega )`$, where $`q`$ is the Edwards Anderson parameter. The low energy excitations producing the $`\delta (\omega )`$ response are due to two different contributions. On one hand, the freezing of slow collective excitations that we identified as groups of spins that build unfrustrated clusters and coexist in a given typical realization of disorder. The nature of these excitations is the slow quantum tunneling between (cluster) states related by a time-reversal operation. On the other hand, there is also a second contribution to the $`\delta (\omega )`$ response due to a 0-mode like excitation originated in the spin-rotation invariance of the model and the fact that there is a finite mean total spin. Interestingly, while the local magnetic moment of the system, i.e. the mean magnetic moment per site, does vanish in the thermodynamic limit, we found that the total spin is $`S0`$ and of $`𝒪(\sqrt{N})`$.

The regular part $`\chi _{reg}^{\prime \prime }(\omega )`$ collects the contribution of incoherent excitations due to quantum fluctuations plus the effect of high energy excitations due to the flipping of some few unbinded spins. Both kinds of excitations seem to be very well described in terms of heuristic mean-field theories that provide a qualitative insight.

When the temperature $`T`$ increases the response of the low energy excitations is found to acquire a finite width of order $`T`$ due to the melting of frozen spins. As these degrees of freedom gradually unbind and become part of the regular contribution, the weight of the high energy excitations gets also gradually reduced. However, the shape of this part of the response is preserved within a wide temperature range as long as there is a separation of energy scales between the slow and fast excitations.

We have also studied the behavior of the specific heat as a function of $`T`$. We found compelling evidence that in the thermodynamic limit $`C_v`$ has a smooth maximum at a temperature $`T_M`$, significantly larger than the transition temperature $`T_g`$. This behavior is in excellent qualitative agreement with that experimentally observed in real spin glasses and is in marked contrast with the predictions based on classical models such as the Sherrington-Kirpatrick hamiltonian. We have shown that the origin of the unusual behavior can be traced to the excitations due to quantum disorder which dominate the behavior of the specific heat.

In conclusion, based in our numerical analysis, the picture we propose for the spin-glass state is that a fraction of the system forms coexisting unfrustrated clusters which are frozen or undergoing a slow collective motion. The rest of the system remains in an uncorrelated incoherent state mostly described by local spin dynamics. Such a picture manifests itself in the different contributions to the dynamical response and particularly in $`C_v(T)`$, since the incoherent excitations produce a smooth variation in the specific heat and this quantity does not show any signature of the paramagnetic to spin-glass transition. We find that this picture can be quite general and that it provides a possible explanation for the behavior of the specific heat in real materials which remained, so far, a standing question in the physics of spin glass systems.

Another interesting observation is that the solutions of the SU(M) generalization of the model treated in the large M limit produces several qualitatively different phases depending on the value of the ”size” of the spin, $`S`$ . In particular, for large $`S`$ the system is in a spin-glass state with a dynamic response that has a $`\delta `$-function at $`\omega =0`$, while for $`S0`$ the system goes into a spin-liquid state with a response function that remains regular and goes to a constant at $`\omega 0`$. Therefore, under the light of our results, the large M limit of the model seems to capture, in its different regimes, several aspects present in the M=2 solution.

From the technical point of view, we can say that the present approach offers a very appealing alternative for the study of disordered systems. Without the need of approximations, it enables the treatment of paramagnetic and disordered phases on equal footing and provides a deep insight on the nature of the relevant excitations. At least for the present model, we were able to find reliable extrapolations to the thermodynamic limit for most of the studied quantities, thus overcoming the obstacle of the finite-size effects of the systems under study.

Finally, many interesting issues are left for future work, such as a systematic study of the role of dimensionality or lattice connectivity and that of anisotropy.

## VI Acknowledgements

L. A. thanks G. Santoro for his comments on this manuscript. We acknowledge support from CONICET and from grants of ANCyPT, Fundación Antorchas and ECOS-SECyT.


format: text

text:
# Point source detection and extraction from simulated Planck TOD using optimal adaptive filters

## 1 Introduction

One of the most thrilling challenges in the study of the Cosmic Microwave Background (CMB) is to deal with the problem of separating the cosmological signal from the different foregrounds and noises that appear in CMB experiments. This problem will be specially relevant in future high-resolution experiments such as MAP (Bennett et al. 1996) and Planck (Mandolesi et al. 1998, Puget et al. 1998). From the point of view of determining constraints on fundamental cosmological parameters, the different foregrounds (synchrotron emission, galactic dust, free-free radiation, thermal and kinetic Sunyaev-Zel’dovich (SZ) emission from galaxy clusters and extragalactic point sources) are considered as contaminants and therefore must be removed together with the noise in order to extract the cosmological signal. Besides, knowledge about each of these ’contaminants’ has a great scientific relevance on itself. Therefore, it is of exceptional importance to be provided with good techniques of denoising and foreground separation.

There are several methods already available to perform the separation of the different components in CMB observations, such as the ones based on Wiener filter (Tegmark and Efstathiou 1996, Bouchet and Gispert 1999) and on Maximum Entropy Methods (MEM, Hobson et al. 1998, 1999; Stolyarov et al. 2001). All these methods take advantage of the different statistical properties of the CMB and the foregrounds (e.g. different angular power spectra) as well as the knowledge of their distinctive frequency dependence. When neither the spatial distribution of the foregrounds nor its frequency dependence are well known the performance of these separation methods dramatically drops. Such is the case of point sources, whose spatial distribution and abundance remains uncertain and whose frequency dependence is not well known. Besides, they can show temporal and even spectral variability. Extragalactic point sources should be removed from the maps before any analysis should be performed. Maximum Entropy Methods can produce a catalogue of point sources as a residual (noise) from the separation process. However, only the faintest sources are recovered, and the brightest ones can still be observed in the residuals. The Independent Component Analysis technique (Baccigalupi et al 2000, Maino et al 2001) has also been applied to this problem with promising results: it has the advantage that, unlike the previous methods, it does not need any prior knowledge of the components to be separated, but its weakest point is actually the separation of point sources.

Filtering techniques have been successful in both denoising and extracting the brightest point sources from CMB maps. Tegmark and Oliveira-Costa 1998 presented a filtering scheme that minimizes the variance of the map. However, their filter does not produce a maximum gain (considered as the ratio between the $`\sigma `$-level of the source after filtering and the $`\sigma `$-level of the source before filtering) at the scale of the source. All the point sources in a CMB map will have the same profile and size (the convolution of a $`\delta `$-Dirac source with the antenna beam) and the filtering process should take advantage of this fact using a filter optimal for that particular profile and size (scale). Wavelet techniques are specially suitable to deal with scales and spatial location at the same time. The application of the Mexican Hat Wavelet (MHW) to realistic simulations was presented in Cayón et al. (2000) and extended in Vielva et al. (2000a). The principal advantage of this method is that no specific assumptions about the statistical properties of the point source population or the underlying emission from the CMB and other foregrounds is required. The MHW and MEM techniques have successfully been combined to extract both faint and bright point sources at the same time (Vielva et al. 2000b).

It is desirable to exploit to the maximum the mentioned scale distinctiveness of point sources. In Sanz et al. 2001 (S01, hereafter), a method was presented to, for a given source brightness profile and given the statistical properties of the background (CMB plus foregrounds and noise), derive the analytic shape of the optimal filter for the particular considered case. That kind of filter was called ’optimal pseudo-filter’. S01 defines as an optimal pseudo-filter one that a) is unbiased, b) gives a maximum at the position and scale of the source and c) gives the minimum variance of pseudo-filter coefficients (i.e. is an efficient estimator of the amplitude of the sources) under assumptions a) and b). They considered one, two and three-dimensional cases and showed that the MHW is optimal (in the mentioned sense) or nearly optimal for a two-dimensional Gaussian profile and reasonable foreground conditions. After carefully thinking about the terminology used in that previous work, we have decided to substitute the rather confusing ’pseudo-filter’ term by ’adaptive filter’. This is so because a filter that satisfies conditions a) to c) adapts itself to the characteristics of the signal and the noise. In the following we will use ’adaptive filter’ in the same sense as ’pseudo-filter’ in S01.

Although two-dimensional maps are the most useful and extended form to show CMB data, they are only available after an exhaustive process of analysis and reduction of the raw data from CMB experiments. These experiments scan different patches of the sky in a sequence producing a unidimensional set of time ordered data (TOD). TODs suffer from lower signal-to-noise ratios than final two-dimensional maps, but on the other hand are less likely to show artifacts coming from the data reduction, such as pixel-to-pixel noise correlations. Moreover, TODs can be analyzed in (almost) real time during the observations in order to produce early (preliminary) catalogs of sources.

In this paper we apply the techniques developed in S01 to realistic simulations of the TOD coming from one of the 30GHz Planck LFI’s channels. In section 2 we summarize some of the conclusions of S01 and present the semi-analytic adaptive filter that should be used in a realistic case. Section 3 describes the simulations used in this work. In section 4 we describe the analysis of the simulated data. In section 5 we describe the performance of the optimal adaptive filter, comparing it with other filtering schemes such as Gaussian filter and MHW. Finally, in section 6 we discuss our conclusions and give an outline of future work in this field.

## 2 One-dimensional adaptive filter

The general formalism of adaptive filters in an n-dimensional image was presented in S01. TODs can be considered as a particular case of a one-dimensional image (a spectrum is another interesting case). The image data values can be expressed as

$$y(t)=s(t)+n(t),$$
(1)
where $`t`$ is the time, $`s(t)`$ represents a symmetric source and $`n(t)`$ is a homogeneous and isotropic background with mean value $`<n(t)>=0`$ and characterized by the power spectrum $`P(q)`$ ($`q`$ is the absolute value of the ’wave vector’ associated to $`t`$). If $`A=s(0)`$ is the amplitude of the source we can introduce the profile $`\tau (t)`$ as $`s(t)=A\tau (t)`$. Let us introduce a centrally symmetric adaptive filter $`\mathrm{\Psi }(t;b,R)`$

$$\mathrm{\Psi }(t,b,R)=\frac{1}{R}\psi (\frac{|tb|}{R})$$
(2)
where $`b`$ defines a translation whereas $`R`$ defines a scaling. The adaptive filtered field is defined as:

$$w(R,b)=𝑑ty(t)\mathrm{\Psi }(t;b,R)𝑑qe^{iqt}y(q)\psi (Rq)$$
(3)
In the last equivalence we have expressed the convolution as a product in Fourier space.

Following S01, it is possible to derive the analytic form of the optimal adaptive filter for a given source profile, once we have defined what we mean by ’optimal’. Let us define an optimal adaptive filter as one that satisfies the following conditions:

1. There exists a scale $`R_o`$ such that at the point source position ($`b=0`$) $`<w(R,0)>`$ has a maximum at that scale.
2. $`<w(R_o,0)>=s(0)`$, i.e. $`w(R_o,0)`$ is an unbiased estimator of the amplitude of the source.
3. The variance of $`w(R,b)`$ has a minimum at the scale $`R_o`$, i.e. we have an efficient estimator.

Given these three conditions, the solution (adaptive filter) is found to be:

$$\stackrel{~}{\psi }(q)\psi (R_oq)=\frac{\tau (q)}{2P(q)\mathrm{\Delta }}[b+c(a+b)\frac{d\mathrm{ln}\tau }{d\mathrm{ln}q}]$$
(4)
$`a{\displaystyle 𝑑q\frac{\tau ^2}{P}},b{\displaystyle 𝑑q\frac{\tau }{P}\frac{d\tau }{d\mathrm{ln}q}},`$
$`c{\displaystyle 𝑑q\frac{1}{P}[\frac{d\tau }{d\mathrm{ln}q}]^2},\mathrm{\Delta }=acb^2`$ (5)
The limits of the integrals go from $`0`$ to $`\mathrm{}`$.

In S01 analytic expressions for Gaussian sources and backgrounds of the type $`P(q)q^\gamma `$ were derived. In a more realistic case, the background can not be modelled in such a simple way, and integrals in (2) should be numerically estimated. When dealing with real data (or realistic simulations such as those used in this work) we must perform the following steps: first, determine the power spectrum of the background directly from the image. Second, evaluate integrals (2). Third, build the adaptive filter (4) and make the convolution (3). Finally, we can proceed to detect the sources, for example looking for peaks above a certain $`\sigma `$-level in the coefficient (filtered) image.

When determining the power spectrum of an image we obtain the power spectrum of both the background and the sources together. In the following we consider that the contribution of the point sources to the total power spectrum is negligible. This is a reasonable assumption in a realistic case, specially at medium and high wavelengths where the emission of IR and radio sources is weak. Another problem related to power spectrum determination is the variance of the power spectrum estimator. The variance of one-dimensional power spectrum estimators tends to be larger than that of two-dimensional cases because the number of sample points is usually smaller. A possible solution to this problem is to estimate low resolution power spectra. However, a large amount of information is lost at scales that are crucial for the determination of the adaptive filters. For this work we chose to average the estimated power spectra of adjacent rings of the TOD (corresponding to regions of the sky separated by only a few arcminutes and therefore possessing very similar underlying power spectra).

A typical profile of the optimal adaptive filter (in Fourier space) for a section of the simulated TOD used in this paper is shown in figure 1. The profile (solid line) is irregular due to the roughness of the estimated power spectrum. These irregularities reflect the particularities of the data and define the scales where the adaptive filter is more or less efficient. For comparison, a Mexican Hat Wavelet (dashed line) and a Gaussian (dotted line), both of them with a width equal to the width of the source, are also shown.

## 3 Data

For this work we chose the 30 GHz LFI28 channel of Planck because of the relatively small size of its resulting TOD. Higher frequencies, in spite of being expected to show more contribution from point sources due to their higher resolution and the contribution of IR sources at $`\nu 300`$ GHz, lead to huge TOD sets and require many hours of computation to be simulated and analysed on an average computer and will be analyzed in a future work.

Our data simulates a 6 months run of the 30 GHz LFI28 channel of Planck, covering the whole sky except for two circles of $`1\mathrm{°}82`$ around the ecliptic poles. The sky simulation consist of the addition (in flux) of a CMB simulation and three templates, each for the Galactic synchrotron radiation, Galactic dust and thermal SZ. These templates have been generated using proper models that are freely distributed among the Planck collaboration222http://planck.MPA-Garching.MPG.DE/SimData/. This sky is then combined with a map of point sources, generated from a catalog performed following the model of Toffolatti et al. 1998 (that is not public for the Planck collaboration) and then ’observed’ by the Planck Pipeline Simulator for one of the LFI 30 GHz radiometers; resulting in a TOD for the observed sky (CMB, foregrounds and point sources) plus instrumental noise. The data are ordered in circular rings centered on points situated on the Ecliptic. The angle between the pointing axis and the rotation axis for the LFI28 instrument is $`88\mathrm{°}18`$. The maximum separation between two consecutive rings is $`2.5\mathrm{}`$ (at the intersection with the Ecliptic plane). Each ring results from the average of 60 revolutions of the detector around the rotation axis, corresponding to one hour of integration time, and contains 1950 measures of antenna temperature. There are a total of 4383 of such rings, leading to 8546850 measures of temperature. The antenna has a FWHM of $`33\mathrm{}`$ and its response slightly differs from a circular Gaussian one.

The simulation contains CMB emission, different extended foregrounds (Galactic synchrotron, dust and free-free, thermal and kinetic SZ emission from clusters, etc), point sources and instrumental noise. Both white noise and $`1/f`$ noise are present. The knee frequency $`f_\nu `$ is set to be $`<20`$ mHz (less than the frequency of rotation). In the lower panel of figure 2 we show a segment of one of the rings of the simulation. There is a bright point source near pixel 400 (in fact, it is the brightest source in the simulation). Apart from this extraordinarily bright source, the features that dominate the image are Galactic emission (the large peaks around pixels 150 and 940) and noise.

## 4 Data analysis

The complete set of simulated TOD was filtered using the optimal adaptive filter described in section 3. Each individual ring was filtered separately. The power spectrum that appears in equations (4) and (2) was obtained by averaging the estimated power spectra of twenty-one consecutive rings (the ring that is being filtered, the ten previous rings and the ten subsequent rings).

In order to detect the sources from the filtered image one can set a certain threshold over the dispersion of the filtered ring and look for connected regions (peaks) above that threshold. That would be the most direct detection method if the different rings were independent of their neighbours. This is not the case, since adjacent rings scan very close regions and each source is expected to appear in more than one ring. The use of information coming from neighbouring rings allows to increase the effective signal to noise ratio of the detections and to discard spurious ’sources’ due to noise fluctuations. The most straightforward way to detect sources is then to perform a kind of two-dimensional thresholding, looking for connected pixels at equal latitude that appear in several adjacent rings. In a first approximation, the position of the source will be the position of the maximum of that region of connected pixels. We will show that this approximation is good enough for the purpose of locating the sources with an error comparable to $`1/3`$ of the antenna FWHM.

Two different regimes of ’noise’ have to be removed in order to optimize the detection of the sources. The large scale features due to Galaxy foregrounds as well as CMB fluctuations are strongly correlated between a given ring of data and its neighbours. On the other hand, the small scale noise, dominated by white instrumental noise, is expected to be nearly independent from one ring to another. This suggests a further step in the idea of combining information from nearby rings to increase the signal to noise ratio of the sources. By averaging each ring data (before filtering) we can construct a ’synthetic TOD’ in which the large scale fluctuations are almost the same than in the original TOD but the noise at the scale of the pixel has been greatly diminished. The $`1/f`$ will still be present in the averaged TOD. The point sources of the original TOD are replaced in the averaged TOD by Gaussians of amplitude

$$A_a=\frac{1}{N_a}\underset{i=1}{\overset{N_a}{}}Ae^{\frac{\stackrel{}{x_o}\stackrel{}{x_i}^2}{2\sigma _b^2}}$$
(6)
where $`A`$ is the true amplitude of the source, $`N_a`$ is the number of rings that are averaged, $`\sigma _b`$ is the beam width and the distance $`\stackrel{}{x_o}\stackrel{}{x_i}^2`$ is the geodesic (spherical) angular distance between the pixel corresponding to the position of the source and the pixel that is being added to calculate the average. This is true when the beam is a perfect Gaussian because in that case the average is a weighted sum of Gaussians of equal width, that is, a new Gaussian of the same width and amplitude given by eq. (6). Therefore, we can filter the rather denoised, ’synthetic TOD’, instead of the original TOD and correct the amplitude of the sources that we will detect using eq. (6) in order to recover an estimate $`\widehat{A}`$ of the true amplitude $`A`$.

Although the Planck Pipeline Simulator uses a realistic beam, we assumed for this work that the beam is a perfect Gaussian. In section 5 it will be shown that this approximation works reasonably well. However, future work should take into account real beams; the 30 GHz detector is expected to be slightly elliptic. This asymmetry will introduce complicate effects in the map making and analysis. In the case of TOD the effect of ellipticity, when projected over the scan trajectory, is a change in the ‘effective width’ of the beam. Therefore, the scale-adaptive filter should be calculated for such ‘effective width’. More complicated asymmetries will require more careful treatment. Future work will deal with this issue.

## 5 Results

In order to determine the goodness of the averaging of close rings we performed two different filterings of the TOD. In the first case, the raw data of the TOD were filtered ring by ring as described in the section 4. In the second case, a ’synthetic TOD’ in which each individual ring was the average of the original ring (at that position) with the twelve neighbours was constructed. The number of rings averaged is such that at the maximum ring spatial separation region on the sky (the ecliptic Equator) the separation between the two most separated rings is approximately equal to the FWHM of the beam. After filtering, the detection and extraction was performed looking for sets of 5 or more connected pixels (in the 2-dimensional sense). A lesser number of connected pixels required would lead to a much higher number of spurious detections. To check if the detections correspond to real sources we compare with the catalogue of the 1000 brightest sources present in the simulations. This catalogue is complete and its flux limit is clearly lower than the expected flux limit of the detected sources.

In the lower panel of figure 3 the number of detections (defined as the number of peaks encountered above a certain threshold that correspond to a real source in the reference catalog) above several $`\sigma `$ thresholds is shown for the case where the data have been averaged (filled circles and solid line) and the case where they have not (open circles and dashed line). The remaining quantities of interest (such as estimation of the amplitudes, fluxes, etc.) concerning the detection/extraction of sources over the averaged and then filtered TOD are shown in table 1 and will be discussed later. As expected, the number of sources detected when analysing the synthetic TOD is 4 to 5 times higher to the number of sources detected over the original, filtered TOD for every $`\sigma `$ level.

Apart from having the greatest number of detections, it is also of great importance to reduce as much as possible the number of spurious detections. In the lower panel of figure 4 (labeled as ’simple detection’) the ratio between spurious and ’true’ detections is represented for the case where the data have been averaged (filled circles and solid line) and the case where they have not (open circles and dashed line). The ratio is lower for the case in which the data have not been averaged. It is evident that a kind of compromise has to be reached between gain and reliability. If we arbitrarily set a maximum proportion of spurious sources versus true ones, say a $`10\%`$ (represented in fig. 4 with an horizontal dashed line), we can determine the minimum $`\sigma `$ level that satisfies this condition. In this example, for the case of non-averaged rings, we can reach the $`3\sigma `$ and find 80 sources (to a minimum flux of 4.33 Jy) with a $`7.5\%`$ of spurious detections. For the case of averaged rings, we must go to the $`4\sigma `$ level, where we find 224 sources (to a minimum flux of 0.885 Jy) with a $`9.8\%`$ of spurious detections. We conclude that the averaging of neighbouring rings is a valid strategy to reduce pixel-scale noise. In the following, all the results will refer to filtering of averaged rings.

The number of detections and spurious sources found with the optimal adaptive filter applied to a synthetic TOD is shown in table 1. The determination of the position of the source, the mean relative absolute error in the determination of the amplitude (defined as $`m.r.a.e.=A_0A_e/A_0`$, where $`A_e`$ is the estimated amplitude and $`A_0`$ is the real amplitude), the mean bias (defined in the same way as the m.r.a.e, but without the absolute value), the minimum flux reached and the completeness flux are also included in table 1. In each case the mean error in the position of the sources is comparable with the size of the ’pixel’ of the TOD, $`11\mathrm{}`$. The determination of the amplitude using eq. (6) has relative errors ranging from $`24.32\%`$ at $`3\sigma `$ threshold to $`19.82\%`$ at $`5.0\sigma `$. The error decreases as the detection threshold increases. This indicates that the estimation of the amplitude of weak sources is less accurate than the estimation of the amplitude of bright ones. Also, the determination of the amplitude is biased to higher values at low $`\sigma `$ levels and to lower values (negative bias) at high $`\sigma `$ levels. The positive bias for weak sources arises due to the peak finding algorithm: it finds preferently the maxima in pixels where the noise contribution is positive.

Let us consider for a while the possible causes for the mentioned bias and how to deal with it. Scale-adaptive filter is designed to be an unbiased estimator of the amplitude of the sources. Thus, *the method is unbiased*. However, in practise it is found a small bias. The origin of this bias can be found in two different kind of effects. On the one hand, the bias introduced by the detection method: this bias was described in the last paragraph and is related with the well-known ‘detection bias’. On the other hand, bias may arise due to the non-ideality of the data. In first place, the profiles considered in the design of the filters are continuous whereas real data is pixelised. Therefore, the correlation between the (continuous) filter and the (pixelised) source profile is not perfect, this resulting in a non-ideal performance of the filter. In second place, the limits in the integrals (2) are from $`0`$ to $`\mathrm{}`$, whereas in real data the frequencies are limited by the sizes of the data set and the pixel. Therefore, integrals (2) can only be approximately calculated. This leads to an inaccuracy in the determination of the shape of the filter and, more important for the bias, the *normalization* (the $`\mathrm{\Delta }`$ in the mentioned integrals). Another source of bias appears when the assumption about the profile of the sources is wrong. For example, we have assumed that the detector beam is Gaussian, while in fact it is not; this could explain the bias in our results. In SO1 it was found using simple simulations that the bias due to the non-ideality of the data was negative. Similarly, the negative bias that appear here is interpreted in the way described above. In practise, it is difficult to distinguish from the data the different contributions to the bias among the ones above mentioned. A possibility to overcome the bias is to calibrate it using a great number of simulations with the same background and artificial sources of known amplitude. Once removed the ‘systematic’ bias, the remaining error will be statistical. This kind of study will be carried on in a future work.

In the case we are considering the bias is never greater than $`16\%`$. At intermediate thresholds, the (positive) effects of detection bias and the non-ideality of the data (negative) tend to cancel and, on average, the amplitude estimates are unbiased. To compare with other ’classical’ filters we repeated the process using a Gaussian filter and a Mexican Hat Wavelet (MHW, hereafter), both of them with a width equal to the beam width of the antenna ($`33\mathrm{}`$). The normalization of both filters was chosen so that the coefficient at the position of the source is equal to the amplitude of the source (that is, the filtering process does not change the amplitude of the sources). The adaptive filter automatically satisfies this condition (condition 2 for an optimal adaptive filter). While the MHW and the optimal adaptive filters are both band-pass filters, the Gaussian is a low-pass filter, so the comparison with the Gaussian is a bit unfair: the Gaussian is expected to perform significantly worse than the other two filters. The number of detections above several thresholds for the two ’classical’ filters together with the optimal adaptive filter are shown in table 1. The number of detections is similar for the MHW and the optimal adaptive filter, yet are slightly higher for the optimal adaptive filter. The lowest number of detections corresponds to the Gaussian filter. In the lower panel of figure 5 the number of detections with the three different filters is compared. Detections with the optimal adaptive filter are shown with open circles and solid line. The open boxes and dashed line corresponds to MHW detections and the triangles and dashed line corresponds to Gaussian filter detections. The ratio between spurious and true detections for the three filters is shown in the top panel of figure 5. Except for the $`2.5\sigma `$ level, the Gaussian filter produces the worst ratio. The optimal adaptive filter gives spurious to detected ($`e/d`$) ratios that quickly decline with increasing $`\sigma `$ thresholds. The $`e/d`$ ratio for the MHW remains almost constant with $`\sigma `$ in the considered cases and clearly exceeds the ratio obtained with the optimal adaptive filter. The m.r.a.e. and the bias in the determination of the amplitude are huge in the case of the Gaussian filter. Both have very similar values. That means that the main source of error is systematic (the filter is biased). Considering the m.r.a.e, the MHW seems to give amplitude estimates a few percent better than the optimal adaptive filter. Flux limits are similar in the MHW and optimal adaptive filter cases. The Gaussian filter leads to higher inferior flux and completeness limits.

Figure 2 provides a useful insight into what is happening with the different filters. The Gaussian filter smoothes the image, removing very efficiently the small scale noise but allowing the large structures (the Galaxy and others extended fluctuations) to remain in the image. The naive $`\sigma `$ thresholding counts these bright, large structures as sources, leading to a big relative number of spurious detections. Besides, the sources that by chance lie on ’valleys’ of the background can not be detected. On the other hand, sources that lie on areas of positive background are enhanced and can be more easily detected. This can explain the large and positive bias in the detections with the Gaussian filter.

Both MHW and optimal adaptive filter are better prepared to deal with this problem than the Gaussian window. Their profiles in Fourier space drop to zero at low frequencies and thus they are efficient in removing large scale structures. Images filtered with the MHW and the optimal adaptive filter in figure 2 are similar. A visual inspection reveals that MHW smoothes better the high-frequency fluctuations. The optimal adaptive filter is more efficient in removing medium and large structures. The fluctuations around pixel 150 (corresponding to one of the two observations of the Galaxy in the ring) are better removed with the optimal adaptive filter than with the MHW. This is more apparent in figure 1 where the profiles of the different filters in Fourier space clearly indicate the faster drop of the adaptive filter at low frequencies (large scales) and its slower drop at high frequencies (small scales). The MHW can have the same problem as the Gaussian filter. However, the probability of this problem should be smaller due to its better efficiency in removing the Galaxy and other large-scale structures. This effect explains the higher $`e/d`$ ratio and the trend to positive bias.

We conclude that the optimal adaptive filter detects point sources better than the MHW and the Gaussian window. The number of detections is comparable to the number of detections with the MHW and clearly higher than with the Gaussian window. The relative number of spurious detections with the optimal adaptive filter is lower, except for very low detection thresholds, than with the MHW and the Gaussian window. Over $`4\sigma `$ the contamination of spurious detections is lower than $`10\%`$. At this level ($`4\sigma `$) the number of expected sources in all the sky is of a few hundreds (224 in our simulation) above fluxes of around 0.9 Jy.

Some tests can be performed in order to discover if the number of spurious sources, the most unpleasant effect of the filtering and detection process, can be reduced. Where do these spurious detections come from? One possibility is that the peak finding algorithm is detecting the Galaxy or other large-scale features. Such structure will appear in several adjacent rings, as sources do, and therefore a possibility of confusion exists. To test this potential source of contamination, we repeated the analysis excluding a band of $`5\mathrm{deg}`$ centered on the Galactic plane (corresponding to a $`4.36\%`$ of the sky area). The top panels of figures 3 and 4 show the number of detections and the ratio of spurious/true detections, respectively, for the optimal adaptive filter in the cases where the rings have been averaged (as explained before, filled circles) and where they have not (open circles). In the first section of table 2 are the results for the case of averaged rings (the tabulated quantities are the same than in table 1) are shown. The decrease in the number of detections corresponds to the one expected for a uniform distribution of sources in the sky (around $`5\%`$). This indicates that the density of detections around the Galactic plane is not substantially different from the density of detections in other regions, less ’contaminated’, of the sky. This can be seen in figure 6, where the $`4\sigma `$ detections have been represented in Galactic coordinates (the Galactic Plane being represented by a dashed line).

The ratio between spurious and true detections remains almost untouched for low $`\sigma `$. For higher $`\sigma `$ levels, the proportion of spurious detections that are due to the Galactic plane increases dramatically (comparing tables 1 and 2 we see that this proportion increases from a $`7\%`$ at $`2.5\sigma `$ to a $`89\%`$ at $`5\sigma `$). This indicates that the contamination by the Galaxy dominates at high signal to noise ratios, whereas low-intensity contamination is dominated by noise fluctuations. In all cases, the number of peaks that correspond to the Galaxy is much smaller than the number of detected sources, thus implying that the optimal adaptive filter deals efficiently with large scale backgrounds.

We stated before that the higher number of spurious detections of MHW could be due to its non-optimal performance at large scales. To further test this hypothesis we repeated the test for the MHW, now considering only the peaks found outside a $`5\mathrm{°}`$ band centered in the Galactic plane, as we did with the optimal adaptive filter before. The results are shown in the last section of table 2. Most of the spurious sources (specially at high $`\sigma `$) lied near the Galactic plane, as we expected. And yet the remaining number of spurious sources is still greater than in the equivalent optimal adaptive filter case. The number of detections and the flux limits remain similar to the optimal adaptive filter case. The m.r.a.e. is also similar in the two cases. Finally, the mean bias in the determination of the amplitude is negative for high $`\sigma `$ thresholds, as happens with the optimal adaptive filter. Note that the MHW used in this work has the same scale as the antenna. In fact, that scale is not the optimal for detection (Vielva et al. 2000a). The optimal scale of the MHW for a particular case has to be determined from the power spectrum of the data. This in a certain way mimics the determination of the optimal scale $`R_0`$ that is automatically included in the optimal adaptive filter via eqs. (4) and (2). The MHW at its optimal scales resembles the shape of the optimal adaptive filter in Fourier space (in the sense that its maximum is located near the maximum of the optimal adaptive filter) and thus the effectiveness of both filters should be similar.

In order to further decrease the number of spurious detections, we could take advantage of the fact that, in many realistic cases, due to the sky coverage of the experiment and its scanning strategy, many positions of the sky can be measured more than once (that is, at different epochs). For example, a source of $`0\mathrm{deg}`$ latitude will be detected once when the center of the ring is located on longitude $`\varphi =\varphi _{source}R`$, being $`R`$ the radius of the ring, and once again several months later, when the center of the ring is located on longitude $`\varphi =\varphi _{source}+R`$. When this occurs, it would be possible to almost duplicate the amount of information in some areas of the sky and therefore improve both the sensitivity and the reliability of the detection. However, such a refined detection can not be done when the data do not cover the whole sky, and therefore is useless for the construction of an early catalogue of sources during the mission’s flight.

## 6 Conclusions

We have simulated and analysed a sequence of time ordered data such as the one that the future Planck 30 GHz LFI28 channel will produce after the first 6 months of flight. The data include all the main foregrounds as well as CMB fluctuations, point sources and instrumental noise. The resulting TOD has been ring-averaged in order to remove pixel-scale noise and then filtered with an optimal adaptive filter that includes the spectral properties of the data in order to maximize the detection of sources of a particular shape (Gaussian) and scale (the scale of the antenna). The optimal adaptive filter was designed to produce an unbiased, efficient estimator of the amplitude of the sources at their position and to give a maximum of detections at the characteristic scale of the sources. The detection of the sources was performed by thresholding the filtered TOD and looking for connected sets of peaks belonging to adjacent rings. At a $`4\sigma `$ detection level 224 sources over a flux of 0.88 Jy are detected with a mean relative error (in absolute value) of $`20.98\%`$ and a systematic bias of $`7.69\%`$. The position of the sources in the sky is determined with errors inferior to the size of the antenna. The catalogue of detected sources is complete at fluxes $``$ 4.337 Jy. The number of spurious detections is 22.

The performance of the optimal adaptive filter has been compared with the performances under the same conditions of a Gaussian window and a Mexican Hat Wavelet (MHW) of width equal to the beam width. The number of sources detected with the optimal adaptive filter is comparable to the number of sources detected with the MHW and much higher than the number of sources detected with the Gaussian. However, the optimal adaptive filter finds a significantly lower number of spurious sources than the other two filters. This is due to the fact that the optimal adaptive filter removes better the large scale fluctuations (e.g. the Galaxy) than the other two filters. To test this hypothesis a $`5\mathrm{deg}`$ band around the Galactic Plane was removed and the analysis was repeated with the optimal adaptive filter and the MHW. This test shows that most of the spurious detections that were found with the MHW were located on the Galactic Plane, whereas spurious detections that were found with the optimal adaptive filter are uniformly distributed in the sky.

The mean absolute error in the determination of amplitudes is slightly higher in the case of the optimal adaptive filter than in the MHW case. The estimation of the amplitudes with the Gaussian window suffer from very large errors. In most cases, the Gaussian detects only the Galaxy.

The number of spurious detections can be reduced by means of further analysis after the filtering. For example, when the sky coverage of the scan is big enough some areas of the sky can be observed twice or more times, increasing the signal to noise ratio of the sources that lie in such areas. This point is out of the scope of the present work, in which we only present the filter as a first step in the data reduction.

In conclusion, the optimal adaptive filter is an efficient, unbiased and reliable tool for the detection and extraction of punctual sources from TOD. A few hundred of sources over 1 Jy will be detected in the 30 GHz channels of the future Planck mission with $`10\%`$ of spurious detections. The catalogue will be complete at fluxes $`4`$ Jy. One possible application of this technique would be the elaboration of early catalogues of sources. A greater number of detections is expected at higher frequencies. The simulation and analysis of such frequencies will be performed in a future work.

## Acknowledgments

The authors strongly appreciate all the comments, suggestions and careful reading done by Matthias Bartelmann and Fabio Pasian.

The TOD used in this work for point source detection was generated by the Planck pipeline simulator of the DPC Level-S. The sky simulations, pointings and other data is available for the Planck consortia at http://www.mpa-garching.mpg.de/SimData. We specially thank Matthias Bartelmann and Klaus Dolag for their help in the generation of the simulations and their support during this work.

The work presented here has used extensively the software package HEALPix (Hierarchical, Equal Area and iso-Latitude Pixelisation of the sphere); we thank the developing group: Krzysztof M. Gorski, E.F. Hivon, Benjamin D. Wandelt, Antony J. Banday, F.K. Hansen and Matthias Bartelmann. The HEALPix package is available at this address: http://www.eso.org/science/healpix

We thank Laura Cayón and Patricio Vielva for useful suggestions and comments. DH acknowledges support from a Spanish MEC FPU fellowship. We thank FEDER Project 1FD97-1769-c04-01, Spanish DGESIC Project PB98-0531-c02-01 and INTAS Project INTAS-OPEN-97-1192 for partial financial support.


text:
# A first test of Wigner function local realistic model

Abstract

We perform a first experimental test of the local realistic model of Ref. . Our results disfavour the model and confirm standard quantum mechanics predictions.

PACS: 03.65.Bz

Keywords: Bell inequalities, non-locality, hidden variable theories

In 1935 Einstein-Podolsky-Rosen  suggested that Quantum Mechanics (QM) could be an incomplete theory, representing a statistical approximation of a complete deterministic theory where observable values are fixed by some hidden variable. A fundamental progress in discussing possible extensions of QM was the discovery of Bell  that any realistic Local Hidden Variable (LHV) theory must satisfy certain inequalities, which can be violated in QM leading to the possibility of an experimental test of the validity of QM as compared to LHV. Since then many experiments have been devoted to a test of Bell inequalities , leading to a substantial agreement with quantum mechanics and disfavouring realistic local hidden variable theories. However, due to the low total detection efficiency, no experiment has yet been able to exclude definitively realistic local hidden variable theories  It must be noticed that a recent experiment  performed using Be ions has reached very high efficiencies, around 98 %, but in this case the two subsystems (the two ions) are not really separated systems during measurement and the test cannot be considered a real implementation of a detection loophole free test of Bell inequalities ., for it is necessary a further additional hypothesis , stating that the observed sample of particles pairs is a faithful subsample of the whole. This problem is known as detection or efficiency loophole. Because of this loophole a local realistic model can still be built in order to reproduce present experimental data.

Very recently Casado et al.  have presented a local realistic model addressed to be compatible with the available experimental data on Bell inequalities obtained with entangled photons, which are the most severe test of local realism realised up to now. This model represents the completion of series of papers where this scheme has been developed . The main idea is that the probability distribution for the hidden variable is given by the Wigner function, which is positive for photons experiments. Furthermore a model of photodetection, which departs from quantum theory, is built in order to reproduce available experimental results.

A great merit of this model is that it gives a number of constraints, which do not follow from the quantum theory and are experimentally testable.

In particular, there is a minimal light signal level which may be reliably detected: a difference from quantum theory is predicted at low detection rates, namely when the single detection rate $`R_S`$ is lower than

$$R_S<\frac{\eta F^2R_c^2}{2Ld^2\lambda \sqrt{\tau T}}$$
(1)
where $`\eta `$ is the detection quantum efficiency, F is the focal distance of the lens in front of detectors, $`R_c`$ is the radius of the active area of the non-linear medium where entangled photons are generated, $`\tau `$ is the coherence time of incident photons, d is the distance between the non-linear medium and the photo-detectors, $`\lambda `$ the average wavelength of detected photons. L and T are two free parameters which are less well determined by the theory : L can be interpreted as the active depth of the detector, while T is the time needed for the photon to be absorbed and should be approximately less than 10 ns , being, in a first approximation, the length of the wave packet divided for the velocity of light.

How the departure from quantum mechanics would happen is not completely clear in Ref., anyway a lowering of the visibility should be expected .

In order to perform a first experimental test of this model we have reconsidered our set-up for measuring Bell inequalities based on the generation of non-maximally entangled photon states by the superposition of the parametric down conversion (PDC) light of two type I crystals . The use of non-maximally entangled state was suggested by the theoretical discovery that for non-maximally entangled pairs a total efficiency larger than 0.67  (in the limit of no background) is required to obtain an efficiency-loophole free experiment (even if the violation turns to be smaller), whilst for maximally entangled pairs this limit rises to 0.81.

More in details, our experiment follows and develops  an idea by Hardy  and contemplates the creation of a polarisation (non maximally-) entangled states of the form

$$|\psi =\frac{|H|H+f|V|V}{\sqrt{(1+|f|^2)}}$$
(2)
(where $`H`$ and $`V`$ indicate horizontal and vertical polarisations respectively). The state of Eq. 2 would be maximally entangled when $`f=1`$.

With reference to the sketch of figure 1, our set up is made of two crystals of LiIO3 (10x10x10 mm) 250 mm apart, a distance smaller than the coherence length of the pumping laser. This guarantees undistinguishability in the creation of a couple of photons in the first or in the second crystal. A couple of planoconvex lenses of 120 mm focal length centred in between focalises the spontaneous emission from the first crystal into the second one maintaining the angular spread. A hole of 4 mm diameter is drilled into the centre of the lenses in order to allow transmission of the pump radiation without absorption and, even more important, without adding stray-light, because of fluorescence and diffusion of the UV radiation. This configuration, which performs as a so-called ”optical condenser”, was chosen among others as a compromise between minimisation of aberrations (mainly spherical and chromatic) and losses due to the number of optical components. The pumping beam at the exit of the first crystal is displaced from its input direction by birefringence: the small quartz plate (5 x5 x5 mm) in front of the first lens of the condensers compensates this displacement, so that the input conditions are prepared to be the same for the two crystals. Finally, a half-wavelength plate immediately after the condenser rotates the polarisation of the Argon beam and excites in the second crystal a spontaneous emission cross-polarised with respect to the first one. The dimensions and positions of both plates are carefully chosen in order not to intersect the spontaneous emissions at 633 and 789 nm. With a phase matching angle of $`51^o`$ they are emitted at $`3.5^o`$ and $`4^o`$ respectively.

We use as photo-detectors two avalanche photodiodes with active quenching (EG&G SPCM-AQ) with a sensitive area of 0.025 $`mm^2`$ and dark count below 50 counts/s. PDC signal is coupled to an optical fibre (carrying the light on the detectors) by means of a microscope objective with magnification 20, preceded by a polariser.

The output signals from the detectors are routed to a two channel counter, in order to have the number of events on single channel, and to a Time to Amplitude Converter circuit, followed by a single channel analyser, for selecting and counting coincidence events.

In our experiment the single detection rate is at the most of $`R_S=10^5`$ counts per second. Referring to the parameters of Eq. 1 we have $`\eta =0.51\pm 0.02`$ (a value which we have directly measured by using PDC detector calibration ), F= 0.9 cm, $`R_c=1`$ mm, d= 0.75 m and $`\tau =4.210^{13}s`$ (due to spectral selection by an interferential filter). L can be estimated of $`310^5`$ m . This leads to $`T\text{ }>1s`$, extremely higher than the limit of 10 ns suggested in the model . Thus, we are strictly in the condition where quantum mechanics predictions are expected to be violated and, in particular, a strong reduction of visibility is expected.

Nevertheless, our experimental data on coincidences of the two photodetectors show, when the angle of a polariser is varied, a high visibility, $`V=0.98\pm 0.01`$. This result is in strong disagreement with the Wigner function local realistic model of ref. .

With the purpose of performing a further test of model let us consider our results about the Clauser-Horne sum:

$$CH=N(\theta _1,\theta _2)N(\theta _1,\theta _2^{})+N(\theta _1^{},\theta _2)+N(\theta _1^{},\theta _2^{})N(\theta _1^{})N(\theta _2)$$
(3)
which is strictly negative for every local realistic theory. In (3), $`N(\theta _1,\theta _2)`$ is the number of coincidences between channels 1 and 2 when the two polarisers are rotated to an angle $`\theta _1`$ and $`\theta _2`$ respectively, whilst $`N(\theta _1^{})`$ and $`N(\theta _2)`$ denote single counts on detectors 1 and 2 respectively.

For quantum mechanics $`CH`$ can be larger than zero. For a maximally entangled state the largest value is obtained for $`\theta _1=67^o.5`$ , $`\theta _2=45^o`$, $`\theta _1^{}=22^o.5`$ , $`\theta _2^{}=0^o`$.

For non-maximally entangled states the angles for which CH is maximal are somehow different and the largest value smaller. The angles corresponding to the maximum can be evaluated maximising Eq. 3 with

$`\begin{array}{c}N[\theta _1,\theta _2]=[ϵ_1^{||}ϵ_2^{||}(Sin[\theta _1]^2Sin[\theta _2]^2)+\hfill \\ ϵ_1^{}ϵ_2^{}(Cos[\theta _1]^2Cos[\theta _2]^2)\hfill \\ (ϵ_1^{}ϵ_2^{||}Sin[\theta _1]^2Cos[\theta _2]^2+ϵ_1^{||}ϵ_2^{}Cos[\theta _1]^2Sin[\theta _2]^2)\hfill \\ +|f|^2(ϵ_1^{}ϵ_2^{}(Sin[\theta _1]^2Sin[\theta _2]^2)+ϵ_1^{||}ϵ_2^{||}(Cos[\theta _1]^2Cos[\theta _2]^2)+\hfill \\ (ϵ_1^{||}ϵ_2^{}Sin[\theta _1]^2Cos[\theta _2]^2+\hfill \\ ϵ_1^{}ϵ_2^{||}Cos[\theta _1]^2Sin[\theta _2]^2)\hfill \\ +(f+f^{})((ϵ_1^{||}ϵ_2^{||}+ϵ_1^{}ϵ_2^{}ϵ_1^{||}ϵ_2^{}ϵ_1^{}ϵ_2^{||})(Sin[\theta _1]Sin[\theta _2]Cos[\theta _1]Cos[\theta _2])]/(1+|f|^2)\hfill \end{array}.`$ (11)
where (for the case of non-ideal polariser) $`ϵ_i^{||}`$ and $`ϵ_i^{}`$ correspond to the transmission when the polariser (on the branch $`i`$) axis is aligned or normal the polarisation axis respectively.

In our experiment we have created a state with $`f0.4`$ which gives the largest violation of Clauser-Horn inequality for $`\theta _1=72^o.24`$, $`\theta _2=45^o`$, $`\theta _1^{}=17^o.76`$ and $`\theta _2^{}=0^o`$.

However, because of low detection efficiency we must, as in any experiment performed up to now, substitute in Eq. 3 single counts $`N(\theta _1^{})`$ and $`N(\theta _2)`$ with coincidence counts $`N(\theta _1^{},\mathrm{})`$ and $`N(\mathrm{},\theta _2)`$, where $`\mathrm{}`$ denotes the absence of selection of polarisation for that channel. This is one of the form in which detection loophole manifests itself. In the following we will denote Clauser-Horn sum where the former substitution is done with $`CH_{exp}`$.

Our experimental result is $`CH_{exp}=513\pm 25`$ coincidences per second, which is more than 20 standard deviations from zero and compatible with the theoretical value predicted by quantum mechanics. Thus, our results show a strong violation of Clauser-Horn inequality, in agreement with standard quantum mechanics, even in a regime where the bound of Eq. 1 is strongly violated and therefore strongly disfavour the model of Ref. .

When more strict bounds on the parameters of the model and more precise experimental predictions will be available  an experiment, specifically addressed to testing it, will be able to give a conclusive test of the model itself.

Acknowledgements We would like to acknowledge support of MIUR.

Figure caption

\- Sketch of the source of polarisation entangled photons. NLC1 and NLC2 are two $`LiIO_3`$ crystals cut at the phase-matching angle of $`51^o`$. L1 and L2 are two identical piano-convex lenses with a hole of 4 mm in the centre. C is a 5 x 5 x 5 mm quartz plate for birefringence compensation and $`\lambda /2`$ is a first order half wave-length plate at 351 nm. UV identifies the pumping radiation at 351 nm.