text:
# 1 Introduction

## 1 Introduction

Classical methods for exploring the mechanisms of brain function involve recording the electrical activity of single nerve cells. Much has been learned from this approach, but there are several reasons to go beyond single neuron recording \[Rieke et al., 1997\]. First, multineuron recording greatly increases the efficiency of studying the properties of single neurons. Second, recording simultaneously from many neurons allows access to the precise temporal relations among action potentials in multiple neurons \[Usrey et al., 1998\]. This can provide a testing ground for the hypothesis that these temporal relationships carry significant information. Finally, multineuron recording experiments might give us a glimpse into the collective dynamics in neural networks, e.g. the existence of multiple stable states and possibility of switching between them \[Amit, 1989, Abeles and Gerstein, 1988\].

This paper is concerned with one of the many technical problems that arise in trying to adddress the above questions, namely the problem of sorting out signals from multiple neurons as they appear on multiple electrodes \[McNaughton et al., 1983, Reece and O’Keefe, 1989\]. In a multineuron recording each cell can appear on multiple electrodes and multiple cells can appear on each electrode. Finding the spike times for each cell is difficult because spikes from different cells are very similar in shape, making it hard to distinguish among them. In addition, the events which might be most interesting, synchronous spiking of nearby cells, are among the most difficult to disentangle.

## 2 The spike sorting problem

Throughout the nervous system, cells generate stereotyped electrical pulses termed action potentials or spikes. Thus, in the absence in the absence of noise or variability, it is expected that each neuron in the system would generate a characteristic signal at each electrode. This waveform can be viewed as a single point in signal space. In the presence of noise, however, these discrete points spread into clusters of points, each point representing a single spike. One of the objectives of the spike sorting method described here is to make a model of these clusters that is accurate enough to allow the assignment of any given voltage waveform to one of the clusters. Even though it is unlikely that all the details of the cluster model are correct, our hope is that any such errors will not hamper spike classification.

If $`𝒱`$ is the set of voltages measured on all of the electrodes during a short time segment, the probability distribution of these voltages may be decomposed as

$$P(𝒱)=\underset{c}{}P(𝒱|c)P(c),$$
(1)
where $`P(c)`$ is the total probability of observing a spike from cell $`c`$ and $`P(𝒱|c)`$ is the distribution of voltage waveforms that arise from this cell. Conversely, if a set of voltage signals $`𝒱`$ is observed, the probability that it comes from cell $`c`$ can be derived from Equation (1) by applying Bayes’ rule,

$$P(c|𝒱)=\frac{P(𝒱|c)P(c)}{_c^{}P(𝒱|c^{})P(c^{})}.$$
(2)
Ideally one would like assignments which are certain. Statements of the form “probably a spike from cell 2, but maybe from cell 1,” “probably a spike from cell 3, but maybe overlapping spikes from cells 4 and 5,” defeat our purpose. To avoid this type of problem the distributions $`P(𝒱|c)`$ and $`P(𝒱|c^{})`$ should not overlap when $`cc^{}`$. Furthermore, it is hoped that if these distributions overlap only slightly, a precise model for the distribution of each cluster may not be necessary and a small number of parameters may suffice to make reliable distinctions among the clusters. If a reliable model of this form can be made, then the assignment problem is solved. Note also that, in this limit, our prior assumptions about the likelihood of different neural responses plays no role.

The approach to the clustering problem described here is as follows. First it is assumed that individual clusters $`P(𝒱|c)`$ have a Gaussian form in which each frequency component fluctuates independently. Next, the best possible clusters are found and the mean and variances at each frequency are calculated for every cluster. If this Gaussian model were correct, the probability that a given cluster could generate a particular waveform would depend only on the $`\chi ^2`$ distance between that waveform and the cluster mean, with appropriate weighting by the variances. These $`\chi ^2`$ values provide a set of new dimensions along which distributions should be nonoverlapping if certainty in the assignments has been achieved. In addition, it is possible to check if the clustering is self consistent because every member waveform of a given cluster will yield a distance of $`\chi ^21`$ for that cluster and $`\chi ^2>>1`$ for all of the other clusters.

This combination of $`\chi ^2`$ below a threshold for one cluster and above threshold for all other clusters is the signature of unambiguous assignment. The reader should note that we achieve this result despite the fact that the real clusters need not obey the assumptions of our simple model. This can be seen by looking closely at the form of the resultant $`\chi ^2`$ distributions. The implementation of these ideas takes on the following form. First, a statistical model of the spikes from different cells is built by clustering. Second, in a completely independent stage we detect the occurrence of these spikes in the data. This is initially done for single spikes. Finally, by superimposing these statistical models, overlaping spikes are treated as well.

We emphasize once more that we are not proposing the Gaussian model as an exact model of the relevant probability distributions. On the contrary, the goal is to show that the Gaussian model suggests dimensions along which clusters are discriminable, and that the real clusters are almost perfectly discriminable along these dimensions. Once this is established, the precise form of the model is irrelevant. In cases where nonGaussian behavior is known to be crucial, one might start with different assumptions but the general strategy would be the same: we want to exhibit explictly the discriminability of clusters along some small number of relevant directions in waveform space, and the starting model is just an aid to finding these dimensions. Similar ideas using a $`\chi ^2`$ test were proposed in recent work by Pouzat et al. \[Pouzat, 2002\].

## 3 The spike sorting algorithm

The first step is to build the set of clusters. This is done by identifying the different spike shapes in the data and constructing a statistical model for each recognized spike type. As with many other spike sorting algorithms the work described here is based on the following assumptions: a) different spikes from the same cell are very similar, b) spikes from different cells have different waveforms on at least one electrode, c) if a cell fires once, it fires many times and d) overlaps are fairly rare.

### 3.1 Clustering

Before the actual clustering can be done a large number of each of the different spike types is needed. To do this the data is broken into short frames. The content of each frame is then examined to see if it contains “clean” spikes which are described below. A frame refers to a set of data snippets, one from each electrode at a given time. We emphasize that during this initial pass through the data we are interested only in collecting clean representatives of every spike type and no attempt is made to deal with frames that did not obviously contain only one single spike.

The object of the clustering is to group similar spikes together. This must work in spite of the fact that two similar spikes may be shifted slightly in time \[Lewicki, 1998\]. To deal with this problem great care is given to accurately aligning the spikes during the clustering process.

Objects to be clustered- While we measure time domain signals, from now on each frame, $`f`$, is represented by a vector composed of the concatenated Fourier transforms of the voltage waveforms of the data from each of the $`N_e`$ electrodes, i.e.,

$$S_f(\omega )=[S_{f,1}(\omega ),S_{f,2}(\omega ),\mathrm{}S_{f,N_e}(\omega )].$$
(3)
The alignment of spikes in the frequency domain is achieved simply by multiplying the Fourier components by $`e^{i\omega \tau }`$ where $`\tau `$ is the necessary time delay (see below). It should be noted that although the work presented here deals with sorting concatenated spectra, other objects can be sorted as well. Extensive experiments on sorting different objects were done while developing the methods described here but for our data $`S_f(\omega )`$ proved to be the most useful. For example, reference \[Rinberg et al., 1999\] describes the sorting of transfer functions between electrodes which is independent of the spike shape. This may be of interest in cases where the spike shape can change \[Fee et al., 1996b\] (e.g. in bursting cells). In certain cases power spectra, which are invariant to time shifts, can be used as well.

Description of the algorithm- The clustering algorithm is outlined in Box(1). Its various steps are described below.

* Line 1: Initialization- First all frames are averaged yielding an average signal,

$`\overline{S}_{0,e}(\omega )`$ $`={\displaystyle \frac{1}{N_f}}{\displaystyle \underset{f=1}{\overset{N_f}{}}}S_{f,e}(\omega )`$ (4)
and variance,
$`\sigma _{0,e}^2(\omega )`$ $`={\displaystyle \frac{1}{N_f}}{\displaystyle \underset{f=1}{\overset{N_f}{}}}\left|S_{f,e}(\omega )\overline{S}_{c,e}(\omega )\right|^2`$ (5)
where $`N_f`$ is the number of frames, $`\omega `$ is the discrete frequency index and the indices $`f`$, $`c`$, and $`e`$ refer to the frame, cluster and electrode respectively. For this initial averaging the time shifts, $`\tau _f`$, described below, are set to zero for all frames.
* Loop 2: Split clusters- In every pass, except the first, each cluster is split in two using a small random vector $`\delta _{c,e}`$, i.e., $`\overline{S}_{c,e}\overline{S}_{c,e}\pm \overline{\delta }_{c,e}`$. This is repeated until some criterion is met. Establishing this criterion proved to be a difficult problem to solve generally. An a priori knowledge of the expected number of cells was found to be the most reliable criterion for stopping the cluster splitting.
* Loop 3: Reassign frames- This loop executes an expectation maximization algorithm. The average and variance at each frequency component of every cluster are calculated using the appropriate time delay for every member frame calculated against each cluster. The time delay $`\tau _f`$ that minimizes $`\chi _{f,c}^2(\tau _f)`$ is calculated against all clusters. The frame is then reassigned to the cluster that matches it most closely, i.e., the one that yields the minimum

$`\chi _{f,c}^2(\tau _f)`$. This is done until the frames are distributed in such a way that these parameters no longer change \[Papoulis, 1965\].
* Line 4: Finalize clusters- After the clustering ends some clusters may actually be identical within statistical error. These are merged into a single cluster. Others might clearly contain frames of different types and are split into two clusters. Still others might contain frames that clearly do not contain spikes or perhaps contain very few spikes. These are discarded.

### 3.2 Detection

The end result of the clustering phase is a statistical template for all of the spike types found in the data. In the next and final phase the data is scanned to find all occurrences of each spike type. The basic idea of the detection is to cut the data into short frames and determine which cluster best describes the data in that frame and its precise timing. Spikes that are well centered in the frames are detected and subtracted from the data. The data is then reframed and the process is repeated until no single spikes remain in the data. Finally, this process is repeated for overlaps. This process is shown schematically in Figure (1).

Single spike detection- Single spike detection is done by finding the $`(c,\tau _f)`$ pair that minimizes

$$\chi _f^2(c,\tau _f)=\underset{\omega ,e}{}\frac{1}{\sigma _{c,e}^2(\omega )}\left|S_{f,e}(\omega )\overline{S}_{c,e}(\omega )e^{i\omega \tau _f}\right|^2$$
(6)
for every frame. In practice Equation (6) is expanded into three terms,

$`\chi _f^2(c,\tau _f)`$ $`=A(c)+B(c)+2C(c,\tau _f),`$ (7)
$`A(c)`$ $`=\underset{\omega ,e}{{\displaystyle }}{\displaystyle \frac{1}{\sigma _{c,e}^2(\omega )}}|\overline{S}_{c,e}(\omega )|^2,`$ (8)
$`B(c)`$ $`=\underset{\omega ,e}{{\displaystyle }}{\displaystyle \frac{1}{\sigma _{c,e}^2(\omega )}}|S_{f,e}(\omega )|^2,`$ (9)
$`C(c,\tau _f)`$ $`=\underset{\omega ,e}{{\displaystyle }}{\displaystyle \frac{1}{\sigma _{c,e}^2(\omega )}}\mathrm{Re}\left(S_{f,e}(\omega )\overline{S}_{c,e}^{}(\omega )e^{i\omega \tau _f}\right).`$ (10)
$`A(c)`$ can be calculated in advance once the frame clustering has been done. $`B(c,\tau _f)`$ and $`C(c,\tau _f)`$ are calculated in the detection algorithm, outlined in pseudo-code in Box (2). The main ideas are described below.

* Line 1: Frame data- The complete data set is broken up into nonoverlapping equal sized frames.
* Loop 2: Process each frame- It is assumed that a frame can contain either noise (see next section), a spike, part of a spike or an overlap of 2 spikes. Frames containing noise are discarded leaving only those containing single or multiple spike events.
* Loop 3: Check fit to each cluster- Here the chosen frame is compared to each cluster. Since the spike in the frame may not be centered it is necessary to align the frame to the cluster. For each frame the time delay, $`\tau _f`$, that maximizes the cross correlation term, $`C(c,\tau _f)`$, is found for each cluster.
* Line 4: Calculate $`\chi _{𝐟,𝐜}^\mathrm{𝟐}(\tau _𝐟,𝐜)`$ for all clusters- If the spike is not near the edge of the frame, i.e., $`\left|\tau _f\right|<\tau _{th}`$, then $`B(c,\tau _f)`$ is calculated yielding the final term in Equation (6). All of the terms of $`\chi _{f,c}^2(\tau _f,c)`$ have now been calculated yielding an estimate of the similarity of the frame to each of the clusters.
* Line 5: Finalize spike detection- The cluster that yields the smallest value for $`\chi _i^2(\tau ,c)`$ is the cluster most similar to the frame being tested. It is not enough to find the cluster for which $`\chi ^2`$ is smallest. A fit is accepted only if $`\chi _i^2(\tau ,c)<\chi _{th}^2`$. Frames that were not good matches to any cluster most likely contained an overlap, noise or a partial spike that will most likely be found when the frames are shifted in a later pass through the data.

Multiple spike (overlap) detection- Once all of the single spikes have been detected and removed, overlapping spikes are detected and removed as well. Here we generalize the single spike case to two spikes, thus, we look for the set of $`(\tau _1,\tau _2,c_1,c_2)`$ that minimizes

$$\chi ^2(\tau _{f,1},\tau _{f,2},c_1,c_2)=\underset{\omega ,e}{}\frac{1}{\sigma _{c_1,c_2,e}^2}\left|S_{f,e}(\omega )S_{c_1,e}(\omega )e^{i\omega \tau _{f,1}}S_{c_2,e}(\omega )e^{i\omega \tau _{f,2}}\right|^2.$$
(11)
The algorithm used to find the clusters and time delays is very similar to that used for the single spike case with some differences. First, the minimization executed on line 3 of Box (2) is over two variables $`\tau _{f,1}`$ and $`\tau _{f,2}`$. This is computationally intensive but since there are a finite number of possible time delays it remains manageable on a personal computer. As before many of the calculations can be performed once the cluster centers are known.

Another issue which needs to be addressed when dealing with two spikes is the calculation of $`\sigma _{c_1,c_2,e}^2(\omega )`$. Here we assume that there is an inherent noise in the spike shape and an additive background noise. Thus the variances are given by

$$\sigma _{c_1,c_2,e}^2\left(\omega \right)=\sigma _{c_1,e}^2\left(\omega \right)+\sigma _{c_2,e}^2\left(\omega \right)\sigma _n^2\left(\omega \right),$$
(12)
where $`\sigma _n^2(\omega )`$ is the background noise computed from regions devoid of spikes. While this assumption is not true in general it is reasonable since the variances of the different clusters are assumed to be independent. It proved to be a reliable working model. This measure of the two-spike variance takes into account the contribution of each cluster but does not overcount the background noise.

## 4 Application to multi-electrode data

The techniques described above were applied to data recorded from the escape response system of the American cockroach Periplaneta americana \[Camhi and Levy, 1989, Kolton and Camhi, 1995, Liebenthal et al., 1994, Westin et al., 1977, Camhi, 1984\]. Neural activity was recorded using 8 hook electrodes attached to the two bilateral abdominal connectives. In this arrangment each electrode measures a weighted sum of the activity from the different neurons in the connective. A more detailed description of the experimental setup is given in reference \[Rinberg and Davidowitz, 2002\]. Typical experiments lasted for several hours and yielded about 8 GB of raw data (see Figure (2)). The aim of the work described here is to unravel the individual spike times from this multi-electrode data.

### 4.1 Statistical model of the clusters

Here we describe the application of the clustering algorithm described in Section (3.1) to the multielectrode cockroach data.

Frame selection- Only frames that fulfill the following criteria were collected to produce a model of the clusters: a) the signal in the middle of the frame is above some threshold, $`v_{th}^{mid}`$, on at least one electrode and b) the signal at the frame edges is less than some other threshold, $`v_{th}^{edge}`$, on all of the electrodes. This idea is illustrated in Figure (3).

The thresholds were proportional to the background noise levels. This background signal, $`\overline{v}_b`$, was the average of several averages,

$$v_b=(\frac{1}{N_p}\underset{i=1}{\overset{N_p}{}}\left|v_i\overline{v}\right|^2)^{\frac{1}{2}},$$
(13)
calculated from regions devoid of spikes at the beginning, middle and end of the experiment (usually silences between trials when no stimulus was presented). Here $`v_i`$ are the data points while the number of data points was typically about $`N_p=510^4`$. The thresholds described above were defined as $`v_{th}^{mid}=4\overline{v}_b`$ and $`v_{th}^{edge}=1.5\overline{v}_b`$. These thresholds worked well for our data but will likely be different for data from other experiments.

While the choice of thresholds is to some extent arbitrary, there are several guidelines that can be followed. If the threshold is set too high low energy spikes will be missed. If it is set too low, a large cluster containing noise will appear. We have checked that lowering the thresholds does not influence the larger amplitude clusters identified by our algorithm. The setting of the threshold is thus a compromise between identifying all the small amplitude spikes and minimizing computing time.

To avoid missing spikes because of voltage drifts or overlapping tails of nearby spikes, the frames are first detrended by subtracting a linear fit to the first and last 0.5 ms of each frame.

Time shift between electrodes- Spikes appear on each electrode at different times because of the finite propagation time of the action potentials along the neurons in the connective. This is evident in the data shown in Figure (3). Frames are thus defined with time delays between the different channels. A further complication arises from the fact that the time delays on the different electrodes may be different for different cells. An average delay was found by calculating the time of the maximum of the cross correlation of voltage traces from the different electrodes. Because these delays can change during the course of an experiment these inter-channel delays were in turn averaged from widely separated segments of data. Typically, these delays were between 0.1 and 0.8 ms depending on the positioning of the electrodes.

Frame size- To keep calculation times short, a small frame size is desirable. On the other hand, the frames have to be big enough to account for the different propagation times of the different cells. For the experiments described here frames were 3.2 ms long (64 data points).

The clustering begins after 10,000 frames containing candidate spikes have been collected. While each frame is represented by the vector defined in Equation (3) only the low frequency components (the first 16 complex numbers) of each Fourier transform were used because the higher frequencies were found to be indistinguishable from noise. Thus $`S_{f,e}(\omega )`$ was of dimension $`\frac{1}{2}32N_e`$ complex numbers, where $`N_e`$ is the number of electrodes.

Since about 7 cells are expected on each side the clustering is stopped once 16 clusters have been found on a side. Up until this point the clustering is automatic. Some results of this automatic clustering are shown in Figure (4). Intervention is now needed to refine the clustering. To do this the $`\chi ^2`$ distributions of the clusters are examined manually. Clusters are then merged, split or discarded as necessary. This process could be automated, but was not deemed worth the effort. Results of the clustering after the final manual intervention are shown in Figure (5) and Figure (6).

Track slow changes- Thus far, the first $`10^4`$ frames have been used to establish the statistics of each cluster. The accuracy of the spike detection can be improved if the change in spike shapes are tracked during the course of a long experiment. For the cockroach experiments described here, the trial period of 100 s was chosen as the time scale upon which changes are tracked. To do this the cluster statistics $`\overline{S}_{c,e}(\omega )`$ and $`\sigma _{c,e}(\omega )`$ are updated with consecutive 100 s segments of data on a first-in, first-out basis. Once the cluster statistics have been computed, the next $`N`$ “clean” frames from the data are found and are each assigned to the cluster that is closest based on the $`\chi ^2`$ distance defined in Box (1). These frames are then added to these clusters while the same number of the earliest frames are removed. As usual, frames that exceed a certain $`\chi ^2`$ for all clusters are discarded. The cluster statistics are then recalculated, yielding dynamic clusters of time dependent spike statistics $`\overline{S}_{c,e}(\omega ,t)`$ and $`\sigma _{c,e}(\omega ,t)`$ that are used as the templates for the spike detection described in the next section. Results of spike shape tracking are shown in Figure (7).

### 4.2 Detection

For the detection phase the complete data set is broken into nonoverlapping frames containing 64 data points (3.2 ms). In practice, $`B(c)`$ is calculated only after zero padding the frame on the outside 2 ms. This is done to eliminate the possibility of a nearby spike, with an overlapping tail, causing $`\chi _{f,c}^2(\tau _f,c)`$ to be too large. $`C(c,\tau _f)`$ does not need to be recalculated because the clusters themselves are averages of many clean spikes and therefore inherently zero padded.

A fit is accepted only if $`\chi _i^2(\tau ,c)<\chi _{th}^2`$ which was determined by the shape of the distribution. Typically $`\chi _{th}^22`$, though this threshold would likely be different for other data. The probability of detecting a spike is smaller if it overlaps with the tails of nearby spikes. This problem can be greatly reduced if, once spikes have been classified, the average spike from the corresponding cluster is subtracted from the data (see Figure (1)). Every pass through the data leaves a smaller number of unclassified spikes. After each pass the data is reframed with a time shift of $`\frac{1}{4}`$ frame and the process is repeated. Ideally, after 4 passes all single spike events have been detected.

Multiple spike (overlap)detection- Once all of the single spikes have been detected and removed, overlapping spikes are detected and removed as well. Figures (8 and 9) show the results of describing an overlap as single and double spike events. Notice that the automatic recognition of an event as being a two spike event is dependent on being able to reject the “best” single spike description of the event. Thus, in Figure (8), the best single spike description is quite good, and the identification of a spike from cell (b) is correct, but the value of $`\chi ^2=3.5`$ is unacceptably high, as can be seen from the $`\chi ^2`$ distributions. Once we explore the space of two spike events we find a unique description with $`\chi ^21`$ (see Figure (9)).

## 5 Discussion

The spike sorting method described here has several key features. First, the experimental design is such that the full waveforms from all electrodes are recorded continuously during the course of an experiment. This is an advantage over many spike sorting techniques that are based on the idea of feature clustering (see \[Lewicki, 1998\] for a recent review of spike sorting techniques as well as an older review in \[Schmidt, 1984\]). In feature clustering one or more features of the spikes (spike height, peak to peak amplitude, rise times, etc.) are extracted and clustered. Since only a few features of the spikes are used, subtle differences between spikes from different cells can be lost. In addition, it is not a priori clear which, and how many, features to use—are there better (optimal) ones? When applied properly this technique can work well for a small number of cells but does not work well with many cells.

Another important feature of the present technique is that the construction of models for the spike shapes from individual cells (clustering) is separated from the problem of finding those spikes in the data (detection). The advantage here is that one has the leisure to look for clean examples of each spike type and thus to build a good statistical model of each spike type. Only then does one look for the occurrence of these spikes in the data. These models are time dependent in the sense that they track the change in spike shape during the course of an experiment. This has traditionally been a problem in template-matching spike sorting techniques in which a model of the spike shapes are constructed \[Bergman and DeLong, 1992\]. In this technique, these models are compared to a given spike and a decision is made as to whether it belongs to the class defined by this template or not. One of the more advanced implementations of this technique \[Lewicki, 1994\] works well in many cases but relies heavily on a Gaussian model of the noise. Another technique \[Fee et al., 1996a\] that makes no assuptions about the waveform noise uses the spike shapes as well as refractory period statistics to classify cells. It works well with bursting cells but in its present form this technique does not treat overlaps well. This is another advantage of the separation of clustering from detection: it simplifies the problem of overlap detection which has traditionally been one of the most difficult parts of the spike sorting problem.

Many comparisons and computations that go into the sorting process are much easier when working in the frequency domain. The ease in which temporal alignment is achieved is one of the advantages of performing the detection in the frequency domain since sub-bin time shifts in the time domain would require resampling \[Lewicki, 1998\]. In the frequency domain this is done simply by multiplying the Fourier components by $`e^{i\omega \tau }`$ which is equivalent to a time shift of $`\tau `$ in the time domain. In addition, it is even possible to sort spikes from bursting cells (in which spike shapes can change drastically over short time scales) by ignoring the spike shapes entirely and sorting on transfer function ratios \[Rinberg et al., 1999\]. Since this is an independent type of information it can be used to check the validity of sorting results. Sorting on the full Fourier transforms of the voltage waveforms yielded excellent results for our data but other investigators will likely find other combinations that work better for their data.

The statistical methods used in this sorting program allow us to decide whether the clustered spikes really are discriminable. In addition, this statistical analysis allows us to develop more rigorous criteria for accepting or rejecting the possible detection of a spike.

Finally we note that in a typical experiment $`400,000`$ spikes are found. Of these 90% were found with single spike detection, 8% were found with overlap detection only 2% were events that could not be classified.

## 6 Summary

In this paper the problem of unraveling multielectrode neural data has been addressed. Special attention has been paid to the detection of overlapping spikes which poses obvious difficulties for any sorting method. The problem of spike sorting has been separated into two independent parts. First, a statistical model of all possible spikes found in the data is constructed. Only then are these spikes detected in the data using strict statistical criteria to quantify the quality of this detection. Overlaps are dealt with after all possible single spike events have been detected. These techniques were applied to multielectrode data from the American cockroach with good results.

## 7 Acknowledgements

The authors would like to thank N. Brenner and R. de Ruyter van Steveninck for many helpful suggestions and discussions.

1 set $`\tau _f:=0`$ and calculate initial $`\overline{S}_{c,e}(\omega )`$ and $`\sigma _{c,e}(\omega )`$
2 repeat split clusters
3 repeat update frame assignment
find average and variance for each cluster using $`\tau _f`$     $`\overline{S}_{c,e}(\omega )`$ $`:={\displaystyle \frac{1}{N_c}}\underset{fc}{{\displaystyle }}S_{f,e}(\omega )e^{i\omega \tau _f}`$    $`\sigma _{c,e}(\omega )`$ $`:={\displaystyle \frac{1}{N_c}}\underset{fc}{{\displaystyle }}\left|S_{f,e}(\omega )e^{i\omega \tau _f}\overline{S}_{c,e}(\omega )\right|^2`$   for every frame, $`f`$, find $`\tau _f`$ that minimizes   $`\chi _{f,c}^2(\tau _f)`$ $`:=\underset{\omega ,e}{{\displaystyle }}{\displaystyle \frac{1}{\sigma _e^2(\omega )}}\left|S_{f,e}(\omega )e^{i\omega \tau _f}\overline{S}_{c,e}(\omega )\right|^2`$     assign frame to the cluster that yields smallest $`\chi _{f,c}^2(\tau _f)`$
until clusters are stationary
until number of clusters $`>`$ expected number of cells
4 finalize clustering by merging, splitting or removing clusters
5 repeat track slow changes in spike shapes
find spikes in next time chunk of data and update cluster statistics
until end of data is reached

Box 1: Pseudo code outline of the frame clustering algorithm. See text for details.

1 Frame whole data set with no overlaps
2 repeat get next frame, $`f`$
if frame has signal above background then
3 repeat for each cluster, $`c`$
find time delay, $`\tau _f`$, that maximizes
$$C(c,\tau _f)=\underset{\omega ,e}{}\frac{\mathrm{Re}\left(S_{f,e}(\omega )\overline{S}_{c,e}^{}(\omega )e^{i\omega \tau _f}\right)}{\sigma _{c,e}^2(\omega )}$$
4 if $`\left|\tau _f\right|<\tau _{th}`$ then zero pad frame and calculate $`\chi _f^2(\tau _f,c)`$
until all clusters checked
5 find cluster, $`c`$, that yields smallest $`\chi _{f,c}^2(\tau _f,c)`$
if $`\chi _f^2(\tau _f,c)<\chi _{th}^2`$ keep this spike and subtract it from data
until end of data reached

Box 2: Pseudo code outline of the single spike detection algorithm. See text for details.

Figure captions


format: text

text:
# A Potential Energy Landscape Study of the Amorphous-Amorphous Transformation in H2O

## Acknowledgments

We thank NSF Chemistry Program, MIUR Cofin 2002 and Firb and INFM Pra GenFdt for support and SHARCNET and BU Computation Center for a generous allocation of CPU time. FS thanks UWO for its hospitality.


text:
# Circle actions and 𝐙/𝑘-manifolds

## Abstract

We establish an $`S^1`$-equivariant index theorem for Dirac operators on $`𝐙/k`$-manifolds. As an application, we generalize the Atiyah-Hirzebruch vanishing theorem for $`S^1`$-actions on closed spin manifolds to the case of $`𝐙/k`$-manifolds.

Résumé français On établit un théorème d’indice $`S^1`$-équivariant pour les opérateurs de Dirac sur des $`𝐙/k`$ variétés. On donne une application de ce résultat, qui généralise le théorème d’Atiyah-Hirzebruch sur les actions de $`S^1`$ aux $`𝐙/k`$ variétés.

Titre français Actions du cercle et $`𝐙/k`$ variétés.

§1. $`S^1`$-actions and the vanishing theorem

Let $`X`$ be a closed connected smooth spin manifold admitting a non-trivial circle action. A classical theorem of Atiyah and Hirzebruch \[AH\] states that $`\widehat{A}(X)=0`$, where $`\widehat{A}(X)`$ is the Hirzebruch $`\widehat{A}`$-genus of $`X`$. In this Note we present an extension of the above result to the case of $`𝐙/k`$-manifolds, which were introduced by Sullivan in his studies of geometric topology. We recall the basic definition for completeness (cf. \[F\]).

Definition 1.1 A compact connected $`𝐙/k`$-manifold is a compact manifold $`X`$ with boundary $`X`$, which admits a decomposition $`X=_{i=1}^k(X)_i`$ into $`k`$ disjoint manifolds and $`k`$ diffeomorphism $`\pi _i:(X)_iY`$ to a closed manifold $`Y`$.

Let $`\pi :XY`$ be the induced map. In what follows, we will call an object $`\alpha `$ (e.g., metrics, connections, etc.) of $`X`$ a $`𝐙/k`$-object if there will be a corresponding object $`\beta `$ on $`Y`$ such that $`\alpha |_X=\pi ^{}\beta `$. We make the assumption that $`X`$ is $`𝐙/k`$ oriented, $`𝐙/k`$ spin and is of even dimension.

Let $`g^{TX}`$ be a $`𝐙/k`$ Riemannian metric of $`X`$ which is of product structure near $`X`$. Let $`R^{TX}`$ be the curvature of the Levi-Civita connection associated to $`g^{TX}`$. Let $`E`$ be a $`𝐙/k`$ complex vector bundle over $`X`$. Let $`g^E`$ be a $`𝐙/k`$ Hermitian metric on $`E`$ which is a product metric near $`X`$. Let $`^E`$ be a $`𝐙/k`$ connection on $`E`$ preserving $`g^E`$ such that $`^E`$ is of product structure near $`X`$. Let $`R^E`$ be the curvature of $`^E`$. Let $`D_+^E:\mathrm{\Gamma }(S_+(TX)E)\mathrm{\Gamma }(S_{}(TX)E)`$ be the associated Dirac operator on $`X`$ and $`D_{+,X}^E`$ (and then $`D_Y^E`$) be its induced Dirac operator on $`X`$ (and then on $`Y`$). Let $`\overline{\eta }(D_Y^E)`$ be the reduced $`\eta `$-invariant of $`D_Y^E`$ in the sense of \[APS\]. Then

$$\widehat{A}_{(k)}(X,E)=_X\stackrel{1/2}{det}\left(\frac{\sqrt{1}R^{TX}/4\pi }{\mathrm{sinh}(\sqrt{1}R^{TX}/4\pi )}\right)\mathrm{tr}\left[e^{\frac{\sqrt{1}}{2\pi }R^E}\right]k\overline{\eta }(D_Y^E)\mathrm{mod}k𝐙$$
$`(1.1)`$
does not depend on ($`g^{TX}`$, $`g^E`$, $`^E`$) and determines a topological invariant in $`𝐙/k𝐙`$ (cf. \[APS\] and \[F\]). Moreover, Freed and Melrose \[FM\] have proved a mod $`k`$ index theorem, giving $`\widehat{A}_{(k)}(X,E)𝐙/k𝐙`$ a purely topological interpretation. When $`E=𝐂`$ is the trivial vector bundle over $`X`$, we usually omit the superscript $`E`$.

Theorem 1.2 If $`X`$ admits a nontrivial $`𝐙/k`$ circle action preserving the orientation and the Spin structure on $`TX`$, then $`\widehat{A}_{(k)}(X)=0`$. Moreover, the equivariant mod $`k`$ index in the sense of Freed and Melrose vanishes.

It turns out that the original method in \[AH\] is difficult to extend to the case of manifolds with boundary to prove Theorem 1.2. Thus we will instead make use of an extension of the method of Witten \[W\]. Analytic localization techniques developed by Bismut-Lebeau \[BL, Sect. 9\] and their extensions to manifolds with boundary developed in \[DZ\] play important roles in our proof.

§2. A mod $`k`$ localization formula for circle actions

We make the assumption that the $`𝐙/k`$ circle action on $`X`$ lifts to a $`𝐙/k`$ circle action on $`E`$. Without loss of generality, we may and we will assume that this $`𝐙/k`$ circle action preserves $`g^{TX}`$, $`g^E`$ and $`^E`$. Let $`D_{+,APS}^E:\mathrm{\Gamma }(S_+(TX)E)\mathrm{\Gamma }(S_{}(TX)E)`$ be the elliptic operator obtained by imposing the standard Atiyah-Patodi-Singer boundary condition \[APS\] on $`D_+^E`$.

Let $`H`$ be the Killing vector field on $`X`$ generated by the $`S^1`$ action on $`X`$. Then $`H|_XX`$ induces a Killing vector field $`H_Y`$ on $`Y`$. Let $`_H`$ denote the corresponding Lie derivative acting on $`\mathrm{\Gamma }(S_\pm (TX)E)`$. Then $`_H`$ commutes with $`D_{+,APS}^E`$.

For any $`n𝐙`$, let $`F_\pm ^n`$ be the eigenspaces of $`\mathrm{\Gamma }(S_\pm (TX)E)`$ with respect to the eigenvalue $`2\pi n`$ of $`\frac{1}{\sqrt{1}}_H`$. Let $`D_{+,APS}^E(n):F_+^nF_{}^n`$ be the restriction of $`D_{+,APS}^E`$ on $`F_+^n`$. Then $`D_{+,APS}^E(n)`$ is Fredholm. We denote its index by $`\mathrm{ind}(D_{+,APS}^E(n))𝐙`$.

Let $`X_H`$ (resp. $`Y_H`$) be the zero set of $`H`$ (resp. $`H_Y`$) on $`X`$ (resp. $`Y`$). Then $`X_H`$ is a $`𝐙/k`$-manifold and there is a canonical map $`\pi _{X_H}:X_HY_H`$ induced from $`\pi `$. We fix a connected component $`X_{H,\alpha }`$ of $`X_H`$, and we omit the subscript $`\alpha `$ if there is no confusion.

We identify the normal bundle to $`X_H`$ in $`X`$ to the orthogonal complement of $`TX_H`$ in $`TX|_{X_H}`$. Then $`TX|_{X_H}`$ admits an $`S^1`$-invariant orthogonal decomposition $`TX|_{X_H}=N_{m_1}\mathrm{}N_{m_l}TX_H`$, where each $`N_\gamma `$, $`\gamma 𝐙`$, is a complex vector bundle on which $`gS^1𝐂`$ acts by multiplication by $`g^\gamma `$. By using the same notation as in \[LMZ, (1.8)\], we simply write that $`TX|_{X_H}=_{v0}N_vTX_H.`$ Similarly, let $`E|_{X_H}`$ admits the $`S^1`$-invariant decomposition $`E|_{X_H}=_vE_v.`$

Let $`S(TX_H,(detN)^1)`$ be the complex spinor bundle over $`X_H`$ associated to the canonically induced Spinc structure on $`TX_H`$. It is a $`𝐙/k`$ Hermitian vector bundle and carries a canonically induced $`𝐙/k`$ Hermitian connection.

Recall that by \[AH, 2.4\], one has $`_vvdimN_v0\mathrm{mod}2𝐙`$. Following \[LMZ, (1.15)\], set

$$R(q)=q^{\frac{1}{2}_v|v|dimN_v}_{v>0}\left(\mathrm{Sym}_{q^v}(N_v)detN_v\right)_{v<0}\mathrm{Sym}_{q^v}\left(\overline{N}_v\right)\underset{v}{}q^vE_v=_nR_nq^n,$$

$$R^{}(q)=q^{\frac{1}{2}_v|v|dimN_v}_{v>0}\mathrm{Sym}_{q^v}\left(\overline{N}_v\right)_{v<0}\left(\mathrm{Sym}_{q^v}(N_v)detN_v\right)\underset{v}{}q^vE_v=_nR_n^{}q^n.$$

Then each $`R_n`$ (resp. $`R_n^{}`$) is a $`𝐙/k`$ Hermitian vector bundle over $`X_H`$ carrying a canonically induced $`𝐙/k`$ Hermitian connection. For any $`n𝐙`$, let

$$D_{X_H,+}^{R_n}:\mathrm{\Gamma }(S_+(TX_H,(detN)^1)R_n)\mathrm{\Gamma }(S_{}(TX_H,(detN)^1)R_n)$$

be the canonical twisted Spinc Dirac operator on $`X_H`$. Let $`D_{X_H,+,APS}^{R_n}`$ be the corresponding elliptic operator associated to the Atiyah-Patodi-Singer boundary condition \[APS\]. We will use similar notation for $`R_n^{}`$.

Theorem 2.1 For any integer $`n𝐙`$, the following identities hold,

$$\mathrm{ind}D_{+,APS}^E(n)\underset{\alpha }{}(1)^{_{0<v}dimN_v}\mathrm{ind}D_{X_{H,\alpha },+,APS}^{R_n}\mathrm{mod}k𝐙,$$
$`(2.1)`$
$$\mathrm{ind}D_{+,APS}^E(n)\underset{\alpha }{}(1)^{_{v<0}dimN_v}\mathrm{ind}D_{X_{H,\alpha }+,APS}^{R_n^{}}\mathrm{mod}k𝐙.$$
$`(2.2)`$
Proof. For any $`T𝐑`$, following Witten \[W\], let $`D_{T,+}^E:\mathrm{\Gamma }(S_+(TX)E)\mathrm{\Gamma }(S_{}(TX)E)`$ be the Dirac type operator defined by $`D_{T,+}^E=D_+^E+\sqrt{1}Tc(H)`$. Let $`D_{T,+,APS}^E`$ be the corresponding elliptic operator associated to the Atiyah-Patodi-Singer boundary condition \[APS\]. Clearly, $`D_{T,+,APS}^E`$ also commutes with the $`S^1`$-action. For any integer $`n`$, let $`D_{T,+,APS}^E(n)`$ be the restriction of $`D_{T,+,APS}^E`$ on $`F_+^n`$. Then $`D_{T,+,APS}^E(n)`$ is still Fredholm. By an easy extension of \[DZ, Theorem 1.2\] to the current equivariant and $`𝐙/k`$ situation, one sees that $`\mathrm{ind}(D_{T,+,APS}^E(n))\mathrm{mod}k𝐙`$ does not depend on $`T𝐑`$ (compare with \[TZ, Theorem 4.2\]).

Let $`D_{T,+,X}^E:\mathrm{\Gamma }((S_+(TX)E)|_X)\mathrm{\Gamma }((S_+(TX)E)|_X)`$ be the induced Dirac type operator of $`D_{T,+}^E`$ on $`X`$. For any integer $`n`$, let $`D_{T,+,X}^E(n):F_+^n|_XF_+^n|_X`$ be the restriction of $`D_{T,+,X}^E`$ on $`F_+^n|_X`$. Also, the induced Dirac operators $`D_{+,X_H}^{R_n}`$ and $`D_{Y_H}^{R_n}`$ can be defined in the same way as in Section 1.

Let $`a_n>0`$ be such that $`\mathrm{Spec}(D_{Y_H}^{R_n})[2a_n,2a_n]\{0\}.`$ By combining the techniques in \[BL, Sect. 9\], \[BZ, Sect. 4b)\] and \[LMZ, Sect. 1.2\], one can prove the following analogue of \[BZ, Theorem 3.9\], stating that there exists $`T_1>0`$ such that for any $`TT_1`$,

$$\mathrm{\#}\{\lambda \mathrm{Spect}(D_{T,+,X}^E(n)):a_n\lambda a_n\}=dim(\mathrm{ker}D_{+,X_H}^{R_n})=kdim(\mathrm{ker}D_{Y_H}^{R_n}).$$
$`(2.3)`$
If $`dim(\mathrm{ker}D_{Y_H}^{R_n})=0`$, then by (2.3), one sees that when $`TT_1`$, $`D_{T,+,X}^E(n)`$ is invertible. Then $`\mathrm{ind}(D_{T,+,APS}^E(n))`$ itself does not depend on $`TT_1`$. Moreover, by combining the techniques in \[LMZ, Sect. 1.2\] and \[DZ, Sect. 3\], one can further prove that there exists $`T_2>0`$ such that when $`TT_2`$,

$$\mathrm{ind}(D_{T,+,APS}^E(n))=\underset{\alpha }{}(1)^{_{0<v}dimN_v}\mathrm{ind}D_{X_{H,\alpha },+,APS}^{R_n}$$
$`(2.4)`$
(compare with \[DZ, (2.13)\]). From (2.4) and the mod $`k`$ invariance of $`\mathrm{ind}(D_{T,+,APS}^E(n))`$ with respect to $`T𝐑`$, one gets (2.1).

In general, $`dim(\mathrm{ker}D_{Y_H}^{R_n})`$ need not be zero, and the eigenvalues of $`D_{T,+,X}^E(n)`$ lying in $`[a_n,a_n]`$ are not easy to control. Thus the above arguments no longer apply directly. Instead, we observe that $`dim(\mathrm{ker}(D_{Y_H}^{R_n}a_n))=0`$, and we use the method in \[DZ\] to perturb the Dirac type operators under consideration.

To do this, let $`\epsilon >0`$ be sufficiently small so that $`g^{TX}`$, $`g^E`$ and $`^E`$ are of product structure on $`[0,\epsilon ]\times XX`$. Let $`f:X𝐑`$ be an $`S^1`$-invariant smooth function such that $`f1`$ on $`[0,\epsilon /3]\times X`$ and $`f0`$ outside of $`[0,2\epsilon /3]\times X`$. Let $`r`$ denote the parameter in $`[0,\epsilon ]`$. Let $`D_{X_H,a_n,+}^{R_n}`$ be the Dirac type operator acting on $`\mathrm{\Gamma }(S_+(TX_H,(detN)^1)R_n)`$ defined by $`D_{X_H,a_n,+}^{R_n}=D_{X_H,+}^{R_n}a_nfc(\frac{}{r}).`$ Let $`D_{X_H,a_n,+,APS}^{R_n}`$ be the corresponding elliptic operator associated to the Atiyah-Patodi-Singer boundary condition \[APS\]. By an easy extension of \[DZ, Theorem 1.2\] (compare with \[TZ, Theorem 4.2\]), we see that,

$$\underset{\alpha }{}(1)^{_{0<v}dimN_v}\mathrm{ind}D_{X_{H,\alpha },a_n,+,APS}^{R_n}\underset{\alpha }{}(1)^{_{0<v}dimN_v}\mathrm{ind}D_{X_{H,\alpha },+,APS}^{R_n}\mathrm{mod}k𝐙.$$
$`(2.5)`$
For any $`T𝐑`$, let $`D_{T,a_n,+}^E:\mathrm{\Gamma }(S_+(TX)E)\mathrm{\Gamma }(S_{}(TX)E)`$ be the Dirac type operator defined by $`D_{T,a_n,+}^E=D_{T,+}^Ea_nfc(\frac{}{r})`$. Let $`D_{T,a_n,+,APS}^E`$ be the corresponding elliptic operator associated to the Atiyah-Patodi-Singer boundary condition. Let $`D_{T,a_n,+,APS}^E(n)`$ be its restriction on $`F_+^n`$. Then $`D_{T,a_n,+,APS}^E(n)`$ is still Fredholm. By another extension of \[DZ, Theorem 1.2\], one has

$$\mathrm{ind}D_{T,a_n,+,APS}^E(n)\mathrm{ind}D_{T,+,APS}^E(n)\mathrm{mod}k𝐙.$$
$`(2.6)`$
Moreover, since $`D_{Y_H}^{R_n}a_n`$, which is the induced Dirac type operator from $`D_{X_H,a_n,+}^{R_n}`$ through $`\pi _{X_H}`$, is invertible, by combining the arguments in \[LMZ, Sect. 1.2\] with those in \[DZ, Sect. 3\], one deduces that there exists $`T_3>0`$ such that for any $`TT_3`$, the following analogue of (2.4) holds,

$$\mathrm{ind}D_{T,a_n,+,APS}^E(n)=\underset{\alpha }{}(1)^{_{0<v}dimN_v}\mathrm{ind}D_{X_{H,\alpha },a_n,+,APS}^{R_n}.$$
$`(2.7)`$
From (2.5)-(2.7) and the mod $`k`$ invariance of $`\mathrm{ind}(D_{T,+,APS}^E(n))`$ with respect to $`T𝐑`$, one gets (2.1).

Similarly, by taking $`T\mathrm{}`$, one gets (2.2). $`\mathrm{}`$

§3. Proof of Theorem 1.2

We apply Theorem 2.1 to the case $`E=𝐂`$.

First, if $`X_H=\mathrm{}`$, by Theorem 2.1, it is obvious that for each $`n𝐙`$,

$$\mathrm{ind}(D_{+,APS}(n))0\mathrm{mod}k𝐙.$$
$`(3.1)`$
When $`X_H\mathrm{}`$, we see that $`_v|v|dimN_v>0`$ (i.e., at least one of the $`N_v`$’s is nonzero) on each connected component of $`X_H`$. Then by (2.1) and by the definition of the $`R_n`$’s, we deduce that for any integer $`n0`$, (3.1) holds. Similarly, by (2.2) and by the definition of the $`R_n^{}`$’s, one deduces that (3.1) holds for any integer $`n0`$.

In summary, for any $`n𝐙`$, (3.1) holds.

From (1.1), (3.1), by the Atiyah-Patodi-Singer index theorem \[APS\], and using the obvious fact that $`\mathrm{ind}(D_{+,APS})=_n\mathrm{ind}(D_{+,APS}(n)),`$ one gets $`\widehat{A}_{(k)}(X)=0`$. $`\mathrm{}`$

Remark 3.1 By combining Theorem 2.1 with the arguments in \[LMZ, Sects. 2-4\], one should be able prove an extension of the Witten rigidity theorem, of which a $`K`$-theoretic version has been worked out in \[LMZ\], to $`𝐙/k`$-manifolds. This, together with some other consequences of Theorem 1.2, will be carried out elsewhere.

Acknowledgements The author would like to thank Xiaonan Ma for helpful conversations.

References

\[AH\] M. F Atiyah and F. Hirzebruch, Spin manifolds and groups actions, Essays on Topology and Related Topics, Mémoireé dédié à Georges de Rham. (Ed. A. Haefliger and R. Narasimhan), Springer-Verlag, (1970), 18-28.

\[APS\] M. F. Atiyah, V. K. Patodi and I. M. Singer, Spectral asymmetry and Riemannian geometry I. Proc. Cambridge Philos. Soc. 77 (1975), 43-69.

\[BL\] J.-M. Bismut and G. Lebeau, Complex immersions and Quillen metrics, Publ. Math. IHES. 74 (1991), 1-297.

\[BZ\] J.-M. Bismut and W. Zhang, Real embeddings and eta invariant. Math. Ann. 295 (1993), 661-684.

\[DZ\] X. Dai and W. Zhang, Real embeddings and the Atiyah-Patodi-Singer index theorem for Dirac operators. Asian J. Math. 4 (2000), 775-794.

\[F\] D. S. Freed, $`𝐙/k`$-manifolds and families of Dirac operators. Invent. Math. 92 (1988), 243-254.

\[FM\] D. S. Freed and R. B. Melrose, A mod $`k`$ index theorem. Invent. Math. 107 (1992), 283-299.

\[LMZ\] K. Liu, X. Ma and W. Zhang, Rigidity and vanishing theorems in $`K`$-theory. Commun. Anal. Geom. 11 (2003), 121-180.

\[TZ\] Y. Tian and W. Zhang, Quantization formula for symplectic manifolds with boundary. Geom. Funct. Anal. 9 (1999), 596-640.

\[W\] E. Witten, Fermion quantum numbers in Kaluza-Klein theory. Shelter Island II: Proceedings of the 1983 Shelter Island Conference on Quantum Field theory and the Fundamental Problems of Physics. (Ed. R. Jackiw et. al.), The MIT Press, (1985), 227-277.


text:
# On the gravitational collapse in anti-de Sitter space-time

## Abstract

We study the semiclassical evolution of a self-gravitating thick shell in Anti-de Sitter space-time. We treat the matter on the shell as made of quantized bosons and evaluate the back-reaction of the loss of gravitational energy which is radiated away as a non-adiabatic effect. A peculiar feature of anti-de Sitter is that such an emission also occurs for large shell radius, contrary to the asymptotically flat case.

thanks: E-mail: alberghi@bo.infn.itthanks: E-mail: casadio@bo.infn.it

One of the main issues arising in quantum gravity is the search for a unitary description of the process of black hole formation and evaporation . The conjectured AdS-CFT correspondence may lead to a way of solving this puzzle, as it provides a gauge theory (hence unitary) description for processes occurring in string theory and asymptotically anti-de Sitter (AdS) space-times . So far, however, it has been difficult to find a concrete description of such a process within this framework. For this reason, one is led to examine simplified models, for black hole formation, such as the collapse of a spherical shell of matter, in order to get an insight of what a full quantum gravity description would be. This is the path followed in Refs. . The main unsolved issue in these analysis is the description of how the initial stages of the evolution, for a shell with a radius much greater then its horizon, may be related to the late stages of the collapse, when a black hole is about to be formed. In Ref. , it is suggested that a solution may be provided by taking into account the radiation emitted by the shell as it collapses. This is just what we are going to describe in the following.

In a series of papers , we have studied the semiclassical dynamics of a self-gravitating shell made of bosonic matter (modelled as a set of $`N1`$ “microshells”) which collapses in an asymptotically flat (Schwarzschild) space-time. The (average) radius of the shell is taken to be a function of (proper) time $`r=R(\tau )`$, and its thickness $`\delta `$ is small but finite ($`0<\delta R`$) throughout the evolution. The main result of Ref.  was that tidal forces acting inside the shell induce non-adiabatic (i.e., proportional to the contraction velocity) changes in the quantum mechanical states of the microshells and give rise to a probability of excitation which becomes appreciable as the shell approaches its own gravitational radius $`R_H=2M_s`$, where $`M_s`$ is the Arnowitt-Deser-Misner (ADM) mass of the shell 111With use units with $`c=1`$ and $`\mathrm{}_p^2=\mathrm{}G_\mathrm{N}`$ is the Planck length squared.. If the bosonic microshells are coupled to a radiation field, they will then emit the excess energy and the ADM mass $`M_s`$ will steadily decrease. This induces an effective (quantum) tension which will modify the trajectory of the shell radius. We explicitly determined the evolution as long as $`R`$ is not too close to $`R_H`$, where some of the approximations employed break down, and showed that the shell contracts at constant (terminal) velocity.

In the present paper we want to generalize our previous results considering a shell embedded in asymptotically AdS space-time with cosmological constant $`\mathrm{\Lambda }=3/\mathrm{}^2`$ (for the thermodynamics of such a case see Ref. ). As we shall show below, one expects strong non-adiabatic effects not just at small shell radius (when $`R2M_s`$ and the self-gravity of the shell dominates), but also for $`R\mathrm{}`$ (when it is the cosmological constant that mainly determines the shell motion). In particular, we will focus on $`R`$ large, in order to highlight the differences with respect to the case already analyzed in Refs. .

The Schwarzschild-AdS metrics inside and outside the shell at $`r=R(\tau )`$ can be written as

$`ds^2=f_{\mathrm{i}/\mathrm{o}}dt^2+f_{\mathrm{i}/\mathrm{o}}^1dr^2+r^2d\mathrm{\Omega }^2,`$ (1)
where $`d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$, $`\tau `$ is the shell proper time,

$`d\tau \sqrt{f_\mathrm{o}(R)}dt=\sqrt{f_\mathrm{i}(R)}dt,`$ (2)
and

$`f_{\mathrm{i}/\mathrm{o}}(r)=1{\displaystyle \frac{2M_{\mathrm{i}/\mathrm{o}}}{r}}+{\displaystyle \frac{r^2}{\mathrm{}^2}}.`$ (3)
with $`M_{\mathrm{i}/\mathrm{o}}`$ the ADM masses for $`r<R`$ and $`r>R`$ respectively. A straightforward application of the junction conditions  yields the equation of motion for the shell,

$`\dot{R}^2=\left({\displaystyle \frac{M_\mathrm{o}M_\mathrm{i}}{M}}\right)^2+{\displaystyle \frac{M_\mathrm{o}+M_\mathrm{i}}{R}}+{\displaystyle \frac{M^2}{4R^2}}{\displaystyle \frac{R^2}{\mathrm{}^2}}1,`$ (4)
where a dot denotes derivative with respect to $`\tau `$. The proper mass $`M`$, as well as $`M_{\mathrm{i}/\mathrm{o}}`$, may depend on $`\tau `$ . In particular, we shall consider $`M`$ constant 222As in Ref. , we are interested in the loss of gravitational energy from the shell. Allowing for a decrease of the proper mass $`M`$ would further strengthen the effect, and will be the subject of future works. and initial conditions at $`\tau =0`$ such that $`\dot{R}(0)=0`$ \[$`R(0)`$ equals the classical turning point\]. In order to compare with the asymptotically flat case (obtained as $`\mathrm{}\mathrm{}`$), we shall also set $`M_\mathrm{i}=0`$ and consider the range of parameters

$`\delta MM_s(0)\mathrm{}R(0),`$ (5)
where $`M_sM_\mathrm{o}(\tau )`$. Note also that $`R(0)\mathrm{}`$ requires $`M_s(0)/MR(0)/\mathrm{}`$.

Let us now note that the proper acceleration of the shell radius given by

$`\ddot{R}={\displaystyle \frac{M_s}{2R^2}}{\displaystyle \frac{M^2}{4R^3}}{\displaystyle \frac{R}{\mathrm{}^2}},`$ (6)
is dominated by the last term for large $`R`$. One then obtains that, for $`R\mathrm{}`$ and near the turning point, the acceleration of the shell is much larger than it would be in asymptotically flat space ($`\mathrm{}\mathrm{}`$), so that the contraction velocity increases faster and non-adiabatic effects might occur before $`R`$ approaches $`R_H`$.

The motion of a microshell in the system of $`N1`$ microshells can be determined from the junction equation (4). In particular, one can obtain a mean field Hamiltonian equation  for the relative displacement $`\overline{r}`$ of each microshell of proper mass $`m`$ with respect to the average radius $`R`$ of the shell, which takes the form

$`H{\displaystyle \frac{1}{2}}(M+m)\dot{R}^2+{\displaystyle \frac{1}{2}}\mu \dot{\overline{r}}^2+V_0+V_m(\overline{r})=0,`$ (7)
where $`\mu m(Mm)/M`$ is the effective mass of the microshell and $`V_0`$ contains all the terms which do not depend on $`\overline{r}`$. Since $`mM=Nm`$, it is consistent just to retain terms linear in $`m`$ and approximate $`\mu m`$. Further, we are interested in the case when the shell radius is large and its thickness small \[see Eq. (5)\]. The time evolution of $`R`$ is thus determined by

$`{\displaystyle \frac{1}{2}}M\dot{R}^2V_0,`$ (8)
and that of $`\overline{r}`$ by

$`{\displaystyle \frac{1}{2}}m\dot{\overline{r}}^2V_m(\overline{r}).`$ (9)
The potential, to first order in $`\overline{r}`$, is given by

$`V_m={\displaystyle \frac{M_sm}{2R}}{\displaystyle \frac{\overline{r}}{R}}\times \{\begin{array}{cc}\left(1{\displaystyle \frac{M^2}{2RM_s}}\right)\hfill & \overline{r}>+\frac{\delta }{2}\hfill \\ & \\ \left(1{\displaystyle \frac{M^2}{2RM_s}}\right)\hfill & \overline{r}>+\frac{\delta }{2}.\hfill \end{array}`$ (13)
Up to corrections of order $`\overline{r}^2/\mathrm{}^2`$, which are negligible in the range of parameters (5), the above expression is identical to the asymptotically flat case treated in Ref.  and represents the tidal force responsible for the microshell confinement. This is expected in a proper time formalism, since the binding potential $`V_m`$ just depends on the local geometry.

The Hamiltonian constraint (9) can be quantized and, via the Born-Oppenheimer approximation, on assuming negligible quantum fluctuations of the metric, it leads to the Hamilton-Jacobi equation (8) for the average radius $`R`$ and to a Schrödinger equation for the microshell 333Such a treatment is rather lengthy and we omit the details for the sake of brevity. For the general formalism see Ref. ; its application to the shell model can be found in Ref. .,

$`i\mathrm{}_\tau n=\left[{\displaystyle \frac{\mathrm{}^2}{2m}}_{\overline{r}}^2+V_m\right]n,`$ (14)
in which the potential $`V_m`$ is time-dependent (recall that $`R`$ and $`M_s`$ are in general functions of $`\tau `$). Such equation was analyzed in Ref.  and the low lying modes in the spectrum,

$`\tau ,\overline{r}n=e^{in\mathrm{\Omega }\tau }\mathrm{\Phi }_n(\overline{r}),`$ (15)
were obtained by means of quantum invariants . The low energy levels are given by

$`E_nE_0n\mathrm{}\mathrm{\Omega },`$ (16)
where the fundamental frequency is

$`\mathrm{\Omega }={\displaystyle \frac{1}{R}}\sqrt{{\displaystyle \frac{M_s}{\delta }}},`$ (17)
and the spatial width of the levels is

$`\delta \mathrm{}_m^{2/3}R^{1/3}\left({\displaystyle \frac{R}{M_s}}\right)^{1/3},`$ (18)
$`\mathrm{}_m=\mathrm{}/m=N\mathrm{}_p^2/M`$ being the Compton wavelength of the microshells.

Due to the time-dependence of $`V_m`$, the amplitude of excitation from the ground state $`\mathrm{\Phi }_0`$ to states of even quantum number $`\mathrm{\Phi }_{2n}`$ at the time $`\tau >0`$ is not zero. To leading order in $`\dot{R}`$, it is given by

$`A_{02n}(\tau )(i)^n{\displaystyle \frac{\sqrt{(2n)!}}{6^nn!}}\left({\displaystyle \frac{\delta }{M_s}}\right)^{n/2}\dot{R}^n.`$ (19)
This expression rapidly vanishes for large $`R`$ ($`M_s`$) in the asymptotically flat case. However, in AdS, $`|\dot{R}|`$ grows faster for $`R\mathrm{}`$ \[as can be inferred from Eq. (6)\] and leads to a non negligible amplitude in such a regime.

The microshells can be coupled to an external scalar field $`\phi `$ by an interaction Lagrangian density of the form

$`\widehat{L}_{\mathrm{int}}=e{\displaystyle \underset{n,n^{}n}{}}e^{i(nn^{})\mathrm{\Omega }\tau }\mathrm{\Phi }_n(\overline{r})\mathrm{\Phi }_n^{}(\overline{r})\widehat{\phi }(t,R+\overline{r}),`$ (20)
where $`e`$ is the coupling constant and $`\tau =\tau (t)`$ \[see Eq. (2)\]. The field $`\phi `$ satisfies the Klein-Gordon equation $`\mathrm{}\phi =0`$ everywhere in AdS. However, since we are interested in the flux of radiation outgoing from the shell, we just need solutions outside the shell and neglect the time-dependence of $`R`$ (consistently with the approximation in which we keep the leading term for small shell velocity). The scalar field for $`r>R`$ is given by a sum over the normal modes

$`t,r,\theta ,\varphi n,l,m={\displaystyle \frac{e^{i\omega _{nl}t}}{r}}u_{nlm}(r)Y_l^m(\theta ,\varphi ),`$ (21)
where the $`Y`$’s are standard spherical harmonics. The radial functions are more easily expressed in terms of the turtle-like coordinate $`dr_{}dr/f_\mathrm{o}`$, and read

$`u_{nlm}=c_{nl}\mathrm{cos}^2(r_{})\mathrm{sin}^{1+l}(r_{})P_n^{l+\frac{1}{2},\frac{3}{2}}\left(\mathrm{cos}(2r_{})\right),`$ (22)
where the $`c_{nl}`$’s are normalization constants, the $`P`$’s Jacobi polynomials and

$`\omega _{nl}=\sqrt{{\displaystyle \frac{f_\mathrm{o}(R)}{1+R^2/\mathrm{}^2}}}{\displaystyle \frac{2n+l+3}{\mathrm{}}}.`$ (23)
We may now compute the transition amplitude for the entire process during which a microshell is excited from the ground state $`\mathrm{\Phi }_0`$ to an excited state $`\mathrm{\Phi }_n`$ and then decays back to $`\mathrm{\Phi }_0`$ by emitting scalar quanta $`n,l,m`$,

$`_n(t,0){\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _0^t}𝑑t^{}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle 𝑑\overline{r}_i^{}\mathrm{\Phi }_0^{}(\overline{r}_i^{})\mathrm{\Phi }_0^{}(\overline{r}_i^{})n,0,0\widehat{L}_{\mathrm{int}}(t^{},\overline{r}_i^{})0,0,0},`$ (24)
where we enforced the spherical symmetry (so that only states with $`l=m=0`$ are included) and summed over the $`N`$ microshells located at $`R+\overline{r}_i`$. We further assume that $`A_{02}(t)`$ is the only non-negligible excitation amplitude (that is, we neglect the possibility that each microshell can be excited more than once) and consider a symmetrized state for the $`N`$ microshells. The probability of emitting an energy equal to $`2\mathrm{\Omega }`$ during the interval between $`t=0`$ (when $`\dot{R}=0`$) and $`t`$,

$`P(2\mathrm{\Omega };t){\displaystyle \underset{n^{}n^{\prime \prime }}{}}_n^{}^{}(t,0)_{n^{\prime \prime }}^{}(t,0),`$ (25)
contains the projection of the Wightman function onto states with $`l=m=0`$,

$`\widehat{\phi }(R^{})\widehat{\phi }(R^{\prime \prime })`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{e^{i\omega _{n0}(t^{}t^{\prime \prime })}}{2\omega _{n0}R^{}R^{\prime \prime }}}u_{n00}^{}(R_{}^{})u_{n00}^{}(R_{}^{\prime \prime })`$ (26)
$``$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{e^{i\omega _{n0}(t^{}t^{\prime \prime })}}{2\omega _{n0}R^{}R^{\prime \prime }}}s_n^2(R^{},R^{\prime \prime }),`$
and is finally given by

$`P(2\mathrm{\Omega };t){\displaystyle \frac{e^2N^2}{\mathrm{}^2}}|A_{02}|^2{\displaystyle \underset{n}{}}s_n^2(R,R){\displaystyle \frac{\mathrm{sin}^2\left[\left(2\mathrm{\Omega }\sqrt{f_\mathrm{o}}\omega _{n0}\right)t/2\right]}{\omega _{n0}^2\left(2\mathrm{\Omega }\sqrt{f_\mathrm{o}}\omega _{n0}\right)^2}},`$ (27)
In the above expression, all functions must be evaluated at the time $`t`$ and we also used the fact that the states $`\mathrm{\Phi }_n`$ have a very narrow width $`\delta `$ \[according to Eq. (5)\] and can be approximated as $`\mathrm{\Phi }_n\delta (\overline{r})`$. Note that the lapse $`t`$ should be chosen short enough so that $`R`$ remains approximately constant during the corresponding time interval, but long enough to average over transient effects 444This is the approximation which yields the Planckian spectrum at the Hawking temperature for $`R2M_s`$ ..

For sufficiently large $`t`$ (or $`\tau `$), one can define a probability per unit proper time

$`𝒫_{2\mathrm{\Omega }}(R,\dot{R}){\displaystyle \frac{P(2\mathrm{\Omega };\tau )}{\tau }}e^2\dot{R}^2F(R,M_s)\delta \left(2\mathrm{\Omega }\omega _{n0}/\sqrt{f_\mathrm{o}}\right),`$ (28)
where the function $`F`$ can be straightforwardly deduced from Eq. (27). The above expression shows that emissions occur when some frequency $`\omega _{n\mathrm{\hspace{0.17em}0}}/\sqrt{f_\mathrm{o}}`$ (as measured in the shell frame) of an AdS mode is close to the proper energy gap $`2\mathrm{\Omega }`$, that is

$`n^{3/2}{\displaystyle \frac{MM_s}{N\mathrm{}_pR}}.`$ (29)
From the probability (28) the energy emitted per unit time is straightforwardly obtained as

$`\dot{M}_s=8\pi R^2\mathrm{\Omega }𝒫_{2\mathrm{\Omega }}(R,\dot{R}),`$ (30)
which, together with the dynamical equation (4) constitutes the system of coupled equations that we solved numerically. In Fig. 1 we plot an example of the radius, velocity and ADM mass of the shell for a choice of the parameters particularly suited to obtain neat plots. It is clear that the velocity approaches a terminal value, although our approximations are not entirely reliable for too long times. For more realistic values of the parameters, it suffices to estimate some relevant quantities at $`\tau =0`$ that we give in Table 1. It is clear from those figures that the larger $`\mathrm{}`$, the smaller is the effect (at small times).

To summarize, we have analyzed the semiclassical dynamics of a shell made of bosons in AdS, thus generalizing the results previously obtained in asymptotically flat space-time . As in the latter case, we find that the shell spontaneously emits gravitational energy in the form of radiation when its velocity of contraction is not negligible and the adiabatic approximation fails. In AdS this occurs not only for the shell radius $`R`$ approaching the gravitational (Schwarzschild) radius $`R_H2M_s`$ but also for relatively large radii, $`R\mathrm{}`$, thus increasing the importance of such an effect. Another peculiarity of AdS is that the spectrum of massless radiation is not continuous, and one therefore finds that the emission occurs as a resonance between the shell inner energy gaps and the AdS levels. Our approach therefore suggests that the quantum nature of the collapsing matter implies the existence of radiation and possibly strong backreaction. Hence, the choice of vacua in Ref.  for a shell radius much larger than the horizon radius (Boulware vacuum) and for a shell approacing the horizon (Unruh vacuum) may actually be connected by the evolution of the system when the backreaction is properly considered. Since the radiation is present from the very initial stages of the collapse, this situation is analogous to the “topped-up” Boulware state used in Ref.  with a radiation of dynamical origin.


text:
# A principal positioning system for the Earth 11footnote 1Talk given at the Journées Systèmes de Référence Spatio-Temporels JSR-2002, 25-28 September, 2002, Bucharest (Romania)

## 1 Introduction

The current conception of the global navigation satellite systems (GNSS), like GPS, is based on a Newtonian model corrected numerically of some ”relativistic effects”. The direct relativistic theory suggests not only amelioration in accuracy, but also new functions for such GNSS systems.

The project SYPOR (Système de Positionnement Relativiste) proposes the ideas and instruments needed to carry out these new possibilities. Particularly it aims to endow the constellation of satellites of GALILEO of the necessary elements to constitute, by itself, a primary, autonomous positioning system for the Earth and its neighbors. The word ”autonomous” refers here to the capability of the constellation to provide complete relativistic metric information, i.e. to describe the kinematics and the dynamics both, of the constellation itself and of the users (possibility of gravimetry).

For this goal, the project SYPOR envisages, for the first time in physics and astronomy, to construct in the neighbors of the Earth a relativity-compatible physical coordinate system (relativistic positioning system).

In this paper we sketch the general lines of the project (Section 2), as well as the underlying physical concepts, namely that of relativistic positioning systems (Section 3) and the way to make such systems autonomous, which is illustrated in the two-dimensional case (Section 4).

Some of the ideas and results on this subject are the fruit of a long and friendly collaboration with Lluís Bel (Univ. Pais Vasco, Bilbao, Spain), Joan Josep Ferrando and Juan Antonio Morales (Univ. Valencia, Burjassot, Spain), and Albert Tarantola (Inst. Physique du Globe, Paris, France).

## 2 Sketch of the Project

For most of the needs of geodesy and positioning, the Earth may be considered as a Newtonian system, for which classical mechanics is enough to explain its essential properties. But a constellation of satellites around the Earth, endowed with clocks exchanging their proper times, is a relativistic system in its own right (mainly due to Doppler and gravitational potential ”relativistic effects”). Consequently, the natural conceptual frame to study GNSS is relativity theory.

At present, the GNSS involve the Earth and the constellation of satellites as a sole, coupled system. They start from a terrestrial, non relativistic coordinate system (Coll, 2001) and use the satellites of the constellation as moving beacons to indicate to the users their position with respect to this system.

The project SYPOR offers this result in two steps, that have different levels of conceptual precision and practical accuracy:

1. At the first level, SYPOR proposes the concepts and means to use the sole constellation of satellite-borne clocks as the most accurate, primary, autonomous, relativistically valid, positioning system for the neighbors of the Earth. At this level, any user may know its coordinates with respect to the satellites, its dynamical state (acceleration, rotation), the exact internal configuration of the constellation, and their situation with respect to the ICRS (i.e. all what an user may hope to obtain from a primary system), and any two users may know their relative position, distance and relative orientation.

Such a positioning system is a non usual one, with light-like (rapidly variable) coordinate surfaces, but any conventional system may be defined with respect to it. It is to be noted that, at this level, no synchronization of clocks is at all necessary.

2. At the second level, the usual data of the control segment on the trajectories of the satellites, are ”read” as the data defining the coordinate change to (secondary, non relativistic) terrestrial coordinates (WGS 84 or ITRF classical eference systems).

At this level, any user may know its position (with terrestrial precision) with respect to the Earth, as with the current GPS conception.

The possibility appears for a space agency to concentrate its interest in the first level, the autonomous positioning system, and delegate to global and local Earth agencies the control of the terrestrial coordinates.

To realize these performances, every satellite must be endowed with the following kinematic devices:

\- a device, on every satellite, allowing to exchange proper times with its neighbors (internal control of the parts of the system) (Hammesfahr, 1999),

\- a device, on four at least of the satellites, pointing to the ICRS (International Celestial Reference System) in order to define virtual local charts ”at rest” with respect to the ICRS (external control of the system as a whole),

\- a device, on every satellite, broadcasting over the Earth, beside its proper time, those of their neighbors (strong integrity: control by the users segment).

## 3 Positioning Systems

Relativity theory may be used:

* as a wise algorithm to sprinkle Newtonian expressions with terms corresponding to the ”relativistic effects” necessary for the obtention of the correct numerical values, or

* as the adequate starting frame to rise and to approach the physical situations with the most recent concepts and developments on the space-time.

The first use is undoubtedly correct in some particular scenarios, as may be in approximate numerical computations or for the abstract comparison of the equations of the two theories. The analysis shows that in fact this first use may be correct in the situations in which the physical determination of the coordinates either does not matter, or may be numerically identified to their Newtonian geometric determinations. But it is obviously useless to take advantage of the progress and specific developments of relativity theory in its proper domain, as it is the case in advanced GNSS.

The basic arena of relativity theory is its space-time. Relativistic space-time differs from Newtonian space-time in the following essential point: ”the space” and ”the present” are not now ”physical objects”, but inessential local arbitrary conventions. The three-dimensional Newtonian space has as much physical reality as have the Ptolemaic cristal spheres, and ”past”, ”present” and ”futur” are not exhaustive complementary parts of the space-time. Consequently, ”objects” in their usual sense, like galaxies, stars, planets, mountains do not exist in the relativistic space-time. What one can find in it are rather the ”absolute invariants” that they generate, that is to say, their ”histories”. It is with the histories of the satellites that a relativistic theory of GNSS must be constructed.

A locationing system is a detailed description for the physical construction of a coordinate system. To construct a coordinate system is either to construct its coordinate lines, or to construct its coordinate (hyper-)surfaces. But the choice of one or the other, and the protocol of their physical construction, give rise to very different physical properties (Coll, 2001).

Thus, when the protocol of their physical construction allows a particular observer (generally situated at the origin) to attach to every point of his neighbor a set of coordinates, one has what is currently called a reference system. But if this protocol allows every point of a neighbor to know its proper coordinates, then one has a positionning system.

In Newtonian physics, these functions are exchangeable, and perhaps this is why they are frequently mistaken. But in relativistic physics these two functions have very different physical properties:

* they are always incompatible for a sole coordinate system. So that the ”reference” or ”positioning” character of a locationing system must be previously chosen,

* it is always impossible to construct a positioning system starting from a reference system,

* it is always possible (and very easily) to construct a reference system starting from a positioning one.

As a consequence, the first element to be conceived in a relativistic GNSS must imperatively be its positioning system, and not its terrestrial reference system (WG 84 or IERS) as it is currently the case.

Remember (Coll, 2001) that relativistic positioning systems are generic, free, and immediate, three very important physical properties that no other locationing system may simultaneously offer. Thus, for example, Cartesian coordinates are not generic, the standard synchronization (two-way signals) usual in the construction of reference systems does not give rise to immediate systems; harmonic conditions are generic, but not free, etc.

A general analysis of rigorous mathematical results, physical possibilities and present technical developments leads to the important epistemic result that, among all the relativistic locationing systems, the set of relativistic positioning systems exists but constitute a very little class. And the simplest representative of this class is the one formed by electromagnetic signals broadcasting the proper time $`\tau _i`$ of four independent clocks $`S_i`$ ($`i=1,\mathrm{},4`$).

From now on, we suppose these clocks carried by (not necessarily geodetic) satellites.

In the space-time, the above wave fronts signals, parameterized by the proper time of the clocks, draw four families of physical hypersurfaces moving at the velocity of light, realizing a covariantly null coordinate system $`\{\tau _i\}`$.

Coordinate systems of this class are very unusual. Very different from the current relativistic ones, and still more from the Newtonian ones, they have been studied by a very restricted number of specialists in relativity. For an almost exhaustive bibliography, see (Blagojević et al., 2001).

In such covariantly null coordinate systems, instead of $`\eta _{\alpha \beta }=diag\{1,1,1,1\},`$ the Minkowski metric of the space-time adopts the complementary form:

$$\eta _{\alpha \beta }=\left(\begin{array}{cccc}0& f& g& h\\ f& 0& \mathrm{}& m\\ g& \mathrm{}& 0& n\\ h& m& n& 0\end{array}\right)$$
(1)
where $`f,g,h,\mathrm{},m,n,`$ are strictly positive functions of the coordinates $`\tau _i.`$

We see that in such coordinates there is no time-like asymmetry, the four coordinate surfaces playing exactly the same role. The nullity of all the diagonal terms in the above real expression seems to have erroneously suggested in the past that such coordinate systems would be ”somewhere degenerate”; this uncorrect intuitive feeling is perhaps the cause of the absence of studies on them and of their slow re-discovery by some authors.

The four proper times $`\{\tau _i\}`$ read at a space-time event by a receptor constitute its (covariantly null) coordinates with respect to the four satellites. But such a system can not be considered as primary (with respect to the space-time structure) if we have not sufficiently information to relate to it any other coordinate system (Cartesian, harmonic, etc). And for this task, we need to know the (dynamical) space-time trajectories of the satellites. In principle, there are many ways to do that, one simple one being to force satellites to follow prescribed trajectories, for example geodesics. But the most complete one is that in which these information is generated and broadcasted at every instant by the system of satellites itself, whatever be their trajectories. When this happens, we call the primary system autonomous.

## 4 Autonomous Positioning Systems

How to make autonomous such a system of embarked clocks , arbitrarily synchronized, broadcasting their proper times? The answer is very simple: broadcasting, not only their proper time, but also the proper time of their neighboring satellite’s clocks.

In other words: let $`\tau _{ij},`$ $`ij,`$ be the proper time of the satellite’s clock $`j`$ received by the satellite $`i`$ at its proper time instant $`\tau _i.`$ Then, the broadcasting of the data $`\{\tau _i,\tau _{ij}\}`$ allow to make autonomous the covariantly null space-time positioning system $`\{\tau _i\}`$.

Observe that the sixteen data $`\{\tau _i,\tau _{ij}\}`$ received by an observer contains, of course, the coordinates $`\{\tau _i\}`$ ($`i=1,\mathrm{},4`$), of this observer but also the coordinates $`\{\tau _i,\tau _{ij}\}`$ of every satellite $`i`$ in the totally covariant null coordinate system that the four satellites are generating.

In a four-dimensional Cartesian grid of axis $`\tau _i`$ ($`i=1,\mathrm{},4`$), the data $`\{\tau _i,\tau _{ij}\}`$ received by an user represent at every instant the proper position of the user as well as the positions of the satellites. During an arbitrary interval, these data allow the user to draw in this grid its proper space-time trajectory as well as the trajectory of the four satellites. Of course, two different users will found, in their corresponding grids, that their personal trajectories are different, but in the covering domain of proper times, the trajectories of the four satellites will be necessarily the same.

Such grids are diffeomorphic to the corresponding domain of the space-time, but they are not isometric. The precise deformation that relates them is nothing but the metric tensor $`g_{ij}(\tau _{\mathrm{}})`$ of this domain, which gives the dynamical properties of the trajectories (their inertial and/or gravitational characteristics).

When the grid trajectories of the satellites are particularly simple, their space-time dynamical ones are easy to obtain. For example, in the two-dimensional Minkowski case, if the grid trajectories of the two satellites are non parallel straight lines, one can show that necessarily the satellites follow convergent geodesics, meanwhile if the grid trajectories are parallel straight lines, the two satellites are necessarily uniformly accelerated.

When the grid trajectories are arbitrary, an important general result may be prove: the knowledge of the dynamics of one of the satellites during a relative acausal interval, allows to know the dynamics of both of the satellites as well as that of the user at any later instant.

In other words: if $`\mathrm{\hspace{0.17em}1}`$ and $`\mathrm{\hspace{0.17em}2}`$ denote the two clocks, the relative acausal interval of the satellite $`\mathrm{\hspace{0.17em}1}`$ with respect to the satellite $`\mathrm{\hspace{0.17em}2}`$ at the ”instant” $`\tau _2`$ is the interval of the trajectory of the satellite $`\mathrm{\hspace{0.17em}1}`$ between its proper times $`\tau _1^{}=\tau _{21}`$ and $`\tau _1^{\prime \prime }`$ such that $`\tau _{12}=\tau _2.`$ Then our result states that the dynamical knowledge of the satellite $`\mathrm{\hspace{0.17em}1}`$ during its proper time interval $`(\tau _1^{},\tau _1^{\prime \prime })`$ suffices to know its dynamical trajectory for any $`\tau _1\tau _1^{\prime \prime }`$ and of $`\mathrm{\hspace{0.17em}2}`$ and the user for any $`\tau _2\tau _{12}.`$

The relative acausal interval of, say, the satellite $`\mathrm{\hspace{0.17em}1}`$ with respect to the satellite $`\mathrm{\hspace{0.17em}2},`$ at the ”instant” $`\tau _2`$ is the interval of the trajectory of the satellite $`\mathrm{\hspace{0.17em}1}`$ between its proper times $`\tau _1^{}=\tau _{21}`$ and $`\tau _1^{\prime \prime }`$ such that $`\tau _{12}=\tau _2.`$ Then, in other words, our result states that the dynamical knowledge of the satellite $`\mathrm{\hspace{0.17em}1}`$ during its proper time interval $`(\tau _1^{},\tau _1^{\prime \prime })`$ suffices to know its dynamical trajectory for any $`\tau _1\tau _1^{\prime \prime }`$ and of $`\mathrm{\hspace{0.17em}2}`$ and the user for any $`\tau _2\tau _{12}.`$

Of course, around the Earth more than four satellites will be convenient. Every four neighboring ones will constitute a local chart of the atlas of covariantly null coordinate systems enveloping the Earth. From this atlas, appropriate virtual global conventional coordinate system may be defined.

## 5 References

Blagojević, M., Garecki, J., Hehl, F., Obukhov, Yu., 2001. Real Null Coframes in General Relativity and GPS type coordinates, Phys. Rev. D65 (2002) 044018, e-print: gr-qc/0110078. Almost exhaustive bibliography.

Coll, B., 2001. Physical Relativistic Frames, in Journées 2001 Systèmes de Référence Spatio-Temporels, Brussels, Belgium, 24-26 September 2001, also http://coll.cc . Some properties of classical and relativistic coordinate systems are pointed out.

Hammesfahr, J. et al., 1999. Intersatellite Ranging and Autonomous Ephemeris Determination for future Navigation Systems, Deutsches Zentrum für Luft- und Raumfahrt (DLR), document number ISR-DLR-REP-002. Excellent report on the possibility of intersatellite links (ISL).


text:
# Collision Drag Effect on Two-Fluid Hydrodynamics of Superfluid 3He in Aerogel

## I Introduction

It is well known that there exists a variety of sound modes in liquid 3He. In the normal phase, in addition to the usual hydrodynamic sound mode (the first sound), the so-called zero sound propagates as a well-defined sound mode in the collisionless regime at low temperatures. In the superfluid phase, other two kinds of sound, the second and forth sound, can propagate which reflects the presence of the normal and the superfluid component of liquid. The second sound is the out-of-phase oscillation of the two fluids. The forth sound is the oscillation of only the superfluid component when the normal fluid is clamped owing to its viscosity in a narrow channel.

Recently, there has been increasing interest in the sounds in liquid 3He confined in highly porous aerogel.Nomura ; Gervais ; McKenna ; Golov ; Rainer ; Ichikawa ; Higashitani ; Higa\_lt23 ; Kotera ; Berkutov  The aerogel consists of a random structure of silica strands of approximately 3 nm in diameter. The 3He-aerogel system has been regarded as a model system for studying the impurity scattering effects on the Fermi liquid which undergoes a transition to an unconventional Cooper pairing state. Nomura et al.Nomura  and Gervais et al.Gervais  have observed that the longitudinal sound attenuation in normal liquid 3He in aerogel does not show the characteristic temperature dependence of the Fermi liquid, i.e., the first-to-zero sound transition. In addition, the normal phase attenuation is of the same order of magnitude as the pure 3He case in spite of the presence of random scatterers. They also found the absence of the attenuation peak due to the order parameter collective excitation in the superfluid phase of 3He in aerogel. To study the attenuation theoretically, the homogeneous scattering modelThuneberg  can be a good starting point, because the sound wavelength ($``$ 20 $`\mu `$m) is much longer than the inter-strands spacing ($``$ 30 nm),Nomura ; Gervais  so that the detailed structure of aerogel is expected not to play a significant role. The homogeneous scattering model, however, predicts a significantly large attenuationRainer ; Ichikawa  (which is two order of magnitude larger than the observed value) though it can reproduce the qualitative temperature dependence.

To resolve the discrepancy between the theory and the experiment, Ichikawa et al.Ichikawa  have studied the effect of simultaneous aerogel motion caused by collisions with 3He-quasiparticles. Such a collision drag effect reduces drastically the longitudinal sound attenuation and the experimental results can be reproduced by the collision drag model.Ichikawa ; Higa\_lt23  This model has recently been applied to the transverse sound in the 3He-aerogel system.Higashitani  It was shown that the collision drag effect modifies strongly the nature of the transverse sound propagation; for example, the transverse sound in the hydrodynamic regime, which is an overdamped mode in pure 3He, can propagate over a long distance.

In this paper, we discuss the collision drag effect in superfluid 3He-aerogel system. McKenna et al.McKenna  have already proposed a theory of the sound in superfluid-aerogel system. They modified a set of two-fluid hydrodynamic equations to take into account the interlocking motion of the normal component and the aerogel. This theory was used by Golov et al.Golov  to extract the superfluid density $`\rho _s`$ in the 3He-aerogel system from the sound velocity measurements.

In deriving the modified two-fluid hydrodynamic equations, McKenna et al.McKenna  assumed the normal fluid to be completely locked to the aerogel matrix, i.e., the system is assumed to be in local equilibrium. The resulting equations can determine only the sound velocity and cannot give any information about the attenuation. In this paper, we discuss the longitudinal sound attenuation in superfluid 3He in aerogel when the system deviates slightly from local equilibrium. In addition to the bulk 3He-aerogel system, we shall consider a restricted geometry such that the aerogel is built in a pore formed by sintered silver powders, as used in the fourth sound experimentKotera .

## II Modification of Two-fluid Hydrodynamic Equations

Now we consider the drag force exerted on aerogel molecules. The drag force can be defined as the net momentum transfer per unit time from superfluid 3He to the aerogel during impurity scattering processes. Thus the drag force can be determined from the collision integral in the quasiparticle transport equation.Ichikawa ; Higashitani ; Higa\_lt23  In the normal state, the drag force density $`𝐅`$ is given by $`𝐅=\frac{1}{\tau _{\mathrm{tr}}}(1+F_1^s/3)\rho (𝐯𝐯_a)`$, where $`\tau _{\mathrm{tr}}`$ is the transport relaxation time for the impurity scattering, $`\rho `$ the mass density of liquid 3He, $`F_1^s`$ the Landau parameter and $`𝐯`$ and $`𝐯_a`$ the local velocities of the liquid and the aerogel, respectively. Since $`𝐅`$ is proportional to the relative velocity $`𝐯𝐯_a`$, it can be also interpreted as friction between the liquid and the aerogel. In the superfluid phase, the motion of the superfluid component is expected to uncouple with the aerogel. It is, therefore, plausible to generalize the normal-state result of $`𝐅`$ into

$$𝐅=\frac{1}{\tau _\mathrm{f}}\rho _n(𝐯_n𝐯_a),$$
(1)
where $`\rho _n`$ is the normal-fluid density, $`𝐯_n`$ is the normal-fluid velocity and $`\tau _\mathrm{f}`$ denotes a relaxation time characterizing the frictional effect between the normal fluid and the aerogel. The relaxation time $`\tau _\mathrm{f}`$ is expected to be of the same order as the impurity scattering time $`\tau _\mathrm{i}`$. In fact in the normal phase it is given by $`\tau _{\mathrm{tr}}/(1+F_1^s/3)`$. In superfluid phase, involved microscopic calculations are necessary to obtain an explicit expression for $`\tau _\mathrm{f}`$. We shall, therefore, treat $`\tau _\mathrm{f}`$ as a phenomenological parameter in this paper.

We can now obtain two-fluid hydrodynamic equations taking into account the collision drag effect. In the set of the conventional two-fluid equations, the momentum conservation law is modified because of the momentum loss due to the drag force, while the equation of motion of the superfluid velocity $`𝐯_s`$ remains unchanged. The two-fluid equations for longitudinal sound are thus given by

$`{\displaystyle \frac{\rho }{t}}+\mathrm{div}𝐉=0,`$ (2)
$`{\displaystyle \frac{\mathrm{div}𝐉}{t}}=^2P+{\displaystyle \frac{4}{3}}\eta ^2\mathrm{div}𝐯_n\mathrm{div}𝐅,`$ (3)
$`{\displaystyle \frac{𝐯_s}{t}}={\displaystyle \frac{1}{\rho }}P+ST,`$ (4)
$`{\displaystyle \frac{\rho S}{t}}+\mathrm{div}(\rho S𝐯_n)=0,`$ (5)
where $`\rho =\rho _n+\rho _s`$, $`𝐉=\rho _n𝐯_n+\rho _s𝐯_s`$, $`\eta `$ is the viscosity and $`P`$, $`S`$ and $`T`$ are pressure, entropy and temperature, respectively. The above equations are different from those by McKenna et al.McKenna  in that the drag force $`𝐅`$ and the viscosity $`\eta `$ are included. The former enables us to study the case when $`𝐯_n𝐯_a`$. The latter is necessary to reproduce the normal-state dispersion relation derived previously.Ichikawa ; Higa\_lt23  The two-fluid equations must be supplemented by the equations of motion of aerogel:

$`\rho _a{\displaystyle \frac{𝐯_a}{t}}=P_a+𝐅,`$ (6)
$`{\displaystyle \frac{\rho _a}{t}}+\mathrm{div}(\rho _a𝐯_a)=0,`$ (7)
where $`\rho _a`$ is the aerogel mass density and the restoring force due to the elasticity of the aerogel is denoted as $`P_a`$. The pressure $`P_a`$ is related to the longitudinal sound velocity $`C_{al}`$ of the skeleton aerogel byMcKenna

$$C_{al}^2=P_a/\rho _a.$$
(8)
## III Sound Dispersion Relation

From the above set of equations, we can derive the dispersion relation of the longitudinal sound in superfluid 3He-aerogel system. To do that, we shall assume, as usual, that all the non-equilibrium quantities have the space and time dependence of $`\mathrm{exp}(i𝐪𝐫i\omega t)`$. Then, using eq. (1) and the equations of motion of aerogel, we obtain the following relation between $`𝐯_n`$ and $`𝐯_a`$:

$$(1\frac{\omega _q^2}{\omega ^2})\rho _a𝐯_a=\frac{i}{\omega \tau _\mathrm{f}}\rho _n(𝐯_n𝐯_a)$$
(9)
where $`\omega _q=C_{al}q`$ is the frequency of skeleton aerogel sound. It follows that the local equilibrium ($`𝐯_n=𝐯_a`$) is achieved by taking the limit $`\omega \tau _\mathrm{f}0`$. Combining eq. (9) with our two-fluid equations, we find that the dispersion relation can be written as

$$(z^2C_1^2)(z^2C_2^2)+(i\frac{4\eta \omega }{3\rho _n}+\frac{i}{\omega \tau _\mathrm{f}}\frac{z^2}{1+\frac{i}{\omega \tau _\mathrm{f}}\frac{\rho _n}{\rho _a}\chi _a})(z^2C_4^2)=0$$
(10)
with $`z=\omega /q`$ and

$$\chi _a=\frac{\omega ^2}{\omega ^2\omega _q^2}=\frac{z^2}{z^2C_{al}^2}.$$
(11)
In eq. (10), $`C_1`$, $`C_2`$ and $`C_4=[(\rho _s/\rho )C_1^2+(\rho _n/\rho )C_2^2]^{1/2}`$ are the first, second and fourth sound velocities as defined in the conventional manner. It is easy to show that eq. (10) is reduced to the previous resultIchikawa ; Higa\_lt23  in the normal state. Moreover, eq. (10) with $`\eta =0`$ and in the limit $`\omega \tau _\mathrm{f}0`$ coincides, as expected, with the dispersion relation of McKenna et al.McKenna

$$(z^2C_1^2)(z^2C_2^2)+\frac{\rho _a}{\rho _n}(z^2C_{al}^2)(z^2C_4^2)=0.$$
(12)
Equation (12) determines the velocities of two sound modes (fast and slow modes) in local equilibrium.McKenna ; Golov  The two sound modes in the 3He-aerogel system have been observed by Golov et al.Golov  According to this experiment, the velocities of the fast and slow modes, $`C_f`$ and $`C_s`$, satisfy the conditions such that $`C_f^2C_2^2`$, $`C_f^2C_{al}^2`$ and $`C_1^2C_s^2C_2^2`$. Taking into these conditions, we find from eq. (12) that $`C_f`$ and $`C_s`$ are approximately given by

$`C_f^2={\displaystyle \frac{C_1^2+(\rho _a/\rho _n)C_4^2}{1+\rho _a/\rho _n}},`$ (13)
$`C_s^2={\displaystyle \frac{C_{al}^2}{1+\rho \rho _n/\rho _a\rho _s}}.`$ (14)
When $`\tau _\mathrm{f}`$ and $`\eta `$ are finite, the complex velocities $`z`$’s have corrections which give the attenuation $`\alpha =\mathrm{Im}q`$ of the two hydrodynamic sound modes. Using eq. (10) and keeping only the leading order correction terms, we obtain the attenuations of the fast and slow modes, $`\alpha _f`$ and $`\alpha _s`$, as

$`\alpha _f={\displaystyle \frac{\omega }{2C_f^3}}{\displaystyle \frac{C_1^2C_4^2}{(1+\rho _a/\rho _n)^2}}({\displaystyle \frac{4\eta \omega }{3\rho _nC_f^2}}+{\displaystyle \frac{\rho _a^2}{\rho _n^2}}\omega \tau _\mathrm{f})`$ (15)
$`\alpha _s={\displaystyle \frac{\omega }{2C_s}}{\displaystyle \frac{\rho _n/\rho _a}{1+\rho \rho _n/\rho _a\rho _s}}({\displaystyle \frac{4\eta \omega }{3\rho _nC_s^2}}+{\displaystyle \frac{\rho ^2}{\rho _s^2}}\omega \tau _\mathrm{f})`$ (16)
The expressions for the attenuations include terms $`\eta `$ and $`\tau _\mathrm{f}`$, reflecting the viscous effect and the frictional effect, respectively. The information on the phenomenological parameters $`\eta `$ and $`\tau _\mathrm{f}`$ can thus be obtained from attenuation measurements. We shall show in the next section that a more useful and direct information on the frictional relaxation time $`\tau _\mathrm{f}`$ can be provided by the forth sound attenuation measurement.

## IV Fourth Sound in Narrow Pore

Now we discuss the sound propagation in a superfluid 3He-aerogel system built in a narrow pore. Such a system has recently been prepared by Kotera et al.Kotera  using sintered silver powders of $``$ 70 $`\mu `$m in diameter to study the fourth sound propagation in aerogel. The sintered silver powders provide pores with the average diameter of $`L10`$ $`\mu `$m.

In such a system, the sintered silver will play a role to clamp the aerogel strand ends. The longitudinal wave of the aerogel along the pore channel is, therefore, accompanied by a transverse standing wave so that the fixed end boundary condition should be satisfied. It follows that the eigenfrequency $`\omega _q`$ of the aerogel is replaced by

$$\omega _q^2=C_{al}^2q^2+\left(\frac{\pi C_{at}}{L}\right)^2,$$
(17)
where $`C_{at}`$ is the transverse sound velocity of aerogel. The second term gives an energy gap due to the standing wave motion. Using $`C_{at}=178.5`$ m/s for 91 % porous aerogelBerkutov , the gap at $`q=0`$ is estimated to be 56 MHz for $`L=10`$ $`\mu `$m. We see that a condition $`\omega _q^2\omega ^2`$ is fully satisfied in low-frequency acoustic experiments in which the frequencies are typically about 10 kHz or less.

This observation shows that the aerogel motion in the pore is hardly excited by the low frequency sounds. In the limit $`\omega ^2/\omega _q^20`$, the aerogel cannot move and then the dispersion relation (10) takes the form

$$(z^2C_1^2)(z^2C_2^2)+(i\frac{4\eta \omega }{3\rho _n}+\frac{i}{\omega \tau _\mathrm{f}}z^2)(z^2C_4^2)=0.$$
(18)
The fourth sound can propagate in the limit $`\omega \tau _\mathrm{f}0`$ because the normal fluid is clamped by the “hard” aerogel. The attenuation $`\alpha _4`$ of the fourth sound is obtained from eq. (18). The result is given, up to leading order in $`\omega \tau _\mathrm{f}`$, by

$$\alpha _4=\frac{\omega }{2C_4^3}(C_1^2C_4^2)\omega \tau _\mathrm{f}.$$
(19)
Note that eq. (19) does not involve the viscosity $`\eta `$ because its contribution to $`\alpha _4`$ occurs only in higher order terms in the $`\omega \tau _\mathrm{f}`$ expansion. The attenuation of the fourth sound in the 3He-aerogel system is thus dominated by the frictional relaxation time $`\tau _\mathrm{f}`$. Using eq. (19), the parameter $`\tau _\mathrm{f}`$ can be determined from simultaneous measurements of $`C_4`$ and $`\alpha _4`$.

The reason why the frictional effect dominates the fourth sound damping is the following. In the pore channel with the aerogel, the normal flow no longer obeys the Hagen-Poiseuille law but is governed by the Drude law.Einzel  In the presence of aerogel, the viscous penetration depth $`\delta `$ for $`\omega \tau _\mathrm{f}1`$ is given byEinzel

$$\delta =\sqrt{\frac{\eta \tau _\mathrm{f}}{\rho _n}}$$
(20)
and the viscosity in the superfluid phase is

$$\eta \rho _n\overline{v}_\mathrm{F}^2\tau _\eta ,$$
(21)
where $`\overline{v}_F`$ is an appropriate average of the quasi-particle group velocity. At low temperatures, the viscous relaxation time $`\tau _\eta `$ is also limited by the impurity scattering time $`\tau _\mathrm{i}\tau _\mathrm{f}`$. In fact, in case of the sound experiment by Nomura et al.Nomura  (98% aerogel, 16 bar) one can see that $`\tau _\eta \tau _\mathrm{i}`$ happens at 10 mK. We find therefore that the viscous penetration depth in the presence of the aerogel is given by

$$\delta \overline{v}_\mathrm{F}\tau _\mathrm{f}$$
(22)
and is as small as the mean free path.

## V Concluding Remarks

We have discussed hydrodynamic sound propagation in superfluid 3He-aerogel system. In a bulk system, the aerogel oscillates together with the normal liquid by the collision drag effect. As a result, two kinds of sound (the fast mode and the slow mode)McKenna  can propagate as was observed by Golov et al.Golov . The deviation in the local velocity between the normal component and the aerogel leads to the damping of the sounds. We have given explicit expressions for the sound absorption. On the contrary in a narrow pore, the aerogel cannot move because the aerogel is clamped by the pore surface. In that case, the fourth sound can propagate along the pore channel. It is to be emphasized that the inter-locking between the normal component and the aerogel occurs not by the viscous effect as so far considered but by the collision drag effect. In other words, the normal component is in the Drude regime.Einzel

The modified two-fluid theory presented in this paper is still a phenomenological theory and includes many unknown parameters. In addition to the superfluid mass density, microscopic information on the relaxation times in dirty superfluid 3He is necessary for more detailed comparison with experiments. This problem and also the justification of our two-fluid model shall be discussed elsewhere.

###### Acknowledgements.

This work is supported by a Grant-in-Aid for COE Research (No. 13CE2002) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.