text:
# Zipf’s law in human heartbeat dynamics

## Abstract

It is shown that the distribution of low variability periods in the activity of human heart rate typically follows a multi-scaling Zipf’s law. The presence or failure of a power law, as well as the values of the scaling exponents, are personal characteristics depending on the daily habits of the subjects. Meanwhile, the distribution function of the low-variability periods as a whole discriminates efficiently between various heart pathologies. This new technique is also applicable to other non-linear time-series and reflects these aspects of the underlying intermittent dynamics, which are not covered by other methods of linear- and nonlinear analysis.

The nonlinear and scale-invariant aspects of the heart rate variability (HRV) have been studied intensively during the last decades. This continuous interest to the HRV can be attributed to the controversial state of affairs: on the one hand, the nonlinear and scale-invariant analysis of HRV has resulted in many methods of very high prognostic performance (at least on test groups) Amaral ; Ivanov ; Thurner ; Saermark ; on the other hand, practical medicine is still confident to the traditional “linear” methods. The situation is quite different from what has been observed three decades ago, when the “linear” measures of HRV became widely used as important noninvasive diagnostic and prognostic tools, soon after the pioneering paper Hon . Apparently, there is a need for further evidences for the superiority of new methods and for the resolution of the existing ambiguities.

During recent years the main attention of studies has been focused on the analysis of the scale-invariant methods. It has been argued that measures related to a certain time-scale (e.g. 5 min) are less reliable, because the characteristic time-scales of physiological processes are patient-specific. The scale-invariant measures are often believed to be more universal and sensitive to life-threatening pathologies Amaral ; Ivanov . However, carefully designed time-scale-related measures can be also highly successful, because certain physiological processes are related to a specific time scale Thurner .

The scale invariance has been exclusively seen in the heart rhythm following the (multi)fractional Brownian motion (fBm) Peng . It has been understood that the heart rhythm fluctuates in a very complex manner and reflects the activities of the subject (sleeping, watching TV, walking etc.) Poon ; Struzik  and cannot be adequately described by a single Hurst exponent of a simple fBm. In order to reflect the complex behavior of the heart rhythm, the multi-affine generalization of the fBm has been invoked Amaral ; Ivanov ; it has been claimed that the multifractal scaling exponents are of a significant prognostic value.

The approach based on fBm addresses long-time dynamics of the heart rhythm while completely neglecting the short-scale dynamics on time scales less than one minute (the respective frequencies are typically filtered out Peng ). The short-time variability has been described only by the so called linear measures, such as $`p_{NN50}`$ (the probability that two adjacent normal heart beat intervals differ more than 50 milliseconds). Meanwhile, the level of the short-time variability of the human heart rate varies in a very complex manner, the high- and low-variability periods are deeply intertwined Poon . This is a very important aspect, because the low-variability periods are the periods when the heart is in a stressed state, with high level of signals arriving from the autonomous nervous system. The conventional linear measures are not appropriate for describing such a complex behavior. Thus, there is a clear need for suitable nonlinear methods.

In this paper we present a new scale-invariant description of the short-time variability of the heart rate. It is shown that the distribution of low-variability periods in the activity of a normal heart follows the Zipf’s law. It is also shown that the distribution function of the low-variability periods contains a considerable amount of diagnostically valuable information. This new technique is also applicable to other non-linear time-series, such as EEG signals and financial data KSK .

Our analysis is based on ambulatory Holter-monitoring data (recorded at Tallinn Diagnostic Centre) of 218 patients with various diagnoses. The groups of patients are shown in Table 1. The sampling rate of ECG was 180 Hz. The patients were monitored during 24 hour under normal daily activities. The preliminary analysis of the ECG recordings was performed using the commercial software; this resulted in the sequence of the normal-to-normal (NN) intervals $`t_{NN}`$ (measured in milliseconds), which are defined as the intervals between two subsequent normal heartbeats (i.e. normal QRS complexes).

Originally, the Zipf’s law addressed the distribution of words in a language Zipf : every word has assigned a rank, according to its “size” $`f`$, defined as the relative number of occurrences in some long text (the most frequent word obtains rank $`r=1`$, the second frequent — $`r=2`$, etc). The empirical size-rank distribution law $`f(r)r^\alpha `$ is surprisingly universal: in addition to all the tested natural languages, it applies also to many other phenomena. The scaling exponent $`\alpha `$ is often close to one (e.g. for the distribution of words). Typically, the Zipf’s law is applicable to a dynamical system at statistical equilibrium, when the following conditions are satisfied: (a) the system consists of elements of different size; (b) the element size has upper and lower bounds; (c) there is no intermediate intrinsic size for the elements. As already mentioned, the human heart rhythm has a complex structure, where the duration $`\tau `$ of the low-variability periods varies in a wide range of scales, from few to several hundreds of heart beats. Thus, one can expect that the distribution of the low-variability periods follows the Zipf’s law

$$r\tau ^\gamma .$$
(1)
However, the scaling behavior should not be expected to be perfect. Indeed, the heart rate is a non-stationary signal affected by the non-reproducible daily activities of the subjects. The non-stationary pattern of these activities, together with their time-scales, is directly reflected in the above mentioned distribution law. This distribution law can also have a fingerprint of the characteristic time-scale (around ten to twenty seconds) of the blood pressure oscillations. Finally, there is a generic reason why the Zipf’s law fails (or is non-perfect) at small rank numbers. The Zipf’s law is a statistical law; meanwhile, each rank-length curve is based on a single measurement. Particularly, there is only one longest low-variability period (and likewise, only one most-frequent word), the length of which is just as long as it happens to be, there is no averaging whatsoever. For large rank values $`r`$, the relative statistical uncertainty can be estimated as $`1/\sqrt{r}`$.

To begin with, we define the local variability for each ($`i`$-th) interbeat interval as the deviation of the heart rate from the local average,

$$\delta (i)=\frac{|t_{NN}(i)t_{NN}(i)|}{t_{NN}(i)}.$$
(2)
The angular braces denote the local average, calculated using a narrow (5 beats wide) Gaussian weight function. Further, we introduce a threshold value $`\delta _0`$; $`i`$-th interbeat interval is said to have a low variability, if the condition

$$\delta (i)\delta _0$$
(3)
is satisfied. A low-variability period is defined as a set of consecutive low-variability intervals; its length $`\tau `$ is measured in the number of heartbeats. Finally, all the low-variability periods are arranged according to their lengths and associated with ranks. The rank of a period is plotted versus its length in a logarithmic graph, see Fig. 1; Zipf’s law would correspond to a straight descending line.

For a very low threshold parameter $`\delta _0`$, all the low-variability periods are very short, because it is difficult to satisfy the stringent condition (3). In that case, the inertial range of scales is too short for a meaningful scaling law. On the other hand, for a very high value of $`\delta _0`$, there is a single low-variability period occupying the entire HRV-recording. Between these two cases, there is such a range of the values of $`\delta _0`$, which leads to a non-trivial rank-length law. For a typical healthy patient, the $`r(\tau )`$-curve is reasonably close to a straight line, and the scaling exponent $`\gamma `$ is a function of the threshold parameter $`\delta _0`$. Thus, unlike all the other well-known applications of the Zipf’s law, we are dealing with a multi-scaling law.

Recently, Ivanov et. al. Ivanov  have reported that anomalous multifractal spectra of the HRV signal indicate an increased risk of sudden cardiac death. Therefore, it is natural to ask, does the presence or failure of the multiscaling behavior indicate the healthiness of the patient? In what follows we discuss a somewhat more general question: what is the relationship between the properties of the distribution function of the low variability periods and the diagnosis of the patient. Testing the prognostic significance for predicting sudden cardiac death, which is also of a great importance, has been postponed due to the nature of our test groups.

First, let us analyze the correlation between the diagnosis of a patient and the scaling exponent $`\gamma `$. To begin with, we have to determine the optimal value for the threshold parameter $`\delta _0`$. For a meaningful analysis, the scaling behavior should be as good as possible. It turned out that for a typical patient, the best approximation of the function $`r(\tau )`$ with a power law is achieved for $`\delta _00.05`$ (see Fig. 2a); in what follows, all the values of the exponent $`\gamma `$ are calculated for $`\delta _0=0.05`$. It should be noted that for some patients, the length-rank distribution is still far from a power law (see Fig. 2b).

The slope of a curve on the logarithmic plot is calculated using root-mean-square (rms) fit for such a range of lengths $`[\tau _{\text{start}},\tau _{\text{end}}]`$, for which the $`r(\tau )`$-curve is nearly a power law, and the scaling range width $`\mathrm{\Delta }=\mathrm{ln}\tau _{\text{end}}\mathrm{ln}\tau _{\text{start}}`$ is as large as possible. Bearing in mind the statistical nature of the Zipf’s law and non-stationarity of the underlying signal, we have chosen a not very stringent definition of what is “nearly a power law”, see Fig. 3. Around the rms-fit-line, two limit lines are drawn; $`\tau _{\text{start}}`$ and $`\tau _{\text{end}}`$ correspond to the points, where the $`r(\tau )`$-curve crosses the limit lines.

Note that the precise placement and shape of the limit lines is arbitrary, i.e. small variations do not lead to qualitative effects. Here, the distance of the limit lines from the central line has been chosen to be $`\mathrm{ln}2`$ at $`\tau =\tau _{\mathrm{max}}`$, and zero at $`\tau =1`$, where $`\tau _{\mathrm{max}}`$ is the length of the longest low-variability period. Admitting mismatch $`\mathrm{ln}2`$ at $`\tau =\tau _{\mathrm{max}}`$ is motivated by the observation that due to the lack of any statistics, the longest low-variability period could have been easily twice as long as we measured it to be. For large rank values, the statistical uncertainty is assessed as $`\sqrt{r}`$; in logarithmic graph, this would correspond to an exponentially decreasing (with increasing $`r`$) distance between a limit curve and the central fit-line. However, the above mentioned effect of the non-stationary pattern of the subjects daily activities makes the situation more complicated. There is no easy way to quantify this effect and therefore, we opted for the simplest possible solution, simple straight limit lines.

The scaling exponent $`\gamma `$ has been calculated for all the patients and Student test was applied to every pair of groups. In most cases, the significance was quite low; two best distinguishable groups were RR and FSK, the result of Student test being 5.7%. Therefore, one can argue that the slopes of linear parts are highly personal characteristics depending also on the daily habits of the subjects, which are weakly correlated with diagnosis.

Further we tested, how is the failure of the power law correlated with the diagnosis. The width of the scaling range $`\mathrm{\Delta }`$ was used as a measure of how well the curve is corresponding to a power law. The Student test results for the parameter $`\mathrm{\Delta }`$ turned out to be similar to what has been observed for the parameter $`\gamma `$: the correlation between the failure of the power law and diagnosis was weak. Thus, a rank-length curve resembling the one depicted by a dashed line in Fig. 2, does not hint to heart pathology. It should be also noted that the dashed curve in Fig. 2 can be considered as a generalized form of scale-invariance with scale-dependent differential scaling exponent. Such a behavior seems to be universal; for instance, certain forest fire models Chen  lead to the differential fractal dimension depending on the local scale.

Finally, we analyzed the diagnostic significance of the parameters $`\mathrm{ln}\tau _{\text{end}}`$ and $`\mathrm{ln}\tau _{\text{start}}`$. This analysis does make sense, because typically, the start- and end-points of the scaling range correspond to certain physiological time-scales. The parameter $`\mathrm{ln}\tau _{\text{end}}`$ provided, indeed, a remarkable resolution between the groups of patients, see Table 2. According to the Student test, the healthy patients, were distinct from five heart pathology groups with probability $`p<1.6`$%. The parameter $`\mathrm{ln}\tau _{\text{start}}`$ was diagnostically less significant.

Unfortunately, the calculation of the parameter $`\tau _{\text{end}}`$ is technically quite a complicated task, not suited for clinical practice. Therefore, we aimed to find a simpler alternative to it. Basically, the strategy was to find a simple parameter reflecting the behavior of the rightmost (large-$`\tau `$) part of the $`r(\tau )`$-curve. An easy option is $`\mathrm{ln}\tau _{\mathrm{max}}`$, which has been already analyzed JKMS . This parameter has indeed a considerable diagnostic value, but its reliability is decreased by the above discussed statistical fluctuations. Better alternatives are provided by (a) the overall number of low-variability periods $`r_{\mathrm{max}}`$ (which is small, if there are lot of long low-variability periods); (b) the coordinates of specific points of the rank-length curve. Here we chose a set of critical ranks $`R=10`$, $`20`$ or $`40`$, and determined the respective lengths $`\tau _R`$ so that $`r(\tau _R)=R`$. We also fixed a set of critical length values, $`T=50`$, $`100`$, or $`200`$, and determined the respective rank numbers $`r_T=r(T)`$, see Fig. 4. Both techniques turned out to be of high diagnostic performance; illustrative $`p`$-values are given in Table 2. Parameters $`\tau _{10}`$ and $`\tau _{20}`$ performed less well than $`\tau _{40}`$ (for instance, the $`p`$-values for the healthy and VES-subject groups were 0.60%, 0.58% and 0.34%, respectively), and are not presented in tabular data. Similarly, $`r_{100}`$ turned out to be more efficient than $`r_{50}`$ and $`r_{200}`$ (the respective healthy and VES-group $`p`$-values being $`0.02`$%, $`0.01`$%, and $`0.09`$%). It also outperforms $`\tau _{40}`$, but is sometimes less efficient than $`r_{\mathrm{max}}`$ or $`\tau _{\text{end}}`$ (see Table 2). Hence, various heart pathologies seem to affect the heart rate dynamics at the time scale around 100 heart beats (one to two minutes).

In conclusion, new aspect of non-linear time-series has been discovered, the scale-invariance of low-variability periods. We have shown that the distribution of low variability periods in the activity of human heart rate typically follows a multi-scaling Zipf’s law. The presence or failure of a power law, as well as the values of the scaling exponents, are personal characteristics depending on the daily habits of the subjects. Meanwhile, the distribution function of the low-variability periods as a whole contains also a significant amount of diagnostically valuable information, the most part of which is reflected by the parameters $`r_{100}`$, $`r_{\mathrm{max}}`$, and $`\tau _{\text{end}}`$, see Table 2. These quantities characterize the complex structure of HRV signal, where the low- and high variability periods are deeply intertwined, aspect which is not covered by the other methods of heart rate variability analysis (such as fractional Brownian motion based multifractal analysis). As a future development, it would be of great importance to analyze the prognostic value of the above mentioned parameters for patients with sudden cardiac death.

The support of ESF grant No. 4151 is acknowledged.


format: text

text:
# On the Global Anisotropy of Cosmic Ray Data above 4×10¹⁹ eV

## Introduction

A puzzle has existed for more than 30 years regarding cosmic ray events with energies exceeding $`4\times 10^{19}`$ eV, a value in the range of the Greisen, Zatsepin and Kuzmin (GZK) bound . The nature of the primary particles causing these events is controversial (see e.g. Refs.  for recent review of the field). If the primaries are protons from cosmological distances, then their sources should be isotropically distributed. The scrambling effects of intervening magnetic fields are difficult to assess, but on very basic physics the magnitudes of the fields (or associated correlation lengths) would have to be exceedingly small to avoid isotropization. Diffuse propagation of cosmic rays from a few prominent sources may result in isotropic arrival directions as suggested in Refs. . Cosmic rays from our own Galaxy core are not considered a viable explanation for the events approaching the GZK limit . If our Galaxy substantially modifies propagation or contributes to production of the primaries above $`4\times 10^{19}`$ eV, then any consequent anisotropy should be correlated with known source(s) of the Galaxy. On several bases, and in particular on the basis of symmetry, one asks whether the highest energy events are isotropic, up to the cuts imposed by observational limitations.

Several studies  at much lower energies have demonstrated a significant isotropy in arrival directions. The AGASA group has presented an analysis of a sizable (4$`\%`$) anisotropy of several hundred thousand events with energies above $`10^{17}`$ eV. The study uses right ascension, declination $`(RA,DEC)`$ coordinates  and finds a statistically significant first Fourier moment in $`\mathrm{sin}(RA)`$. The definition of right ascension however, involves historical human conventions and a choice of axis based on the Earth’s orbit. The coordinate system axis defining $`RA`$ requires two parameters for specification, which are implicitly used in the analysis, potentially affecting the claimed statistical significance of effects. In addition, $`\mathrm{sin}\varphi `$ is not a particular sensitive or informative statistic, potentially diminishing the impact of 216,000 events.

Here we examine all available data from the AGASA experiment with energy above $`4\times 10^{19}`$ eV. We take care to define certain quantities of the analysis which previously have been used in a coordinate dependent way. Since the coordinate system is a human artifact, this is to be avoided, unless there is a preferred coordinate system based on other information. Our methods are taken from Refs.  and a large literature on invariant angular statistics. Recently Sommers  has advocated an approach using the same basic concepts.

Statistical comparisons require formulation of a well defined null hypothesis. We take particular care on this point. The term “null” often indicates conditions of “no signal”, which would be appropriate if one considers anisotropy a signal. More deeply, the scientific method proceeds by ruling out possibilities, rather than attempting to prove particular notions. The symmetries of data are the simplest and purest characterizations capable of being tested, or ruled out. Thus a solid null hypothesis is the foundation of any further claims.

The elimination of coordinate dependent artifacts, the testing of a symmetry-based null, and the objectivity of maximum likelihood analysis create very powerful tools. If this method is valid, then it should be able to determine the systematic bias in real data objectively. The main purpose of the paper is to demonstrate the power of the tools for future use. However we also test our method to find significant evidence for anisotropy in the distribution of the highest energy events, which has been a question of intrinsic interest for many years.

## Methodology

The method of maximum likelihood (see Refs.  for example) allows one to test objectively between different statistical models. It is good practice to test models which are continuously related to the null when parameters are varied. Continuity in the likelihood test tends to assure that the same features of the data are tested by the null and competing models.

Let $`\stackrel{}{x}`$ denote an element of a data set, for example the coordinates of an observed track. Given a normalized distribution $`f(\stackrel{}{x})`$, the likelihood of the data in the distribution is defined to be the product of the distribution evaluated at each data point over the data set. Maximum likelihood occurs at maximum log likelihood $``$, defined as

$$=\underset{j}{\overset{N}{}}log[f(\stackrel{}{x}_j)],$$
(1)
where $`N`$ is the number of data points.

Hypothesis testing is done by comparing the difference of maximized log-likelihoods ($`T`$) for a given data set evaluated with two different distributions $`f(\stackrel{}{x}),f_{\mathrm{null}}(\stackrel{}{x})`$, where $`f_{\mathrm{null}}(\stackrel{}{x})`$ is the null distribution. For large $`N1`$, the statistic $`2T=2(_{\mathrm{model}}^{\mathrm{max}}_{\mathrm{null}}^{\mathrm{max}})`$ has a chi-squared distribution for $`p`$ parameters ($`\chi _p^2`$). We will verify this and also determine the distribution independently by Monte Carlo. Note that it is not possible and not our objective to find the “exact” distribution behind the observed data. Instead, emphasis lies on testing the null distribution. It is sufficient to show that a trial distribution fits the data sufficiently well to rule out the null; one does not imply, or conclude from this that the trial distribution is the last word.

## Covariant and Invariant Statistics

For covariant quantities we map the astronomical coordinates $`(DEC,RA)(\pi /2\theta ,\varphi )\widehat{x}`$, where $`\widehat{x}`$ is a unit vector on the dome of the sky,

$$\widehat{x}=(\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{cos}\theta ).$$

Naturally $`\widehat{x}`$ transforms like a vector when the coordinate system is changed. The expansion in spherical harmonics are the same thing: vectors are the $`l=1,m`$ representations of spherical harmonics. By contracting vector indices it is straightforward to make trial distributions and statistically valid quantities that are scalars under rotations, and therefore independent of the coordinate system.

The procedure is very natural, and so simple that we should expand on alternatives which do not have the same features. It is very common in the subject of “circular statistics” to see quantities such as $`<\theta >=_i^N\theta _i/N,\mathrm{\Delta }\theta =\sqrt{<\theta ^2><\theta >^2}`$. Such quantities can be computed but they are so faulty as to be nearly meaningless. The fault is immediately seen by calculating $`<\theta >`$ for an isotropic distribution on a circle. The distribution has no preferred orientation, yet the average angle $`<\theta >\pi `$ yields a preferred point, which is unacceptable.

## Analysis

We apply our method to the set of 58 published tracks from the AGASA group with energies not less than $`4\times 10^{19}`$ eV (Table 1, Refs. ). The average angular resolution of the track directions is stated to be 1.8 degrees  containing one sigma ($`68\%`$) of the events. The tracks come from the region of polar angle $`100^{}<\theta <10^{}`$. This data set could be, in principle, anisotropic in any direction.

The world’s published data consists of some 200 points. However not all are available, while resolution, cuts and data quality varies, and the dangers of combining data from different groups suggest that restricting the study to a homogeneous set is the logical first step. We also want to make it clear that objective and highly significant conclusions can be extracted from carefully constructed tests without needing particularly large data sets. Our methods, of course, extend readily to much larger data sets anticipated from new arrays including HIRES and AUGER.

The most basic model  to analyze unimodal spherical data {$`\widehat{x}(\theta ,\varphi )`$} is the Fisher distribution:

$$f_{\mathrm{Fisher}}(\theta ,\varphi )=\frac{\kappa }{4\pi \mathrm{sinh}\kappa }\sqrt{1(\widehat{n}\widehat{x})^2}e^{\kappa (\widehat{n}\widehat{x})}$$
(2)
where $`\widehat{n}(\alpha ,\beta )`$ is the direction of symmetry-breaking axis. This distribution has a long history as the spherical generalization of the venerated von-Mises distribution on the circle. Both distributions have the analytic properties of the Gaussian distribution, as seen by writing

$$e^{\kappa (\widehat{n}\widehat{x})^2/2}e^{\kappa (\widehat{n}\widehat{x})},\widehat{n}^2=\widehat{x}^2=1,$$

suppressing normalization factors. Parameter $`\kappa `$ is called the concentration parameter, and determines the extent that data is concentrated along the anisotropy axis $`\widehat{n}`$.

The Fisher distribution has the minimal number of parameters possible on the sphere, and exhibits cylindrical symmetry about the axis $`\widehat{n}`$. The pre-factor $`\sqrt{1(\widehat{n}\widehat{x})^2}`$ is explained as follows. Choose axis direction $`\widehat{n}=`$ (0, 0, 1) or along the north pole . This hides 2 parameters, and reduces the distribution to the simpler form:

$$f_{\mathrm{Fisher}}(\theta ,\varphi )=\frac{\kappa }{4\pi \mathrm{sinh}\kappa }\mathrm{sin}\theta e^{\kappa \mathrm{cos}\theta }.$$
(3)
The pre-factor $`\mathrm{sin}\theta `$ takes into account the Jacobian of solid angle on the sphere. Cylindrical symmetry about $`\widehat{n}`$ is also obvious now, since $`\varphi `$ does not appear anywhere in (3). For the limit of small anisotropy $`\kappa 1`$, the exponential can be expanded

$$e^{\kappa \mathrm{cos}\theta }1\kappa \mathrm{cos}\theta ,$$

which is one of the distributions discussed by Sommers , and which needs the additional Jacobian factor to test anisotropy.

The data set is subject to a cut in declination, ($`10^{}\theta 100^{}`$). We modify the Fisher distribution to take this cut into account as:

$$f_{\mathrm{Fisher}}^{\mathrm{cut}}(\theta ,\varphi )=\frac{f_{\mathrm{Fisher}}(\theta ,\varphi )g(\theta )}{f_{\mathrm{Fisher}}(\theta ,\varphi )g(\theta )𝑑\theta 𝑑\varphi };$$
(4)
$$g(\theta )=[\eta (\theta 10^{})\eta (\theta 100^{})],$$
(5)
where $`\eta `$ is the Heaviside step function. The step functions have an invariant representation we do not bother to write here. Now (4) is the normalized Fisher distribution in the cut region and serves as a trial model of anisotropy. A variant of the Fisher distribution known as Watson distribution is useful in cases where the data is in either bipolar or girdle form and is given by

$$f_{\mathrm{Watson}}(\theta ,\varphi )=\frac{\sqrt{1(\widehat{n}\widehat{x})^2}e^{\kappa (\widehat{n}\widehat{x})^2}}{4\pi _0^1e^{\kappa u^2}𝑑u}$$
(6)
which can be modified to take in to account the cuts in the data in the same way in (5).

We also need a null distribution on the cut region to test our model using the method of maximum likelihood. The null isotropic distribution on the sphere with the cut in angle $`\theta `$ is given by,

$$f_{\mathrm{null}}^{\mathrm{cut}}(\theta ,\varphi )=\frac{\mathrm{sin}\theta g(\theta )}{\mathrm{sin}\theta g(\theta )𝑑\theta 𝑑\varphi }.$$
(7)
The null distribution coincides with $`f_{\mathrm{Fisher}}^{\mathrm{cut}}(\theta ,\varphi )`$ in the limit $`\kappa 0`$. Thus when likelihood is evaluated, the variation of parameter $`\kappa `$ allows the analysis to be continuously connected with the null. However, the detectors usually do not have uniform acceptance in declination. AGASA  gives a distribution in declination of all events above $`>10^{19}`$ eV. We fit the distribution by a function of the form ($`a\mathrm{sin}b\theta +c`$) which is related to the isotropic distribution in $`\theta `$, namely $`d(\mathrm{cos}\theta )`$. While the isotropic distribution peaks at $`\theta =90^{}`$, the detector response drops as one move away from $`\theta =57^{}`$. The distributions are plotted in Fig. 1.

## Results

Evaluation of likelihood for the AGASA data set is straightforward. One finds a maximum log likelihood in the Fisher distribution to be

$$_{\mathrm{Fisher}}=128.10$$

units at

$$\kappa =1.35\pm 0.41;\widehat{n}=(4.22^{},88.34^{}).$$

In equatorial coordinates ($`RA`$, $`DEC`$), the axis of concentration lies along the direction ($`5^h\mathrm{\hspace{0.33em}53.36}^m,\mathrm{\hspace{0.33em}85.78}^{}`$). We have also evaluated the maximum log likelihood for Watson distribution (6) and the result is: $`_{\mathrm{Watson}}=129.71`$ units.

The null hypothesis is evaluated by running the same analysis with our fit function to the detector response in $`\theta `$ shown in Fig. 1. The log likelihood is found to be

$$_{\mathrm{null}}=128.81$$

units. The quantity $`2T^{\mathrm{max}}=2(_{\mathrm{model}}_{\mathrm{null}})`$ is 1.42 units in case of Fisher distribution ($`2T_{\mathrm{Fisher}}^{\mathrm{max}}`$) and 1.81 units in case of Watson distribution ($`2T_{\mathrm{Watson}}^{\mathrm{max}}`$). However, if the detector response is not explicitly known or is not well-defined, one uses the isotropic null distribution (7). The log likelihood is $`_{}^{}{}_{\mathrm{null}}{}^{}=134.39`$ units for isotropic null case. Correspondingly $`2T_{\mathrm{Fisher}}^{{}_{}{}^{}\mathrm{max}}=12.57`$ and $`2T_{\mathrm{Watson}}^{{}_{}{}^{}\mathrm{max}}=9.36`$.

For large $`N1`$ it is known that $`2T`$ is distributed like $`\chi _p^2(2T)`$ of $`p`$ parameters. In practice 30-40 points is sufficient for large $`N`$ to apply, but the question also depends on the dimensionality of the data. We have $`3`$ parameters in our case, so one should expect a $`\chi _3^2`$ distribution. This is a prediction of pure statistical theory, without detail about the problem at hand, and something that may be questioned.

Real data sets often contain correlations or irregularities of myriad possible origin. Such features may upset analytic estimates. We therefore made an independent Monte Carlo study of the statistic $`2T`$ based on the features of the problem. We took ten thousand random samples of $`58`$ data points inside the cut region. Randomness was implemented by selecting points from a distribution flat in $`\varphi `$, and flat in $`\mathrm{cos}\theta `$. For each sample of 58 we varied $`\kappa `$ and $`\widehat{n}`$ to determine the maximum likelihood, and we calculated the likelihood for $`\kappa =0`$; the result was a trial value of $`2T`$. The distribution of $`2T`$ values for 10,000 runs is shown in Fig 2 along with the $`\chi _3^2`$ distribution. Agreement is excellent.

Statistical significance is determined by calculating $`P`$-values, also known as confidence levels, which are defined as the integrated probability of $`2T`$ to fluctuate in the null distribution above the value determined for the data. A comparison of $`P`$-values is also a method to determine whether the distributions have any unexpected long tails. The $`P`$-value for the $`\chi _3^2`$ distribution to give $`2T1.42`$ or $`2T1.81`$ is very large. It means the AGASA data set is rather isotropic and do not strongly point towards a unique source direction or do not suggest that the data points are concentrated in a girdle form. This is the conclusion from using detector response in declination reported by AGASA group. On the other hand, from isotropic null distribution, the $`P`$-value for the $`\chi _3^2`$ distribution to give $`2T12.57`$ is $`3.9\times 10^4`$, or 0.039%. In terms of 2-sided Gaussian statistics, which are often used for comparisons, $`P=3.9\times 10^4`$ corresponds to a “3.55 $`\sigma `$” effect.

## Observations

We plotted the data in several coordinate systems to visually examine signals of anisotropy, and to check the fit to the Fisher distribution. Results in Galactic Coordinates and Aitoff-Hammer equal-area projection are shown in Fig 3. The advantage of an equal-area projection is that an isotropic distribution will appear uniform, and anisotropy is not unduly distorted by the coordinates. Galactic coordinates are used so that any correlation with the galaxy plane might be easily seen. The equatorial plane is shown by the solid curve.

The figure shows the Fisher “hot-spot” extracted from the fit, namely the best-fit anisotropy axis $`\widehat{n}`$. As can be seen, the Fisher axis is close to the equatorial north pole and its error cone actually encloses the north pole. Also shown are the locations of some other features of obvious interest: the hot-spot of the lower energy AGASA analysis , and the approximate location of the galaxy core. We added the direction of our system’s motion relative to the cosmological background radiation, as deduced from the dipole component of the the $`3^{}`$ cosmic microwave background, for another comparison.

## Conclusion

Unresolved questions regarding the origin of cosmic rays with energies above $`4\times 10^{19}`$ eV suggests testing the degree of isotropy of the data. The method we have developed is capable of ruling out isotropy and eliminating any coordinate dependent artifacts. Given the detector bias reported by AGASA, the group’s cosmic ray tracks above $`4\times 10^{19}`$ eV are consistent with an isotropic distribution. The question of anisotropy probably hinges on a finer understanding of angular bias than currently available in published sources: for example, error bars on the angular acceptance would be useful.

Acknowledgements: This study springs out of collaboration with Pankaj Jain, who we thank for productive comments and criticsm. Work was supported in part under the Department of Energy grant number DE-85ER40212 by the University of Kansas General Research Fund, and the Kansas Institute for Theoretical and Computational Science/K\*STAR program.

| Table 1: AGASA Events Above $`4\times 10^{19}`$ eV (J2000 coordinates) | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Date | Energy | $`RA`$ | $`DEC`$ | Date | Energy | $`RA`$ | $`DEC`$ |
|  | ($`\times 10^{19}`$ eV) |  |  |  | ($`\times 10^{19}`$ eV) |  |  |
| 84/12/12 | 6.81 | 22h 21m | 38.4 | 84/12/17 | 9.79 | 18h 29m | 35.3 |
| 86/01/05 | 5.47 | 4h 38m | 30.1 | 86/10/23 | 6.22 | 14h 02m | 49.9 |
| 87/11/26 | 4.82 | 21h 57m | 27.6 | 89/03/14 | 5.27 | 13h 48m | 34.7 |
| 89/08/16 | 4.07 | 5h 51m | 58.5 | 90/11/25 | 4.51 | 16h 17m | $``$7.2 |
| 91/04/03 | 5.09 | 15h 47m | 41.0 | 91/04/20 | 4.35 | 18h 59m | 47.8 |
| 91/05/31 | 5.53 | 3h 37m | 69.5 | 91/11/29 | 9.10 | 19h 06m | 77.2 |
| 91/12/10 | 4.24 | 0h 12m | 78.6 | 92/01/07 | 4.51 | 9h 36m | 38.6 |
| 92/01/24 | 4.88 | 17h 52m | 47.9 | 92/02/01 | 5.53 | 0h 34m | 17.7 |
| 92/03/30 | 4.47 | 17h 03m | 31.4 | 92/08/01 | 5.50 | 11h 29m | 57.1 |
| 92/09/13 | 9.25 | 6h 44m | 34.9 | 93/01/12 | 10.1 | 8h 17m | 16.8 |
| 93/01/21 | 4.46 | 13h 55m | 59.8 | 93/04/22 | 4.42 | 1h 56m | 29.0 |
| 93/06/12 | 6.49 | 1h 16m | 50.0 | 93/12/03 | 21.3 | 1h 15m | 21.1 |
| 94/07/06 | 13.4 | 18h 45m | 48.3 | 94/07/28 | 4.08 | 4h 56m | 18.0 |
| 95/01/26 | 7.76 | 11h 14m | 57.6 | 95/03/29 | 4.27 | 17h 37m | $``$1.6 |
| 95/04/04 | 5.79 | 12h 52m | 30.6 | 95/10/29 | 5.07 | 1h 14m | 20.0 |
| 95/11/15 | 4.89 | 4h 41m | 29.9 | 96/01/11 | 14.4 | 16h 06m | 23.0 |
| 96/01/19 | 4.80 | 3h 52m | 27.1 | 96/05/13 | 4.78 | 17h 56m | 74.1 |
| 96/10/06 | 5.68 | 13h 18m | 52.9 | 96/10/22 | 10.5 | 19h 54m | 18.7 |
| 96/11/12 | 7.46 | 21h 37m | 8.1 | 96/12/08 | 4.30 | 16h 31m | 34.6 |
| 96/12/24 | 4.97 | 14h 17m | 37.7 | 97/03/03 | 4.39 | 19h 37m | 71.1 |
| 97/03/30 | 15.0 | 19h 38m | $``$5.8 | 97/04/10 | 3.89 | 15h 58m | 23.7 |
| 97/04/28 | 4.20 | 2h 18m | 13.8 | 97/11/20 | 7.21 | 11h 09m | 41.8 |
| 98/02/06 | 4.11 | 9h 47m | 23.7 | 98/03/30 | 6.93 | 17h 16m | 56.3 |
| 98/04/04 | 5.35 | 11h 13m | 56.0 | 98/06/12 | 12.0 | 23h16m | 12.3 |
| 98/09/03 | 4.69 | 19h36m | 50.7 | 98/10/27 | 6.11 | 3h45m | 44.9 |
| 99/01/22 | 7.53 | 19h11m | 5.3 | 99/07/22 | 4.09 | 7h39m | 32.2 |
| 99/07/28 | 7.16 | 3h46m | 49.5 | 99/09/22 | 10.4 | 23h03m | 33.9 |
| 99/09/25 | 4.95 | 22h40m | 42.6 | 99/10/20 | 6.19 | 4h37m | 5.1 |
| 99/10/20 | 4.29 | 4h02m | 51.7 | 00/05/26 | 4.98 | 14h08m | 37.1 |


text:
# Lorentz violations in canonical quantum gravity

## I Introduction

It has been recognized for some time that applying the rules of quantum mechanics to general relativity is problematic . The fact that general relativity has a dimensionful coupling constant immediately implies that a straightforward perturbative approach will lead to a non-renormalizable theory. The current majority-held point of view towards the problem of quantum gravity is that general relativity is pathological and it has to be replaced by a more general theory to be quantized. Current candidates include “M-Theory”, a non-perturbative superset of string theories .

A group of researchers has over the last decade devoted themselves to reexamine the original question: is general relativity impossible to quantize? The main rationale for this is that there is a feeling that the question has not been examined in a detailed enough way. Yes, apparent non-renormalizability suggests a problem, but it is not clear if this is a problem of the theory itself, or of the way in which perturbation theory is being implemented. A striking example of this distinction is given by gravity in $`2+1`$ dimensions . In $`2+1`$ dimensions the Einstein equations imply space-time is locally flat. The only possible degrees of freedom of gravity have to do with global properties of the space-time. The resulting theory is not even a field theory proper: its reduced phase space is finite dimensional. Nevertheless, if one does a straight application of perturbation theory to the quantization of general relativity in $`2+1`$ dimensions, the theory appears non-renormalizable for the same reasons as gravity in $`3+1`$ dimensions. Witten was the first to point out that a non-perturbative quantization should be possible, and after that was accomplished, other researchers showed that one could indeed study the theory perturbatively. The overriding question is: could we be in the same situation with gravity in $`3+1`$ dimensions? Other examples of theories that do not appear to exist perturbatively but do non-perturbatively are known .

It is not surprising that quantizing general relativity should be tricky. The theory is invariant under space-time diffeomorphisms, reflecting the coordinate invariance of general relativity. Most traditional quantum field theory techniques are not well suited for diffeomorphism invariant theories, since they break such invariance. We have only comparatively recently started to understand how to quantize diffeomorphism invariant theories, but most of the results have been for “topological field theories” (again, theories with a finite number of degrees of freedom) . With almost no exception, in general relativity is the first time we confront the quantization of a theory that is both diffeomorphism invariant and also has an infinite number of degrees of freedom.

Non-perturbative quantizations of general relativity (both using canonical and path integral techniques) were first attempted in the 60’s and 70’s. Most of these attempts ended unsuccessfully. The expressions encountered either for the quantum Hamiltonian constraint or the path integral were ill-defined and there was no idea of how to regularize them. A revival of interest took place when Ashtekar  discovered a new set of variables in terms of which general relativity resembled a Yang–Mills theory. This opened the possibility of applying techniques of non-perturbative quantization that were successfully applied in the Yang–Mills context. The new variables introduced by Ashtekar consist of replacing the spatial metric (in canonical quantization, the theory is based on studying a slice of space as it evolves in time) by a set of three vector fields (triads), $`E_i^a`$ (the reader can pretend one is dealing with the three electric fields of an $`SU(2)`$ Yang–Mills theory), and the conjugate momentum behaves like an $`SU(2)`$ connection $`A_a^i`$. In terms of these variables the constraints become polynomial and include among them a Gauss law. That is, one can view the phase space of general relativity as a submanifold of that of Yang–Mills theory.

Having the theory cast in terms of a Yang–Mills-like connection, allows one to expand the wavefunctions in terms of traces of holonomies (Wilson loops). The coefficients of the expansion constitute the “loop representation”, which had been studied for Yang–Mills theories by Gambini and Trias  in the 80’s and was introduced into the gravitational context by Rovelli and Smolin . Expressing things in terms of loops is quite attractive in the gravitational context, since one can embody diffeomorphism invariance in a remarkably clear way: one simply considers functions of loops that are invariant under deformations of the loops. Such functions have been studied by mathematicians since the time of Gauss, and are called knot invariants .

Ashtekar and Lewandowski  introduced a set of mathematical tools to deal with gauge-invariant diffeomorphism invariant functions of a connection. These tools consist of a set of simple functions (cylindrical functions) that can be rigorously endowed with a measure. This is highly non-trivial, since it is a measure of integration in an infinite-dimensional space. There are few examples of such measures available. One now therefore has a definite, mathematically controlled setting in which to construct a theory of quantum gravity.

The basis of Wilson loops is overcomplete, and for many years this caused problems when working in the loop representation. Rovelli and Smolin  showed that one could find an elegant way of labeling the independent elements of the basis of loops in terms of “spin networks”. The latter are a mathematical construction first introduced by Penrose consisting of embedded graphs of lines intersecting at multivalent vertices, each line labeled by a representation of $`SU(2)`$. The Ashtekar–Lewandowski measure is particularly simple to understand in terms of spin networks.

Several attractive results have been obtained using these techniques. For instance, the spectrum of the quantum operators associated with the area of a surface and of the volume  of a region were studied in detail, and found to be discrete. Expressions for the Hamiltonian constraint of gravity that are finite and well defined were also proposed . This constitutes the first proposal ever of a non-trivial, well defined theory of quantum gravity. The present efforts are concentrating on demonstrating that the theory constructed has the correct classical and semi-classical limit.

There has also been progress in applying these mathematical tools to define a path-integral quantization of general relativity . Impressive recent results have found a certain sense in which the Lorentzian path integral can be made finite.

## II Lorentz violations

Why should there be Lorentz violations in non-perturbative quantum gravity? The answer can only be tentative, since we do not have a definite answer concerning the semi-classical limit of the theory. Fortunately, progress in this area is happening at a significant pace  and I expect in near future CP meetings a more definitive answer will be present.

To see why one might expect violations, let us consider at a very broad and heuristic sense how would one construct a semi-classical limit of the theory. When one has a quantum theory under control, one should be able to build coherent states that are peaked around a solution of the classical equations of motion, and to therefore study quantum fluctuations. Since gravity is now under control, we can build such states. However, there is an unusual twist when one carries out this straightforward idea in the case of gravity. The twist is that presumably one does not want to consider a coherent state that approximates any classical metric. The reason for this is that one can obviously construct many classical metrics that do not correspond to a “classical” situation at all (for instance a metric with gravitational waves of wavelength shorter than the Planck length). Nothing prevents one from mathematically constructing such coherent states. One should therefore introduce into the semi-classical construction a proviso that one will only consider states that approximate metrics that are “reasonably smooth” at a certain length scale. This length scale is absolute. Therefore, unless one elaborates further, Lorentz violations will occur.

As I stated before, this subject is rapidly evolving. Amelino-Camelia has several articles on this point and I would refer the reader to them .

## III Concrete calculations

It is evidently that more work is needed before we can perform a concrete calculation showing Lorentz violations. Up to present only a handful of preliminary investigations have been carried out. In a piece of work with Gambini  we studied the coupling of gravity to Maxwell theory when gravity is treated quantum mechanically. We took a naïve point of view in which we assumed that the quantum state was given by a spin network state and studied the coupling Hamiltonian pretending that the Maxwell fields were in a semi-classical coherent state and only kept the leading terms. This calculation is inadequate at many levels: first of all we did not consider a quantum state that satisfies all constraints, and considered a semi-classical limit of them. We just considered a generic spin network state and treated the Maxwell theory as a regular field theory. Under these assumptions, we showed that if the spin network state had certain properties (namely it was parity-violating), one got Lorentz violating corrections to Maxwell theory, that implied a birefringence of space-time at the level of $`\mathrm{}_p/\lambda `$ where $`\mathrm{}_P`$ is Planck’s length and $`\lambda `$ is the wavelength of the light. Such effect could influence the light arriving from gamma-ray bursts and be within a few orders of magnitude of observation (very much like the string inspired violations of Amelino-Camelia et al. ). It was later noted that such birefringence is severely constrained by optical observations of polarized sources . Alfaro et al.  have performed a reexamination of these calculations, with a more careful motivation for the quantum states considered, and found certain discrepancies in the results with those of reference. The same authors  have also extended these results to propagations of Fermions, finding possibly observable effects in the time arrival of neutrinos. The Fermionic calculations have a further drawback in that a given expectation value for the connection has to be assumed in the quantum state. The results are therefor highly dependent on this value, which is not currently fixed by any a priori argument. There is rather general consensus that effects at higher orders in $`\mathrm{}_P/\lambda `$ than the one we found are likely be present. Such effects appear to weak to be observed.

## IV Conclusions

At present it appears possible that the canonical quantization of gravity will, in its semi-classical limit, contain violations of Lorentz symmetry. Some preliminary heuristic calculations have exhibited this effect. The encouraging part is that rigorous semi-classical studies are currently under way and it is expected that in the next few years these questions will be given a definitive answer.

###### Acknowledgements.

This work was supported in part by grants NSF-PHY0090091, NSF-PHY-9800973, NSF-INT-9811610, and by funds of the Horace C. Hearne Jr. Institute for Theoretical Physics.