text:
# What’s Wrong with these Observables?

## Abstract

An imprecise measurement of a dynamical variable (such as a spin component) does not, in general, give the value of another dynamical variable (such as a spin component along a slightly different direction). The result of the measurement cannot be interpreted as the value of any observable that has a classical analogue.

1. MEMORABILIA

I was never compelled to learn solid state physics, and my first encounter with David Mermin was his paper  “Is the moon there when nobody looks? Reality and the quantum theory”. That paper gave a wonderful proof of Bell’s theorem. I wrote a comment  saying “Mermin gives a remarkably simple illustration of Bell’s theorem, but leaves the impression that something mysterious is implied. The situation is much simpler: the pair of photons is a single, nonlocal, indivisible entity … It is only because we force upon the photon pair the description of two separate particles that we get the paradox of Einstein, Podolsky and Rosen.” David answered to me, privately, “one person’s mystery is another person’s explanation.”

The following year, we met in a conference that Danny Greenberger had organized in the World Trade Center in New York City  and it was friendship at first sight. I gave a talk “When is a quantum measurement?” after which Abner Shimony came to me and said, with a big friendly smile “Asher, you understand nothing! You speak just like Niels Bohr!” One person’s insult is another person’s compliment.

Over the years, David and I collaborated many times, but I never had the honor of being his co-author. Many of our discussions revolved around the Kochen-Specker (KS) theorem . In 1990, David relayed to me information that Conway and Kochen had found a new proof of that theorem using only 33 rays, instead of 117 in the original proof. The only clue I had was that they used a cubic lattice for constructing their rays. This led me to find a much simpler proof, also with 33 rays, but this one using an irrational lattice. I wrote to Simon Kochen, asking him to inform me when he publishes his proof, so that I could publish mine immediately afterwards. Kochen answered that my proof was indeed very simple, but meanwhile he had found one with only 31 rays. I was both happy and disappointed. However, soon afterwards, I found a proof with only 24 rays, in four dimensions. I wrote to Kochen that $`24<31`$ and that I would publish my results, while mentioning the existence of their proof with 31 rays . Conway and Kochen then never bothered to publish the latter.

In my paper  I also stated that 24 was the minimal number number in four dimensions. David was not convinced by my argument. We exchanged several letters on this issue but could not resolve it. Meanwhile I was curious to see the 31 rays in three dimensions. I guessed they were a subset of the same cubic lattice and I did the search with a computer program. The algorithm appears in my book , where there is also an exercise: write a computer program for the four-dimensional case, and check that 24 rays are the minimum number. I never bothered to solve that boring exercise, but two students took up the challenge and found that it was possible to remove any one of the 24 rays, and still have a KS set. Michael Kernaghan, in Canada, found a KS set with 20 rays  and then Adán Cabello, together with José Manuel Estebaranz and Guillermo García Alcaine in Madrid, found a set of 18 rays . They still hold the world record (probably for ever).

Soon after I found my $`33+24`$ rays, there was the Gulf war and Saddam Hussein sent numerous scud missiles on Israeli targets. Each time there was a raid alert, we had to go into an air-tight room and don gas masks (fortunately they were never needed). Each alert lasted about half an hour, until the scud debris were examined by civil defense experts. To help time pass, I read Boojums  through the goggles of the mask. This was sometimes difficult, because the goggles were fogged by my breath, but this was always enjoyable and instructive.

It is a pleasure to dedicate this article to David Mermin, for his birthday and many more birthdays.

2. WHAT CAN BE OBSERVED?

Standard texbooks on quantum mechanics tell you that observable quantities are represented by Hermitian operators, their possible values are the eigenvalues of these operators, and that the probability of detecting eigenvalue $`\lambda _n`$, corresponding to eigenvector $`u_n`$, is $`|u_n|\psi |^2`$, where $`\psi `$ is the (pure) state of the quantum system that is observed. With a bit more sophistication to include mixed states, the probability can be written in a general way $`u_n|\rho |u_n`$. Really bad books also claim that the state of the physical system after the measurement collapses into the corresponding $`u_n`$. This is sheer nonsense. (Finding appropriate references is left as an exercise for the reader.)

The simple and obvious truth is that quantum phenomena do not occur in a Hilbert space. They occur in a laboratory. If you visit a real laboratory, you will never find there Hermitian operators. All you can see are emitters (lasers, ion guns, synchrotrons and the like) and detectors. The experimenter controls the emission process and observes detection events. The theorist’s problem is to predict the probability of response of this or that detector, for a given emission procedure. Quantum mechanics tells us that whatever comes from the emitter is represented by a state $`\rho `$ (a positive operator, usually normalized to 1). Detectors are represented by positive operators $`E_\mu `$, where $`\mu `$ is an arbitrary label whose sole role is to identify the detector. The probability that detector $`\mu `$ be excited is $`\mathrm{tr}(\rho E_\mu )`$. A complete set of $`E_\mu `$, including the possibility of no detection, sums up to the unit matrix and is called a positive operator valued measure (POVM) .

The various $`E_\mu `$ do not in general commute, and therefore a detection event does not correspond to what is commonly called the “measurement of an observable.” Still, the activation of a particular detector is a macroscopic, objective phenomenon. There is no uncertainty as to which detector actually clicked.

3. IMPRECISE MEASUREMENTS

There has recently been considerable controversy about the possibility of testing the physical implications of the KS theorem or even “nullifying” the latter . The “observables” that are usually considered are spin components of a spin-1 particle, which is the paradigm of the KS theorem. Opinions are varied  but more often than not assume that the result of an imprecise measurement of a spin component $`𝐧𝐉`$ is the value of another spin component, along some direction $`𝐧^{}`$ which is close to the correct n. As shown below, this assumption is generally unfounded (except in the trivial case of spin-$`\frac{1}{2}`$ particles).

When I first learnt of this “nullification,” my gut feeling was that the claim was not even wrong. However, finding a decisive argument that would convince the authors  proved as elusive as giving an experimental proof of a violation of Bell’s inequality. It is always possible to find loopholes. All I can do is to challenge them to find a loophole in the following argument.

Whatever the causes of imprecision are, if we don’t measure exactly $`𝐧𝐉`$, we have to give a physical interpretation to the result of our experimental procedure. Do we measure some other nearby operator,

$$A=𝐧𝐉+\underset{m,n}{}a_{mn}J_mJ_n/\mathrm{},$$

where the unknown dimensionless coefficients $`a_{mn}=a_{nm}^{}`$ are small? If the array $`a_{mn}`$ is antisymmetric, and therefore imaginary, we have $`A=𝐦𝐉`$, where $`𝐦𝐧`$ (but m need not be a unit vector). However, in general, the operator $`A`$ defined above also contains bilinear terms which cannot be reduced to another component of J. These terms have no classical analogue and no name in our vocabulary. Quantum mechanics has a much richer choice of dynamical variables than classical mechanics.

Actually, the result of the measurement procedure may not even be the value of any dynamical variable that looks like $`A`$. Traditional concepts such as “measuring Hermitian operators,” that were borrowed or adapted from classical physics, are not appropriate in the quantum world. In the latter, as explained above, we have emitters and detectors, and calculations are performed by means of POVMs.

Any lack of precision in the experimental setup simply means that the positive operator $`E_\mu `$ which is actually implemented in the laboratory only approximates the original $`E_\mu `$ that the experimenter intended to use. That $`E_\mu `$ may have been associated with a classical interpretation such as the value of a component of angular momentum. The actual $`E_\mu `$ may have no classical interpretation at all. In particular, there is no reason to expect it to correspond to a component of J along a slightly different direction.

What does this imply for attempts to test experimentally the KS theorem? That theorem is only a statement about Euclidean geometry. Any purported experimental test should be analyzed as explained above: particles are emitted according to a specified procedure, and quantum theory is used to predict the probability of excitation of various detectors. After we take into account known imperfections of the emission and detection processes, any discrepancy would imply either a poor understanding of the laboratory equipment, or a failure of quantum theory. The KS theorem itself is not involved. Its only role is to restrict our freedom of concocting realistic non-contextual theories that would replace quantum mechanics.

Acknowledgment

This work was supported by the Gerard Swope Fund and the Fund for Encouragement of Research.


format: text

text:
# 1 Introduction

## 1 Introduction

In general, to understand hadronic uncertainties we need a separation of short $`pQ`$ and long $`p\mathrm{\Lambda }_{\mathrm{QCD}}`$ distance fluctuations. For processes in QCD with momentum transfers $`Q\mathrm{\Lambda }_{\mathrm{QCD}}`$, the short distance part is calculable in terms of Wilson coefficients or hard scattering functions. The long distance contributions can be arranged into universal non-perturbative matrix elements which can be extracted from data or calculated on the lattice. This process of separating short and long distance fluctuations is sometimes referred to as factorization.

In $`B`$ physics a proper understanding of hadronic uncertainties is crucial, partially since the size of typical expansion parameters $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$, $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_c`$, or $`\sqrt{\mathrm{\Lambda }_{\mathrm{QCD}}/m_b}`$ leaves room for power corrections to play a non-negligible role. For many processes Heavy Quark Effective Theory (HQET) provides the conceptual framework for quantifying the factorization between the scales $`m_{b,c}`$ and $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. A well known example is the extraction of the CKM matrix element $`|V_{cb}|`$ from knowledge about the form factors in $`BD^{}\mathrm{}\nu `$ decays. At zero recoil the leading result is fixed by Heavy Quark Symmetry , and the first power corrections vanish . More recently important progress has been made at reducing the dominant model dependence by computing the matrix elements of $`1/m_Q^2`$ operators on the Lattice . For inclusive $`BX_c\mathrm{}\nu `$ and $`BX_u\mathrm{}\nu `$ decays the non-perturbative HQET matrix elements $`\overline{\mathrm{\Lambda }}`$ and $`\lambda _1`$ are actively being extracted from experimental data .

Since the $`B`$ is so heavy, many of its decays produce energetic light hadrons. For these decays the energy of the hadron $`E_H`$ in the $`B`$ rest frame is an additional perturbative scale, and HQET alone does not separate the perturbative and non-perturbative information. Examples of such processes include the decays $`BD\pi `$, $`B\pi \pi `$, $`BK\pi `$, the large recoil region in $`B\pi \mathrm{}\nu `$, $`B\rho \mathrm{}\nu `$, $`BK^{}\gamma `$, and $`BK\mathrm{}^+\mathrm{}^{}`$, and the endpoint spectra of the inclusive decays $`BX_u\mathrm{}\nu `$ and $`BX_s\gamma `$. The nature of factorization in these decays shares features in common with many exclusive and inclusive hard scattering processes. Examples are $`\gamma ^{}\gamma \pi ^0`$ at large $`q^2`$ and the $`x1`$ endpoint region of deep inelastic scattering.

In this talk I discuss an effective field theory that has been developed for processes with energetic hadrons, which is referred to as the Soft-Collinear Effective Theory (SCET) . SCET can be used for both hard scattering processes  and B decays. This theory makes symmetries relevant in the large energy limit explicit at the level of the Lagrangian and operators (such as the reduction of spin structures, helicity constraints, and collinear gauge invariance). Furthermore, SCET has a transparent power counting in $`\lambda =\mathrm{\Lambda }_{\mathrm{QCD}}/Q`$ (or $`\lambda =\sqrt{\mathrm{\Lambda }_{\mathrm{QCD}}/Q}`$) so that power corrections can be investigated in a systematic way . This includes processes not amenable to an operator product expansion such as exclusive decays. The renormalization group improvement of operators in the effective theory sums single infrared logs, as well as double Sudakov logarithms when they appear . Finally, SCET allows proofs of factorization theorems to be simplified and carried out in a gauge invariant way.

## 2 Formalism

The factorization of scales in the effective theory is carried out by describing long distance fluctuations with $`p^2\stackrel{<}{}Q^2\lambda ^2`$ using effective theory fields, and those with $`p^2Q^2\lambda ^2`$ by computable short distance Wilson coefficients. Typical processes require collinear fields and in addition either soft or usoft fields. Examples are $`BX_s\gamma `$ which needs collinear and usoft fields, and $`BD\pi `$ which needs collinear and soft fields. This field content is summarized in Table 1 together with the scaling of the momenta and fields with the expansion parameter $`\lambda `$ . The momenta scales $`Q`$, $`Q\lambda `$, and $`Q\lambda ^2`$ are separated by making phase redefinitions to pull out the larger momenta, $`\varphi _n(x)=_pe^{ipx}\varphi _{n,p}(x)`$. Derivatives on the new fields then always pick out the small scale, $`^\mu \varphi _{n,p}(x)(Q\lambda ^2)\varphi _{n,p}(x)`$, while the large momenta are picked out by introducing label operators, for example $`\overline{𝒫}\xi _{n,p}=(\overline{n}p)\xi _{n,p}`$. Since $`\overline{𝒫}\lambda ^0`$ in the power counting the hard coefficients $`C(\overline{𝒫})`$ are arbitrary functions of this operator , which can be determined by matching. More generally we have functions $`C(\omega _i)_i\delta (\omega _i\overline{𝒫})`$ where the delta functions are inserted inside collinear operators in the most general locations allowed by gauge and reparameterization invariance.

Furthermore, there are gluon fields which are order $`\lambda ^0`$ in the power counting, namely $`\overline{n}A_{n,q}1`$. Integrating out offshell fluctuations builds up Wilson lines in these fields, such as in the example  of matching the full theory heavy-to-light current $`\overline{u}\mathrm{\Gamma }b`$ onto the SCET current, $`J_0=C(\overline{𝒫})\overline{\xi }_{n,p}W\mathrm{\Gamma }h_v`$. In fact the $`\overline{n}A_{n,q}`$ field can be traded for the Wilson line

$`W`$ $`=`$ $`\left[{\displaystyle \underset{\mathrm{perms}}{}}\mathrm{exp}\left({\displaystyle \frac{g}{\overline{𝒫}}}\overline{n}A_{n,q}(x)\right)\right],`$ (1)
since the covariant derivative $`i\overline{n}D_c=\overline{𝒫}+g\overline{n}A_{n,q}=W\overline{𝒫}W^{}`$.

Soft Wilson lines $`S_n[nA_s]`$ are also built up by integrating out offshell fluctuations . Beyond tree level the structure of operators containing factors of $`W`$ or $`S_n`$ is protected by collinear and soft gauge transformations .

In general we have three types of gluon fields, collinear, soft, and usoft. These are the fields associated with gauge transformations $`U(x)`$ which have support over collinear, soft, and usoft momenta respectively . The usoft fields $`A_{us}^\mu `$ are dynamical quantum fields which due to their slow variation appear as background fields to the soft and collinear quarks and gluons. For a gauge transformation with support over collinear momenta it is convenient to factor out the large momentum components, $`U(x)=_Re^{iRx}𝒰_R(x)`$. The collinear gluon field then transforms such that

$`i𝒟^\mu +A_{n,p}^\mu 𝒰_QA_{n,R}^\mu 𝒰_{Q+Rp}^{}+{\displaystyle \frac{1}{g}}𝒰_Qi𝒟^\mu 𝒰_{Qp}^{},`$ (2)
where $`i𝒟^\mu =n^\mu \overline{𝒫}/2+𝒫_{}^\mu +\overline{n}^\mu inD/2`$. On the other hand an usoft gauge transformation has $`i^\mu U(x)(Q\lambda ^2)U(x)`$. For a fixed $`R`$ the transformation $`𝒰_R(x)`$ itself has an usoft component (having no large phase), so there is an overlap between what is meant by collinear and usoft gauge transformations. Effectively therefore one also has

$`A_{us}^\mu 𝒰_RA_{us}^\mu 𝒰_R^{}+𝒰_R{\displaystyle \frac{i}{g}}_\mu 𝒰_R^{}.`$ (3)
This transformation is usoft – it does not induce large momentum in the usoft gauge field because the large momenta labels in $`𝒰_R`$ and $`𝒰_R^{}`$ cancel. Gauge invariance constrains the form of the Lagrangian at leading and subleading orders. Beyond leading order the fact that the transformations can not be completely separated is crucial to correctly constrain operators as pointed out in Ref. . For instance, the collinear gauge invariance of the subleading term $`\overline{\xi }_{n,p}iD/_{us}^{}W(1/\overline{𝒫})W^{}iD/_c^{}\overline{n}/\xi _{n,p^{}}`$ relies on the transformation in Eq. (3). For the explicit form of the leading order gluon and quark actions $`_c^{(0)}`$ we refer the reader to Ref. , for higher order terms in the collinear quark action to Refs. , and for the mixed collinear-usoft quark action to Ref. .

Besides the constraints from gauge invariance on collinear operators there are addition constraints from the way in which Lorentz invariance is realized in the effective theory. The collinear fields are defined by introducing two auxillary light-like vectors, $`n`$ and $`\overline{n}`$, such that $`n\overline{n}=2`$. Naively the presence of these vectors breaks Lorentz invariance. However, in practice Lorentz invariance is restored order by order in the power counting by a reparameterization invariance (RPI) . For the collinear theory the study of RPI was initiated in Ref.  and generalized to the three most general classes of allowed transformations in Ref. . For Lagrangians and operators with collinear fields, RPI gives non-trivial constraints between the Wilson coefficients which appear different orders in the power expansion. In general there is no way of deducing these constraints using only the full theory.

Finally, it is worth discussing why the proof of factorization theorems is simplified by using the effective theory. The factorization between hard and collinear fluctuations or soft and collinear fluctuations is simplified by the fact that it takes place at the level of matching onto the effective theory. The resulting structures are constrained by the symmetries of the low energy theory as already discussed. The factorization between collinear and usoft interactions is simplified by the fact that many cancellations occur in a universal way at the level of the effective Lagrangian. For instance, at lowest order the actions for usoft and collinear particles can be factorized by a simple field redefinition on the collinear fields, $`\xi _{n,p}=Y_n\xi _{n,p}^{(0)}`$ and $`A_{n,p}=Y_nA_{n,p}^{(0)}Y_n^{}`$, where $`Y_n=P\mathrm{exp}[ig_{\mathrm{}}^x𝑑snA_{us}(sn)]`$. This transformation moves all leading order usoft interactions from the collinear Lagrangian $`_c^{(0)}`$ into the external operators and currents, after which cancellations due to the unitarity of the usoft Wilson line, $`Y_n^{}Y_n=1`$, are more readily seen .

## 3 Results

I briefly discuss two B-physics applications of SCET, the exclusive decay $`BD\pi `$ and the inclusive process $`BX_s\gamma `$.

For the decay $`BD\pi `$ the energy of the outgoing pion in the rest frame of the $`B`$ is $`E_\pi =2.31\mathrm{GeV}`$. Since this energy is large it is useful to consider this decay as being in the situation where $`Q\mathrm{\Lambda }_{\mathrm{QCD}}`$ with $`Q=m_b`$, $`m_c`$, or $`E_\pi `$. In this limit, SCET has been used to prove the following factorization theorem

$`D^{()}\pi |H_w|B=NF^{BD^{()}}(0){\displaystyle _0^1}𝑑\xi T(\xi ,Q,\mu )\varphi _\pi (\xi ,\mu )+\mathrm{},`$ (4)
where the ellipses denote terms that vanish faster than the leading term as $`Q\mathrm{}`$.

The relevant terms in the electroweak Hamiltonian are $`H_W=C_1O_1+C_8O_8`$, where $`O_1=\overline{d}_L\gamma ^\mu u_L\overline{c}_L\gamma _\mu b_L`$ and $`O_8=\overline{d}_LT^A\gamma ^\mu u_L\overline{c}_LT^A\gamma _\mu b_L`$. Eq. (4) was proposed in Ref. , proven to two-loops in Ref. , and proven to all orders in $`\alpha _s`$ in Ref. . The idea behind the proof is shown in Fig. 2. After integrating out offshell fluctuations the leading order interactions involve soft gluons exchanged between quarks in the $`B`$ and $`D`$ which build up the $`BD`$ form factor $`F^{BD}`$, and collinear gluons exchanged between quarks in the pion building up the light-cone pion wavefunction $`\varphi _\pi (x)`$. It should be remarked that the effective theory analysis is carried out at an operator level so it does not rely on perturbation theory.

For the proccess $`BX_s\gamma `$ in the region of photon energies $`E_\gamma \stackrel{>}{}m_B/2\mathrm{\Lambda }_{\mathrm{QCD}}2.2\mathrm{GeV}`$ the particles in $`X_s`$ are collimated in a collinear jet with offshellness $`p_X^2m_B\mathrm{\Lambda }_{\mathrm{QCD}}`$. In this process the $`B`$ meson is dominated by usoft dynamics and the leading power prediction from factorization is

$`{\displaystyle \frac{1}{\mathrm{\Gamma }_0}}{\displaystyle \frac{d\mathrm{\Gamma }}{dE_\gamma }}=H(m_b,\mu ){\displaystyle _{2E_\gamma m_b}^{\overline{\mathrm{\Lambda }}}}𝑑k^+S(k^+,\mu )J(k^++m_b2E_\gamma ,\mu ).`$ (5)
The SCET has been used to give a simple direct proof of this result . After the collinear field redefinitions, the usoft interactions rearrange themselves to leave only diagrams such as the one shown in Fig. 2. In Eq. (5) $`H`$ encodes calculable $`m_b`$ scale contributions, $`J`$ encodes $`\sqrt{m_b\mathrm{\Lambda }_{\mathrm{QCD}}}`$ scale contributions, and $`S`$ is the non-perturbative shape function . This factorization formula is required to describe the CLEO data on the photon energy spectrum .

To summarize, the soft-collinear effective theory allows factorization proofs to be simplified and formulated at the level of an effective Lagrangian and operators. The proofs involve essentially the same steps for many exclusive and inclusive processes. Finally, SCET provides us with a new framework for investigating power corrections, and in the future should be used to classify subleading non-perturbative matrix elements in a similar way to HQET.

This work was supported in part by the Department of Energy under the grant DE-FG03-00-ER-41132. I would like to thank Christian Bauer and Dan Pirjol for enjoyable collaboration on the results presented here.


text:
# The ANTARES project

## 1 Introduction

A promising challenge for exploring the Universe is the detection of high-energy ($`1`$ GeV) neutrinos. Such neutrinos could be produced by powerful cosmic accelerators, for example those in supernova remnants, active galactic nuclei, compact binaries, such as micro-quasars, and those producing gamma-ray bursts. Neutrino astronomy complements high-energy gamma astronomy since the early Universe cannot be probed with high-energy photons due to photon-matter and photon-photon interactions - gamma rays of a few hundred TeV from the Galactic Centre cannot survive their journey to the Earth. The weakly interacting nature of neutrinos, coupled with the fact that they point to their source of origin without deviation makes them unique ‘probes’ with which to investigate regions at distances larger than 50 Mpc.

The production of neutrinos in cosmic accelerators requires hadronic mechanisms to be active in the source. Whilst most of the observed high energy emissions to date can be explained by electromagnetic processes, recent observations of gamma-rays up to 5 TeV from SNR RXJ1713-39 by the Cherenkov telescope CANGAROO  suggest that electromagnetic mechanisms are inconsistent with the data but $`\pi ^0`$ decay models can readily explain the measured energy spectrum. There is, however, some debate on this result .

The technique employed by neutrino telescopes is dictated by the small neutrino cross section and the large background due to atmospheric muons. Natural Cherenkov radiators, such as water or ice, provide a large active volume at reasonable cost and the indirect detection of neutrinos through muons produced in charged current interactions increase the ’target’ region by the muon range. Detectors are constructed at large depths where the atmospheric muon flux is significantly reduced compared to that at the surface. Upward-going muons produced by neutrinos having crossed the Earth, are recognised as the products of neutrino interactions in or close to the instrumented region.

The Cherenkov light emitted by charged particles in deep water or ice is detected using an array of photomultiplier tubes (PMTs) which are housed, together with some associated electronic components, in a high pressure-resistant glass sphere known as an optical module (OM). The muon direction and energy are measured using the arrival times and amplitudes of the PMT pulses. The detector sensitivity increases with energy due to the increase in the $`\nu `$-N cross section, the longer muon range and the increase in the amount of emitted Cherenkov light through secondary particles.

## 2 Status of the ANTARES project

The ANTARES (Astronomy with a Neutrino Telescope and Abyss environmental RESearch) project  started in 1996 and involves physicists, sea science experts and engineers from France, Germany, Italy, Russia, Spain, The Netherlands and the United Kingdom. The ANTARES location (42 50’N, 6 10’E) coupled with the Earth’s rotation gives an annual sky coverage of about 3.5$`\pi `$ sr. The instantaneous overlap with the AMANDA experiment, located at the South Pole , will be approximatively 0.5$`\pi `$ sr giving an integrated common coverage of 1.5$`\pi `$ sr. The Galactic Centre, an important potential source of high-energy neutrinos, will be visible for 67% of the day.

The detector will be installed at 2400 m depth in the Mediterranean Sea, 37 km off-shore of La Seyne sur Mer, near Toulon (France). An extensive programme of site evaluation has provided the relevant environmental parameters of the detector site. Surveys have confirmed the suitability of the sea bed for detector deployment and sea currents have been measured, they average $`6`$ cm/s with a maximum observed value of 19 cm/s. The average loss of light transmission through OM glass spheres due to bio-fouling and sedimentation has been determined by long-term measurements to be only $`2\%`$ at the equator of a glass sphere one year after deployment. Moreover, this loss is a decreasing function of the zenith angle . In order to be sure that the biofouling effect will be negligible for long-term measurements, ANTARES PMTs will be oriented downwards.

Water transparency affects the light detection efficiency and thus determines the natural scale for spacing the OMs. The amount of light scattering predominantly affects track reconstruction and is a factor having a relevant effect on the angular resolution. For blue light the absorption length is $`55`$ m while the effective scattering length (the scattering length divided by $`1\mathrm{cos}\theta `$, where $`\theta `$ is the scattering angle) is above the 100 m range. At the ANTARES site the light scattering contribution is small compared to Antarctic ice. In fact, 95% of the photons emitted by a pulsed isotropic LED source at a distance of 24 m from a PMT are collected within 10 ns (30 ns) for blue (UV) light.

The detector, illustrated in Fig. 1, consists of a 3-dimensional array of OMs  arranged in strings made of mechanically resistent electro-optical cables. They will be anchored at the sea bed and held taut by buoys at the top of the string. At least 10 strings each equipped with 90 OMs arranged in 30 storeys are foreseen. Strings will be separated on the sea bed by approximately 60 m and storeys will be vertically separated by 12 m. Each storey is equipped with 3 OMs oriented at $`45^{}`$ to the downward vertical. Readout electronics for each group of 3 OMs are located in a titanium container mounted just above the OMs. The 10” Hamamatsu R7081-20 PMTs in the OMs are sensitive to single photons . The front-end electronics of ANTARES OMs includes the Analogue Ring Sampler (ARS) ASIC  which has been developed to process PMT signals and measures the arrival time and charge for single photoelectrons (99% of the events) and the full pulse shape for larger photoelectron signals. The ARS has been shown to have an intrinsic time resolution of $`\sigma =0.3`$ ns. The overall time resolution, including the PMT transit time spread, is about $`\sigma =1.5`$ ns. Signals are digitized, then transmitted via optical fibres in cables. The cables from individual strings will be connected (using a submarine) to a junction box at the end of an electro-optical cable which sends signals to a shore station where data are recorded. LED beacons (one per six storeys) will be able to emit precise short pulses of blue light in order to calibrate the detector’s timing. An extra string, called the ”Instrumentation String”, will be devoted to the measurement of environmental parameters and will include, at its base, a laser beacon that will also contribute to the timing calibration.

Trigger logic and electronics, track reconstruction and background rejection must take into account the presence of background light. The background rate is seen to consist of two constant components: the rate due to the $`\beta `$-decay of 40K present in sea water and a continuum rate that varies on a time scale of several hours due to bioluminescence. The background light from these components produces an almost constant counting rate of about 60 kHz on a 10” PMT. Furthermore, short bursts of bioluminescence have been measured with a rise time of a few milliseconds and a few seconds duration. These bursts induce a MHz level counting rate in a PMT which leads, on average, to less than 5% dead time on any individual PMT with negligible correlation between storeys.

The ANTARES collaboration has successfully achieved the following:

* from Nov. 1999 to June 2000, a so-called ”demonstrator string” (of a different mechanical design to the current design) was deployed at a depth of 1100 m and connected to shore with a 37-km electro-optical cable. This string was instrumented with 7 PMTs and more than 50000 7-fold coincidences were recorded, zenith angles were reconstructed using times and amplitudes. The shape of the zenith angular distribution of atmospheric muons was reproduced reasonably well, despite there being only a small number of PMTs and a single string . After reconstruction, the data were seen to include a $`50\%`$ contribution of multiple muons - as expected from such a shallow site. The demonstrator string allowed the ANTARES relative and absolute positioning to be tested with a system of acoustic rangemeters, compasses and tiltmeters. Relative distances between 2 elements were measured with 5 cm accuracy and absolute positioning was obtained with $`1`$ m accuracy. These values garantee a sufficiently accurate positioning for astrophysics searches.
* A test of OM sphere implosion was performed in June 2000 using a mechanical structure consisting of 2 mechanically complete but un-instrumented storeys. One OM was artificially weakened and imploded at a depth of 2600 m. After the implosion, the other spheres and the electronics box in the same storey had also imploded, but the components in the adjacent storey and the mechanical cable remained intact;
* the electro-optical cable to carry electrical power and data between the ANTARES site and shore was successfully deployed in October 2001;
* a mechanical test of a prototype string deployment has been performed in November 2001;
* in April 2002 the industrial production of 900 Optical Modules started.

The planning of the experiment includes the following key dates:

* Autumn 2002: deployment of a prototype string containing 5 storeys (called sector string) and of the instrumentation string;
* December 2002: connection to the junction box and start of data taking;
* June 2003: deployment of first two strings;
* December 2004: operation of the 10 string detector.

## 3 Expected performance

The ANTARES scientific programme is mainly devoted to the detection of neutrinos of astrophysical origin produced in point-like sources or coming from the whole sky as a diffuse flux. The neutrino signal can be detected above the background due to atmospheric neutrinos at energies greater than 10 to 100 TeV as a result of the harder neutrino spectrum expected from cosmic accelerators compared to neutrinos from cosmic ray interactions in the atmosphere. Unfortunately, estimates of neutrino fluxes at large energies are affected by an uncertainty of about 2 orders of magnitude on the prompt neutrino contribution coming from charmed meson decays. Since their spectrum is harder than those of the pion/kaon component of the atmospheric neutrino spectrum, prompt neutrino fluxes could constitute a relevant source of background for the cosmic neutrino diffuse flux search, while for point-like and transient sources the knowledge of direction and time will help to suppress this background.

An intensive study to simulate the 10 string detector response to different neutrino fluxes and to reconstruct events has been undertaken. The muon production in the neutrino interaction has been simulated by using deep inelastic scattering with CTEQ3-DIS parton functions . The Earth shadowing effect due to increasing neutrino cross-sections is properly taken into account. Muons are propagated to the detector using MUM code  taking account of different media (rock/water). The water properties are simulated including the scattering of photons and a 60 kHz background count rate due to 40K decays. A simulation of the detector response including the behaviour of the PMT and front-end electronics is also performed.

The important parameters which characterise a neutrino telescope are its effective area (which includes reconstruction and selection efficiencies), the angular resolution and the energy resolution. The muon effective area, a ratio of the selected event rate to the flux of incident muons, depends on the selection criteria used for any specific analysis. For point-like source searches, particularly when a signal should be observed across the entire sky and not necessarily associated to any known object, strict selection criteria are mandatory to keep only well reconstructed events. Thus, good pointing accuracy and a good rejection of the background of atmospheric muons is necessary. Other searches, for instance for known transient sources, such as gamma-ray bursts, which are almost background free on timescales of hundreds of seconds, can be performed by replacing stringent cuts on the quality of the reconstruction by an angular cut around the presumed source. In this case, the effective area calculation is made by requiring that the angle between the direction of the reconstructed muon and the neutrino direction is less than a defined angular cut.

The effective area for the ANTARES detector when using these two approaches is illustrated below. Fig. 2 depicts the effective area of the 10-string detector for muon events as a function of neutrino energy assuming an isotropic distribution of neutrinos inducing muons. The solid line is obtained applying strict selections in order to have high-quality reconstructed events (see below for the corresponding angular resolution plot in Fig. 4). The upper dashed curve, which exceeds the geometrical area at high energies, is obtained for reconstructed muon events within an angular error lower than $`1^{}`$ from the neutrino direction after relaxing quality cuts (e.g. for transient source searches).

Fig. 3 illustrates how the effective area depends on the intrinsic angular resolution for events selected by applying strict reconstruction quality cuts. Together with the curve including all selected events, the curves for a reconstruction error of the muon (angle between the simulated muon and the reconstructed one) less than $`1^{}`$ and than $`0.3^{}`$ are shown demonstrating that most of the events are well reconstructed. For a typical $`E^2`$ cosmic accelerator neutrino spectrum  96% (72%) of the events are reconstructed with an error less than $`1^{}`$ ($`0.3^{}`$).

The intrinsic angular resolution of the telescope, defined as the median angular separation between the ’true’ and the reconstructed muon track, has been estimated. In Fig. 4 the median value of the distribution of the angle between the reconstructed and simulated muon is shown as a solid line, while the dashed line is the median angle between the reconstructed muon and the parent neutrino as a function of the neutrino energy. The reconstruction quality cuts applied to obtain this plot are the same as those used for the solid line in Fig. 2. It should be noted that the improvement of the angular resolution with increasing energy results in an increase of the signal to noise ratio for neutrino astrophysics studies at energies greater than 10 TeV. At these energies, the pointing accuracy is only limited by the intrinsic angular resolution with a limiting value of $`0.12^{}`$. This estimate includes light scattering effects. Good angular resolution facilitates powerful background rejection when searching for point-like sources. After an optimisation of the signal ($`E^2`$) to noise ratio the optimal search cone is seen to have a half-width of $`0.7^{}`$ around the known source direction or around any measured neutrino event when looking for clusters around it. Furthermore, it is found that the ’optimal’ grid in the sky should have $`70\times 70`$ bins in right ascension and declination.111While the background is expected to be constant in right ascension, it depends on declination; hence the binning in declination is variable in order to have constant background. These results can be compared to those from AMANDA B-10 which report an optimal angular binning of $`12^{}\times 12^{}`$ .

### 3.1 Sensitivity to point-like sources and to diffuse fluxes

The effective area in Fig. 2 is directly correlated to the properties of the detector. To understand the response to a possible signal, the effective area has to be convoluted with a model-specific neutrino flux. Taking as a reference an $`E^2`$ signal spectrum, the response curve (differential rate of events versus energy) has a broad maximum around neutrino energies of $`10`$ to $`10^2`$ TeV and around $`3`$ to $`30`$ TeV of the corresponding muon energies. The maximum of the response curve moves to lower energies for softer spectra. The effective area averaged over an $`E^2`$ spectrum is of the order of 22000 m2 and decreases for softer spectra. For point-like sources it is also important to consider the dependence of the effective area on declination. This, together with an estimate of the background, allows the upper limits of the muon flux for any given type of source to be calculated. This study is performed considering physical backgrounds, including atmospheric muons and neutrinos. The resulting upper limits in units cm-2 s-1 of muon flux for 1 year of data taking and for two neutrino models are shown in Fig. 5 where the integration over neutrino energies is made assuming $`E_{min}=10`$ GeV.

To illustrate the expected sensitivity of ANTARES to point-like sources the following example can be considered. Recently a new calculation on neutrino fluxes from galactic micro-quasars has been published  which predicts some promising expected fluxes. In particular, ANTARES expects to measure 6.5 and 4.3 events per year respectively for GX339-4 and SS433 compared with a background of 0.3 events per year in a $`1^{}`$ cone.

Diffuse fluxes from astrophysical models are expected to exceed the atmospheric neutrino background at energies above 10 to 100 TeV. In Fig. 6 the expected differential event rate versus neutrino energy is shown for atmospheric neutrinos (Bartol flux ) and for the Waxman-Bahcall upper limit of $`4.5\times 10^8E^2`$ GeV-1 cm-2 s-1 sr-1  which is valid for extra-galactic sources transparent to nucleons. The shaded band at high energy is derived from a set of prompt atmospheric neutrino models (references in ). Currently the spread between the predictions covers 2 orders of magnitude.

Fig. 6 indicates that the atmospheric neutrino background can be suppressed using the information on the reconstructed energy. Below 100 GeV, the neutrino energy is calculated from the range of the muon in the detector (for neutrino oscillation studies), while at energies above 100 GeV an energy estimator based on the PMT charge amplitudes currently allows muon energies to be reconstructed within a factor of 4 at 1 TeV decreasing to a factor of 3 at 10 TeV and reaching the value of 2 in the region 10 to $`10^7`$ TeV. The ANTARES sensitivity to a diffuse neutrino flux can thus be estimated by imposing a cut on the reconstructed energy which minimizes the expected background with respect to the signal. The predicted sensitivity is $`E^2d\mathrm{\Phi }/dE=10^7`$ GeV cm-2 s-1 sr-1 above $`10^5`$ GeV, one order of magnitude better than AMANDA B-10  and Baikal  but still above the Waxman-Bahcall upper limit .

The expected event rates for the atmospheric neutrino background, for a flux equal to the Waxman-Bahcall  limit and for the more optimistic upper limit in reference  are given in Tab. 1 for reconstructed muon energies greater than 10 TeV and 100 TeV. The ’true’ event rates for the same threshold on the true muon energy are also quoted as references. The integrated event rates are illustrated in Fig. 7 where solid lines are the expected spectra and dashed lines are the reconstructed ones. The present ANTARES energy reconstruction reproduces harder spectra better than those of atmospheric neutrinos. For harder energy spectrum a small fraction of low energy events with over-estimated energy cannot contaminate significantly high energy events as it is the case for atmospheric neutrinos, which have the steepest energy spectrum. This effect is enhanced by the fact that the energy resolution is a slightly decreasing function of the energy.

### 3.2 Sensitivity to WIMPS

There is strong evidence in favour of a flat Universe with matter density $`\mathrm{\Omega }_{matter}0.35`$. This result, combined with Big Bang nucleosynthesis (baryon density $`\mathrm{\Omega }_{baryon}0.05`$) and with simulations of structure formation, indicates that most of the dark matter in the Universe could be cold (non-relativistic at decoupling time) and non-baryonic. Weakly Interacting Massive Particles (WIMPs) are considered among the most promising candidates, since a typical weak interaction cross-section would provide the relic density to close the Universe. In a supersymmetric framework, assuming R-parity conservation, a theoretically-favoured candidate is the lightest supersymmetric particle, which in many regions of the SUSY parameter space is the lightest neutralino $`\chi `$. In minimal supergravity inspired models (mSUGRA), using renormalization group equations and requiring radiative electroweak symmetry breaking, 5 parameters222These are: the common gaugino mass $`m_0`$, the common scalar mass $`m_{1/2}`$, a common trilinear coupling A, the higgsino mass parameter sign sign($`\mu `$) and the ratio of the Higgs vacuum expectation values $`\mathrm{tan}\beta `$. determine the MSSM physics including the neutralino sector.

Neutralinos in the galactic halo with typical velocities of few hundreds of kilometres per second could be slowed down sufficiently by scattering inside celestial bodies to become gravitationally trapped in their core. Indirect dark matter searches aim to detect high-energy neutrinos from WIMP annihilation with subsequent heavy quark and gauge boson decays in the core of the Sun, the Earth and the Galactic Centre. Indirect detection is complementary to direct techniques which observe the nuclear recoil from spin-independent WIMP-nucleon scattering in scintillation detectors and in cryogenic phonon Ge and Si detectors. These two techniques have different backgrounds and systematic errors and they have better sensitivity in different regions of parameter space. For axially-coupled neutralinos indirect detection is favoured because the Sun is largely made of hydrogen. In addition, dark matter in the Galactic Centre is not accessible to direct techniques.

Indirect searches aim to detect a statistically significant excess against the background of atmospheric neutrinos pointing back towards the Sun, the Galactic Centre and the nadir respectively. Recent simulation results on ANTARES sensitivity have been presented in Ref. . In Fig. 8 the sensitivity curve (90% c.l.) for 3 years of data taking as a function of neutralino mass is shown. Models are described in . The search cone for neutrino events around the core of the Sun has been optimized considering the angular distribution of the expected signal. For comparison, limits from other experiments are shown.

### 3.3 Atmospheric neutrino oscillations

Three different neutrino event topologies can be used to study neutrino oscillations: contained events with the vertex and the induced muon inside the implemented region, semi-contained events with the vertex inside the detector but with the muon leaving it and through-going events with the vertex outside the detector. Work is in progress to understand the systematic effects and to improve the precision on the oscillation parameter measurement. The measurement of the ratio, E/L, of the neutrino energy to the baseline length through the Earth for atmospheric neutrinos can provide evidence of neutrino oscillations. The best estimator of this quantity is $`E_\mu /\mathrm{cos}\theta `$, where $`E_\mu `$ is connected to neutrino energy even though the hadronic component takes away part of the energy in the interaction, and $`\mathrm{cos}\theta _\mu `$ is proportional to L, even though at neutrino energies less than about 10 GeV the average angle between the muon and the neutrino is not negligible. The visible muon energy is reconstructed from the muon range for contained events selected by a containment cut based on PMTs along the track that have not been hit; for partially-contained events an estimated correction is applied. An analysis using the 720 single-string and 2100 multi-string contained and semi-contained events expected in 4 yrs of data taking indicates that the oscillation pattern can be reconstructed in the region allowed by Super-Kamiokande (see Fig. 9).

## 4 Conclusions

The ANTARES neutrino telescope is foreseen to be fully deployed by the end of 2004. The R&D phase of the project, which started in 1996, has finished. During this phase a detailed assessment of the main requirements for an undersea neutrino telescope was made. The phase culminated in 1999 with the successful operation for several months of a demonstrator string which allowed the reconstruction of atmospheric muons. The project has now entered the construction and deployment phase for a 0.1 km2 scale detector. An electro-optical cable has been successfully deployed and a prototype string and a string dedicated to environmental parameter measurements will be deployed during this Autumn.

Simulations tools are providing relevant information concerning the scientific programme of the experiment, namely neutrino astrophysics, dark matter and neutrino oscillations. The predicted ANTARES sensitivity to fluxes from cosmic accelerators are about an order of magnitude better than results presented to date by other neutrino telescopes.


text:
# 1 The mean spatial separation between orbits generated in frozen-𝑁 systems from initial conditions separated in configuration space by a distance 𝛿𝑟(0)=10⁻⁵. Each curve was generated by averaging over 100 pairs of initial conditions evolved in 𝑁-body realisations of the chaotic triaxial ellipsoid plus black hole potential. The four curves correspond (from top to bottom) to 𝑁=10^4.5, 10⁵, 10^5.5, and 10⁶. The solid line corresponds to a slope of 0.025, generated as a least squares fit to the 𝑁=10⁶ data over the interval 32<𝑡<128. The triple dot-dashed line has a slope 0.022, equal to the mean value of the smooth potential Lyapunov exponent 𝜒_𝑆. The dashed line has a slope 0.75, equal to the mean value of the 𝑁-body Lyapunov exponent 𝜒_𝑁. The dot-dashed curve overlaying the data for 𝑁=10^4.5 represents the function 𝛿𝑟=𝐴(𝑡-𝑡₀) for 𝐴=0.008 and 𝑡₀=12.0.

SMOOTH POTENTIAL CHAOS AND N-BODY SIMULATIONS

HENRY E. KANDRUP

Department of Astronomy, Department of Physics, and

Institute for Fundamental Theory, University of Florida, Gainesville, Florida 32611

kandrup@astro.ufl.edu

IOANNIS V. SIDERIS

Department of Astronomy, University of Florida, Gainesville, Florida 32611

sideris@astro.ufl.edu

ABSTRACT

Integrations in fixed N-body realisations of smooth density distributions corresponding to a chaotic galactic potential can be used to derive reliable estimates of the largest (finite time) Lyapunov exponent $`\chi _S`$ associated with an orbit in the smooth potential generated from the same initial condition, even though the $`N`$-body orbit is typically characterised by an $`N`$-body exponent $`\chi _N\chi _S`$. This can be accomplished either by comparing initially nearby orbits in a single $`N`$-body system or by tracking orbits with the same initial condition evolved in two different $`N`$-body realisations of the same smooth density.

1. INTRODUCTION AND MOTIVATION

At the present time, there are two general approaches which can be used to model the structure and evolution of systems like elliptical galaxies. On the one hand, one can perform detailed $`N`$-body simulations, solving the coupled equations of motion either exactly or in some approximation. On the other, one can try to construct equilibrium models as solutions to the collisionless Boltzmann equation, either analytically (e.g., using the approach developed by Hunter & Qian ) or numerically (e.g., by implementing some version of Schwarzschild’s  method).

Both approaches would appear to be extremely useful, but each has its limitations. Despite significant advances in hardware, direct $`N`$-body simulations still cannot be performed for particle number as large as $`N\mathrm{\hspace{0.33em}10}^{11}`$, so that one cannot consider point masses in the integrations as corresponding to individual stars in real galaxies. Integrations in smooth potentials have the advantage that one can consider characteristics which, presumably, correspond to the orbits of individual stars. However, the very assumption that the system can be described by a smooth potential constitutes an idealisation which, albeit generally accepted, has never been proven rigorously.

Obvious questions to be answered thus include the following: To what extent is it true that, for sufficiently large $`N`$, solutions to the full gravitational $`N`$-body problem can be mimicked by motions in a smooth potential? And, to the extent that the smooth potential approximation is not completely sufficient, to what extent can discreteness effects really be modeled by friction and noise in the context of a Fokker-Planck description (cf. Rosenbluth, MacDonald, & Judd 1957)? In particular, to what extent are ‘real’ $`N`$-body orbits well mimicked by solutions to a (time-dependent) Langevin equation (cf. Chandrasekhar 1943) which incorporates dynamical friction and Gaussian white noise? Fokker-Planck descriptions were formulated originally to extract statistical information about long time behaviour, assuming implicitly that the bulk potential is integrable or nearly integrable. However, recent years have seen a growing recognition that galactic potentials may admit a fair amount of chaos; and analyses of the short time behaviour of individual orbits in chaotic potentials have provided compelling evidence that friction and noise can dramatically accelerate phase space transport (Lieberman & Lichtenberg 1972, Lichtenberg & Wood 1989, Kandrup, Pogorelov, & Sideris 2000, Siopis & Kandrup 2000). Does this really mean that discreteness effects can be important already on time scales much shorter than the relaxation time $`t_R`$?

Closely related to these issues is the nature of the continuum limit. In what sense is it true that, as $`N\mathrm{}`$, orbits in the $`N`$-body potential converge towards characteristics in some smooth potential? Superficially, at least, it might seem that such a convergence is impossible. Dating back to Miller (1964), it has been recognised that the $`N`$-body problem is chaotic in the sense that individual orbits exhibit exponential sensitivity towards small changes in initial conditions; and it seems generally accepted today that, when expressed in units of inverse dynamical times $`t_D^1`$, the largest $`N`$-body Lyapunov exponent $`\chi _N`$ does not converge towards zero as $`N\mathrm{}`$ (cf. Kandrup & Smith 1991, Goodman, Heggie & Hut 1994), even if the $`N`$-body system samples an integrable density distribution. Indeed, recent work, both numerical (Hemsendorf & Merritt 2002) and analytic (Pogorelov 2001), suggests that, even for a density distribution corresponding to an integrable potential, the largest Lyapunov exponent may actually increase with increasing $`N`$. As probed in the usual way, the $`N`$-body problem may become more chaotic as $`N`$ increases!

A complete resolution to this apparent conundrum will require long time integrations of systems with very large $`N`$, which is impractical using current hardware. However, considerable insight into the continuum limit can be, and has been, obtained by studying the properties of orbits and orbit ensembles evolved in frozen-$`N`$ systems, i.e., fixed (in time and space) $`N`$-body realisations of specified smooth density distributions (Kandrup & Sideris 2001, Sideris & Kandrup 2002). In particular, that work led to several significant conclusions: (1) The largest $`N`$-body Lyapunov exponent $`\chi _N`$ does not decrease with increasing $`N`$, even for an integrable density distribution. However, there is still a clear, quantifiable sense in which, as $`N`$ increases, the $`N`$-body orbits become progressively more similar to smooth potential characteristics. (2) As $`N`$ increases the Fourier spectra associated with $`N`$-body orbits more closely resemble the spectra associated with characteristics in the smooth potential, be these either regular or chaotic. (3) Alternatively, viewed macroscopically $`N`$-body orbits and smooth characteristics with the same initial condition typically diverge as a power law in time on a time scale $`t_G(N)`$ that increases with increasing $`N`$. For the case of regular characteristics, $`t_GN^{1/2}t_D`$; for chaotic characteristics, $`t_G(\mathrm{ln}N)t_D`$. It follows that, for sufficiently large $`N`$, $`N`$-body orbits and smooth potential characteristics remain close for comparatively long times.

This would seem a result of some significance, but it still begs an important issue. If $`N`$-body orbits converge towards smooth potential characteristics, it should be possible, at least for sufficiently large $`N`$, to extract information about any chaos that may be associated with the bulk potential. In particular, it must be possible to extract estimates of finite time (cf. Grassberger, Badii, & Politi 1988) Lyapunov exponents $`\chi _S`$ for motion in the smooth potential, even though the $`N`$-body orbits themselves are characterised by exponents $`\chi _N`$ which, typically, are much larger than $`\chi _S`$. (Chaotic orbits in generic smooth potentials typically have a largest Lyapunov exponent $`\chi _St_D^1`$. For both interacting \[Hemsendorf & Merritt 2002\] and frozen-$`N`$ \[Kandrup & Sideris 2000\] Plummer systems, $`\chi _N\mathrm{\hspace{0.33em}20}t_D^1`$ for $`N\mathrm{\hspace{0.33em}10}^5.`$)

The aim of this note is to demonstrate that this can in fact be done. In particular, it will be shown that there are at least two different ways in which estimates of $`\chi _S`$ can be extracted from frozen-$`N`$ systems, one involving a comparison of orbits generated from nearby initial conditions and the other involving integrations of the same initial condition in two different $`N`$-body realisations of the same smooth density distribution. Section 2 describes the algorithms and then exhibits representative results for two simple model potentials. Section 3 interprets the success of these algorithms by postulating the existence of two ‘types’ of chaos, microscopic chaos, or microchaos, which is generic to the $`N`$-body problem, and macroscopic chaos, or macrochaos, possibly associated with the bulk potential which, if present, will also be manifested in the $`N`$-body problem.

2. NUMERICAL ESTIMATION OF SMOOTH POTENTIAL LYAPUNOV EXPONENTS

2.1. Models Considered

The algorithms described in this Section will be applied to representative orbits evolved in two different density distributions:

1. A spherically symmetric Plummer sphere, for which

$$\rho _P(r)=\left(\frac{3M}{4\pi b^3}\right)\left(1+\frac{r^2}{b^2}\right)^{5/2}.$$
(1)
This corresponds via Poisson’s equation to a spherically symmetric, and hence integrable, potential

$$\mathrm{\Phi }_P(r)=\frac{GM}{\sqrt{r^2+b^2}}.$$
(2)
Units were so chosen that $`G=M=b=1`$.

2. A constant density triaxial ellipsoid, for which

$$\rho _E(𝐫)=\frac{3M}{4\pi abc}\times \{\begin{array}{cc}m^2\hfill & \text{if }m^2\mathrm{\hspace{0.33em}1}\text{,}\hfill \\ 0\hfill & \text{if }m^2>1\text{,}\hfill \end{array}$$
(3)
with

$$m^2=\left(\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}\right),$$
(4)
perturbed by a spherically symmetric density spike (black hole) of mass $`M_{BH}`$, this corresponding to a potential

$$\mathrm{\Phi }_E(𝐫)=\mathrm{\Phi }_0+\frac{1}{2}\left(\omega _a^2x^2+\omega _b^2y^2+\omega _c^2z^2\right)\frac{GM_{BH}}{\sqrt{r^2+ϵ^2}},$$
(5)
with $`ϵ=10^3`$. Attention focused on the case $`M=1.0`$ and $`M_{BH}=10^{1.5}M\mathrm{\hspace{0.33em}0.0316228}`$, and units were again chosen so that $`G=1`$. The axis ratios were taken as $`a=1.95`$, $`b=1.50`$, and $`c=1.05`$, which yield (cf. Bertin 2000) $`\mathrm{\Phi }_01.00608`$, $`\omega _a\mathrm{\hspace{0.33em}0.4663}`$, $`\omega _b\mathrm{\hspace{0.33em}0.5508}`$, and $`\omega _c\mathrm{\hspace{0.33em}0.6753}`$. For energies sufficiently small that orbits are restricted to $`m<1`$, the phase space is almost completely chaotic (Kandrup & Sideris 2002). This implies that one need not worry about transitions between regular and chaotic behaviour which can be induced by discreteness effects in more complex potentials like the triaxial generalisations of the Dehnen (1993) potential, which admit a complex coexistence of both regular and chaotic orbits.

Frozen-$`N`$ orbits in these systems were integrated with a particle-particle numerical scheme using a variable timestep integrator with accuracy parameter $`10^8`$ which conserved energy to at least one part in $`10^6`$. The $`1/r`$ kernels for the individual masses were regulated through the introduction of a softening parameter $`ϵ=10^5`$.

3.2. The Numerical Algorithms

Algorithm 1. Involving a comparison of orbits generated from nearby initial conditions in the same frozen-$`N`$ density distribution.

The key realisation here is that, even though orbits generated from nearby initial conditions diverge initially at a rate $`\mathrm{\Lambda }\chi _N`$, this divergence quickly saturates once the growing separation between the orbits becomes large compared with the typical distance between neighboring point masses. If the smooth potential characteristic is regular, this initial exponential divergence is replaced immediately by a more modest power law divergence which proceeds on a time scale $`t_GN^{1/2}t_D`$. If instead the characteristic is chaotic, the initial exponential divergence at a rate $`\mathrm{\Lambda }\chi _N`$ is replaced by a slower exponential divergence at a rate $`\mathrm{\Lambda }\chi _S`$ which, in turn, typically proceeds until the separation becomes ‘macroscopic’, i.e., comparable to the size of the entire accessible phase space region. At this point, this second exponential divergence saturates and is replaced by a power law divergence which proceeds on a time scale $`t_G(\mathrm{ln}N)t_D`$.

These statements can be corroborated straightforwardly by tracking the separation

$$\delta r(t)=|𝐫_1𝐫_2|$$
(6)
and plotting $`\mathrm{ln}\delta r(t)`$ as a function of $`t`$. By so doing, one discovers that, just as for the smooth potential, nearby initial conditions can exhibit significantly different values of $`\chi _S`$. (Nearby initial conditions evolved in a smooth potential can have significantly different finite time Lyapunov exponents even if, for late times, these exponents converge towards the same asymptotic $`\chi _{\mathrm{}}`$.) If, alternatively, one wishes to estimate a ‘typical’ $`\chi _S`$ for orbits in some given phase space region, one can select an ensemble of initial conditions from that region, evolve each initial condition into the future, and then extract a mean $`\chi _S`$ from the time-dependent mean separation

$$\delta r=\frac{1}{k}\underset{i=1}{\overset{k}{}}\delta r_i.$$
(7)
Figure 1 exhibits the results of such a computation for an ensemble of 100 chaotic initial conditions evolved in a frozen-$`N`$ realisation of the the ellipsoid plus black hole system. The different curves in the Figure represent frozen-$`N`$ backgrounds with $`N`$ varying between $`N=10^{4.5}`$ and $`N=10^6`$. Figure 2 exhibits analogous data generated for regular initial conditions in the Plummer potential. In each case, the initial conditions sampled a phase space region of size $`\mathrm{\Delta }r\mathrm{\Delta }v\mathrm{\hspace{0.33em}10}^3`$, with the perturbed orbits being generated from initial conditions displaced from the original initial conditions by a distance $`\delta r(0)=10^5`$ in a randomly chosen direction.

It is evident that, for both the ellipsoid and Plummer potentials, the mean separation $`\delta r`$ begins by diverging at a rate that is comparable to the value of the largest $`N`$-body Lyapunov exponent, $`\chi _N`$. For the case of the ellipsoid potential, this exponential divergence is (at least for sufficiently large $`N`$) eventually replaced by a slower exponential divergence at a rate $`\chi _S`$ which persists until $`\delta r`$ becomes ‘macroscopic,’ i.e., comparable to the size of the accessible phase space region. At this point the exponential divergence is replaced by a slower power law separation which proceeds until $`\delta r`$ saturates. For the case of the Plummer potential, the second exponential phase is absent, the initial exponential phase being replaced immediately by a power law growth. This is especially evident from Figure 3, which plots the same data as Figure 2 on a linear scale.

As noted already, for both regular and chaotic potentials the first exponential phase typically stops once the separation $`\delta r`$ becomes large compared with the mean interparticle spacing $`n^{1/3}`$, with $`n`$ a typical number density. For the constant density ellipsoid, $`n^{1/3}\mathrm{\hspace{0.33em}0.427}N^{1/3}`$. For $`N=10^5`$ this corresponds to $`n^{1/3}\mathrm{\hspace{0.33em}0.00919}`$ and $`\mathrm{ln}n^{1/3}4.69`$; for $`N=10^6`$, $`n^{1/3}\mathrm{\hspace{0.33em}0.00427}`$ and $`\mathrm{ln}n^{1/3}5.45`$. Alternatively, the second, slower exponential phase exhibited by chaotic potentials typically ceases once the separation has become ‘macroscopic’, i.e., comparable to the size of the accessible phase space region. This macroscopic scale appears to coincide with the scale on which two initially nearby orbits evolved in the smooth potential will cease their exponential divergence.

As discussed more carefully in Section 3, the comparatively sharp break between the first and second phases indicates that one can implement a de facto distinction between microscopic chaos associated primarily with close encounters and macroscopic chaos associated with the bulk potential, even though the effects of these sources of chaos are not completely decoupled.

It is evident that this prescription for estimating $`\chi _S`$ can only work for comparatively large values of $`N`$, where the typical interparticle spacing is much smaller than the total size of the accessible configuration space region. If this condition is not satisfied, $`\delta r(t)`$ will become ‘macroscopic’ almost as soon as it becomes large compared with $`n^{1/3}`$, so that the intermediate second stage disappears. For the particular model exhibited in Figure 1, one requires $`N>10^5`$ or so in order to obtain a reasonable estimate of $`\chi _S`$. Indeed, it is evident from the dot-dashed curve in Figure 1 that, after the initial exponential divergence and before the final saturation, the data for $`N=10^{4.5}`$ can be well fit by a linear growth law $`\delta r=A(tt_0)`$.

Algorithm 2. Involving a comparison of orbits generated from a single initial condition evolved in two different frozen-$`N`$ density distributions that sample the same smooth density.

The smooth exponent $`\chi _S`$ also provides information about the rate of divergence associated with orbits generated in two different frozen-$`N`$ simulations. Specifically, if one computes $`\delta r`$ for a collection of orbits evolved in two different frozen-$`N`$ density distributions associated with the same chaotic potential, one again finds an evolution manifesting the same three stages. This is illustrated in Figure 4, which was generated for the ellipsoid potential for the same initial conditions as Figure 1. Figure 5 exhibits analogous data generated for the Plummer potential. It is evident once again that, for the Plummer potential, the intermediate stage is absent.

Given that $`\chi _S`$ provides information about the rate of divergence of orbits in two different frozen-$`N`$ backgrounds, each of which can be viewed intuitively as a ‘perturbation’ of the smooth density distribution, one might also expect that $`\chi _S`$ provides information about the rate at which orbits in a single frozen-$`N`$ simulation diverge from smooth characteristics with the same initial condition. As illustrated in Figure 6, this expectation is in fact correct. The macroscopic power law divergence between the frozen-$`N`$ orbit and the corresponding smooth potential characteristic, which parallels the behaviour observed if Algorithm 1 or 2 is implemented, has been discussed extensively elsewhere (cf. Figure 8 in Kandrup & Sideris 2001 and Figure 2 in Sideris & Kandrup 2002 and the accompanying discussion).

3.3. Why these algorithms work

Given the assumption that orbits in an $`N`$-body system ‘feel’ two different sorts of chaos which act on different scales, it is not surprising that one can derive estimates of both the $`N`$-body Lyapunov exponent $`\chi _N`$ and the smooth potential $`\chi _S`$ from a comparison of initially proximate orbits in a single frozen-$`N`$ system. However, the fact that $`\chi _S`$ is also related to the rate of divergence of orbits in different frozen-$`N`$ simulations is, perhaps, less obvious. The key to understanding this phenomenon is the fact that discreteness effects really can be mimicked by dynamical friction and Gaussian white noise in the context of a Langevin description.

Specifically (Sideris and Kandrup 2002), at both the level of individual orbits, as probed, e.g., by Fourier spectra, and at the level of orbit ensembles, as probed, e.g., by the efficiency of phase mixing for both regular and chaotic orbit ensembles, discreteness effects associated with an $`N`$-body density distribution are extremely well reproduced by Gaussian white noise with a ‘temperature’ $`\mathrm{\Theta }|E|`$, with $`E`$ the orbital energy, and a coefficient of dynamical friction $`\eta \mathrm{\hspace{0.33em}1}/N`$. This dependence on $`N`$ is of course very similar to the scaling $`t_R\eta ^1N/(\mathrm{ln}N)`$ predicted in a conventional Fokker-Planck description; and, indeed, given the limited range in $`N`$ that can be probed numerically (the notion of a smooth potential appears to break down for $`N10^3`$; simulations with $`N\mathrm{\hspace{0.33em}10}^6`$ become prohibitively expensive computationally!), the numerical simulations are completely consistent with this scaling.

But why does this explain the observed divergence of orbits in different frozen-$`N`$ distributions? The crucial point here, as described, e.g., in Habib, Kandrup & Mahon (1997) or Kandrup and Novotny (2002), is that, viewed mesoscopically, an ensemble of noisy orbits with fixed $`\mathrm{\Theta }`$ and $`\eta `$, each generated from the same chaotic initial condition or from a set of very nearby initial conditions, will typically disperse in such as a fashion that

$$\delta r(\mathrm{\Theta }\eta )^{1/2}\mathrm{exp}(\chi _St).$$
(8)
Given the assumed scaling $`\eta \mathrm{\hspace{0.33em}1}/N`$, it then follows that

$$\mathrm{ln}\delta r=\mathrm{const}+\frac{1}{2}\mathrm{ln}\eta +\chi _St=\mathrm{const}\frac{1}{2}\mathrm{ln}N+\chi _St.$$
(9)
Numerical simulations demonstrate that noisy orbits diverge at the same rate $`\chi _S`$ observed for orbits in frozen-$`N`$ simulations; and the connection between $`N`$ and $`\eta `$ implicit in a Fokker-Planck description makes a specific prediction as to the $`N`$-dependence of the exponential prefactor. If discreteness effects really can be modeled as Gaussian noise, $`\delta r`$ should satisfy eq. (9). To the extent that the ‘mean’ trajectory associated with the noisy ensemble coincides, at least approximately, with the smooth potential characteristic, the same scaling should also be observed when comparing noisy orbits and smooth potential characteristics.

This scaling implies (a) that $`\mathrm{ln}\delta r`$ should grow linearly at a rate $`\chi _S`$, independent of $`N`$, but (b) that, for fixed $`t`$, an increase in $`N`$ or a decrease in $`\eta `$ by an order of magnitude should decrease $`\mathrm{ln}\delta r`$ by $`(1/2)\mathrm{ln}10\mathrm{\hspace{0.33em}1.15}`$. That this scaling is in fact realised for both noisy orbits and orbits in frozen-$`N`$ backgrounds is evident from Figure 7. Here the solid curves exhibit the results of noisy integrations of the same initial conditions used to generate Figure 1 and 4, all computed with $`\mathrm{\Theta }=1.0`$ but allowing for values of $`\eta `$ extending from $`10^4`$ to $`10^7`$. The dots accompanying the upper two curves represent results from frozen-$`N`$ integrations for $`N=10^{4.5}`$ and $`10^{5.5}`$ (the same data plotted in Figure 6). The results from Sideris and Kandrup (2002) suggest a best fit correspondence $`\mathrm{log}_{10}\eta =\mathrm{log}_{10}N+p`$, with $`p\mathrm{\hspace{0.33em}0.5}`$; and indeed, it is apparent visually that the noisy and frozen-$`N`$ curves for, e.g., $`\eta =10^5`$ and $`N=10^{5.5}`$ are extremely similar.

Presuming that this scaling holds for smaller $`\eta `$ as well, the lowest curve in Figure 7 should correspond, at least approximately, to $`N=10^{7.5}`$.

III. DISCUSSION

This paper has described two algorithms which can be used to obtain estimates of the largest (finite time) Lyapunov exponents $`\chi _S`$ associated with a smooth density distribution from frozen-$`N`$ realisations of that density. The first involves comparing two orbits in a single frozen-$`N`$ system generated from nearby initial conditions. The second involves comparing orbits with the same initial condition evolved in two different frozen-$`N`$ systems, each sampling the same smooth density distribution. The success of the first algorithm emphasises the fact that detailed information about the bulk potential really is buried in an $`N`$-body simulation. The success of the latter emphasises another important point, namely that smooth potential Lyapunov exponents $`\chi _S`$ also provide information about the divergence of the same initial condition in different $`N`$-body systems, i.e., information about the extent to which, viewed mesoscopically, discreteness effects limit the intrinsic reliability of orbits in a pointwise sense. As stressed already, the success of this alternative algorithm reflects the fact that discreteness effects really can be well mimicked by Gaussian white noise in the context of a Fokker-Planck description.

The key point in all this is that, for the case of a chaotic bulk potential, two nearby initial conditions will, when evolved into the future, exhibit a three stage evolution, reflecting the effects of both microscopic and macroscopic chaos, i.e., microchaos and macrochaos.

(1) For early times and small separations, the orbits will diverge exponentially at a rate comparable to a typical $`N`$-body Lyapunov exponent $`\chi _N`$.

(2) However, once the separation between the orbits becomes large compared with the typical interparticle spacing, this divergence ceases and is replaced by a slower divergence at a rate $`\chi _S`$, which proceeds until the separation becomes macroscopic.

(3) At still later times the orbits exhibit a more modest power law divergence.

The computations described here only yielded estimates of finite time Lyapunov exponents, not the true Lyapunov exponent as defined in a late time limit, which, given a complex phase space, can be much larger or smaller. This, however, is not necessarily bad. Although old in physical time, galaxies are young objects when expressed in terms of the dynamical time $`t_D`$ – typically no more than $`\mathrm{\hspace{0.33em}100}200t_D`$ in age – so that such asympotic limits are not well motivated physically. However, one might argue that, even though an asymptotic limit is not justified for individual orbits, the true Lyapunov exponent is important in that (cf. Kandrup & Mahon 1994) it characterises the average instability associated with the invariant measure, i.e., a uniform sampling of the chaotic portions of the constant energy hypersurface. The obvious point, then, is that to obtain an estimate of the true Lyapunov exponent it suffices to repeat the calculations described here for an ensemble of initial conditions sampling the constant energy hypersurface which, as described elsewhere (cf. Kandrup, Sideris, & Bohn 2001) is straightforward numerically. Alternatively, one can actually compute estimates of $`\chi _S`$ in the usual way (Benettin, Galgani, and Strelcyn 1976) by tracking the evolution of a small perturbation which is periodically renormalised, provided only that one makes sure that the perturbation always remains large compared with the scale on which the microscopic chaos saturates, i.e., very large compared with a typical interparticle spacing.

As a practical matter, it would appear that the algorithm can work for any $`N`$-body system in which the typical interparticle spacing is sufficiently small compared to the size of the system. If $`N`$ is not sufficiently large, the first exponential phase will not saturate until the separation of the originally proximate orbits has become macroscopic, so that the second exponential phase is lost. For the models considered here, one requires $`N>10^5`$ or so. For systems manifesting very high density contrasts, e.g., triaxial generalisations of the cuspy Dehnen potentials, one may require much larger $`N`$ to obtain an adequate sampling of the central region. In point of fact, however, the requirement of large $`N`$ is more than a practical consideration: it would appear that, if $`N`$ is too small, the very notion of a bulk potential becomes suspect. One finds, e.g., that, for the ellipsoid plus black hole potential, discreteness effects can be reasonably well modeled by Gaussian white noise for $`N\mathrm{\hspace{0.33em}10}^4`$, but that this model fails for substantially smaller $`N`$ (Sideris & Kandrup 2002).

Viewed from the standpoint of nonlinear dynamics, the $`N`$-body problem – or at least the frozen-$`N`$ model considered here – constitutes an interesting example of a system in which chaos can arise for different reasons on different scales. Because gravity is strong on short scales, one finds generically that close encounters between nearby particles trigger chaos on a time scale that is typically short compared with $`t_D`$. As probed by $`N`$-body Lyapunov exponents computed in the usual way, this microchaos does not decrease with increasing $`N`$; if anything, $`\chi _N`$ increases with increasing $`N`$. However, the ‘range’ of the chaos, expressed relative to the total size of the system, does decrease with increasing $`N`$ since the typical interparticle spacing scales as $`N^{1/3}`$. Alternatively, because of the long range character of the gravitational interaction, one also encounters the possibility of macrochaos which will arise if the bulk density distribution corresponds to a bulk potential that admits global stochasticity. Both forms of chaos can play a nontrivial role in the $`N`$-body problem and, as has been shown here, information about $`\chi _N`$ and $`\chi _S`$ can both be extracted from a judicious analysis of numerical data.

In this sense, it would appear that although, strictly speaking, the rate at which $`N`$-body orbits diverge is set by the Lyapunov exponent $`\chi _N`$, it may be misleading to assert (Heggie 1991) that “the approximation of a smooth potential is useful for studying orbits, but not for studying their divergence.” The smooth potential Lyapunov exponent $`\chi _S`$ does indeed provide useful information regarding the divergence of orbits on mesoscopic scales large compared with the interparticle spacing but small compared with the size of the system.

Portions of this paper were written while HEK was a visitor at the Aspen Center for Physics, the hospitality of which is acknowledged gratefully. This work was supported in part by the NSF through AST0070809.

REFERENCES

Bertin, G. 2000, Dynamics of Galaxies (Cambridge: Cambridge University Press)

Chandrasekhar, S. 1943, Rev. Mod. Phys. 15, 1

Dehnen, W. 1993, MNRAS, 265, 250

Bennetin, G., Galgani, L., & Strelcyn, J.-M. 1976, Phys. Rev. A 14, 2338

Goodman, J., Heggie, D. C., & Hut, P. 1994, ApJ, 415, 715

Grassberger, P., Badii, R., & Politi, A. 1988, J. Stat. Phys. 51, 135

Habib, S., Kandrup, H. E., Mahon, M. E. 1997, ApJ, 480, 155

Heggie, D. 1991, in Predictability, Stability, and Chaos in $`N`$-Body Dynamical Systems, ed. A.

Roy (New York: Plenum), 47

Hemsendorf, M. & Merritt, D. 2002 (astro-ph/0205538)

Hunter, C. & Qian, E. 1993, MNRAS, 262, 401

Kandrup, H. E. & Mahon, M. E. 1994, A&A, 290, 762

Kandrup, H. E. & Novotny, S. 2002 (astro-ph/0204019)

Kandrup, H. E., Pogorelov, I. V., & Sideris, I. V. 2000, MNRAS, 311, 719

Kandrup, H. E. & Sideris, I. V. 2001, Phys. Rev. E 64, 056209

Kandrup, H. E. & Sideris, I. V. 2002, Celestial Mechanics 82, 61

Kandrup, H. E., Bohn, C. L., & Sideris, I. V. 2001, Phys. Rev. E 65, 016214

Kandrup, H. E. & Smith, H. 1991, ApJ, 374, 255

Lichtenberg, A. J. & Wood, B. P. 1989, Phys. Rev. Lett. 62, 2213

Lieberman, M. A. & Lichtenberg, A. J. 1972, Phys. Rev. A 5, 1852

Miller, R. H. 1964, ApJ, 140, 250

Pogorelov, I. V. 2001, University of Florida Ph. D. dissertation

Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957, Phys. Rev. 107, 1

Schwarzschild, M. 1979, ApJ, 471, 82

Sideris, I. V. & Kandrup, H. E. 2002, Phys. Rev. E 66, 066203

Siopis, C. & Kandrup, H. E. 2000, MNRAS, 319, 43