text:
# How big is a quantum computer?

Accounting for resources is the central issue in computational efficiency. We point out physical constraints implicit in information readout that have been overlooked in classical computing. The basic particle-counting mode of read-out sets a lower bound on the resources needed to implement a quantum computer. As a consequence, computers based on classical waves are as efficient as those based on single quantum particles.

Every quantum information processor must be coupled to a decoder, a device that condenses the contents of the quantum register into the classical information that makes up the computational output. The details of the decoder itself can have a profound effect on the power of a given quantum information processor. In the laboratory, decoding is usually done by particle counting, using a detector with finite spatial, temporal, spectral, etc., resolution – and this has some significant implications for the resources that may be needed to implement a quantum processor.

The formal theory of processor operation (the quantum state-space algebra) does not sharply distinguish between a quantum information processor consisting of a few particles occupying many states or many particles occupying a few states. But obviously the decoder (particle counter) can make the distinction easily. The distinction is closely related to the different kinds of entanglement central to advantages accrued by working in the quantum domain. For a single particle, entanglement is usually defined in terms of states belonging with different degrees of freedom of the particle (and associated with commuting observables). For many-particle systems entanglement is possible in addition between identical degrees of freedom of different particles. In any case a particle-counting register readout simply provides a dichotomic answer to the question as to whether there was a particle in a particular mode of the system as a whole.

Information processing by photons provides an illustration. A typical apparatus contains photodetectors whose function is to count by photoionization the number of photons in a certain spatial and spectral range. The photon-counting detector can therefore be understood to give a direct answer to the question: How many photons are in detector mode $`\lambda `$? Note that the mode label $`\lambda `$ refers to a set of numbers specifying the modal character. For a photon, these might be the frequency, polarization and wavevector. A different and equally uncomplicated example is provided by the measurement of the electronic state of a Rydberg atom. Standard ramped-field ionization detection gives an answer to the question: Is a particular electronic state, representing an eigenmode of the Schrödinger equation, occupied or not? In fact, because of the particle-counting register readout, a quantum information processor that is implemented by means of the modal entanglement of a single particle can be implemented with equal efficiency by classical wave interference. It may be that detection schemes that do not rely on particle counting could show improvements beyond those of classical wave interferometers. However, we are aware of no laboratory schemes for implementing such measurements.

Any von-Neumann type measurement requires for each of $`N`$ eigenvalues of interest a distinct measurement capability. In the simplest example one can envisage these might be the $`N`$ channels in a Zeeman or Stern-Gerlach apparatus. For such a case, one says that a “resource” of order $`N`$ must be on hand. With this notion of a measurement and its resources, we can quantify the “resource space” needed to implement a combination of processor and decoder. It is simply the total spatio-temporal volume of the modes defining the measurement capabilities. This aspect of the resource requirement is usually ignored in classical computer science, because it is assumed that such a volume may be infinitesimal. Yet this is completely unphysical: quantum mechanics requires that the phase-space volume associated with each mode must be at least $`\mathrm{}`$, and consequently the spatio-temporal volume of a mode cannot be zero for a system with finite energy. In addition to counting the resources required for the processor, quantum mechanics therefore forces us to reassess the resources required for readout.

In the specific case of a single-particle processor, the quantum register that is processed during the computation is expressed by a quantum state in which the particle occupies a superposition of $`N`$ different modes. The readout of such a register then requires at least $`\mathrm{log}_2N`$ distinct detectors or generalized resources. The classical outputs of these detectors then have to be polled to locate the sequence of detectors that fired. This is simply the classical search problem for a sorted database, and it is evident that such a search can be done most efficiently in $`\mathrm{log}_2N`$ steps by means of a binary search. This method of resource counting can be adopted to show that any single-particle processor can be efficiently implemented using classical wave interference, even when the system relies exclusively on the entanglement of different degrees of freedom of the particle.

A description of the quantum processor in terms of the system state vector alone makes it difficult to distinguish single-particle from many-particle situations, let alone to understand the optimal combinations of particles and degrees of freedom that might be needed to accomplish a certain processing task. The essential roles of interference and entanglement can be examined unambiguously most simply in the language of second quantization. The field operator $`\widehat{\psi }(x)`$ describes the annihilation of atoms, photons or other types of particles at position $`x`$. If we introduce a complete set of orthonormal spatial modes $`u_\lambda (x)`$, where $`\lambda `$ labels the different modes including the possibility of having several degrees of freedom, we may introduce mode annihilation operators $`\widehat{a}_\lambda `$ via $`\widehat{\psi }(x)=_\lambda \widehat{a}_\lambda u_\lambda (x)`$, where $`[\widehat{a}_\lambda ,\widehat{a}_\mu ^{}]_\pm =\delta _{\lambda \mu }`$ depending on whether the field describes fermions or bosons.

If the unitary operator $`\widehat{U}`$ defines the evolution of a quantum processor, then the output register is given by

$$\widehat{\psi }(x,t)=\widehat{U}(t)\widehat{\psi }(x,0)\widehat{U}^{}(t),$$
(1)
where $`\widehat{\psi }(x,0)`$ represents the input register. The only restriction on the Hamiltonian generating $`\widehat{U}`$ is that it preserves the particle number of the processor register: If this register contains only a single particle, so that the state is given by $`|\psi =_\lambda \psi _\lambda \widehat{a}_\lambda ^{}|\mathrm{vac}`$, the only terms in the normally ordered Hamiltonian contributing to the dynamical evolution are those bilinear in $`\widehat{a}_\lambda ^{}`$ and $`\widehat{a}_\lambda `$. Thus the unitary evolution (1) reduces from the operators to the mode functions themselves:

$`\widehat{\psi }(x,t)={\displaystyle \underset{\lambda }{}}\widehat{a}_\lambda u_\lambda (x,t),`$ (2)
where $`u_\lambda (x,t)`$ are the time-propagated mode functions of wave mechanics.

The operation of an elementary measurement can be formulated in the language of fields as the action of a linear filter on the underlying mode structure of the system. The register field operator after the filter is given by

$$\widehat{\psi }_\mu (x,t)=𝑑x^{}𝑑t^{}\mathrm{\Gamma }_\mu (x,t|x^{},t^{})\widehat{\psi }(x^{},t^{}),$$
(3)
where $`\mu `$ denotes the mode that is to be probed for particles. The observable corresponding to an elementary measurement is given by the number operator for this mode

$$\widehat{N}_\mu =𝑑x𝑑t\widehat{\psi }_\mu ^{}(x,t)\widehat{\psi }_\mu (x,t).$$
(4)
In each measurement one determines the number of particles in the specified detector mode. This can be described by the projector $`|N;\mu N;\mu |=\widehat{\delta }(\widehat{N}_\mu N)`$, where $`|N;\mu `$ denotes the state with $`N`$ particles in the detected mode $`\mu `$ and $`\delta (N)`$ is understood as the Kronecker delta function. In the case of a single particle, $`N`$ may take the values $`(0,1)`$, so that the projector simplifies to

$$\widehat{\delta }(\widehat{N}_\mu N)=\delta (N)(\widehat{1}_\mu \widehat{N}_\mu )+\delta (N1)\widehat{N}_\mu ,$$
(5)
showing that the measurement projector can be entirely expressed in terms of the bilinear number operator (4).

As indicated previously, the minimum number of detectors required to read out the $`N`$-channel register is $`\mathrm{log}_2N`$. This can be achieved for any quantum system by representing each of the mode labels $`\mu `$ by a set of $`\mathrm{log}_2N`$ bits $`\{b_i\}`$ and grouping the number operators of the modes, $`\widehat{N}_{\{b_i\}}`$, into a specific set of detectable number operators $`\widehat{M}_\alpha `$ according to their common bits. For example, in a system of $`N=8`$ different modes (i.e. with a Hilbert space of eight dimensions) the three number operators are

$`\widehat{M}_{1xx}`$ $`=`$ $`\widehat{N}_{\{1,0,0\}}+\widehat{N}_{\{1,0,1\}}+\widehat{N}_{\{1,1,0\}}+\widehat{N}_{\{1,1,1\}},`$ (6)
$`\widehat{M}_{x1x}`$ $`=`$ $`\widehat{N}_{\{0,1,0\}}+\widehat{N}_{\{0,1,1\}}+\widehat{N}_{\{1,1,0\}}+\widehat{N}_{\{1,1,1\}},`$ (7)
$`\widehat{M}_{xx1}`$ $`=`$ $`\widehat{N}_{\{0,0,1\}}+\widehat{N}_{\{0,1,1\}}+\widehat{N}_{\{1,0,1\}}+\widehat{N}_{\{1,1,1\}}.`$ (8)
The readout would then occur via a cascaded sequence of particle counters. The result of measuring the three observables (6) is the set of bits specifying the mode. Here too the projectors corresponding to the operators $`\widehat{M}_\alpha `$ take the form (5). Note that it is not necessary that each of the bits correspond to a different degree of freedom of the particle, i.e. $`\mu =\{\mu _1,\mu _2,\mathrm{}\}`$. The binary readout scheme is possible even for a single particle excited in one degree of freedom. In practice it may prove simpler to use different degrees of freedom (such as the principle quantum number, and the two angular momentum quantum numbers for a Rydberg electron), though no such scheme has been analyzed in detail as yet. Though the number of detectors is minimal, the number of resources in terms of space-time volume of required modes is not, as can be seen from the fact that each measured number operator involves more than one mode.

If the result of the computation leaves the register in one of the detected modes, then the readout provides directly a classical result for the computation. In general, however, the computation will place the register in a superposition state, so that the readout will be probabilistic. That is, for many runs of the same computation the classical outputs will be different. Therefore the whole computational process has to be repeated many times, either sequentially or in parallel. Thus in both cases one must determine for each mode the expectation value of the projector (5).

For a single particle it is straightforward to show using Eqs (3) – (5) that this expectation value is proportional to the correlation function

$$C(x,t;x^{},t^{})=\widehat{\psi }^{}(x,t)\widehat{\psi }(x^{},t^{}).$$
(9)
Using the general single-particle state $`|\psi =_\lambda \psi _\lambda \widehat{a}_\lambda ^{}|\mathrm{vac}`$ and the dynamics according to Eq. (2) we obtain

$$C(x,t;x^{},t^{})=\psi ^{}(x,t)\psi (x^{},t^{}),$$
(10)
where we have defined the single-particle wavefunction as $`\psi (x,t)=_\lambda \psi _\lambda u_\lambda (x,t)`$. This result holds for both the single-atom and single-photon case and is equivalent to a classical interference pattern obtained by replacing the expectation value (10) by the classical electric-field correlation function

$$C_{\mathrm{class}}(x,t;x^{},t^{})=E^{}(x,t)E(x^{},t^{}).$$
(11)
where $`E(x,t)`$ is an analytic signal. A comparison of Eqs (10) and (11) shows a complete one-to-one correspondence between classical interference and single-particle quantum interference, not only in the dynamics but also in the measurement. We conclude therefore that any quantum computer based on a single-particle quantum register can be implemented with equal efficiency entirely by a classical interferometer. This is because the concept of the entanglement of degrees of freedom of a single particle cannot be attributed an inherently quantum character: it is perfectly understandable in terms of classical wave interference.

This formulation also provides some insight into the issue of the speedup of a particular algorithm using quantum interference. Consider, for example, a quantum search algorithm that is based on single-particle entanglement, that is, on the correlation of different degrees of freedom of a single particle at the level of probability amplitudes. In this case the quantum register may be realized by $`N1`$ modes but the modes can only be occupied by a total population of 1. In the optical case, these might be the modes associated with the polarization and the wavevector degrees of freedom. Alternatively the energy and angular momentum degrees of freedom of a Rydberg electron may be used.

The readout requires $`O(N)`$ resources and at minimum $`\mathrm{log}_2N`$ detectors, whose classical outputs can be searched in not less than $`O(\mathrm{log}_2N)`$ steps. In fact, it is evident that a classical search is always required in the final readout, since the values of the output quantum register only obtain physical reality after having been measured. Therefore, any quantum computation with a single particle cannot be faster than $`O(\mathrm{log}_2N)`$. Moreover, by virtue of the one-to-one correspondence described above, it can be performed using classical interference with the same speed and the same number of resources. This is in striking contrast to the belief that a quantum-search algorithm, even without using multi-particle entanglement, i.e. with only a single particle, can be advantageous for solving the search problem.

The quantum computation itself consists of preparing the input register and performing a unitary $`N\times N`$ transformation $`\widehat{U}`$. According to Grover, this unitary transformation may perhaps consist of a series of operations like passing the single particle through a so-called Oracle, whose function is to alter the phase of one of the modes, followed by an $`N\times N`$-port beamsplitter. These steps are usually referred to as “querying the Oracle” and “inversion about the mean”. Using these operations it has been shown that an unsorted database can be searched by performing only $`O(\sqrt{N})`$ “queries” or even only a single “query” of the Oracle. This is in stark contrast to the corresponding classical search which requires $`O(N)`$ queries.

The caveat, as Steane has pointed out, is that one should be careful of comparing an inefficient classical algorithm with an efficient quantum one. The experimental implementations of the Grover search to date have not in fact implemented a search of an unsorted database, but rather a database with one single marked item. Therefore experimental evidence for the claimed speedup has yet not been achieved, since a sorted database can be searched classically in $`\mathrm{log}_2N`$ steps, which is identical to the readout limit for a quantum computer.

The superiority of the quantum search algorithm becomes apparent only when one examines carefully the notion of a query. It is evident that there is no information gain in the sort of oracle query described above, which involves only unitary transformations. It is in the encoding of such oracles, which must be done by a quantum computer, that the real speedup occurs. Once a properly encoded oracle is available, the number of steps required to perform the information processing is then limited by the final register readout to $`O(\mathrm{log}_2N)`$. This is exactly the same as for a classical information processor, which shows that it is the realization of actual information by the readout, as opposed to predictive information that is contained in the quantum state, that is the ultimate limiting procedure in quantum information processing.

We have shown here that if the readout of the register is performed by particle counting then there exists a one-to-one correspondence of single-particle quantum interference and classical interference. Therefore we conclude that any enhancements in processing power that can be ascribed to quantum interference can also be found in classical wave processors, and this includes systems based on modal entanglement. Multi-particle entanglement, on the other hand, may provide enhancements that cannot be efficiently transcribed to classical interferometers, even when particle counting is used to realize the output information.

Acknowledgments
This work was supported by an ARO-administered MURI Grant No. DAAG-19-99-1-0125 (IAW and JHE), and by NSF Grant No. PHY-9415583 (JHE). SW acknowledges support by the Studienstiftung des deutschen Volkes. We acknowledge several stimulating interactions with P. Bucksbaum, P.L. Knight, M.G. Raymer, C.R. Stroud, Jr. and K. Wodkiewicz.

Correspondence should be addressed to I. A. Walmsley (email: walmsley@optics.rochester.edu).


format: text

text:
# Correlations and fluctuations of a confined electron gas

## Abstract

The grand potential $`\mathrm{\Omega }`$ and the response $`R=\mathrm{\Omega }/x`$ of a phase-coherent confined noninteracting electron gas depend sensitively on chemical potential $`\mu `$ or external parameter $`x`$. We compute their autocorrelation as a function of $`\mu `$, $`x`$ and temperature. The result is related to the short-time dynamics of the corresponding classical system, implying in general the absence of a universal regime. Chaotic, diffusive and integrable motions are investigated, and illustrated numerically. The autocorrelation of the persistent current of a disordered mesoscopic ring is also computed.

PACS numbers: 73.23.-b, 73.20.Dx, 05.45.Mt

A major advance in mesoscopic physics has been the observation of single-electron interference or more generally of quantum mechanical effects in mesoscopic objects. Among these we can mention the supershell structure in atomic clusters, the persistent currents in mesoscopic rings, or the quantization of the longitudinal conductance in microjunctions .

In all these experiments the spatial confinement of the electron gas plays a significant role. When a parameter controlling the system is varied – like the chemical potential, a magnetic field, or the shape of the confining potential – the response of the gas shows irregular fluctuations as a function of the parameter. While the average value of the variation is related to a macroscopic property of the system (like its volume), the wild fluctuations are the imprint of quantum mechanical interference.

Neglecting the interactions between the electrons, our purpose is to provide a general description of the fluctuating properties of the many-particle ground state electron gas when different parameters – like temperature $`T`$, chemical potential $`\mu `$, or external parameter $`x`$ – are varied. Taking as reference parameters $`(x_0,\mu _0,T_0)`$, the system is described in the grand-canonical ensemble by the grand potential

$$\mathrm{\Omega }(x_0,\mu _0,T_0)=𝑑E𝒩(E,x_0)f(E,\mu _0,T_0),$$
(1)
with $`𝒩(E,x_0)=g_s_k\mathrm{\Theta }(EE_k(x_0))`$ the counting function of the single particle states $`E_k(x_0)`$, $`\mathrm{\Theta }`$ the Heavyside step function, $`f(E,\mu _0,T_0)=1/[1+exp((E\mu _0)/kT_0)]`$ the Fermi distribution and $`g_s=2`$ for spin degeneracy.

The reaction of the system to a variation of the external parameter is characterized by the response function

$$R(x_0,\mu _0,T_0)=\mathrm{\Omega }/x_0|_{\mu _0,T_0}.$$

The physical interpretation of $`R`$ depends on the geometry of the sample and the nature of the parameter $`x`$. It is a magnetic moment if $`x`$ is a magnetic field, a persistent current when the sample geometry is annular and $`x`$ is an Aharonov-Bohm flux, or a force when $`x`$ controls the shape of the confining potential (like when deforming a metallic cluster for example).

We characterize the fluctuating properties of the electron gas by computing the autocorrelation function

$$C_\mathrm{\Omega }(x,\mu ,T)=\stackrel{~}{\mathrm{\Omega }}(𝐪_{})\stackrel{~}{\mathrm{\Omega }}(𝐪_+)_{x_0,\mu _0},$$
(2)
with $`𝐪_\pm =(x_0\pm x/2,\mu _0\pm \mu /2,T_0\pm T/2)`$. An analogous definition holds for $`C_R(x,\mu ,T)`$. We denote the fluctuating part by an upper tilde symbol (see below). The brackets denote some energy and parameter averages.

Parametric correlations like $`C_\mathrm{\Omega }`$ or $`C_R`$ have been studied in the past for chaotic and disordered systems for “individual” (single-particle) levels, rather than for a many-body ground state. The single-particle analog of $`C_R(x,0,0)`$ is the velocity-velocity correlation . Although for the latter no global expression is known, its tail $`x\mathrm{}`$ has been shown to be universal . In contrast to this, we will see below that for a Fermi gas the parametric correlators are explicitly computable and are expressed semiclassically, even at zero temperature, in terms of the classical short time dynamics of the system. Because the latter is system-specific, the correlators are in general system-dependent with no universal regime whatsoever. This difference with the universal local properties occurs mainly because the typical time scale of a perturbed ground state Fermi gas is $`h/\mu _0`$, which is much smaller than the Heisenberg time.

We express the counting function in Eq.(1) as a sum of smooth plus oscillatory terms $`𝒩=\overline{𝒩}+\stackrel{~}{𝒩}`$. The smooth part is given by the usual Weyl or Thomas-Fermi approximation. Semiclassically, the oscillatory part is $`\stackrel{~}{𝒩}=2\mathrm{}g_s_p(A_p/\tau _p)\mathrm{sin}(S_p/\mathrm{})`$. The sum is over the periodic orbits $`p`$ with action $`S_p`$ (including the Maslov index), period $`\tau _p=S_p/E`$, and stability amplitude $`A_p`$ . This sum is inserted in Eq.(1) to compute the oscillatory part of the grand potential. To leading order in $`\mathrm{}`$ and for $`kT_0\mu _0`$ ,

$$\stackrel{~}{\mathrm{\Omega }}(x_0,\mu _0,T_0)2\mathrm{}^2g_s\underset{p}{}\stackrel{~}{A}_p\mathrm{cos}(S_p/\mathrm{}),$$
(3)
with $`\stackrel{~}{A}_p=A_p\kappa (\tau _p,T_0)/\tau _p^2`$. The factor $`\kappa (\tau _p,T_0)=\frac{\tau _p/\tau _c}{\mathrm{sinh}(\tau _p/\tau _c)}`$, with $`\tau _c=\mathrm{}/\pi kT_0`$ takes into account finite temperature effects and acts as an exponential cut off for long periodic orbits. All the classical quantities entering this expression are evaluated at $`(x_0,\mu _0)`$. The main dependence on $`x`$ comes from the oscillatory terms via the dependence of the actions, $`S_p=S_p(x_0,\mu _0)`$. Taking this into account, the oscillating part of the response function $`\stackrel{~}{R}=\stackrel{~}{\mathrm{\Omega }}/x_0|_{\mu _0,T_0}`$ is

$$\stackrel{~}{R}(x_0,\mu _0,T_0)2\mathrm{}g_s\underset{p}{}\stackrel{~}{A}_pQ_p\mathrm{sin}(S_p/\mathrm{}),$$
(4)
with $`Q_p=S_p/x_0|_{\mu _0,T_0}`$.

Replacing Eq.(3) in (2) we obtain a double sum over periodic orbits with oscillations depending on the difference of their actions. The averages in Eq.(2) are done over some range $`\mathrm{\Delta }\mu `$ and $`\mathrm{\Delta }x`$ which contain several oscillations of $`\stackrel{~}{\mathrm{\Omega }}`$ but are on the other hand small in comparison to the scale over which the classical functions vary. As mentioned before, the main dependence on $`(x,\mu )`$ comes from the actions, which may be linearized with respect to these parameters. The average restricts the sum to orbits having approximately the same $`Q_p`$ and $`\tau _p`$. Ordering the orbits by their period, taking into account these restrictions and using the semiclassical definition of the parametric form factor $`K(\tau ,x)`$ (i.e., the Fourier transform of the two-point correlation function ) we obtain

$$C_\mathrm{\Omega }=\frac{\mathrm{}^2g_s^2}{2\pi ^2}_0^{\mathrm{}}\frac{d\tau }{\tau ^4}\chi (\tau ,T)\mathrm{cos}(\mu \tau /\mathrm{})K(\tau ,x),$$
(5)
with $`\chi (\tau ,T)=\kappa (\tau ,T_0T/2)\kappa (\tau ,T_0+T/2)`$. The autocorrelation of the response function follows from Eq.(5) by derivation with respect to the external parameter.

Due to the factor $`\tau ^4`$ in the denominator, the integral in (5) is dominated by the short-time dynamics of the system. For chaotic systems a naive replacement of the form factor $`K(\tau ,x)`$ by the random matrix result produces a divergence. In real systems, this divergence is avoided by the cut off introduced by the shortest periodic orbit of period $`\tau _{min}`$ (the form factor vanishes for $`\tau <\tau _{min}`$). Then even at $`T_0=T=0`$ ($`\chi =1`$), the autocorrelation will be described with good accuracy by computing in Eq.(5) only the contribution of the shortest periodic orbits, the longer ones representing a small correction. Because this short-time dynamics is system-specific, then in general we do not expect any universality for $`C_\mathrm{\Omega }`$ and $`C_R`$.

Eq.(5) may therefore be restricted to the diagonal approximation $`p=p^{}`$ for the form factor, $`K_D(\tau )=h^2_pA_p^2\delta (\tau \tau _p)`$, with the result

$$C_\mathrm{\Omega }=\frac{\mathrm{}^2g_s^2}{2\pi ^2}_0^{\mathrm{}}\frac{d\tau }{\tau ^4}\chi \mathrm{cos}(\mu \tau /\mathrm{})\mathrm{cos}(Qx/\mathrm{})K_D(\tau ).$$
(6)
Here the average is evaluated over the $`Q_p`$’s taken among all the periodic orbits having a period between $`\tau `$ and $`\tau +d\tau `$. The expression for the response function is

$$C_R=\frac{g_s^2}{2\pi ^2}_0^{\mathrm{}}\frac{d\tau }{\tau ^4}\chi \mathrm{cos}(\mu \tau /\mathrm{})Q^2\mathrm{cos}(Qx/\mathrm{})K_D(\tau ).$$
(7)
To leading order, Eq.(6) is also valid for systems with a fixed number of particles $`n`$ . The form factor is now evaluated at $`E=\overline{\mu }(n,x_0)`$ obtained by inverting $`n=\overline{𝒩}(\overline{\mu },x_0)`$. On the other hand, Eq.(7) is modified according to $`Q^2Q^2+\tau ^2\left(\overline{\mu }/x_0\right)^2`$. The importance of this new term depends on the nature of the external parameter. If $`x`$ is a flux in an Aharonov-Bohm geometry, $`\overline{𝒩}`$ is independent of $`x`$ and $`\overline{\mu }/x_0=0`$. On the other hand, if $`x`$ is the deformation of a two dimensional billiard, we find that the contribution of this additional term is comparable to the $`Q^2`$ term if the area of the billiard varies, while it is smaller by a factor $`1/n`$ when the area is kept fixed.

Eqs.(6) and (7) are the basic results of this paper. We now proceed to analyze them for different physical situations. First we consider a ballistic motion in a regular cavity (the trajectories of the electrons occur on phase-space tori). As an example of regular dynamics we have studied a Fermi gas contained in a two dimensional rectangular box of fixed area $`𝒜`$ and sides $`a`$ and $`b`$. We choose $`x=a/b`$ as the external parameter. Exact expressions for $`A_p,\tau _p,S_p`$ and $`Q_p`$ are known in this case. We can therefore explicitly evaluate $`C_\mathrm{\Omega }`$ and $`C_R`$ as a function of $`\mu `$ or $`x`$. In Fig. 1 we show the numerical results for $`x_01.3`$ compared to the theoretical curves using only the first $`10`$ orbits as a function of the rescaled variable $`\xi =\sqrt{\mu _0}x/x_0^{3/2}`$. We observe erratic long-range oscillations accurately described by our formulas.

We now turn our attention to chaotic systems. Considering that the form factor is strictly zero for $`\tau <\tau _{min}`$, the simplest approximation one can made for these systems is to assume a random matrix theory  behavior for $`K_D(\tau )`$ starting at $`\tau =\tau _{min}`$

$$K_D(\tau )=\{\begin{array}{cc}0,\hfill & \tau <\tau _{min}\hfill \\ 2\tau /\beta ,\hfill & \tau _{min}<\tau <<\tau _H.\hfill \end{array}$$
(8)
where $`\beta =1(2)`$ for systems with (without) time-reversal symmetry. Because the short-time dynamics usually dominates the correlation functions, this oversimplified approximation that doesn’t consider individually the contribution of other short trajectories is not expected to be very precise in general. However, what is interesting about Eq.(8) is that it leads to universal expressions for $`C_\mathrm{\Omega }`$ and $`C_R`$, since as we will now see the system-dependent features may be re-absorb by a rescaling of the parameters. These universal results then serve as a guide (in particular in experimental situations where the shape of the potential is not well known), and are useful to obtain simple qualitative estimates of the different characteristic scales involved.

Let’s evaluate using (8) the variance of the fluctuations of the grand potential and of the response function. Restricting to zero temperature, we obtain from Eq.(6)

$$C_\mathrm{\Omega }^0=C_\mathrm{\Omega }(0,0,0)=\frac{g_s^2}{2\beta \pi ^2}\frac{\mathrm{}^2}{\tau _{min}^2},$$
(9)
where $`\tau _{min}`$ is computed at $`(x_0,\mu _0)`$. This simple expression relates the typical fluctuations to the period of the shortest orbit. For a billiard, $`\tau _{min}=(m/2\mu _0)^{1/2}L_{min}`$, with $`L_{min}`$ the length of the orbit. Weyl’s law allows to express $`\mu _0`$ in terms of the average number of particles in the gas $`n=\overline{𝒩}(\mu _0,x_0)`$. We find $`C_\mathrm{\Omega }^0n^{2/d}`$, with $`d`$ the space dimensionality. This may be compared to the variance of the sum of $`m`$ energy levels located in a small energy window around $`\mu _0`$, which is universal and proportional to $`m^2\mathrm{log}m`$ . If a finite temperature is considered, we find that $`C_\mathrm{\Omega }^0`$ is damped by a factor $`4(\tau _{min}/\tau _c)\mathrm{exp}(\tau _{min}/\tau _c)`$ when $`\tau _c<<\tau _{min}`$.

For the variance of the response function we must compute the variance of the distribution of the $`Q_p`$’s. We assume for chaotic systems a Gaussian distribution . The variance is $`Q^2=\alpha \tau `$. Explicit expressions for $`\alpha `$ valid for chaotic billiards were given in . Its value depends on energy, with a dependence $`E^{3/2}`$. These results as well as the approximation (8) in Eq.(7) lead to

$$C_R^0=C_R(0,0,0)=\frac{g_s^2}{\beta \pi ^2}\frac{\alpha }{\tau _{min}}.$$
(10)
It is interesting to compare $`C_R^0`$ to the “local” variance of the response of the individual single-particle levels $`E_k`$ located in a small energy window around $`(x_0,\mu _0)`$, $`v^2=(E_k/x)^2=2\alpha g_s^2/\beta \tau _H`$  ($`\tau _H=h\overline{\rho }`$ is the Heisenberg time, conjugate to the mean level spacing). We have

$$C_R^0=(1/2\pi ^2)(\tau _H/\tau _{min})v^2.$$
(11)
Because $`\tau _H>>\tau _{min}`$, the electron gas amplifies the local fluctuations. For billiards the amplification factor is proportional to $`n^{(d1)/d}`$, while the overall dependence on $`n`$ is $`C_R^0n^{4/d}`$.

The autocorrelation functions are now evaluated in the approximation (8). Setting $`\mu =T=T_0=0`$ in (6) and using $`\mathrm{cos}(Qx/\mathrm{})=\mathrm{exp}(\alpha \tau x^2/2\mathrm{}^2)`$ by the Gaussian assumption we get

$$\frac{C_\mathrm{\Omega }(\xi ,0,0)}{C_\mathrm{\Omega }^0}=(1\xi ^2/2)e^{\xi ^2/2}+\xi ^4\mathrm{\Gamma }(0,\xi ^2/2)/4,$$
(12)
where $`\mathrm{\Gamma }(a,z)=_z^{\mathrm{}}t^{a1}e^t𝑑t`$. We have rescaled the parameter according to

$$\xi =\sqrt{\alpha \tau _{min}}x/\mathrm{}=\sqrt{C_R^0/2C_\mathrm{\Omega }^0}x,$$
(13)
which makes, at this level of approximation, $`C_\mathrm{\Omega }/C_\mathrm{\Omega }^0`$ a universal rapidly-decreasing function. The typical parameter scale over which $`C_\mathrm{\Omega }`$ decorrelates is $`\delta x=\mathrm{}/\sqrt{\alpha \tau _{min}}`$. This is consistent with the previous assumptions made, e.g., $`\delta x`$ is small in the semiclassical $`\mathrm{}0`$ limit. It is interesting to compare $`\delta x`$ to the analogous local scale for a single-particle level , $`\delta x_L=1/\overline{\rho }\sqrt{<v^2>}`$. Since $`\delta x=g_s\sqrt{\tau _H/2\beta \pi ^2\tau _{min}}\delta x_L`$, the scale of the Fermi gas is larger by a factor $`n^{(d1)/2d}`$.

Expressed in terms of the same rescaled variable, the autocorrelation function of the response is

$$\frac{C_R(\xi ,0,0)}{C_R^0}=e^{\xi ^2/2}3\xi ^2\mathrm{\Gamma }(0,\xi ^2/2)/2.$$
(14)
Analogous results are obtained when computing the autocorrelation in energy, $`C_\mathrm{\Omega }(0,ϵ,0)`$ and $`C_R(0,ϵ,0)`$, were the rescaled parameter is $`ϵ=\mu \tau _{min}/\mathrm{}`$.

In order to test numerically the goodness of these approximations we have considered an electron gas in a chaotic cavity. The cavity is the two-dimensional Limaçon billiard , whose shape is defined by the conformal transformation $`w=z+bz^2+ce^{i\delta }z^3`$. For parameter values around $`b=c=0.2`$ it was numerically shown that the dynamics is dominated by chaotic trajectories . We take $`x=\delta `$ as the external parameter. Fig. 1 shows the numerical results compared to the theoretical curves. The approximation (8) provides a good qualitative description of the main features. The deviations observed may be due to contributions of periodic orbits other than the shortest one.

Analogous results may be obtained for diffusive systems by using the appropriate form factor in Eqs.(6) and (7). We concentrate on a particular example, the calculation of the autocorrelation of the persistent current of an electron gas in a disordered ring threaded by an Aharonov-Bohm flux $`\varphi `$.

The current at thermodynamic equilibrium is given by $`I=c\frac{\mathrm{\Omega }}{\varphi }`$. Only the fluctuating part of $`\mathrm{\Omega }`$ depends on $`\varphi `$. We replace $`S_p/\mathrm{}S_p/\mathrm{}+2\pi w_p\varphi /\varphi _0`$ in Eq.(3) to include the flux, with $`\varphi _0=hc/e`$ the flux quantum and $`w_p`$ the winding number of the periodic orbit $`p`$ around the ring. The current in the ring is

$$I(\lambda )=\frac{4\pi \mathrm{}^2cg_s}{\varphi _0}\underset{p}{}\stackrel{~}{A}_pw_p\mathrm{sin}\left(S_p/\mathrm{}+2\pi w_p\lambda \right),$$
(15)
where we have introduced $`\lambda =\varphi /\varphi _0`$. The autocorrelation of the current $`C_I(\lambda ,\lambda _R)=I(\lambda _R+\lambda )I(\lambda _R)`$ is computed in the diagonal approximation taking into account the contribution of the primitive orbits and of their time reversal partners. As before, it can be expressed in terms of the form factor, which now takes into account the diffusive motion of the particle . For simplicity we restrict the calculation to a quasi one-dimensional ring and set the temperature to zero. At $`\varphi =0`$ we get the well known variance of the persistent current $`C_I^0`$ . Normalizing $`C_I`$ with this variance and averaging both with respect to the reference flux $`\lambda _R`$, we arrive at

$$\frac{C_I(\lambda )}{C_I^0}=\frac{1}{\zeta (3)}\underset{w=1}{\overset{\mathrm{}}{}}\frac{\mathrm{cos}(2\pi w\lambda )}{w^3},$$
(16)
where $`\zeta (s)`$ is the Riemann zeta function.

Eq.(16) is a universal function, independent of the specific aspects of the dynamics (i.e., independent of the elastic time between collisions or the diffusion constant), but depends on the dimensionality of the system. No rescaling of the parameter was needed. This interesting fact is due to the geometry of the sample and the physical quantity considered. Orbits with zero winding number do not contribute to the current because they have no flux dependence, while those with $`w1`$ have a period $`\tau >>\tau _e`$, the elastic time. As a consequence, all the non-universal shortest orbits gave no contribution, and the autocorrelation is a universal function (at fixed $`d`$).

In conclusion, we have studied the statistical properties of a phase-coherent noninteracting confined electron gas in response to a parameter. The results are applicable to any confined fermionic system. The autocorrelations $`C_\mathrm{\Omega }`$ and $`C_R`$ were computed explicitly. Their behavior is described by the shortest classical periodic orbits, making them system-dependent with in general no universal regime. When the form factor is simplified to take into account only the shortest orbit, it is possible to remove the dependence on specific properties of the system by a rescaling determined by $`\tau _{min}`$. The typical fluctuations of the gas are also expressed by simple formulas in terms of $`\tau _{min}`$. This is different from the rescaling of parametric correlations of single-particle levels, determined by $`\tau _H`$. This latter time is of secondary importance in our context, and long orbits represent a small correction to our results. In the case of persistent currents in a diffusive Aharonov-Bohm geometry, the autocorrelation is independent of the parameters of the system, without any rescaling, due to the suppression by the geometry of the short-time orbits. The analytical results have been verified by numerical simulations on chaotic and regular cavities.

We are especially grateful to H. Bruus, C. H. Lewenkopf and E. R. Mucciolo for providing the program for the numerical simulation of the Limaçon billiard.


text:
# Angular Signatures for Galactic Halo WIMP Scattering in Direct Detectors: Prospects and Challenges

## I Introduction

The effort to directly detect Galactic halo weakly interacting massive particles (WIMPs) which may form the dark matter dominating our Galaxy is now entering a new phase. Numerous experiments are now operating with a sensitivity to scattering cross sections in the range predicted for low energy SUSY neutralinos , and indeed one such experiment has claimed a preliminary detection .

There is a problem however, with existing detection schemes, which the controversy over the recent DAMA results underscores. The signature for elastic scattering of WIMPs on nuclear targets involves the observation of what is essentially excess noise in the detector. Thus, in order to definitively establish a positive detection, one has to convincingly demonstrate that this excess noise is indeed WIMP induced, and not due to some unexpected terrestrial background. The signature focussed on by the DAMA collaboration, an annual modulation due to the Earth’s orbital motion through the WIMP halo, is problematic in this regard. In the first place, the effect is small, at the few percent level at best . In addition, numerous radioactive backgrounds are known to modulate over the course of a year.

The Earth’s motion around the Sun, combined with the Sun’s motion around the Galactic disk, can however provide another, much larger and more unambiguous signature. If detectors can be developed that are sensitive to the direction of recoil of the nucleus impacted by a WIMP, then the expected signal will have an angular dependence characteristic of this velocity relative to the incident WIMP halo phase space distribution.

For a simple isothermal sphere distribution, the resulting angular differential rates can be easily calculated , and it is clear that a large forward-backward asymmetry can be produced. However, it is quite likely that the actual Galactic halo WIMP distribution is not well represented by such a simple analytical approximation. Building detectors with angular sensitivity is a daunting challenge, although testing of at least one such proposed detector, the DRIFT detector  is currently underway. It is therefore appropriate to explore in advance to what extent uncertainties in our knowledge of a halo WIMP distribution will have on the detector parameters required required to separate a WIMP signal from other backgrounds.

Recently we outlined a comprehensive formalism which allows the generation of expected angular event rates in terrestrial detectors with angular sensitivity for any incident WIMP distribution, and made a provisional exploration of the overall event rates required in detectors for simple halo models in order to separate the signal from a flat background . Here, we follow up by providing a detailed derivation and description of our results, and also by incorporating a much broader spectrum of analytical halo models. In addition, we also make a provisional analysis of the possible effects of spin-dependence of WIMP-nucleon cross sections on the resulting angular signatures. Finally, we briefly describe ongoing work aimed at incorporating more realistic WIMP distributions arising from N-body simulations into our analysis, and applying our results to various specific proposed detectors, where energy-dependence of both signals and backgrounds may play a key role. Our present results, however, point out the challenges of discerning an angular signal in advance of specific knowledge of the halo WIMP distribution. In addition, the formalism provided here can be used to generate angular event rates which can in turn be used by experimental groups as an input into their detector Monte Carlo analyses.

The outline of this work is as follows: First, in section II, we present a detailed derivation of the differential angular cross sections in a terrestrial laboratory frame for any incident WIMP distribution. In section III we discuss a broad range of different analytical halo models that have been proposed in the literature, and outline the relevant parameter ranges of importance for influencing the resulting angular signature for WIMP scattering. In this section we also incorporate effects of the detailed motion of the Earth in a Galactic-rest frame, and also examine the effect of incorporating possible spin-dependence of WIMP-nucleus scattering on the resulting angular signals, as well as briefly outlining the generic detector and WIMP parameters we used in generating our results. In section IV we present the results for the expected angular distributions in these model detectors. In section V we then proceed to describe the Monte Carlo analysis we performed to determine the required event rates in order to separate signals from backgrounds for the different halo and detector models. In section VI we summarize our conclusions, and describe ongoing efforts to further refine our analyses and predictions.

## II Differential Event Rates

The WIMP rate as a function of recoil energy for various WIMP models has been carefully studied over the past decade (see  for a review). The angular dependence of the event rate, assuming an incident spherically symmetric isothermal halo WIMP distribution has also been discussed. . In order to explore a more general class of halo distributions, in particular those which might not be spherically these formalisms are not adequate, however, and one must generalize them . We present the derivation of such a generalized formalism here.

The event rate of WIMPs depends on their local density, $`\rho _0`$, and their velocity distribution in the halo, $`f(𝐯)`$. Assume the WIMPS have a uniform spatial density on the scale of the Earth’s motion through the halo over the course of one year. The event rate is then simply given by

$$R\frac{n\sigma }{m_n}v,$$
(1)
where $`m_n`$ is the mass of the nuclear target, $`m_\chi `$ is the mass of the WIMP, $`n=\rho _0/m_\chi `$ is the number density of WIMPs, $`\sigma `$ is the cross section, and $`v`$ is the average velocity of the WIMPs relative to the detector.

Consider a WIMP of mass $`m_\chi `$ moving with velocity $`𝐯=v(\widehat{x}\mathrm{sin}\alpha \mathrm{cos}\beta +\widehat{y}\mathrm{sin}\alpha \mathrm{sin}\beta +\widehat{z}\mathrm{cos}\alpha )`$ in the laboratory frame (figure 1). Let this WIMP scatter off a nucleus of mass $`m_n`$ in the direction ($`\theta ^{}`$, $`\xi `$) in the center-of-mass frame. The recoil energy of the nucleus is

$$Q=\frac{m_nm_\chi ^2}{\left(m_n+m_\chi \right)^2}v^2(1\mu ),$$
(2)
where $`\mu =\mathrm{cos}\theta ^{}`$. Let $`𝐮=u(\widehat{x}\mathrm{sin}\gamma \mathrm{cos}\varphi +\widehat{y}\mathrm{sin}\gamma \mathrm{sin}\varphi +\widehat{z}\mathrm{cos}\gamma )`$ be the recoil velocity of the nucleus in the laboratory frame (figure 1). To determine the relations among these sets of angles we begin by considering the simplified scattering problem of a WIMP incident along the $`z`$-axis. We can then rotate this result for an arbitrary WIMP incident at an angle $`(\alpha ,\beta )`$ to find (in the non-relativistic limit)

$`\mathrm{cos}\gamma `$ $`=`$ $`\sqrt{{\displaystyle \frac{1\mu }{2}}}\mathrm{cos}\alpha \sqrt{{\displaystyle \frac{1+\mu }{2}}}\mathrm{sin}\alpha \mathrm{cos}\xi ,`$ (3)
$`\mathrm{sin}\gamma \mathrm{cos}\varphi `$ $`=`$ $`\sqrt{{\displaystyle \frac{1\mu }{2}}}\mathrm{sin}\alpha \mathrm{cos}\beta \sqrt{{\displaystyle \frac{1+\mu }{2}}}\left(\mathrm{sin}\beta \mathrm{sin}\xi \mathrm{cos}\alpha \mathrm{cos}\beta \mathrm{cos}\xi \right),`$ (4)
$`\mathrm{sin}\gamma \mathrm{sin}\varphi `$ $`=`$ $`\sqrt{{\displaystyle \frac{1\mu }{2}}}\mathrm{sin}\alpha \mathrm{sin}\beta +\sqrt{{\displaystyle \frac{1+\mu }{2}}}\left(\mathrm{cos}\beta \mathrm{sin}\xi +\mathrm{cos}\alpha \mathrm{sin}\beta \mathrm{cos}\xi \right).`$ (5)
Notice in the limit of $`(\alpha ,\beta )(0,0)`$ this reduces to the usual two dimensional result.

To remove the dependence on the center of mass angles we note that

$$\mathrm{sin}\frac{\theta ^{}}{2}=\sqrt{\frac{1\mu }{2}}=\left|\mathrm{cos}\gamma \mathrm{cos}\alpha +\mathrm{sin}\gamma \mathrm{sin}\alpha \mathrm{cos}(\beta \varphi )\right|$$
(6)
and the Jacobian

$$\frac{(\mathrm{cos}\gamma ,\varphi )}{(\mu ,\xi )}=\frac{1}{4}\sqrt{\frac{2}{1\mu }}.$$
(7)
Next, to simplify subsequent formulae, we define a geometrical factor, $`J`$, by

$`J(\alpha ,\beta ;\gamma ,\varphi )`$ $``$ $`{\displaystyle \frac{𝐮𝐯}{|𝐮||𝐯|}}={\displaystyle \frac{1}{4}}\left|{\displaystyle \frac{(\mathrm{cos}\gamma ,\varphi )}{(\mu ,\xi )}}\right|^1`$ (8)
$`=`$ $`\mathrm{cos}\alpha \mathrm{cos}\gamma +\mathrm{sin}\alpha \mathrm{sin}\gamma \mathrm{cos}(\beta \varphi ).`$ (9)
With this information we can now write the event rates solely as a function of the incoming lab angle ($`\alpha `$, $`\beta `$) and the recoiling lab angle ($`\gamma `$, $`\varphi `$).

The event rate per nucleon in the detector is given by

$$dR=f(v,\alpha ,\beta )v^3dvd(\mathrm{cos}\alpha )d\beta \frac{d\sigma }{d\mu }d\mu \frac{d\xi }{2\pi }.$$
(10)
Here we will assume the nucleon-WIMP scattering has the form

$$\frac{d\sigma }{d\mu }=\frac{\sigma _0}{2}F^2(Q),$$
(11)
where $`F^2(Q)`$ is the form factor and $`\sigma _0`$ is energy independent. The number of events in which a nucleus recoils with energy $`Q`$ is

$$\frac{dR}{dQ}=\frac{\sigma _0\rho _0}{2m_r^2m_n}F^2(Q)_{v_{\mathrm{min}}}^{v_{\mathrm{esc}}}v𝑑v𝑑\mathrm{\Omega }_{\alpha ,\beta }f(v,\alpha ,\beta ),$$
(12)
where $`v_{\mathrm{esc}}`$ is the escape velocity of the Galaxy, $`v_{\mathrm{min}}^2=\frac{\left(m_\chi +m_n\right)^2}{2m_\chi ^2m_n}Q`$ is the minimum incident WIMP velocity that can produce a nuclear recoil of energy $`Q`$, and $`m_r=m_\chi m_n/(m_\chi +m_n)`$ is the reduced mass. This formula is well known and, in the case of an isothermal model when we only consider the projected motion of the Earth along the Sun’s velocity, it can be evaluated analytically.

In this work we are most interested in the angular dependent event rate which is found to be

$$\frac{dR}{d\mathrm{\Omega }_{\gamma ,\varphi }}=\frac{\sigma _0\rho _0}{\pi m_nm_\chi }_{v_{\mathrm{min}}/J}^{v_{\mathrm{esc}}}v^3𝑑vF^2\left(Q(v,J)\right)𝑑\mathrm{\Omega }_{\alpha ,\beta }f(v,\alpha ,\beta )J(\alpha ,\beta ;\gamma ,\varphi )\mathrm{\Theta }(J).$$
(13)
Here $`\mathrm{\Theta }(J)`$ is the usual step function and $`J`$ is the geometrical factor discussed above (9) that relates the incident WIMP direction to the recoiled nucleus direction.

We can further derive the event rate as a function of both the deposited energy and the recoil angle,

$$\frac{d^2R}{dQd\mathrm{\Omega }_{\gamma ,\varphi }}=\frac{m_n\sigma _0\rho _0}{8\pi m_r^4m_\chi }QF^2(Q)𝑑\mathrm{\Omega }_{\alpha ,\beta }f(v(Q,J),\alpha ,\beta )\frac{\mathrm{\Theta }(J)}{J^3}.$$
(14)
In principle this would provide the most information allowing for the best separation of signal from background. Note, however, that in order to fully exploit this signature some detector-specific estimation of the energy dependence of various backgrounds must be given.

In this work we will focus on the angular event rate, which does not require a knowledge of the specific energy dependent detector backgrounds in order to explore the separation of signals from backgrounds. Exploiting the full distribution for various proposed detectors will be analyzed in future work.

## III Model Parameters

Armed with the above differential event rates we can now analyze a wide variety of halo models, form factors, WIMP parameters, detector parameters, etc. Here we will discuss the parameters considered in this study.

### A Halo Models

A standard spherically symmetric isothermal sphere reproduces the large scale features of the flat rotation curves around galaxies. However, numerous dynamical arguments suggest that actual halo model may not be well described by such a distribution. As a result, a variety of more complex analytical models have been developed. The models we discuss here can be considered to be of three types: smooth axisymmetric, smooth triaxial, and clumped. The first type includes a standard spherical isothermal WIMP halo, as well as a modified axisymmetric halo, both rotating and non-rotating. Most recently a general triaxial generalization of the spherical isothermal halo has been developed. As an example of the last type of halo model, which ultimately should include realistic halo distributions arising from accretion of sub-systems by the growing galaxy, we consider a caustic model recently developed by Sikivie. While all of these models are analytic approximations, they do represent a broad spectrum of different possibilities, which should give some idea of the range of uncertainty in the predicted angular distributions in detectors.

#### 1 Isothermal Model

A simple spherically symmetric isothermal distribution is given by,

$$f(𝐯)=\frac{1}{\pi ^{3/2}v_0^3}e^{\left|𝐯\right|^2/v_0^2}.$$
(15)
Here $`v^2^{1/2}=\sqrt{3/2}v_0`$ is the dispersion velocity of the dark matter in the halo. Note that we can allow different dispersions in each direction but will not pursue that option further. The standard value chosen for $`v_0`$ is $`220\mathrm{km}/\mathrm{s}`$. We will consider the range of values $`v_0=150\mathrm{km}/\mathrm{s}`$, $`220\mathrm{km}/\mathrm{s}`$, and $`300\mathrm{km}/\mathrm{s}`$ here.

#### 2 Axisymmetric Halo Model

The Evans model  is an axisymmetric halo model that allows for flattening. This model has been studied in the context of annular modulations by Kamionkowski and Kinkhabwala . The distribution function is given by

$$f(𝐯)=\frac{1}{D^2}\left[AR_0\left(v\mathrm{cos}\alpha v_{,x}\right)^2+B\right]e^{2(v\mathrm{\Psi })^2/v_0^2}e^{2v_{}^2\mathrm{\Psi }^2/v_0^2}+\frac{C}{D}e^{(v\mathrm{\Psi })^2/v_0^2}e^{v_{}^2\mathrm{\Psi }^2/v_0^2},$$
(16)
where

$$A=\left(\frac{2}{\pi }\right)^{5/2}\frac{1q^2}{Gq^2v_0^3},B=\sqrt{\frac{2}{\pi ^5}}\frac{R_c^2}{Gq^2v_0},C=\frac{2q^21}{4\pi Gq^2v_0},$$

$$D=R_c^2+R_0^2,\mathrm{\Psi }=v_{,x}\mathrm{cos}\alpha +v_{,y}\mathrm{sin}\alpha \mathrm{cos}\beta +v_{,z}\mathrm{sin}\alpha \mathrm{sin}\beta ,$$
(17)
$`R_c`$ is the core radius, $`R_0`$ is the distance of our solar system from the center of the Galaxy, $`v_0`$ is the circular speed at large radii, and $`q`$ is the flattening parameter ranging from $`q=1`$ to $`q=1/\sqrt{2}`$. In this work we adopt $`R_c=7\mathrm{kpc}`$, $`R_0=8.5\mathrm{kpc}`$, and $`v_0=220\mathrm{km}/\mathrm{s}`$. The most important parameter is the flattening, $`q`$. Here we choose $`q=1`$ (cored isothermal model), $`0.85`$, and $`1/\sqrt{2}`$ (maximum flattening).

#### 3 Triaxial Halo Models

Constructing velocity distributions from general triaxial halos is a difficult problem. Evans, Carollo, and de Zeeuw  have provided such distributions for the logarithmic ellipsoidal potential; the simplest triaxial model and a natural generalization of the isothermal sphere. In this case we can approximate the velocity distribution as

$$f(𝐯)=\frac{1}{\pi ^{3/2}\sigma _x\sigma _y\sigma _z}\mathrm{exp}\left[\frac{v_x^2}{\sigma _x^2}\frac{v_y^2}{\sigma _y^2}\frac{v_z^2}{\sigma _z^2}\right].$$
(18)
When the Sun is on the long axis of the ellipsoid we have the relations

$$\sigma _x^2=\frac{v_0^2}{(2+\gamma )(p^2+q^21)},\sigma _y^2=\frac{v_0^2(2p^21)}{2(p^2+q^21)},\sigma _z^2=\frac{v_0^2(2q^21)}{2(p^2+q^21)}.$$
(19)
When the Sun is on the intermediate axis of the ellipsoid we have the relations

$$\sigma _x^2=\frac{v_0^2p^4}{(2+\gamma )(1+q^2p^2)},\sigma _y^2=\frac{v_0^2(2p^2)}{2(1+q^2p^2)},\sigma _z^2=\frac{v_0^2(2q^2p^2)}{2(1+q^2p^2)}.$$
(20)
Here $`p`$ and $`q`$ describe the axis ratios of the ellipsoid ($`p=q=1`$ is a spheroid) and $`\gamma `$ describes the anisotropy ($`\gamma =1`$ for a sphere). We will consider the cases where $`p=0.9`$ and $`q=0.8`$ for both positions of the Sun. We will also consider $`\gamma =1.78`$ (radially anisotropic) and $`\gamma =16`$ (tangentially anisotropic).

#### 4 Galactic Infall—Sikivie Caustic Model

The caustic model is based on the work of Sikivie and collaborators . It is derived from the assumption that WIMPs continuously fall into our Galaxy from all directions. The WIMPs only get thermalized after many passes through the Galaxy. Thus the WIMP halo of our Galaxy should be made up of an isothermal (or similar) core plus a set of inflowing and outgoing peaks in velocity space with the distribution function given by

$$f(𝐯)=\underset{j}{}\rho _j\delta (𝐯𝐯_j).$$
(21)
Here $`j`$ sums over the velocity flows, $`\rho _j`$ is the density in each flow, and $`𝐯_j`$ is the velocity of each flow. The peaks for one such model are given in table I. Notice that there are two flows for each velocity peak. In this model the total local density of WIMPs is $`\rho _0=0.52\mathrm{GeV}/\mathrm{cm}^3`$ and the first 14 peaks are not thermalized (Sikivie, private communications). From table I we find $`\rho _{\mathrm{caustic}}=0.34\mathrm{GeV}/\mathrm{cm}^3`$. Thus 65% of the WIMPs are in the form of velocity flows and the remaining 35% are in the form of a thermalized distribution. For simplicity we will choose the standard ($`v_0=220\mathrm{km}/\mathrm{s}`$) isothermal model for the thermalized distribution.

Most recently a number of groups have begun to numerically investigate the growth of our Galactic halo via accretion in N-body simulations. While we plan to incorporate these results into a future analysis, it is worth stressing that while the Sikivie model is undoubtedly an idealization of the actual hierarchical accretion process, the resulting distribution, containing a thermalized background plus several different unthermalized velocity flows may share some qualitative features of the N-body distribution, even if it is not likely to agree in detail with these results. It thus provides one tractable example of the complexity that might actually characterize the real galactic halo, and can usefully demonstrate some of the possible implications of phase space flows for signatures in WIMP detectors.

### B Motion of the Earth

In previous analyses it was generally assumed that the Earth orbits the Sun in the same plane that the Sun orbits the center of the Galaxy. In fact this is not true, the Earth’s orbit is tilted by approximately $`60^{}`$ from this plane. When working with an isothermal model and considering either the total number of events or the annular modulation, the effect of the tilt of the Earth’s orbit is insignificant. However, since we are interested here in angular effects and in halo models that need not be spherically symmetric, it is necessary to use the correct velocity of the Earth. In the frame of the Galaxy this is

$`v_{,x}`$ $`=`$ $`0.13v_{}\mathrm{sin}\left(2\pi (tt_p)/\mathrm{year}\right),`$ (22)
$`v_{,y}`$ $`=`$ $`0.11v_{}\mathrm{cos}\left(2\pi (tt_p)/\mathrm{year}\right),`$ (23)
$`v_{,z}`$ $`=`$ $`v_{}\left[1.05+0.06\mathrm{cos}\left(2\pi (tt_p)/\mathrm{year}\right)\right],`$ (24)
where $`v_{}=220\mathrm{km}/\mathrm{s}`$ is the velocity of the Sun in the frame of the Galaxy, $`t`$ is measured in days, and $`t_p153\mathrm{days}`$ (June 2) is the day when the direction of the Earth’s motion around the Sun matches the direction of the Sun’s motion around the Galactic center. The numerical coefficients reflect the $`60^{}`$ tilt of the Earth’s orbit.

The distribution functions above are given in the rest frame of the Galaxy. In the lab frame (rest frame of the Earth) we must shift them as $`f(𝐯+𝐯_{})`$.

### C Form Factor

In WIMP-nucleus scattering, a form factor is incorporated into the cross section, as described above, to account for the loss of coherence over the nucleus for large momentum transfers. The appropriate specific form factor is, of course, nucleus-dependent. However, for the purposes of this analyses we utilize several reasonable analytical approximations which are appropriate in the case of either spin-dependent, or spin-independent scattering, as described below. For a more detailed discussion see .

#### 1 Spin Independent

Here the WIMP couples to various quantum numbers of the entire nucleus. The standard form factor, $`F^2(Q)`$, for such nuclear interactions is the Wood-Saxon form factor. In this work we will consider the simpler exponential form factor

$$F^2(Q)=e^{Q/Q_0}$$
(25)
since it is easier to work with analytically and it produces qualitatively similar results. Here $`Q_0=3/(2m_nr_0^2)`$ and $`r_0=0.3+0.91\sqrt[3]{m_n}`$ is the radius of the nucleus (in femtometers when $`m_n`$ is in GeV).

#### 2 Spin Dependent

The energy dependent component of the form factor for an axial vector interaction is given by

$$F^2(q)S(q)=a_0^2S_{00}(q)+a_1^2S_{11}(q)+a_0a_1S_{01}(q).$$
(26)
Here $`a_0`$ ($`a_1`$) represents the isoscalar (isovector) parameterization of the matrix element. $`S_{00}`$, $`S_{11}`$, and $`S_{01}`$ are obtained from detailed nuclear calculations. In this work we will parameterize this energy dependent part using exponential functions for analytical ease.

For the purposes of our calculation we choose germanium. While it is not likely that this material will in fact be used in direction-sensitive detectors, it has the virtue of being a mid-range nucleus with sensitivity to axial couplings, and has been well studied. Without specific candidates for directional detector targets it thus seems a reasonable first step, in order to get some idea of the possible changes in signature for spin-dependent vs. spin independent interactions. For germanium, we find

$$S_{\mathrm{Ge}}(q)=0.20313\left(1.102a_0^2e^{7.468y}+a_1^2e^{8.856y}2.099a_0a_1e^{8.191y}\right),$$
(27)
where $`y=\left((1\mathrm{fm})|𝐪|/2A^{1/6}\right)^2`$ and is valid for $`y<0.2`$ . For a WIMP that is a pure bino the couplings are $`a_0=0.082\zeta _q`$ and $`a_1=0.276\zeta _q`$ when the European Muon Collaboration values for the spin content of the nucleons are employed. For our purposes in this work $`\zeta _q`$ is an energy independent, overall scale factor. See  for a more thorough discussion of $`S_{\mathrm{Ge}}`$.

### D Rotating Halo

To allow for a rotating halo we employ the standard technique from Lynden-Bell . Let

$$f_+(v,\alpha ,\beta )=\{\begin{array}{cc}f(v,\alpha ,\beta )\hfill & 0<\alpha <\pi /2\hfill \\ 0\hfill & \pi /2<\alpha <\pi \hfill \end{array}$$
(28)
and

$$f_{}(v,\alpha ,\beta )=\{\begin{array}{cc}0\hfill & 0<\alpha <\pi /2\hfill \\ f(v,\alpha ,\beta )\hfill & \pi /2<\alpha <\pi \hfill \end{array}.$$
(29)
Then the rotating distribution function is given by

$$f_R(v,\alpha ,\beta )=(1+\kappa )f_+(v,\alpha ,\beta )+(1\kappa )f_{}(v,\alpha ,\beta ).$$
(30)
Note that this definition does not change the normalization nor dispersion of the distribution function; $`<f_R>=<f>`$ and $`<v^2f_R>=<v^2f>`$, independent of $`\kappa `$. When $`\kappa =0`$ we recover the original distribution function.

To generate a halo with a particular average velocity we note that

$$vf_R=\frac{2v_0}{\sqrt{\pi }}\kappa .$$
(31)
Thus when $`vf_R=v_0`$ we find that $`\kappa =\sqrt{\pi }/2`$.

### E Target Nucleus

As mentioned above we chose $`{}_{}{}^{73}\mathrm{Ge}`$ as our prototypical nuclear target. Germanium has both scalar and axial vector interactions with WIMPs, and thus serves as a good baseline for probing the dependencies in the angular distribution. Throughout, we assume that 25% of the recoil energy is in the detected channel (quenching factor of $`0.25`$). The detector thresholds quoted below have taken this into account. Both argon and xenon have been discusses as likely targets in a TPC such as DRIFT . There are experimental difficulties with using argon (it is naturally radioactive) and little is known about xenon gas (in particular its form factor). Thus we will not examine these targets in detail at this time.

Note that for our purposes the choice of $`{}_{}{}^{73}\mathrm{Ge}`$ as target is somewhat arbitrary. Since the normalization of the overall WIMP-nucleus scattering cross section is not relevant for our discussions, the chief distinction between choosing different target nuclei will be the kinematic dependence on WIMP mass vs target mass, and possible alterations in the spin-dependent form factor parameters. For most materials being used for dark matter detection, quenching factors generally lie in the $`0.10.4`$ range and thus our choice of the measured Ge quenching factor is also relatively generic.

### F WIMP parameters

Throughout this work we focus on the shape of the nuclear recoil distribution. Thus, as emphasized above, we do not focus on the overall event rate normalization. This quantity depends on specifics of the halo density, the WIMP mass and couplings, the target nucleus, detector size, type of interaction, etc. The key point is that when modeling a specific detector this normalization will be essential, but in order to explore generic features of the angular signature of WIMP-nucleus scattering, it is not important.

As noted above, the mass of the WIMP affects both the normalization and shape of the nuclear recoil distribution. Thus, the results remain identical if one appropriate scales both target nucleus and WIMP mass. In this work we will consider $`m_\chi =60\mathrm{GeV}`$ as our standard WIMP mass (as mass matches well with the germanium target, $`m_{\mathrm{Ge}}=73\mathrm{GeV}`$). In some cases we will also show results for $`m_\chi =180\mathrm{GeV}`$ to show how the mass of the WIMP affects the results. The results we present should thus be considered to be appropriate to a WIMP whose mass matches well with the target nucleus mass, or in some cases greatly exceeds it.

## IV Results

The calculations of $`dR/d\mathrm{\Omega }`$ were performed on a $`40\times 40`$ grid in the $`(\mathrm{cos}\gamma ,\varphi )`$ plane. They were calculated in $`5`$ day bins and summed over the year. By summing over the year we average out any annual modulations in both the signal and the background. Furthermore, we are implicitly averaging over the day in which time the detector rotates (along with the Earth) by $`2\pi `$ with respect to the direction of motion. Thus any angular dependent backgrounds, due to, for example, hot spots in the detector, will also be averaged out (at least for the component in the plane of motion of the Earth) and the probability of a background induced nuclear recoil is uniform in angle when averaged over the year, with respect to the direction of the Earth’s motion around the Sun.

### A Angular Distributions

The angular distribution for the isothermal model with $`v_0=220\mathrm{km}/\mathrm{s}`$ is shown in figure 2. Notice the exponential decay in $`\mathrm{cos}\gamma `$ in the forward direction and the independence of $`\varphi `$. The reason for the decay in the forward direction is simple. The Earth is moving through this WIMP halo, and thus more WIMPs are incident from the forward direction than the backward direction, and thus nuclei are preferentially scattered backwards. This exponential decay would lead to an easy statistical separation of signal events from (flat) background events as discussed below. To examine dependence on threshold we plot $`dR/d\mathrm{\Omega }`$ for $`\varphi =0`$ for a number of different threshold energies (figure 3). As we increase the threshold the backward scattering peak in the distribution becomes more pronounced in relation to the rest of the distribution, making it more easily distinguishable from a (flat) background. However, as is also evident from the figure, as we increase the threshold the total number of events decreases. These two effects will be discussed below.

The angular distributions for a number of isothermal, axisymmetric and rotating models are shown in figure 4. Here we plot $`dR/d\mathrm{\Omega }`$ as a function of $`\mathrm{cos}\gamma `$ for $`\varphi =0`$ and $`E_{\mathrm{th}}=0\mathrm{keV}`$. The two solid curves for the $`v_0=220\mathrm{km}/\mathrm{s}`$ isothermal model are for two different WIMP masses, $`m_\chi =60\mathrm{GeV}`$ and $`m_\chi =180\mathrm{GeV}`$. All of the curves are normalized so that $`R=1`$. As can be seen from the figure all of the models are similar in shape (at least at $`E_{\mathrm{th}}=0\mathrm{keV}`$). Thus we do not expect a large variation in the number of events required to identify a WIMP signal as we will discuss later. On the other hand, it also means that it will be difficult to distinguish between different halo models based solely on this information.

The angular distribution for the triaxial model is shown in figure 5. Here there is a panel for each set of parameters. This shows how much the rate can vary even in a single model depending on the parameters. The most difficult combination of parameters to distinguish from background is seen in panel (a) since it looks the most similar to a flat distribution. The rest of the models look more like the isothermal distribution (panels b–d) and will be easier to distinguish from the background.

The angular distribution for the caustic model is shown in figure 6. Here both the pure caustic and a caustic with isothermal component are shown. Notice that the pure caustic model is peaked in the forward direction (opposite the other models) and falls off more slowly. This is because we are located near a caustic and WIMPS are preferentially catching up with the Sun as it moves through the Galaxy, so more are incident from behind than from the forward direction. There are also a number of small features in the $`\varphi `$ direction. When combined with the $`v_0=220\mathrm{km}/\mathrm{s}`$ isothermal model the distribution flattens greatly. This makes is more difficult to distinguish from a flat background than the other models.

### B Signal Identification

To determine the number of signal events necessary to identify a WIMP signal we employ a maximum likelihood analysis along with Monte Carlo generation of sample nuclear recoil distributions. We define a likelihood function

$$=\underset{i=1}{\overset{N_e}{}}P(\mathrm{cos}\gamma _i,\varphi _i),$$
(32)
where $`N_e`$ is the total number of events and $`P(\mathrm{cos}\gamma _i,\varphi _i)`$ is the probability of nuclear recoil in the $`(\gamma _i,\varphi _i)`$ direction for a particular model (e.g. an isothermal model). Here $`P(\mathrm{cos}\gamma _i,\varphi _i)`$ is generated by calculating $`dR/d\mathrm{\Omega }`$ for each recoil direction $`(\gamma _i,\varphi _i)`$ and averaging the result over the year in 5 day bins. At the 95% confidence level when $`\mathrm{log}_{dR/d\mathrm{\Omega }}\mathrm{log}_{\mathrm{flat}}<1`$ the distributions are indistinguishable. We generate 10,000 data sets for each $`N_e`$ and demand that the log-likelihood condition is satisfied less than 5% of the time. The smallest $`N_e`$ for which this occurs is the minimum number of events required to get a 95% detection 95% of the time.

With this procedure we can determine the number of events required to distinguish a signal from a flat background or to distinguish one signal (e.g. an isothermal halo) from another signal (e.g. a caustic halo). For a pure signal we use $`P_{dR/d\mathrm{\Omega }}`$ as determined by the shape of $`dR/d\mathrm{\Omega }`$. For a signal-to-noise ratio, $`S/N`$, we replace $`P`$ with $`P_{\mathrm{total}}=\lambda P_{\mathrm{flat}}+(1\lambda )P_{dR/d\mathrm{\Omega }}`$ where $`\lambda =1/(1+S/N)`$. For the results discussed here we consider a $`40\times 40`$ grid in the $`(\mathrm{cos}\gamma ,\varphi )`$ plane and assume perfect angular knowledge on this grid. This implies an accurate angular recoil detector. We also consider the case where knowledge of just the forward-backward asymmetry is available. In this case $`P`$ is replaced by a binomial distribution.

In figure 7 we show the number of signal events needed to differentiate a WIMP signal from a flat background at the 95% confidence level 95% of the time for a standard isothermal halo. In this figure there are two sets of curves shown. The lower set contains the results for a pure signal (no noise). The upper set contains the results for a signal to noise ratio of one, thus twice as many events as shown would have been detected (roughly half signal and half noise). In each set two curves are shown. The lower curve in each set shows the actual number of events required to distinguish the WIMP signal. For the S/N=1 case this decreases from 32 events at $`E_{\mathrm{th}}=0\mathrm{keV}`$ to 8 events at $`E_{\mathrm{th}}=25\mathrm{keV}`$. As mentioned above, this is due to the fact that the distribution becomes more peaked in the backward direction, thus it is more easily distinguished from a flat distribution. However, as also mentioned above, at a higher threshold energy there are fewer events. Thus, for a fixed-size detector raising the threshold will not in general aid in signal detection, because most of the signal events will be discarded. We attempt to account for this fact in the upper curve in each set. The upper curve is gives the number of events that would have been detected at zero threshold in order to observe the number of events shown in the lower curve for finite thresholds. This curve shows the effect of the decrease in the rate as the threshold is increased. For isothermal halos, at least, it is desirable to choose the lowest possible threshold of the detector for a fixed signal-to-noise ratio in order to improve the ability to detect WIMPs with the fewest number of events.

The effect of changing the width of the isothermal distribution (changing $`v_0`$) is shown in figure 8. As expected from figure 4 at low thresholds the narrower distribution ($`v_0=150\mathrm{km}/\mathrm{s}`$) is more easily distinguished from a flat background than the broader distributions. At higher thresholds the tail of the WIMP distribution becomes important since the number of high energy nuclear recoils falls off more rapidly for narrower distributions. Note that it is still the case that the number of events required for detection is an increasing function of energy.

The number of events required for the axisymmetric isothermal models are shown in figure 9. For reference the $`v_0=220\mathrm{km}/\mathrm{s}`$ isothermal model is also included (solid line). Both the $`q=1`$ and $`q=0.85`$ models are very similar to the $`v_0=220\mathrm{km}/\mathrm{s}`$ isothermal model. The $`q=1/\sqrt{2}`$ (maximal flattening) model is very similar to the $`v_0=150\mathrm{km}/\mathrm{s}`$ isothermal model.

In order to explore the general effects of spin-dependence, WIMP mass, and rotation, we display in figure 10 the predicted event distributions for a $`v_0=220\mathrm{km}/\mathrm{s}`$ isothermal model, while allowing these other parameters to vary. From the figure we see that the spin interaction using the form factor presented earlier leads to a steeper rise in the number of events. This is due to the steeper fall off of the spin-dependent form factor with energy compared to the spin independent form factor, and will thus be nucleus-dependent. A larger mass WIMP leads to more energy transfered to the nucleus, thus the $`m_\chi =180\mathrm{GeV}`$ WIMP requires fewer events at high threshold than the lighter $`m_\chi =60\mathrm{GeV}`$ WIMP does. A (net) corotating halo is more difficult to distinguish than the isothermal sphere as expected. However, it is not as difficult as might be expected since we have cut out some, but not all, of the WIMPs from the forward direction. This shows up more strongly as the threshold is increased.

It is important to note that all of the above distributions require roughly the same number of events in order to distinguish a signal from a flat background. If this were generic, then low-threshold detectors sensitive to roughly 30-50 events would be guaranteed to unambiguously identify a WIMP signal, even if a distinction between halo models might not be easy. Unfortunately, however, two of the models we have considered here, both of which exist in the literature as candidates for our Galactic halo, will require far more events in order to disentangle the signal from a flat background. For the reasons mentioned earlier, it is reasonable to expect that the distributions that arise from N-body simulations may share these features.

First, we consider the caustic distribution. As is expected from figure 6, this most closely resembles a flat distribution and thus requires the most number of events. Also not surprisingly it is a strong function of threshold energy. As the threshold is increased some of the velocity peaks do not contribute detectable nuclear recoils thus greatly changing the shape of the recoil distribution. We note that the number of events required for detection is not a monotonic function. Thus the detection of WIMPs in a caustic halo will depend sensitively on the mass of the WIMP and the threshold of the detector. A range of thresholds would need to be probed. In practice also having energy information about the event will most likely prove to be extremely important to differentiate WIMP events from background events in this case.

Next we consider the triaxial logarithmic ellipsoidal model. As expected from figure 5 most of the models are not too different from the standard isothermal sphere 11. However, when the Sun is on the intermediate axis and we have a radial anisotropy ($`\gamma =1.78`$) approximately 120–150 events are required. This model shows the importance of exploring a broad class of halos and parameters in order to assess the requirements of a a detector with angular resolution.

In the above results we have assumed fairly accurate angular resolution (bins of about 9 degrees in width were used). However a real detector may not be able to achieve this level of angular resolution. To explore this we consider the case of only having knowledge of the forward to backward asymmetry. In this case the events fall into one of two bins. In our likelihood analysis we compare two binomial distributions; one for the data with fraction of events in the forward bin as our probability and the other for the flat background where the probability is a half.

In figure 12 we show the number of signal events needed to differentiate a WIMP signal from a flat background at the 95% confidence level 95% of the time when only two angular bins are employed (forward and backward). This figure is the counterpart of figure 7. Even in this case of relatively little angular information few events are needed. For $`S/N=1`$ if $`E_{\mathrm{th}}<5\mathrm{keV}`$ then $`N_{\mathrm{events}}40\text{}65`$. It is interesting to note that thresholds slightly above zero give a somewhat more discernible signal in this case.

The results for the other models considered in this work are shown in figures 1316. These figures should be compared to figures 811. In most cases we see the slight dip (or at least very slow rise) in the number of events needed around $`E_{\mathrm{th}}2\mathrm{keV}`$. We note that the result for the caustic model is not shown. As previously mentioned the caustic model most closely resembles a flat distribution. This shows up most prominently here where we have limited the angular information available. In this case the number of events needed with $`S/N=1`$ at $`E_{\mathrm{th}}=0\mathrm{keV}`$ is about 5200. Again most of the triaxial models are not that different than the isothermal sphere. However the radially anisotropic model with the Sun on the intermediate axis again requires the most number of events. Here approximately 250–300 events are required for $`E_{\mathrm{th}}<5\mathrm{keV}`$.

The DRIFT detector in its current form will only obtain information about the recoil in two directions. One will be along the axis of the Sun’s motion through the Galaxy ($`z`$-axis in our coordinates) and the other either the $`x`$ or $`y`$-axis. Having information along the $`z`$-axis is necessary to have any ability to identify WIMPs. Since the recoil distribution is essentially independent of $`\varphi `$ (or equivalently the $`x`$, $`y`$-plane, see figure 7) having information on these two axes is almost the same as having the full angular distribution. We have calculated the number of events required to identify a WIMP signature in this case and the results are identical to those discussed above for the isothermal and axisymmetric models. Naturally in the case of the caustic model more events are required since this model has structure in the $`x`$, $`y`$-plane.

This simple study of the DRIFT’s capabilities assumed that we could identify the recoil tracks regardless of their orientation. In practice if we do not have full angular resolution then some tracks will not be resolved. For example, if we are projecting out the $`x`$-axis any event that has a large component along the $`x`$-axis and a small component along the $`y`$-axis will appear as a short track and will thus be difficult to identify. To accurately model DRIFT it would be necessary to take this track identification into account along with other detector characteristics (threshold, etc.) Furthermore it may not be possible to measure the direction of the track, only its axis. In a detector such as DRIFT that relies on ionization to measure the recoil direction, it is difficult to tell where a track begins and ends. To make this distinction we must rely on the ionization (energy deposition) as a function of the energy of the recoiling nucleus. The extent to which this can be used to determine directionality will be crucial for determining the ability of a detector to identify a WIMP signal. We will consider these details in future work.

## V Conclusions

Our results are useful when considering the construction of the next generation of WIMP detectors sensitive to the angular distribution of WIMP scattering events. As we have described, this sensitivity will provide, in principle, a key extra signature needed to unambiguously demonstrate that an observed signal is indeed WIMP-related, and one which will require both far fewer events and also be subject to fewer systematic uncertainties than exist in a search for simple annual modulations.

The formalism we have developed allows one to calculate the angular distribution of recoil events, $`dR/d\mathrm{\Omega }`$, and the distribution in both energy and angle, $`d^2R/dQd\mathrm{\Omega }`$. This formalism accommodates arbitrary velocity space WIMP distributions and correctly employs the full motion of the Earth around the Sun and the Sun through the Galaxy. As such, it can be used as an input to any detector Monte Carlo, allowing a determination of the predicted angular distribution of events for any halo model.

While we hope to utilize our analysis to examine various specific detectors possibilities in the future, we have thus far derived various generic predictions for the range of events that will be required in an idealized detector, given our present uncertainty in the actual Galactic halo WIMP distribution. We feel this is an important first step in order to set the scale of detector parameters that may be required to adequate exploit angular sensitivity.

In particular, we have calculated the recoil distribution from WIMP scattering and the expected number of events necessary to distinguish a WIMP signal from a background with a 95% confidence level for various analytic halo models; isothermal models with $`v_0=150\mathrm{km}/\mathrm{s}`$, $`220\mathrm{km}/\mathrm{s}`$, and $`300\mathrm{km}/\mathrm{s}`$; Evans models with flattening $`q=1`$, $`0.85`$, and $`1/\sqrt{2}`$, and rotating halo models. In all such cases the number of events required is quite low, typically $`50\text{}70`$ for $`E_{\mathrm{th}}<5\mathrm{keV}`$ when the full angular distribution is measured. Even if only the forward to backward ratio is measured relatively few events are required, typically $`70\text{}110`$ for $`E_{\mathrm{th}}<5\mathrm{keV}`$ for these models. This study considered WIMPs with masses $`m_\chi =60\mathrm{GeV}`$ and $`180\mathrm{GeV}`$ and included both spin dependent and independent interactions. We note that spin-dependent interactions may require lower thresholds in order to have sufficient number of events to distinguish a WIMP signal.

The isothermal and Evans models are both axisymmetric. As extreme asymmetric models we also studied the case of caustics (peaks) in the velocity distribution, and a set of triaxial halo models. The former distribution should be considered a first pass at considering more realistic halo distributions arising from hierarchical accretion in N-body simulations. In this case the resulting recoil spectrum is very similar to a flat background and thus is rather difficult to distinguish. Here the number of events required ranges from $`300\text{}600`$ for $`E_{\mathrm{th}}=2\text{}10\mathrm{keV}`$ and depends rather sensitively on the threshold. For the triaxial model most sets of parameters lead to limits similar to the isothermal sphere. Only the case of the Sun on the intermediate axis with a radial anisotropy ($`\gamma =1.78`$) are the requirements raised to 120–150 events for $`E_{\mathrm{th}}<5\mathrm{keV}`$.

We have also performed a preliminary study of a DRIFT-like detector. Here we projected out one of the directions perpendicular to the direction of motion of the Sun since DRIFT will only have resolution in two directions. Here we found that the results can be as robust as having full angular resolution for axisymmetric halo models.

In future work we plan to refine our analysis in several ways. First, we will exploit the output of existing N-body simulations to explore the consequences for the expected angular signals. Next, we plan to carry out an analysis similar to the present analysis for specific proposed detectors. This will allow us to explore the very important issue of utilizing energy as well as angular sensitivity, but it will also require us to adequately model the expected energy dependence of various detector backgrounds. Finally, after completing these analyses, we hope to explore exactly how many events may be needed to distinguish between halo models, in order to determine how an eventual WIMP detection may shed light on the unknown features of our Galactic halo, and thus on important issues associated with galaxy formation and evolution.

If a WIMP signal is ultimately observed in terrestrial detectors it will provide one of the most important observations that has ever been made in particle physics and cosmology. It is thus worthwhile considering in advance not only how one might best unambiguously disentangle such a signal from backgrounds, but also how one might then use this signal as a probe of astrophysics. Our results provide a useful first step in this direction.

We thank P. Sikivie, J. Martoff, and B. Moore for useful discussions about their research results, and J. Heo for his contributions to our initial investigations.