text:
# A model of cardiac tissue as a conductive system with interacting pacemakers and refractory time

## 1 Introduction

One of the remarkable examples of excitable media is the cardiac tissue. Because the stability of its behaviour is essential for living creatures, investigations of processes occurring in the cardiac muscle attract a considerable interest of various scientists. Owing to a great complexity and a wide variety of the processes taking place in the heart, it can be considered with several points of view. One of them relies upon the interpretation of the cardiac tissue as an active conductive system. Then, the cardiac rhythms are described on the basis of the dynamical systems theory (see, for example \[Glass et al., 2002\]).

The excitation wave in the cardiac tissue originates in the sinoatrial node (SA) and spreads rapidly in succession over the right atrium, the left atrium, then the atrioventricular node (AV), bundle of His and Purkinje fibers, and finally to the walls of the right and left ventricles. The normal rhythm of the heart is determined by the activity of the SA node which is called the leading pacemaker (a source of concentric excitation waves) or the first order driver of the rhythm. In addition to the SA node cells, the other parts of the cardiac conductive system reveal an automaticity. So, the second order driver of the rhythm is located in the AV conjunction. The Purkinje fibers are the rhythm driver of the third order. Moreover, for some pathological states of the heart arising of ectopic pacemakers is typical.

Appearance of several excitation sources leads to various disorders in the cardiac rhythm, i.e. arrhythmias \[Schamorth, 1980; Marriot & Conover, 1983; Zipes & Jalife, 1985; Winfree, 1987; Glass & Mackey, 1988\]. Because arrhythmias are dangerous diseases of the heart, their investigations have a great importance. Analysis of the complex cardiac rhythms on the basis of their interpretation as chaotic phenomena can give a clue to the problem of controllability of the complex cardiac dynamics and removing the cardiac tissue to the required regime \[Goldberger & Rigney, 1988; Goldberger, 1990; Garfinkel et al., 1992\].

As known, some arrhythmias can be presented by the interaction of spontaneous nonlinear sources \[Glass et al., 1983; Glass & Mackey, 1988; Kremmydas et al., 1996\]. In some cases, a model of impulse systems describing certain types of such arrhythmias is constructed in the framework of the theory of dynamical systems. This model is represented by coupled low-dimensional maps (circle maps) which can exhibit a complex dynamics. Such an approach is based on the fact that experiments on periodically stimulated cardiac cells are in a close agreement with the dynamics predicted by one-dimensional circle maps \[Courtemanche et al., 1989\].

In our investigations, a quite general model of two nonlinear coupled oscillators describing certain types of cardiac arrhythmias is constructed. The model turns out to be a universal one in the sense that it does not depend on the specific form of interactions, i.e. on the phase response curve (PRC). The experimentally obtained PRC is approximated by a certain polynomial function with plateau. This plateau describes a refractory time when the system does not respond to an external action. Note that the refractory time plays an important role for the normal cardiac functioning. For example, the refractoriness extends over almost the whole period of the cardiac contraction protecting the myocardium from premature heartbeats caused by an external perturbation. The refractoriness provides also the normal sequence of an excitation propagation in the heart tissue and the electrical stability of the myocardium \[Marriot & Conover, 1983; Zipes & Jalife, 1985; Winfree, 1987; Glass & Mackey, 1988\]. In the proposed model, taking into consideration the refractory time possible areas of the phase lockings are investigated. Phenomena of the splitting of resonance tongues and superposition of the synchronization areas are found. Using the obtained results we can define dynamics of the excitable media with two active pacemakers depending on the type and intensity of their interaction and the initial phase difference. Moreover, generalizing the principles of our construction one can develop a quite general theory of excitable media with interacting pacemakers under external actions. This fact has a great practical importance because admits to realize the control of the cardiac rhythms by external stimuli.

## 2 Heart Tissue as a Dynamical System

In certain cases cardiac arrhythmias can be described as an interaction of two spontaneously oscillating nonlinear sources. This interaction can be considered as an influence of some external periodic perturbation on a nonlinear oscillator (with the constant amplitude and frequency). So, for the description of such a situation it is possible to use the well-known circle map \[Glass et al., 1983; Bub & Glass, 1994; Kremmydas et al., 1996\]:

$$x_{n+1}=x_n+f\left(x_n\right)\left(\text{mod  1}\right),$$

where $`x_n`$ is a phase difference in oscillators and the function $`f\left(x\right)`$ determines a change in the phase after action of stimulus. This function is called a phase response curve (PRC).

One of the most important characteristics of the circle map is a rotation number $`\rho `$. It is defined as follows:

$$\rho =\underset{n\mathrm{}}{lim}\frac{x_nx_0}{n}.$$

For stable phase locking $`N:M`$ the rotation number is rational, $`\rho =N/M`$. If it is irrational the system behaviour is quasiperiodic or chaotic.

Analysing dynamics of the constructed model based on the circle map, it is necessary to find a proper analytical approximation of the experimentally obtained phase response curve. This allows us to investigate the basic features of the behaviour of the considered system.

Experiments on the recording of phase shifts have been carried out for a quite large number of various systems. First of all, we are interested in the PRC experimentally obtained from the research of some cardiac tissues. In \[Weidmann, 1961; Jalife & Moe, 1976\] measurements of the cycle durations of the spontaneous beating Purkinje fibres after stimulation by short electric current pulses have been performed. The found phase response curve is shown in Fig.1 (dotted lines). Taking into account this experimental material, it is possible to make the following general conclusions \[Weidmann, 1961; Jalife & Moe, 1976; Glass et al., 1986; Glass & Mackey, 1988\]:

* after perturbation the rhythm is usually restored (after some transient time) with the same frequency and amplitude, but the phase is shifted;
* depending on a phase the single input can lead to either increasing or decreasing of the period of a perturbed cycle;
* at some amplitudes of stimulus the obvious breaks appear.

The basic feature of any approximation of the PRC is the dependence on two physical parameters: on the amplitude of stimulus and the input phase. In the ideal case the other (so-called ”internal”) parameters can be reduced to them.

Taking into account the polynomial function for the approximation of the PRC, we construct a model of two mutually interacting impulse active oscillators.

## 3 Analytical Model with a Mutual Influence of Impulses

Let us consider the system of two nonlinear interacting oscillators (Fig.2). Suppose that the pulse of the first oscillator with period $`T_1`$ beats at $`t_n`$, and the pulse of the second oscillator (with period $`T_2`$) beats at $`\tau _n`$. Then the moments of time of the next appearance of the impulses are defined as follows:

$$\{\begin{array}{cccccccccccccccccccc}t_{n+1}=t_n+T_1,& & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & \\ \tau _{n+1}=\tau _n+T_2.& & & & & & & & & & & & & & & & & & & \end{array}$$

Now, taking into account the change in the period of the first oscillator under the influence of the second impulse by the value of $`\mathrm{\Delta }_1\left(\left(\tau _nt_n\right)/T_1\right)`$, one can get that $`t_{n+1}=t_n+T_1+\mathrm{\Delta }_1\left(\left(\tau _nt_n\right)/T_1\right)`$. By the same manner, for the second oscillator $`\tau _{n+1}=\tau _n+T_2+\mathrm{\Delta }_2\left(\left(t_{n+1}\tau _n\right)/T_2\right)`$. Dividing these expressions by $`T_1`$ we arrive at the corresponding values for the phases:

$$\{\begin{array}{cccccccccccccccccccc}\phi _{n+1}=\phi _n+\frac{1}{T_1}\mathrm{\Delta }_1\left(\delta _n\phi _n\right),& & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & \\ \delta _{n+1}=\delta _n+\frac{T_2}{T_1}+\frac{1}{T_1}\mathrm{\Delta }_2\left(\frac{t_n}{T_2}+\frac{T_1}{T_2}+\frac{1}{T_2}\mathrm{\Delta }_1\left(\delta _n\phi _n\right)\frac{\tau _n}{T_2}\right).& & & & & & & & & & & & & & & & & & & \end{array}$$

Here $`\phi _n=t_n/T_1`$ is a phase of the first perturbed oscillator with respect to the unperturbed one (with period $`T_1`$), and $`\delta _n=\tau _n/T_1`$ is the phase of the second perturbed oscillator with respect to the same first oscillator with period $`T_1`$. Using parameters $`a=T_2/T_1`$ and $`\mathrm{\Delta }_1/T_1=f_1`$, $`\mathrm{\Delta }_2/T_1=f_2`$ one can write:

$$\{\begin{array}{cccccccccccccccccccc}\phi _{n+1}=\phi _n+f_1\left(\delta _n\phi _n\right),& & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & \\ \delta _{n+1}=\delta _n+a+f_2\left(\frac{1}{a}\left(\phi _n+1+f_1\left(\delta _n\phi _n\right)\delta _n\right)\right).& & & & & & & & & & & & & & & & & & & \end{array}$$

Now, omitting intermediate calculations, for the final expression of the phase difference in the oscillators we get

$$x_{n+1}=x_n+a+f_2\left(\frac{1}{a}\left(1+f_1\left(x_n\right)x_n\right)\right)f_1\left(x_n\right)\left(\text{mod  1}\right),$$
(1)
where $`x_n=\delta _n\phi _n`$.

It is obvious, that the PRC changes the form depending on the amplitude of the external stimulus. In the simplest case this dependence can be considered as the multiplicative relation. Then the phase response curves can be written as

$$f_1=\gamma h\left(x\right),f_2=\epsilon h\left(x\right),$$

where $`h\left(x\right)`$ is a periodic function and $`h\left(x+1\right)=h\left(x\right)`$. In such an assumption, the expression (1) takes the form:

$$x_{n+1}=x_n+a+\epsilon h\left(\frac{1}{a}\left(1+\gamma h\left(x_n\right)x_n\right)\right)\gamma h\left(x_n\right)\left(\text{mod  1}\right).$$
(2)
In the present paper we dwell on the investigations of the map (2) with a polynomial function $`h\left(x\right)`$. The obtained results are the continuation of our previous works concerning modeling certain cardiac arrhythmias \[Loskutov, 1994; Loskutov et al., 2002\].

## 4 Phase Diagrams for Unidirectional Coupling of Oscillators

First of all let us analyze the situation when the permanent inputs act on the nonlinear oscillator, i.e. $`f_2\left(x\right)0`$ or $`\epsilon =0`$. As an analytical approximation of the experimental curve in Fig.1, let us consider the following polynomial function:

$$h\left(x\right)=Cx^2\left(\frac{1}{2}x\right)\left(1x\right)^2.$$
(3)
The normalizing factor $`C`$ we choose in such a way that the amplitude of $`h(x)`$ is equal to $`1`$, so that $`C=20\sqrt{5}`$ (see Fig.1, solid line). Then taking into account the refractory time $`\delta `$ the map (2) we can write as follows:

$$x_{n+1}=\{\begin{array}{cccccccccccccccccccc}x_n+a,\hfill & 0x_n\delta ,\hfill & \left(\text{mod  1}\right),\hfill & & & & & & & & & & & & & & & & & \\ x_n+a+C\gamma h\left(\frac{x_n\delta }{1\delta }\right),\hfill & \delta <x_N1,\hfill & \left(\text{mod  1}\right),\hfill & & & & & & & & & & & & & & & & & \end{array}$$
(4)
where $`h()`$ is determined by (3). Now let us compare several cases with different values of the refractory time. First of all consider the case without the refractory period, i.e. $`\delta =0`$. The phase locking regions in the parametric space $`(a,\gamma )`$ obtained by numerical analysis are shown in Fig.3a. Without loss of generality, in this Figure we choose $`a[1,2]`$. Different colours define the phase locking areas with the multiplicity $`N:M`$, where $`N`$ cycles of external stimulus correspond to $`M`$ cycles of nonlinear oscillator. One can see that ”tales” of the main locking regions are slightly splitted and overlap each other at large $`\gamma `$. Note that as follows from the analysis of the system (4) with $`\delta =0.1`$ (Fig.3b), introduction of the refractory time leads to the extension of the phase locking areas and significant splitting and overlapping their ”tales”.

In Fig.4a the numerically constructed phase diagram in the case of $`\delta =0.3`$ is presented. For the comparison, in the given Figure the same $`N:M`$ stable phase lockings as in Fig.4 are shown. One can see that when the value of the refractory time is growing, the $`2:3`$ phase locking area is increasing with simultaneous decreasing of the $`1:1`$ and $`1:2`$ areas.

The phase locking regions in the case of $`\delta =0.5`$ are shown in Fig.4b. This phase diagram is qualitatively different from pictures considered above. The form of $`2:3`$ phase locking area is stretched and looks like an arrow. The forms of $`3:4`$ and $`3:5`$ regions also resemble arrows in the case of $`\delta =0.7`$ (Fig.5a). At $`\delta =0.9`$ all phase lockings are degenerated into vertical lines. This situation is presented in Fig.5b. Note that in the case of $`\delta =1`$ (i.e. the system does not respond to the external action) there is no any dependence on the stimulus amplitude $`\gamma `$.

## 5 Phase Diagrams for Systems with Mutual Interaction

In this section the system (2) at $`\delta =0.1`$ is considered. The analysis is performed in the ($`\gamma `$, $`a`$) and ($`\gamma `$, $`\epsilon `$)–parametric spaces.

### 5.1 Phase locking areas in the ($`\gamma `$, $`a`$)–space

Let us consider the case of a mutual interaction of two impulse systems. Assume that the influence of the first oscillator on the second one is small enough, for example, $`\epsilon =0.1`$. The corresponding phase diagram displaying the possible behaviour regimes of the system for $`\delta =0.1`$ is shown in Fig.6a. One can see that the mutual interaction leads to the deformation and the splitting of the phase locking areas. Note that even for small values of the amplitude of the second stimulus $`\gamma `$, the overlapping of the main phase lockings takes place. Thus, the system dynamics becomes multistable. This corresponds to the situation when the limit state of the map depends on an initial phase difference $`x_n`$. The growth of the refractory time in the model with $`\epsilon =0.1`$ leads to more deep distortion in the forms of the main tongues and disappearance of the splitting areas. If, however, we increase the influence of the first oscillator up to, for instance $`\epsilon =0.5`$, then one can see a very complicated structure with much more deep deformation of the main phase locking areas (see Fig.6b). For example, the $`1:1`$ area is degenerated into a narrow strip, whereas the $`1:2`$ phase locking area increases due to appearance of long narrow tongues.

The numerical analysis shows that at the growth of $`\epsilon `$ up to approximately $`0.5`$ the occupied by the resonance zones area becomes larger. At the same time, the shape of the phase lockings is complexified, and their location is changed. This leads to the almost full mixed picture, such that we can find zones of various multiplicity in a small neighborhood of almost any point ($`\gamma `$, $`a`$). However, at the given values of $`\epsilon `$ the self-similarly structures are clearly observed.

Additionally, we have found that at the further growth of the nonlinearity parameter $`\epsilon `$ the resonance zones decrease, occupying the less space. In this case the mixing of resonance tongues also takes place. Thus, increasing the effect of the influence of oscillators leads to the mixing of initially quite regular structure in ($`\gamma `$, $`a`$)–space.

### 5.2 Phase locking regions in the ($`\gamma `$, $`\epsilon `$)–space

Now we construct the phase diagrams of the interacting oscillators in the space of influence amplitudes, i.e. ($`\gamma `$, $`\epsilon `$). In the first instance, let us consider $`a=2`$ (Fig.7a). This value of the ratio of periods means that for $`\gamma =\epsilon =0`$ the rotation number is rational, and the dynamics of the system is periodic with the $`1:2`$ phase locking. Although at the growth of nonlinearity it is possible to obtain the phase locking areas with another multiplicity, even at the large values of $`\gamma `$ and $`\epsilon `$ the system behaviour is periodic with the $`1:2`$ phase locking.

Another situation is observed at $`a=\pi /2`$. Here the rotation number is irrational at zero stimulus amplitudes, and the system exhibits the property of quasiperiodicity or chaoticity. However, at increasing the nonlinearity the possibility of the appearance of the periodic behaviour exists (Fig.7b). Here, as for a sufficiently large $`\epsilon `$, one can observe decreasing of the area occupying the resonance zones. Therefore, at irrational values of $`a`$ the probability of the complex behaviour of the system (2) is a large enough.

## 6 Analogy with Pathological Heart Rhythms

Summarizing, we shall try to make an analogy between the obtained results and the pathological states of the cardiac tissue. Using the developed models it is possible, for example, to describe the interaction of the sinus and the ectopic pacemakers, the SA and AV nodes and impact of an external perturbation on the sinus pacemaker.

Similar case of mutual interaction of the SA and AV nodes was investigated in \[di Bernardo et al., 1998; Signorini et al., 1998\]. The authors modeled the AV node as a van-der-Pol oscillator and the SA node was considered as a certain modified van-der-Pol oscillator. The obtained waveforms satisfactory replicated the action potentials of the SA and AV node cells. A bifurcation analysis performed in these works showed a possibility to reproduce and classify various types of cardiac pathologies.

Now let us consider some types of arrhythmias one can predict on the basis of our model. If the first pulse oscillator is presented as the SA node and the second one is considered as the AV node, then we come to conclusion that some stable phase lockings correspond to the cardiac pathologies which are observed in a clinical practice. In this case among various obtained lockings one can reveal the normal sinus rhythm ($`1:1`$ phase locking). In addition, in the diagrams we can see the classical rhythms of Wenckebach ($`N:(N1)`$ phase lockings) and $`N:1`$ AV blocks.

When the first pulse system is considered as the AV node and the second one is presented as the SA node, we obtain the known inverted Wenckebach rhythms found in some patients.

Presence of the wide areas of phase lockings (see Fig.37) confirms that in such systems it is possible to observe the various kinds of synchronization of two oscillators qualitatively corresponding to some types of cardiac arrhythmias. The phase diagram allows us to reveal under what conditions of the interaction (i.e. at what values of the parameters $`a`$, $`\gamma `$, $`\epsilon `$ and $`\delta `$) one or another type of synchronization exists. Moreover, the phase pictures indicate that at the increasing the nonlinearity (i.e. at the growth of the parameter $`\gamma `$) areas with various phase lockings are overlapped. The knowledge of such regions and the dynamics permits us to remove the system from an undesirable mode of synchronization to a more appropriated state by the external action.

## 7 Concluding Remarks

In the presented paper a quite general model of two nonlinear interacting impulse oscillatory systems is developed. On the basis of this model it is possible to predict some types of cardiac arrhythmias. The constructed model is a universal one in the sense that it does not depend on the chosen interaction type, i.e. on the form of the phase response curve. Taking into account the refractory time the possible phase locking regions of the polynomial maps which describe a nonlinear oscillator under the permanent inputs, are investigated. It is found that involving the refractory time leads to the extension of the phase locking areas, significant splitting and overlapping their ”tales”. Moreover, the forms of the phase locking areas are stretched and degenerated into vertical lines as the refractory time tends to one.

Detailed analysis of the phase diagram of the system with two mutually interacting oscillators in the ($`\gamma `$, $`a`$)–space shows that besides splitting of the central tongues there is an overlapping of the main regions of synchronization which corresponds to various types of cardiac arrhythmias. This bistability is observed even for small enough values of the interaction. Increasing the refractory time leads to the distortion in the form of the main tongues and disappearance of the splitting areas. For sufficiently large values of the interaction we obtain a very complicated picture, where the phase locking areas are interwoven with each other.

In addition, in the constructed model the phase lockings in the space of the stimulus amplitudes are observed. It is found that the interacting oscillators can be synchronized even if the ratio of their periods is irrational (note, that the probability of this phenomenon is a quite small). However, in the case without coupling this would correspond only to the complex dynamics (quasiperiodic or chaotic).

The obtained results allow us to predict the dynamics of oscillatory systems depending on the initial phase difference, on the type and the intensity of the interaction. Moreover, using the principle of the construction of the model one can develop a quite general theory of the interacting oscillators under periodic perturbation. In this case the knowledge of the multistability areas can help to stabilize the system dynamics and remove the cardiac tissue to the required type of the behaviour.

Bub, G. & Glass, L.  “Bifurcations in a continous circle map: A theory for chaotic cardiac arrhythmia,” Int. J. Bifurcation and Chaos 5(2), 359–371.

Courtemanche, M., Glass, L., Belair, J., Scagliotti, D. & Gordon, D.  “A circle map in a human heart,” Physica D 49, 299–310.

di Bernardo, D. D., Signorini, M. G. & Cerutti, S.  “A model of two nonlinear coupled oscillators for the study of heartbeat dynamics,” Int. J. Bifurcation and Chaos 8(10), 1975–1985.

Garfinkel, A., Spano, M. L. & Ditto, W. L.  “Controlling cardiac chaos,” Science 257, 1230–1235.

Glass, L., Guevara, M. R., Shrier, A. & Perez, R.  “Bifurcation and chaos in a periodic stimulated oscillator,” Physica D 7(1–3), 89–101.

Glass, L., Guevara, M. R. & Shrier, A.  “Phase resetting of spontaneously beating embryonic ventricular heart cell aggregates,” Am. J. Physiol. 251, H1298–H1305.

Glass, L. & Mackey, M.  From clocks to chaos: the rhythms of life (Princeton Univ. Press, Princeton).

Glass, L., Nagai, Yo., Hall, K., Talajie, M. & Nattel, S.  “Predicting the entrainment of reentrant cardiac waves using phase resetting curves,” Phys. Rev. E 65, 021908-1–021908-10.

Goldberger, A. L. & Rigney, D. R.  “ Sudden death is not chaos,” Dynamic patterns in complex systems, eds. Kelso, J. A. S., Mandell, A. J. & Schlesinger, M. F. (World Sci. Pub., Teaneck, NJ), 248–264.

Goldberger, A. L.  “Nonlinear Dynamics, Fractals and Chaos: Applications to Cardiac Electrophysiology,” Ann. Biomed. Eng. 18(2), 195–198.

Jalife, J. & Moe, G. K.  “Effects of electronic potencials on pacemaker activity of canine Purkinje fibers in relation to parasistole,” Circ. Res. 39(6), 801–808.

Kremmydas, G. P., Holden, A. V., Bezerianos, A. & Bountis, T.  “Reprezentation of sino-atrial node dynamics by circle maps,” Int. J. Bifurcation and Chaos 6(10), 1799–1805.

Loskutov, A.  “Nonlinear dynamics and cardiac arrhythmia,” Applied Nonlinear Dynamics 2(3–4), 14–25 (Russian).

Loskutov, A., Rybalko, S. D. & Zhuchkova, E. A.  “Dynamics of excitable media with two interacting pacemakers,” Biophysics 47(5), 892–901.

Marriot, H. J. L. & Conover, M. M.  Advanced Concepts in Cardiac Arrhythmias (C.V. Mosby, St. Louis).

Schamorth, L.  The Disoders of the Cardiac Rhythm (Blackwell, Oxford).

Signorini, M. G., Cerutti, S. & di Bernardo, D. D.  “Simulation of hertbeat dynamics: a nonlinear model,” Int. J. Bifurcation and Chaos 8(8), 1725–1731.

Weidmann, S.  “Effects of current flow on the membrane potential of cardiac muscle,” Journal of Physiology 115, 227.

Winfree, A. T.  When Time Breaks Down: The Three-Dimensional Dynamics of Electrochemical Waves and Cardiac Arrhythmias (Princeton Univ. Press, Princeton).

Zipes, D. P. & Jalife, J.  Cardiac Electrophysiology and Arrhythmias, (Grune and Stratton, Orlando).

## Figure captions

Fig.1. The phase response curves: the experimental curve (dotted line) and its analytical approximation (solid line). The experimentally obtained phase response curve shows a dependence of the duration of perturbed cycle (in $`\%`$) on the phase of input.

Fig.2. A construction of the model of two nonlinear interacting oscillators.

Fig.3. The phase diagram of the map (4): $`a`$) $`\delta =0`$; $`b`$) $`\delta =0.1`$.

Fig.4. The phase locking areas of the map (4): $`a`$) $`\delta =0.3`$; $`b`$) $`\delta =0.5`$.

Fig.5. The phase diagram of the map (4): $`a`$) $`\delta =0.7`$; $`b`$) $`\delta =0.9`$.

Fig.6. The phase locking regions of the system of two mutual interacting oscillators with $`\delta =0.1`$: $`a`$) $`\epsilon =0.1`$; $`b`$) $`\epsilon =0.5`$.

Fig.7. The phase lockings in the space of stimulus amplitudes ($`\delta =0.1`$): $`a`$) $`a=\pi /2`$; $`b`$) $`a=2`$.


format: text

text:
# Landau Ghosts and Anti-Ghosts in Condensed Matter and High Density Hadronic Matter

## 1 Introduction

It was discovered many decades ago, that the dimensionless coupling constant in QED, $`\alpha _0=e_0^2/\mathrm{}c`$, after 1-loop renormalization becomes momentum transfer $`Q^2`$ dependent $`\alpha (Q^2)`$ and that it has a pole at a (space-like $`Q^2>0`$) value which is clearly unphysical. A Landau ghost (or pole) is the name associated with such a singularity. If $`\alpha (Q^2)`$ can turn negative either by passing through a zero or through a pole, it would lead to instability separating a phase where two like charges repel to a different phase where they attract. In the present work, we shall explore a general origin of the ghost (and similar pathologies) via the Dyson equation. While for QED, the phase transition is at unattainable high energies, in certain condensed matter situations the singularity is at very low energies and leads to an observable phase transition called super-radiance. A very practical application is to the dielectric properties of systems for which a one loop renormalization leads to Clausius-Mossotti(CM) and Lorenz-Lorentz(LL) equations. For such phenomenological models, the phase transition computations can be done to completion  and we find quantitatively that water is indeed super-radiant at room temperature, a result predicted earlier by other authors from somewhat different theoretical considerations.

In contrast to QED with a Landau ghost, in the same 1-loop approximation, a non-abelian gauge theory such as QCD for quarks and glue avoids the Landau ghost (provided the number of quark flavors is limited to 16) and is asymptotically free (AF). However, this welcome change in sign of the beta function leading to AF brings with it its own pathology. The Landau pole is now no longer a pole on the real axis in the physical region ($`s=Q^2>0`$) since it acquires an imaginary part, but the sign of the imaginary part is wrong i.e., it is negative. It seems appropriate to call this latter phenomenon an “anti-ghost” since, as we shall see later it also signals an instability. Physically, it leads to the pleasing result that quark and gluon states can only be transient realized under restricted conditions (of large energies and momentum transfers). We shall not consider here how color confinement with normal hadronic phase prevails at low energies and momentum transfers.

The paper is organized as follows. In Sec.2, the genesis of the singularity is considered through a generic discussion of “feedback” which is at the heart of the Dyson equation. In Sec.3, we first develop the CM and LL equations for dielectrics in order to exhibit the resulting phase transition as critical values (of the volume and/or the temperature) leading to super-radiance are reached for the important case of water (and other polar materials). In Sec.4, we briefly discuss QCD and what the anti ghost entails for the resulting physics there. Sec.5 closes the paper with some concluding remarks and directions for future work.

## 2 Dyson equation and the critical feedback

Let $`D`$ be the exact photon propagator, $`D_0`$ the free photon propagator and $`\mathrm{\Pi }`$ its proper self energy part (suppressing all Lorentz indices and momentum labels momentarily). The Dyson equation then reads

$$D=D_0+D_0\mathrm{\Pi }D=(1D_0\mathrm{\Pi })^1D_0,$$
(1)
or, equivalently written in terms of the inverse propagators

$$D^1=D_0^1\mathrm{\Pi }.$$
(2)
For computations to follow, it would be useful to recast the logic behind Eq.(1) in the language of linear feedback. Suppose an external current $`\delta J_{ext}`$ is applied to a system and we wish to find the corresponding change in the potential from $`A`$ to $`A+\delta A`$. The potential change $`\delta A`$ is given in terms of the total change in the current $`\delta J`$ through the free propagator $`D_0`$ by

$$\delta A=D_0\delta J=D_0\left(\delta J_{ext}+\delta J_{ind}\right),$$
(3)
where in Eq.(3), the total change in the current contains a feedback current $`\delta J_{ind}`$ in addition to the applied $`\delta J_{ext}`$. The induced piece is the system’s linear response to the change in the potential through the “self energy” term $`\mathrm{\Pi }`$

$$\delta J_{ind}=\mathrm{\Pi }\delta A.$$
(4)
Now defining the exact propagator through

$$\delta A=D\delta J_{ext},$$
(5)
Eq.(1) follows from Eqs.(3-5). Eq.(1), when depicted pictorially as in Fig., makes abundantly clear the linear feedback mechanism implicit in the Dyson equation.

Restoring explicitly the Lorentz indices, we may write the Feynman gauge, renormalized photon propagator of 4- momentum $`q`$ with

$`s=q^2=(\omega /c)^2|𝐤|^\mathrm{𝟐}`$
($`\omega `$ the photon frequency and $`𝐤`$ its wave-vector), in terms of the “running coupling constant” $`\alpha (s)`$ and the dielectric “tensor” $`ϵ(s)`$ as follows:

$$\alpha _0D^{\mu \nu }(q)=\alpha (s)D_0^{\mu \nu }(q)=\frac{\alpha _0D_0^{\mu \nu }(q)}{ϵ(s)}.$$
(6)
Clearly, both poles and zeroes in $`\alpha `$ (which become respectively, the zeroes and poles of $`ϵ`$) are pathological.

For a renormalization at an arbitrary space-like point ($`s=\mathrm{\Lambda }^2`$),

$$\mathrm{\Pi }_{\mu \nu }(q)=[4\pi \alpha (\mathrm{\Lambda }^2)](q_\mu q_\nu \eta _{\mu \nu }q^2)\pi _c(s).$$
(7)
With $`\pi _c(s=\mathrm{\Lambda }^2)=0`$ we have

$$\alpha (s)=\left[\frac{\alpha (\mathrm{\Lambda }^2)}{1+4\pi \alpha (\mathrm{\Lambda }^2)\pi _c(s)}\right]$$
(8)
and likewise

$$ϵ(s)=1+4\pi \alpha (\mathrm{\Lambda }^2)\pi _c(s).$$
(9)
Now let us consider the allowed singularities in the exact propagator $`D_{\mu \nu }(q)`$. There are branch point singularities in $`D`$ (for time-like momenta at all production thresholds) through unitarity as $`\pi _c(s)`$ develops an imaginary part. In the physical region of time-like momenta, there may even be poles corresponding to resonances on the unphysical Riemann sheet(s). On the other hand, for space-like momenta, $`D`$ is forbidden to have any poles, i.e., on general grounds, we must require that

$$4\pi \alpha (\mathrm{\Lambda }^2)\pi _c(s)>1,\mathrm{all}s<0,$$
(10)
(all space-like momenta for which $`\omega <c|𝐤|`$). Whenever this inequality is not satisfied, we shall have an unphysical Landau pole. The feedback is then beyond its critical value invalidating linear response and it would signal a phase transition to another, super-radiant phase of the theory.

For example, in one fermion loop QED, renormalized on shell at $`s=0`$, the standard vacuum result reads

$$ϵ_{1loop}(s)=1\frac{2\alpha }{\pi }_0^{\mathrm{}}\left(zz^2\right)\mathrm{ln}\left[1\frac{sz(1z)}{\kappa ^2i0^+}\right]𝑑z,$$
(11)
where $`\kappa =mc/\mathrm{}`$, with $`m`$ the fermion mass. In the space-like region $`ϵ`$ decreases to zero at the Landau value

$$s_L=Q_L^2=\kappa ^2e^{3\pi /\alpha },$$
(12)
beyond which it turns negative. Thus, QED is in its normal (stable) phase until an enormously large $`Q^2`$ value above which it would turn super-radiant. In the condensed matter example treated in the next section, a similar phenomenon would occur (this time through a zero in $`\alpha `$) but at physically relevant energies and hence of practical importance.

A discussion of the complementary notion of the anti-ghost is better done via a subtracted dispersion relation for $`ϵ`$

$$ϵ(s)=1+\left(\frac{s+\mathrm{\Lambda }^2}{\pi }\right)_0^{\mathrm{}}\frac{(ds^{})\mathrm{}mϵ(s^{})}{(s^{}+\mathrm{\Lambda }^2)(s^{}si0^+)},$$
(13)
In spinor QED, where only fermion loops contribute, $`\mathrm{}mϵ(s)`$ is positive definite ( for all physical $`s`$). In one loop, it is given by

$$\mathrm{}mϵ(s)=\frac{\vartheta (s4\kappa ^2)\alpha }{3}\left[1+\frac{2\kappa ^2}{s}\right]\sqrt{\left[1\frac{4\kappa ^2}{s}\right]},$$
(14)
thus the system is “dissipative”. Whenever the system is dissipative, the subtraction tells us ( in 1-loop) that there is a Landau ghost in the space like region.

On the other hand, for QCD, along with the fermion loops (quarks in this case) as above, which are dissipative, we also have a larger imaginary part with a negative sign from the gluons

$$\mathrm{}mϵ(s)=\mathrm{}mϵ_{glue}(s)+\mathrm{}mϵ_{quark}(s),$$
(15)
where

$$\mathrm{}mϵ_{glue}(s)=\vartheta (s)\left(\frac{11\alpha _s}{8}\right),$$
(16)
and

$$\mathrm{}mϵ_{quark}(s)=\underset{i=1}{\overset{N_f}{}}\vartheta (s4\kappa _i^2)\left(\frac{\alpha _s}{12}\right)\left[1+\frac{2\kappa _i^2}{s}\right]\sqrt{\left[1\frac{4\kappa _i^2}{s}\right]}$$
(17)
For the number of flavors $`N_f<16`$, the total imaginary part remains negative definite. This cures the Landau ghost problem and as well known we have asymptotic freedom, for in the space-like region $`s=Q^2<\mathrm{\Lambda }^2`$,

$$ϵ(Q^2)=1+\frac{(Q^2\mathrm{\Lambda }^2)}{\pi }_0^{\mathrm{}}\frac{ds\left[\mathrm{}mϵ(s)\right]}{(s+\mathrm{\Lambda }^2)(s+Q^2)}>1.$$
(18)
The price paid for exorcizing the Landau ghost into a benevolent anti-ghost is to turn the system into being non-dissipative i.e., into an “amplifier” which is well known to be transient. We shall return to this problem in Sec.4.

## 3 Landau ghost in condensed matter and super-radiance

As stated earlier, unlike the Landau ghost in QED, in condensed matter we have unphysical singularities at low frequencies which signal phase transitions of practical importance for the physics and the chemistry of substances.

In condensed matter literature, the analogue of $`\alpha \pi _c`$ in Eq.(8) is known as the electric susceptibility $`\chi _e`$ and the object here is to relate it to macroscopic quantities such as the polarizability of the molecules $`\widehat{\alpha }`$ and $`v`$ the volume per molecule.( We have defined the polarizability with a “hat” in order to distinguish it from the coupling constant $`\alpha `$).

The CM equation is an ideal example of a linear feedback paralleling the argument of the last section. To briefly recall what is involved and to set the notation, the electric dipole moment per unit volume or, polarization $`𝐏`$

$$𝐏=\chi _e𝐄=N𝐩_{mol},$$
(19)
where $`N`$ is the number of molecules per unit volume and $`𝐩_{mol}`$ is the mean dipole moment of the molecules. Then, without the “feedback”, justifiable in rare media where molecular separations are large, one assumes that the macroscopic field $`𝐄`$ is the same as that acting on a given molecule

$$𝐩_{mol}=\widehat{\alpha }𝐄,$$
(20)
in terms of the polarizability. So that for rare media

$$\chi _e=N\widehat{\alpha }=\widehat{\alpha }/v,$$
(21)
where $`v`$ is the volume occupied by a single molecule.

In denser media where molecular separations are small, the polarization $`𝐏`$ induces an internal electric field $`𝐄_i`$. If one includes only the contributions from molecules far away, this feedback term is estimated to be

$$𝐄_i=\left(\frac{4\pi }{3}\right)𝐏.$$
(22)
Now, Eq.(19) leads to

$`𝐏`$ $`=`$ $`N𝐩=N\widehat{\alpha }\left(𝐄+𝐄_i\right)`$ (23)
$`=`$ $`N\widehat{\alpha }\left(𝐄+{\displaystyle \frac{4\pi }{3}}𝐏\right)`$
$`=`$ $`\left[{\displaystyle \frac{N\widehat{\alpha }}{1(4\pi /3)N\widehat{\alpha }}}\right]𝐄.`$
It is perhaps not surprising therefore, that the CM equation gives the same expression for $`\chi _e`$ as is obtained from the lowest order mean field theory

$$\chi _e=\left(\frac{(\widehat{\alpha }/v)}{1(4\pi \widehat{\alpha }/3v)}\right),$$
(24)
and for $`ϵ`$ one finds

$$ϵ=\left(\frac{1+2(4\pi \widehat{\alpha }/3v)}{1(4\pi \widehat{\alpha }/3v)}\right).$$
(25)
The normal phase would require for its stability that the static polarizability satisfy

$$4\pi \widehat{\alpha }<3v,$$
(26)
in order for no poles to develop in $`ϵ`$ (or $`\chi _e`$). On the other hand, as the molecular density decreases below a critical value $`v_c`$, i.e., for

$$v<v_c=\left(\frac{4\pi }{3}\right)\widehat{\alpha },$$
(27)
the normal phase would become unstable yielding to a new stable super-radiant phase of Dicke-Lieb- Preparata models.

CM equations become the Lorenz-Lorentz (LL) equations  as static $`\widehat{\alpha },\chi _e\mathrm{and}ϵ`$ get generalized to a finite frequency $`\omega `$. Here one finds giant (macroscopic) Lamb shifts in frequency due to the “feedback”. Let us suppose that a single atom with atomic number $`Z`$ has a strong resonance at a frequency $`\omega _{\mathrm{}}`$, so that

$$\widehat{\alpha }(\omega )\left(\frac{Ze^2}{m}\right)\frac{1}{\omega _{\mathrm{}}^2\omega ^22i\gamma \omega },$$
(28)
where $`m`$ is the electron mass. Eqs(24) and (28) then imply

$$\chi _e(\omega )\left(\frac{Ze^2}{m}\right)\frac{1}{\omega _0^2\omega ^22i\gamma \omega },$$
(29)
where the “renormalized” resonant frequency $`\omega _0`$ is given by

$$\omega _0^2=\omega _{\mathrm{}}^2\left(\frac{4\pi Ze^2}{3mv}\right)=\omega _{\mathrm{}}^2\left[1(\frac{4\pi \widehat{\alpha }}{3v})\right],$$
(30)
where $`\widehat{\alpha }\widehat{\alpha }(\omega =0)`$

At the critical temperature $`T_c`$, the static polarizability stability condition given by Eq.(26) becomes

$$\left(\frac{4\pi \widehat{\alpha }(T_c)}{3}\right)=v(\mathrm{critical}\mathrm{point}).$$
(31)
Dynamically, the renormalized frequency in Eq.(30) is driven to zero at the critical temperature

$$\omega _0\left(TT_c+0^+\right)=0,$$
(32)
which again leads to Eq.(31).

In previous work, we have obtained the phase diagram through (2 and 4 level) model Hamiltonians which allow us to compute $`\widehat{\alpha }(T)`$ for a single molecule. It is to be emphasized that perturbative corrections would shift the exact boundary of the super-radiant phase from its CM value but its occurrence is left intact.

Here we simply estimate the polarizability of water using the Debye model

$$\widehat{\alpha }=\widehat{\alpha }_o+\left(\frac{\mu ^2}{3k_BT}\right).$$
(33)
For water

$$\widehat{\alpha }_0=1.494A^3\mathrm{and}\mu =1.855\times 10^{18}Gausscm^3.$$
(34)
This gives for the polarizability of water at room temperature

$$\widehat{\alpha }_{water}28.6962\AA ^3.$$
(35)
From the molecular density of water, we find $`v30.0137\AA ^3`$, so that the super-radiance condition

$$\frac{v_c}{v}=\left(\frac{4\pi \widehat{\alpha }}{3v}\right)4.0054>1,$$
(36)
is satisfied showing that water is indeed super-radiant at room temperature. Surely we expect numerical corrections to the above result. But the “margin of safety” is reasonable so that the claim about room temperature water being super-radiant is robust.

To show the efficacy of the above criterion, we have carried a similar analysis for a large number of polar and non-polar materials. The results shown in Fig. confirm that at room temperatures, many polar liquids are super-radiant whereas non-polar materials are not.

## 4 Implications of the anti-ghost in QCD

In this section, we resume our discussion of the physical implications arising from the negative sign of $`\mathrm{}mϵ`$ computed in 1-loop QCD. In general, such a negative sign signals a lack of causality, the latter driving the system to become an “amplifier” instead of a “dissipator”. An amplifier (for example a laser with an inverted population entails the notion of a negative temperature) can only be a transient phenomenon since it requires a source to feed energy into it.

We may conclude from the above that the $`J=1`$, color octet channel of QCD represented by a gluon propagator with a negative absorptive part cannot be a part of the physical spectrum. We may conjecture that absorptive parts for all channels carrying color are of the amplifier variety thereby excluding all such states from the physical spectrum. At best then quarks and gluons can occur only transiently as virtual states in a restricted part of the phase space spanned by hadrons, the colorless children of quarks and glue. Only after one is able to produce hadrons as colorless composites of quarks and gluons, would it be possible to check whether these true inhabitants of the “normal” phase of QCD have dissipative absorptive parts. Needless to say that we are nowhere near that goal.

## 5 Conclusions and outlook

We have presented a general discussion of the Landau ghost and related phenomena through the Dyson equation and these pathologies have been found to trigger phase transitions. A convincing case has been made for polar liquids to be in the super-radiant phase at ambient temperatures substantiating previous work on water. Similar arguments for QCD, through a wrong sign absorptive part in the gluon propagator, show that the gluons (and conjectured more generally for any channel carrying color) cannot be a part of the physical spectrum. This generalization needs verification. Also, further work is needed to show that the spectral functions for hadrons (that is the color singlet bound states) in QCD do possess positive absorptive parts as they must to qualify as the true “normal” ground states of the theory.


text:
# Effects of spin and exchange interaction on the Coulomb-blockade peak statistics in quantum dots

## Abstract

We derive a closed expression for the linear conductance through a quantum dot in the Coulomb-blockade regime in the presence of a constant exchange interaction. With this expression we calculate the temperature dependence of the conductance peak-height and peak-spacing statistics. Using a realistic value of the exchange interaction, we find significantly better agreement with experimental data as compared with the statistics obtained in the absence of an exchange interaction.

The conductance through a quantum dot that is weakly coupled to leads displays sharp peaks as an applied gate voltage is varied. Each conductance peak describes the addition of one more electron into the dot, and between peaks the conductance is “blocked” by the Coulomb interaction. The statistics of both peak heights and peak spacings in dots for which the single-electron dynamics is chaotic have been intensively studied in recent years alhassid00 . Some of the experimental observations, e.g. the peak-height distributions at low temperature jalabert92 ; chang96 ; folk96 , could be explained at least qualitatively by the constant-interaction (CI) model, in which the interaction is represented in the simple form of an electrostatic charging energy. Other measured observables, such as the peak-spacing distribution patel98a , have indicated that spin and residual interactions beyond the charging energy should be taken into account. A consistent theoretical approach that provides quantitative agreement with both the measured peak-height and peak-spacing statistics is still lacking.

Recently, a universal Hamiltonian was derived kurland00 ; aleiner02  for a dot with a large Thouless conductance $`g_\mathrm{T}\sqrt{N}`$ ($`N`$ is the number of electrons). An important contribution to the interaction part of this Hamiltonian is a constant exchange interaction in addition to the usual charging-energy term. The remaining interaction terms are suppressed at large $`g_\mathrm{T}`$. Here we study the effect of the exchange interaction on the finite-temperature statistics of both peak heights and spacings. To this end, we derive a closed expression for the conductance in the presence of a constant exchange interaction (in the sequential-tunneling limit). This formula expresses the conductance in terms of quantities that characterize spinless non-interacting electrons. We then calculate the finite-temperature peak-height and peak-spacing statistics and find them both to be sensitive to the exchange interaction. Using an RPA estimate of the exchange interaction for the samples studied experimentally in Refs. patel98a ; patel98b , we obtain very good agreement with the observed temperature dependence of the standard deviation of the peak spacing. We also explain most of the known discrepancies between the experimental peak-height statistics patel98b  and the predictions of the CI model for $`kT0.6\mathrm{\Delta }`$ ($`\mathrm{\Delta }`$ is the mean spacing between spin-degenerate levels).

The universal Hamiltonian of a quantum dot in the limit $`g_\mathrm{T}\mathrm{}`$ is given by kurland00 ; aleiner02

$$\widehat{H}=\underset{\lambda \sigma }{}ϵ_\lambda a_{\lambda \sigma }^{}a_{\lambda \sigma }+\frac{e^2}{2C}\widehat{n}^2J_s\widehat{𝐒}^2,$$
(1)
where $`ϵ_\lambda `$ are spin-degenerate single-particle levels ($`\sigma =\pm 1`$ labels the spin). The second term in Eq. (1), where $`C`$ is the dot’s capacitance and $`\widehat{n}`$ is the total–particle-number operator, accounts for the electrostatic energy of the dot. The third term, in which $`\widehat{𝐒}`$ is the total-spin operator, describes a constant exchange interaction with strength $`J_s`$. The occupation-number operator $`\widehat{n}_\lambda =\widehat{n}_{\lambda +}+\widehat{n}_\lambda `$ of any single-particle orbital $`\lambda `$ commutes with the total spin, $`[\widehat{n}_\lambda ,\widehat{𝐒}]=0`$, and the Hamiltonian $`\widehat{H}`$ is invariant under spin rotations. Thus the eigenstates of $`\widehat{H}`$ are characterized by their particle number $`N`$, the configuration of orbital occupation numbers $`𝐧=\{n_\lambda \}`$ ($`n_\lambda =0,1`$ or $`2`$), the total spin $`S`$, and the spin projection $`S_z=M`$. We label the eigenstates as $`|N𝐧\gamma SM`$ where the quantum number $`\gamma `$ distinguishes between states with the same total spin $`S`$ and particle configuration $`𝐧`$. The eigenenergies are given by $`\epsilon _{𝐧S}^{(N)}=_\lambda ϵ_\lambda n_\lambda +e^2N^2/2CJ_sS(S+1)`$.

In the limit of sequential tunneling (when a typical tunneling width is small compared with $`kT`$ and $`\mathrm{\Delta }`$), the conductance can be calculated using a rate-equations approach. In Ref. alhassid02a , we developed such an approach in the presence of interactions and spin. In particular, an explicit solution exists when the orbital occupation numbers $`n_\lambda `$ are good quantum numbers. Expressing the conductance $`G`$ in a rescaled form $`G=(e^2\overline{\mathrm{\Gamma }}/8\mathrm{}kT)g`$ (where $`\overline{\mathrm{\Gamma }}`$ is an average width of a level), we have, in the vicinity of the $`N+1`$-st Coulomb-blockade peak

$$g=4\underset{\genfrac{}{}{0pt}{}{\lambda 𝐧\gamma S}{𝐧^{}\gamma ^{}S^{}}}{}\stackrel{~}{P}_{𝐧S}^{(N)}f(\epsilon _{S^{}S}^\lambda )|([N+1]𝐧^{}\gamma ^{}S^{}a_\lambda ^{}N𝐧\gamma S)|^2g_\lambda .$$
(2)
Here $`g_\lambda =2\overline{\mathrm{\Gamma }}^1\mathrm{\Gamma }_\lambda ^\mathrm{l}\mathrm{\Gamma }_\lambda ^\mathrm{r}/(\mathrm{\Gamma }_\lambda ^\mathrm{l}+\mathrm{\Gamma }_\lambda ^\mathrm{r})`$ are the single-particle level conductances, where $`\mathrm{\Gamma }_\lambda ^{\mathrm{l},\mathrm{r}}`$ are the partial widths of an electron in orbital $`\lambda `$ to decay to the left or right lead. The equilibrium probability of the dot to be in the state $`|N𝐧\gamma SM`$ is $`\stackrel{~}{P}_{𝐧S}^{(N)}=e^{\beta (ϵ_{nS}^{(N)}\stackrel{~}{ϵ}_\mathrm{F}N)}/Z`$, where the partition function $`Z`$ is a Boltzmann-weighted sum over all possible $`N`$\- and $`(N+1)`$-body states (no other particle numbers contribute because of the charging energy), and $`\stackrel{~}{ϵ}_\mathrm{F}=e\zeta V_\mathrm{g}+ϵ_\mathrm{F}`$ is an effective Fermi energy ($`ϵ_\mathrm{F}`$ is the Fermi energy in the leads, $`V_\mathrm{g}`$ is the gate voltage and $`\zeta =C_\mathrm{g}/C`$ with $`C_\mathrm{g}`$ the dot-gate capacitance). The Fermi-Dirac function $`f(x)=(1+e^{\beta x})^1`$ is evaluated at an electron energy (relative to the Fermi energy) $`\epsilon _{S^{}S}^\lambda =\epsilon _{𝐧^{}S^{}}^{(N+1)}\epsilon _{𝐧S}^{(N)}\stackrel{~}{ϵ}_\mathrm{F}`$ that conserves energy at the transition between states $`|N𝐧\gamma SM`$ and $`|(N+1)𝐧^{}\gamma ^{}S^{}M^{}`$. The corresponding reduced matrix element $`([N+1]𝐧^{}\gamma ^{}S^{}a_\lambda ^{}N𝐧\gamma S)`$ enforces the selection rule $`S^{}=|S\pm 1/2|`$.

Eq. (2) can be rewritten in the form

$$g=\underset{\lambda }{}(w_\lambda ^{(0)}+w_\lambda ^{(1)})g_\lambda ,$$
(3)
where the contributions with $`n_\lambda =0`$ and $`n_\lambda =1`$ are collected in $`w_\lambda ^{(0)}`$ and $`w_\lambda ^{(1)}`$, respectively. For the cases with $`n_\lambda =0`$, the final $`(N+1)`$-particle state is given by $`|(N+1)𝐧^{}\gamma ^{}S^{}M^{}=_{M\sigma }(SM\frac{1}{2}\sigma |S^{}M^{})a_{\lambda \sigma }^{}|N𝐧\gamma SM`$, where $`(SM\frac{1}{2}\sigma |S^{}M^{})`$ is a Clebsch-Gordon coefficient. When $`n_\lambda =1`$ (and hence $`n_\lambda ^{}=2`$), the $`N`$-particle state can be similarly related to the $`(N+1)`$-particle state by changing to a hole representation. This leads to the following reduced matrix elements,

$$(𝐧^{}\gamma ^{}S^{}a_\lambda ^{}𝐧\gamma S)=()^{SS^{}\frac{1}{2}}\{\begin{array}{cc}\sqrt{2S^{}+1}\hfill & \mathrm{if}n_\lambda =0,\hfill \\ \sqrt{2S+1}\hfill & \mathrm{if}n_\lambda ^{}=2.\hfill \end{array}$$
(4)
Using the relation $`\stackrel{~}{P}_{𝐧S}^{(N)}f(ϵ_{S^{}S}^\lambda )=\stackrel{~}{P}_{𝐧^{}S^{}}^{(N+1)}[1f(ϵ_{S^{}S}^\lambda )]`$, we get

$$w_\lambda ^{(0)}=4\underset{S}{}b_{\lambda ,N,S}P_{N,S}\underset{S^{}=S\pm 1/2}{}(2S^{}+1)f(ϵ_{S^{}S}^\lambda ),$$
(5a)
$$w_\lambda ^{(1)}=4\underset{S^{}}{}c_{\lambda ,N+1,S^{}}P_{N+1,S^{}}\underset{S=S^{}\pm 1/2}{}(2S+1)[1f(ϵ_{S^{}S}^\lambda )],$$
(5b)
where the quantities $`b_{\lambda ,N,S}=\frac{1}{2}(\widehat{n}_\lambda 1)(\widehat{n}_\lambda 2)_{N,S}`$ and $`c_{\lambda ,N,S}=\frac{1}{2}\widehat{n}_\lambda (\widehat{n}_\lambda 1)_{N,S}`$ ensure that the sum is only over contributions with $`n_\lambda =0`$ or 1, respectively. They are defined in terms of thermal expectation values at constant particle number $`N`$ and spin $`S`$, i.e. $`\widehat{X}_{N,S}=\mathrm{Tr}_{N,S}[\widehat{X}e^{\beta \widehat{H}}]/\mathrm{Tr}_{N,S}[e^{\beta \widehat{H}}]`$. The quantity $`P_{N,S}`$ is the probability to find the dot with $`N`$ electrons and spin $`S`$,

$$P_{N,S}=e^{\beta [F_{N,S}+U_{N,S}]}/Z,$$
(6)
where $`F_{N,S}=\beta ^1\mathrm{ln}\mathrm{Tr}_{N,S}e^{\beta _{\lambda \sigma }ϵ_\lambda a_{\lambda \sigma }^{}a_{\lambda \sigma }}`$ is the free energy of $`N`$ non-interacting electrons with total spin $`S`$ and $`U_{N,S}=e^2N^2/2CJ_sS(S+1)\stackrel{~}{ϵ}_\mathrm{F}N`$.

The spin-projected trace of a scalar observable can be calculated from traces at fixed spin projection $`M`$ using $`\mathrm{Tr}_{N,S}\widehat{X}=\mathrm{Tr}_{N,M=S}\widehat{X}Tr_{N,M=S+1}\widehat{X}`$. For spin-$`1/2`$ particles, the projection on fixed particle number $`N`$ and spin projection $`M`$ is equivalent to projecting on a fixed number of spin-up and spin-down particles $`n_\pm =N/2\pm M`$. Therefore,

$$\mathrm{Tr}_{N,S}\widehat{X}=\mathrm{tr}_{\frac{N}{2}+S,\frac{N}{2}S}\widehat{X}\mathrm{tr}_{\frac{N}{2}+S+1,\frac{N}{2}(S+1)}\widehat{X},$$
(7)
where the traces “$`\mathrm{tr}`$” on the r.h.s. are evaluated at fixed $`n_+`$ and $`n_{}`$. Using $`\widehat{X}=e^{\beta _{\lambda \sigma }ϵ_\lambda a_{\lambda \sigma }^{}a_{\lambda \sigma }}`$ in Eq. (7) we find that the free energy in Eq. (6) is given by

$$e^{\beta F_{N,S}}=e^{\beta (\stackrel{~}{F}_{\frac{N}{2}+S}+\stackrel{~}{F}_{\frac{N}{2}S})}e^{\beta (\stackrel{~}{F}_{\frac{N}{2}+S+1}+\stackrel{~}{F}_{\frac{N}{2}(S+1)})}.$$
(8)
The free energy $`\stackrel{~}{F}_q`$ in Eq. (8) is defined for $`q`$ spinless particles $`e^{\beta \stackrel{~}{F}_q}=\stackrel{~}{\mathrm{t}}\mathrm{r}_qe^{\beta _\lambda ϵ_\lambda c_\lambda ^{}c_\lambda }`$ where $`c_\lambda ^{}`$ and $`c_\lambda `$ create and annihilate spinless particles in non-degenerate levels with energies $`ϵ_\lambda `$. The quantity $`c_{\lambda ,N,S}`$ from Eq. (5b) can now be expressed as

$$c_{\lambda ,N,S}=\frac{\stackrel{~}{n}_\lambda _{\frac{N}{2}+S}\stackrel{~}{n}_\lambda _{\frac{N}{2}S}e^{\beta (\stackrel{~}{F}_{\frac{N}{2}+S}+\stackrel{~}{F}_{\frac{N}{2}S})}\stackrel{~}{n}_\lambda _{\frac{N}{2}+S+1}\stackrel{~}{n}_\lambda _{\frac{N}{2}(S+1)}e^{\beta (\stackrel{~}{F}_{\frac{N}{2}+S+1}+\stackrel{~}{F}_{\frac{N}{2}(S+1)})}}{e^{\beta (\stackrel{~}{F}_{\frac{N}{2}+S}+\stackrel{~}{F}_{\frac{N}{2}S})}e^{\beta (\stackrel{~}{F}_{\frac{N}{2}+S+1}+\stackrel{~}{F}_{\frac{N}{2}(S+1)})}}$$
(9)
where $`\stackrel{~}{n}_\lambda `$ is the particle-number operator of a non-degenerate orbital $`\lambda `$. The function $`b_{\lambda ,N,S}`$ from Eq. (5a) is expressed by replacing $`\stackrel{~}{n}_\lambda `$ by $`(1\stackrel{~}{n}_\lambda )`$ in Eq. (9). The complete expression for the conductance is then obtained from Eqs. (3), (5), (6), (8), (9) and the relation indicated in the previous sentence. Thus the dot’s conductance in model (1) is determined in terms of the free energy $`\stackrel{~}{F}_q`$ and single-particle occupation numbers $`\stackrel{~}{n}_\lambda _q`$ of $`q`$ non-interacting spinless fermions. Both $`\stackrel{~}{F}_q`$ and $`\stackrel{~}{n}_\lambda _q`$ are familiar from earlier works in the framework of the CI model, and can be expressed in closed form using particle-number projection \[see Eqs. (140) in Ref. alhassid00 \].

In chaotic dots, the single-particle Hamiltonian in (1) is described by random-matrix theory. We have studied the statistics of peak heights and spacings for both the orthogonal and unitary symmetries. The dimension of the configuration sum in Eq. (2) increases combinatorially with the number of single-particle orbitals and a direct use of (2) becomes impractical at higher temperatures. In contrast, the closed expression we derived greatly facilitates the calculation of the conductance for a rather large model system of 50 single-particle orbitals $`\lambda `$. We checked that our results are not affected by the finite size of the system up to temperatures of $`kT3\mathrm{\Delta }`$.

Theoretical calculations of the width $`\sigma (\mathrm{\Delta }_2)`$ of the peak-spacing distribution, based on a spinless CI model alhassid99 , describe qualitatively the observed decrease of this quantity with increasing temperature patel98a . However, a proper modeling of the peak-spacing distribution itself requires the inclusion of spin. When spin is included and in the absence of an exchange interaction, the calculated values of $`\sigma (\mathrm{\Delta }_2)`$ (long-dashed line in Fig. 1) show a large discrepancy with the experimental values (symbols). Fig. 1 also shows $`\sigma (\mathrm{\Delta }_2)`$ for non-zero values of $`J_s`$. For a gas constant of $`r_s1.2`$ (that corresponds to the samples used in the experiments), the RPA estimate is $`J_s0.3\mathrm{\Delta }`$ oreg02 , and we find for this value a very good agreement with the measurements. The results for $`J_s=0.5\mathrm{\Delta }`$ underestimate the experimental widths. We remark that at temperatures $`kT0.4\mathrm{\Delta }`$, the model (1) does not describe well the shape of the peak-spacing distribution, and it is necessary to include the fluctuating part of the universal Hamiltonian to explain the absence of bimodality usaj01 ; alhassid02 . At higher temperatures, the bimodality is absent already in model (1) and the residual interaction has a negligible effect on the width.

Another measured quantity is the ratio between the standard deviation $`\sigma (g_{\mathrm{max}})`$ and the average $`\overline{g}_{\mathrm{max}}`$ of the peak heights $`g_{\mathrm{max}}`$ patel98b . The experimental data for this ratio (symbols in Fig. 2) are seen to be suppressed in comparison with the results of model (1) without an exchange term (long-dashed line in the left panel of Fig. 2). Spin-orbit interaction was proposed as a mechanism for this suppression at low temperatures held02 . It was necessary to assume a spin-orbit coupling that is sufficiently strong to completely decorrelate the spin-up and spin-down levels. However, spin-orbit effects are likely to be suppressed in the small dots used in the experiment. To determine whether an exchange interaction can explain the observed suppression of $`\sigma (g_{\mathrm{max}})/\overline{g}_{\mathrm{max}}`$, we calculated this ratio versus temperature $`kT`$ for different strengths of the exchange interaction (see Fig. 2). In the elastic limit (left panel), a realistic exchange interaction of $`J_s=0.3\mathrm{\Delta }`$ leads to closer agreement with the data. The remaining small discrepancy at temperatures $`kT0.6\mathrm{\Delta }`$ can probably be accounted for by adding a realistic weak spin-orbit interaction. It still remains to explain the discrepancy at higher temperatures, where inelastic scattering may play a role. The calculation of Ref. rupp02,  showed that the suppression of $`\sigma (g_{\mathrm{max}})/\overline{g}_{\mathrm{max}}`$ due to inelastic scattering is small for $`J_s=0`$. In the right panel of Fig. 2 we show results for the rapid-thermalization limit of strong inelastic scattering in the presence of an exchange interaction. While the agreement (for $`J_s=0.3\mathrm{\Delta }`$) is now better at low temperatures, we do not expect inelastic scattering to be important at these temperatures. At higher temperatures, the rapid-thermalization limit does not describe the data, and it would be interesting to determine the effect of an additional weak spin-orbit term.

For $`kT\mathrm{\Delta }`$ and $`J_s=0`$, the peak-height distribution $`P(g_{\mathrm{max}})`$ can be calculated analytically jalabert92  and is shown for the unitary symmetry as a solid line in the left panel of Fig. 3 (compared to the case of spinless electrons, the peak heights are rescaled alhassid00  by $`8(\sqrt{2}1)^21.37`$). Also shown (histogram) is the peak-height distribution calculated at $`kT=0.01\mathrm{\Delta }`$ and $`J_s=0.5\mathrm{\Delta }`$. No significant effect due to exchange is observed except for a small enhancement of the probability at small peak heights.

At finite temperature, the exchange interaction has a stronger effect on the peak-height distribution. The right panel of Fig. 3 compares the histogram (gray shaded) of the experimental data for $`P(g_{\mathrm{max}}/\overline{g}_{\mathrm{max}})`$ at $`kT=0.1\mathrm{\Delta }`$ with the calculated histograms for the cases of no exchange ($`J_s=0`$) and $`J_s=0.3\mathrm{\Delta }`$. This latter realistic value of the exchange interaction explains the observed suppression of the probability at small peak heights.

The weak-localization effect in the average peak height attracted recent attention both in experiment and theory folk01 ; beenakker01 ; eisenberg02 ; rupp02 . Its suppression at higher temperatures was suggested as a signature of inelastic scattering in the dot. The effect is quantified by the parameter $`\alpha =1(\overline{g}_{\mathrm{max}}^{\mathrm{GOE}}/\overline{g}_{\mathrm{max}}^{\mathrm{GUE}})`$. In the rapid-thermalization limit, $`\alpha `$ decreases rapidly with increasing temperature from its value of $`0.25`$ at $`kT\mathrm{\Delta }`$ beenakker01 . In contrast, if inelastic scattering is negligible, $`\alpha `$ was expected to be temperature independent. However, calculations for $`J_s=0`$ showed a slight suppression of the elastic $`\alpha `$ around $`kT0.25\mathrm{\Delta }`$ rupp02 ; eisenberg02 . This was understood by the fact that close lying levels and hence higher conductances are more likely for the orthogonal symmetry. The effect of the exchange interaction on $`\alpha `$ is shown in Fig. 4. We find that the dip in $`\alpha `$ around $`kT0.25\mathrm{\Delta }`$ is flattened out and in the small-temperature limit $`\alpha `$ becomes larger than $`0.25`$. While $`\alpha `$ is seen to be sensitive to the exchange interaction at low temperatures, the experimental uncertainties of Ref. folk01  are too large to observe this effect. In the rapid-thermalization limit, $`\alpha `$ is insensitive to $`J_s`$.

In conclusion, we have derived a closed expression for the conductance in the presence of spin and exchange interaction. Using this formula we studied the dependence of the peak height and spacing statistics on the exchange interaction and found a significantly better quantitative agreement with experiment than in the absence of an exchange interaction.

This work was supported in part by the U.S. DOE grant No. DE-FG-0291-ER-40608. We thank C. M. Marcus for helpful discussions.


text:
# First Results from KamLAND: Evidence for Reactor Anti-Neutrino Disappearance

(December 9, 2002)

## Abstract

KamLAND has been used to measure the flux of $`\overline{\nu }_e`$’s from distant nuclear reactors. In an exposure of 162 ton$``$yr (145.1 days) the ratio of the number of observed inverse $`\beta `$-decay events to the expected number of events without disappearance is $`0.611\pm 0.085(\mathrm{stat})\pm 0.041(\mathrm{syst})`$ for $`\overline{\nu }_e`$ energies $`>`$ 3.4 MeV. The deficit of events is inconsistent with the expected rate for standard $`\overline{\nu }_e`$ propagation at the 99.95% confidence level. In the context of two-flavor neutrino oscillations with CPT invariance, these results exclude all oscillation solutions but the ‘Large Mixing Angle’ solution to the solar neutrino problem using reactor $`\overline{\nu }_e`$ sources.

preprint: Submitted to PRL, Revised Dec. 9

The primary goal of the Kamioka Liquid scintillator Anti-Neutrino Detector (KamLAND) experiment kamland  is a search for the oscillation of $`\overline{\nu }_e`$’s emitted from distant power reactors. The long baseline, typically 180 km, enables KamLAND to address the oscillation solution of the ‘solar neutrino problem’ using reactor anti-neutrinos under laboratory conditions. The inverse $`\beta `$-decay reaction, $`\overline{\nu }_e+pe^++n`$, is utilized to detect $`\overline{\nu }_e`$’s with energies above 1.8 MeV in liquid scintillator (LS) RMP . The detection of the $`e^+`$ and the 2.2 MeV $`\gamma `$-ray from neutron capture on a proton in delayed coincidence is a powerful tool for reducing background. This letter presents the first results from an analysis of 162 ton$``$yr of the reactor $`\overline{\nu }_e`$ data.

KamLAND is located at the site of the earlier Kamiokande Kamioka , with an average rock overburden of 2,700 m.w.e. resulting in 0.34 Hz of cosmic-ray muons in the detector volume. As shown in Fig. 1, the neutrino detector/target is 1 kton of ultra-pure LS contained in a 13-m-diameter spherical balloon made of 135-$`\mu m`$-thick transparent nylon/EVOH (Ethylene vinyl alcohol copolymer) composite film. The balloon is supported and constrained by a network of kevlar ropes. The LS is 80% dodecane, 20% pseudocumene (1,2,4-Trimethylbenzene), and 1.52 g/liter of PPO (2,5-Diphenyloxazole) as a fluor. A buffer of dodecane and isoparaffin oils between the balloon and an 18-m-diameter spherical stainless-steel containment vessel shields the LS from external radiation. During the filling procedure a water extraction and nitrogen bubbling method purify , optimized for KamLAND, was used to purify the LS and buffer oil; PPO prepurification and dust removal were especially important.

The specific gravity of the buffer oil is adjusted to be 0.04% lower than that of the LS. An array of 1,879 photomultiplier tubes (PMTs), mounted on the inner surface of the containment vessel, completes the inner detector (ID) system. This array includes 1,325 specially developed fast PMTs with 17-inch-diameter photocathodes and 554 older Kamiokande 20-inch PMTs Kamiokande . While the total photo-cathode coverage is 34%, only the 17-inch PMTs corresponding to 22% coverage are used for the analysis in this letter. A 3-mm-thick acrylic barrier at 16.6-m diameter helps prevent radon emanating from PMT glass from entering the LS. The containment vessel is surrounded by a 3.2 kton water-Cherenkov detector with 225 20-inch PMTs. This outer detector (OD) absorbs $`\gamma `$-rays and neutrons from the surrounding rock and provides a tag for cosmic-ray muons. The primary ID trigger threshold is set at 200 PMT hits, corresponding to about 0.7 MeV. This threshold is lowered to 120 hits for 1 msec after the primary trigger to observe lower energy delayed activity. The OD trigger threshold is set to provide $`>99\%`$ tagging efficiency.

The inner detector is calibrated with $`\gamma `$-ray sources of 68Ge, 65Zn, 60Co, and Am-Be, deployed at various positions along the vertical axis. These sources provide calibration energies in the 0.5 to 7.6 MeV region. Energy measurements are based on the number of detected photoelectrons (p.e.) with corrections for PMT gain variation, solid angle, the density of PMTs, shadowing by suspension ropes, and light transparencies of the LS and buffer oil. At the center, about 300 p.e. per MeV are observed. Fig. 2 (a) shows the fractional deviation of the reconstructed energies from the source energies. The 68Ge and 60Co sources emit two coincident $`\gamma `$-rays and are plotted at an average energy in Fig. 2 (a). The observed energy resolution is $``$7.5%/$`\sqrt{E(\mathrm{MeV})}`$.

The energy calibration is augmented with studies of 40K and 208Tl (which are contaminants in the detector), Bi-Po sequential decays, the spallation products 12B and 12N, and $`\gamma `$’s from thermal neutron captures on protons and 12C. The reconstructed energy varies by less than 0.5% within the 10-m-diameter fiducial volume; local variations near the chimney region are about 1.6%. The energy scale exhibits less than 0.6% variation in time during the entire data run. Corrections for quenching and Cherenkov light production are included, and contribute to the systematic error in Fig. 2 (a). The estimated systematic error for the energy scale is 1.9% at the 2.6 MeV energy threshold which corresponds to a 2.1% uncertainty in the detected neutrino rate. Correlated decays from Bi-Po and 8He/9Li yield energies between 0.6 MeV and $``$15 MeV and are used to extend the range of the energy calibration.

The source positions are reconstructed from the relative times of PMT hits. Energy-dependent radial adjustments are used to reproduce the known source positions to $``$5 cm; the typical position reconstruction resolution is 25 cm. Vertex reconstruction performance throughout the LS volume is verified by reproducing the uniform distribution of 2.2 MeV capture $`\gamma `$’s from spallation neutrons, as shown in Fig. 2 (b).

The data presented in this letter were obtained during a period from March 4 through October 6, 2002. In total we collected 370 million events in 145.1 days of live time, corresponding to an average trigger rate of $``$30 Hz. Events with fewer than 10,000 p.e. ($``$30 MeV) and no prompt OD veto signal are classified as ‘reactor $`\overline{\nu }_e`$ candidates’; more energetic events are treated as ‘muon candidates.’

The criteria for selection of $`\overline{\nu }_e`$ events are (1) a fiducial volume cut ($`R<5`$ m), (2) a time correlation cut (0.5 $`\mu `$sec $`<\mathrm{\Delta }T<`$ 660 $`\mu `$sec), (3) a vertex correlation cut ($`\mathrm{\Delta }R<1.6`$ m), (4) a delayed energy window cut (1.8 MeV $`<E_{delay}<`$ 2.6 MeV), and (5) a cut on the delayed vertex position requiring it to be more than 1.2 m from the central vertical axis (to eliminate background from thermometers deployed in this region to monitor LS temperature). The overall efficiency associated with criteria (2)-(5) combined with the effect of (1) on the delayedvertex is $`(78.3\pm 1.6)\%`$.

Positrons annihilate in the LS and the detected energy includes the kinetic energy plus twice the electron rest mass; thus $`e^+`$ from $`\overline{\nu }_e`$ events produce $`E_{prompt}=E_{\overline{\nu }_e}\overline{E}_n0.8`$ MeV, where $`\overline{E}_n`$ is the average neutron recoil energy. Anti-neutrinos emitted by 238U and 232Th decays in the Earth, ‘geo-neutrinos’ ($`\overline{\nu }_{geo}`$), contribute low-energy events with $`E_{prompt}<2.49`$ MeV. For example, model Ia of $`\overline{\nu }_{geo}`$ in geonue  predicts about 9 $`\overline{\nu }_{geo}`$ events in our data set. However the abundances of U and Th and their distributions in the Earth are not well known. To avoid ambiguities from $`\overline{\nu }_{geo}`$’s we employ (6) a prompt energy cut, $`E_{prompt}>2.6`$ MeV in the present analysis.

Low-energy $`\gamma `$-rays from 208Tl ($``$3 Hz) entering from outside the balloon are potential sources of background in the range up to 3 MeV and are strongly suppressed by the fiducial volume cut (1). The fiducial volume is estimated using the uniform distribution of spallation neutron capture events shown in Fig. 2 (b). The ratio of these events in the fiducial volume to those in the total volume agrees with the geometric fiducial fraction to within 4.06%. The same method is used for higher energy events, 12N, 12B $`\beta `$’s following muon spallation, and agrees within 3.5%. Accounting for uncertainty in the LS total mass of 2.1%, we estimate the total systematic uncertainty of the fiducial volume to be 4.6%. The density of the LS is 0.780 g/cm3 at $`11.5^{}`$C; the hydrogen-to-carbon ratio is computed from the LS components to be 1.969 and was verified by elemental analysis ($`\pm 2\%`$). The specific gravity of the LS is measured to 0.01% precision using a commercial density meter, and we assign an additional 0.1% systematic error due to the uncertainty in the LS temperature. The 408 ton fiducial mass thus contains $`3.46\times 10^{31}`$ free target protons.

The trigger efficiency was determined to be 99.98% using LED light sources. The combined efficiency of the electronics, data acquisition, and event reconstruction was carefully studied using time distributions of uncorrelated events from calibration $`\gamma `$ sources. We find that this combined efficiency is greater than $`99.98\%`$. The vertex fitter yields $`>99.9\%`$ efficiency within 2 m of known source positions. Using calibrated 60Co and 65Zn sources, the overall efficiency was verified to within the 3% uncertainty in source strengths. The detection efficiency of delayed events from the Am-Be source (4.4 MeV prompt $`\gamma `$ and 2.2 MeV delayed neutron capture $`\gamma `$ within 1.6 m) was verified with 1% uncertainty.

From studies of Bi-Po sequential decays, the equilibrium concentrations of 238U and 232Th in the LS are estimated to be $`(3.5\pm 0.5)\times 10^{18}`$ g/g and $`(5.2\pm 0.8)\times 10^{17}`$ g/g, respectively. The observed background energy spectrum constrains the 40K contamination to be less than $`2.7\times 10^{16}`$ g/g. The accidental background, obtained from the observed flat distribution in the delayed time window 0.020-20 sec, is $`0.0086\pm 0.0005`$ events for the present data set.

At higher energies, the background is dominated by spallation products from energetic muons. We observe $``$3,000 neutron events/day/kton. We also expect $``$1,300 events/day/kton hagner  for various unstable products.

Single neutrons are efficiently suppressed by a 2-msec veto following a muon, but care is required to avoid neutrons which mimic the $`\overline{\nu }_e`$ delayed coincidence signal. Most fast neutrons are produced by energetic muons which pass through both the OD and the surrounding rock. This background is evaluated by detecting delayed coincidence events with a prompt signal associated with a muon detected only by the OD. As expected, a clear concentration of events near the balloon edge is observed. The expected background inside the fiducial volume is estimated by extrapolating the vertex position distribution and considering an OD reconstruction efficiency of 92%. To estimate the number of background events due to neutrons from the surrounding rock, the OD-associated rate is scaled by the relative neutron production and neutron shielding properties of the relevant materials. We estimate that the total fast neutron background is less than 0.5 events for the entire data set.

Most radioactive spallation products simply beta decay, and are effectively suppressed by requiring a delayed neutron signal. Delayed neutron emitters like 8He ($`T_{1/2}=119`$ msec) and 9Li (178 msec) are eliminated by two time/geometry cuts: (a) a 2-sec veto is applied for the entire fiducial volume following a showering muon (those with more than $`10^6`$ p.e., $``$3 GeV, extra energy deposition), (b) for the remaining muons, delayed events within 2 sec and 3 m from a muon track are rejected. The efficiency of these cuts is calculated from the observed correlation of spallation neutrons with muon tracks. The remaining 8He and 9Li background is estimated to be $`0.94\pm 0.85`$. The dead time due to the spallation cuts is 11.4%. This method is checked by exploiting the time distribution of the events after a detected muon to separate the short-lived spallation-produced activities from $`\overline{\nu }_e`$ candidates. The uncorrelated $`\overline{\nu }_e`$ event distribution has a characteristic time constant of $`1/R_\mu 3`$ sec, where $`R_\mu `$ is the incident muon rate. Spallation products have a much shorter time constant ($``$0.2 sec). These methods agree to 3% accuracy. As shown in Table 1 the total number of expected background events is $`0.95\pm 0.99`$, where the fast neutron contribution is included in the error estimate.

Instantaneous thermal power generation, burn-up and fuel exchange records for all Japanese commercial power reactors are provided by the power companies. The time dependence of the thermal power generation data is checked by comparison with the independent records of electric power generation. The fission rate for each fissile element is calculated from these data, resulting in a systematic uncertainty in the $`\overline{\nu }_e`$ flux of less than 1%. Averaged over the present live-time period, the relative fission yields from the various fuel components are 235U : 238U : 239Pu : 241Pu = 0.568 : 0.078 : 0.297 : 0.057. The $`\overline{\nu }_e`$ spectrum per fission and its error (2.48%) are taken from the literature spectra . These neutrino spectra have been tested to a few percent accuracy in previous short-baseline reactor $`\overline{\nu }_e`$ experiments RMP ; BUGEY . The finite $`\beta `$-decay lifetimes of fission products introduce an additional uncertainty of 0.28% to the $`\overline{\nu }_e`$ flux; this is estimated from the difference of the total $`\overline{\nu }_e`$ yield associated with shifting the run time by one day. The contribution to the $`\overline{\nu }_e`$ flux from Korean reactors is estimated to be (2.46 $`\pm `$ 0.25)% from the reported electric power generation rates. Other reactors around the world give an average (0.70 $`\pm `$ 0.35)% contribution, which is estimated by using reactor specifications from the International Nuclear Safety Center INSC . The uncertainties for the event rate calculation are summarized in Table 2. The errors from reactors outside Japan are included in the table under ‘Reactor Power’.

Although the anti-neutrino flux at the location of KamLAND is due to many nuclear reactors at a variety of distances, the $`\overline{\nu }_e`$ flux is actually dominated by a few reactors at an average distance of $``$180 km. More than 79% of the computed flux arises from 26 reactors within the distance range $`138`$-$`214`$ km. One reactor at 88 km contributes an additional 6.7% to the flux and the other reactors are more than 295 km away. This relatively narrow band of distances implies that for some oscillation parameters KamLAND can observe a distortion of the $`\overline{\nu }_e`$ energy spectrum.

The flux of anti-neutrinos from a reactor a distance $`L`$ from KamLAND is approximately proportional to the thermal power flux $`P_{th}/4\pi L^2`$, where $`P_{th}`$ is the reactor thermal power. The integrated total thermal power flux during the measurement live time is $`254\mathrm{Joule}/\mathrm{cm}^2`$. The systematic error assigned to the thermal power is conservatively taken as 2% from the regulatory specification for safe reactor operation. The corresponding expected number of reactor neutrino events (in the absence of neutrino oscillations) in the fiducial volume for this data set is $`86.8\pm 5.6`$.

The distribution of prompt and delayed energies for the present sample before the energy cuts are applied is shown in Fig. 3. A clear cluster of events from the 2.2 MeV capture $`\gamma `$’s is observed. One event with delayed energy around 5 MeV is consistent with a thermal neutron capture $`\gamma `$ on 12C. The space-time correlation of the prompt and delayed events is in good agreement with expectations, and the observed mean neutron capture time is $`188\pm 23\mu `$sec. After applying the prompt and delayed energy cuts, 54 events remain as the final sample. The ratio of the number of observed reactor $`\overline{\nu }_e`$ events to that expected in the absence of neutrino oscillations is

$$\frac{N_{obs}N_{BG}}{N_{expected}}=0.611\pm 0.085(\mathrm{stat})\pm 0.041(\mathrm{syst}).$$

The probability that the KamLAND result is consistent with the no disappearance hypothesis is less than 0.05%. Fig. 4 shows the ratio of measured to expected flux for KamLAND as well as previous reactor experiments as a function of the average distance from the source.

The observed prompt energy spectrum is shown in Fig. 5. The expected positron spectrum with no oscillations and the best fit with two-flavor neutrino oscillations above the 2.6 MeV threshold are shown. A clear deficit of events is observed. The measured spectrum is consistent (93% confidence) with a distorted spectrum shape as expected from neutrino oscillations, but a renormalized no-oscillation shape is also consistent at 53% confidence.

The neutrino oscillation parameter region for two-neutrino mixing is shown in Fig. 6. The dark shaded area is the LMA region at 95% C.L. derived from solar . The shaded region outside the solid line is excluded at 95% C.L. from the rate analysis with $`\mathrm{\Delta }\chi ^23.84`$ and

$$\chi ^2=\frac{(0.611R(\mathrm{sin}^22\theta ,\mathrm{\Delta }m^2))^2}{0.085^2+0.041^2}.$$

Here, $`R(\mathrm{sin}^22\theta ,\mathrm{\Delta }m^2)`$ is the expected ratio with the oscillation parameters.

The spectrum of the final event sample is then analyzed with a maximum likelihood method to obtain the optimum set of oscillation parameters with the following $`\chi ^2`$ definition:

$`\chi ^2`$ $`=`$ $`\chi _{\mathrm{rate}}^2(\mathrm{sin}^22\theta ,\mathrm{\Delta }m^2,N_{\mathrm{BG1}2},\alpha _{14})`$
$`2\mathrm{log}L_{\mathrm{shape}}(\mathrm{sin}^22\theta ,\mathrm{\Delta }m^2,N_{\mathrm{BG1}2},\alpha _{14})`$
$`+\chi _{\mathrm{BG}}^2(N_{\mathrm{BG1}2})+\chi _{\mathrm{distortion}}^2(\alpha _{14}),`$
where $`L_{\mathrm{shape}}`$ is the likelihood function of the spectrum including deformations from various parameters. $`N_{\mathrm{BG1}2}`$ are the estimated number of 9Li and 8He backgrounds and $`\alpha _{14}`$ are the parameters for the shape deformation coming from energy scale, resolution, $`\overline{\nu }_e`$ spectrum and fiducial volume. These parameters are varied to minimize the $`\chi ^2`$ at each pair of $`[\mathrm{\Delta }m^2,\mathrm{sin}^2\theta ]`$ with a bound from $`\chi _{\mathrm{BG}}^2(N_{\mathrm{BG1}2})`$ and $`\chi _{\mathrm{distortion}}^2(\alpha _{14})`$. The best fit to the KamLAND data in the physical region yields $`\mathrm{sin}^22\theta =1.0`$ and $`\mathrm{\Delta }m^2=6.9\times 10^5\mathrm{eV}^2`$ while the global minimum occurs slightly outside of the physical region at $`\mathrm{sin}^22\theta =1.01`$ with the same $`\mathrm{\Delta }m^2`$. These numbers can be compared to the best fit LMA values of $`\mathrm{sin}^22\theta =0.833`$ and $`\mathrm{\Delta }m^2=5.5\times 10^5\mathrm{eV}^2`$ from solar . The 95% C.L. allowed regions from the spectrum shape analysis are shown in Fig. 6. The allowed regions displayed for KamLAND correspond to $`0<\theta <\frac{\pi }{4}`$ consistent with the solar LMA solution, but for KamLAND the allowed regions in $`\frac{\pi }{4}<\theta <\frac{\pi }{2}`$ are identical hitoshi .

Another spectral shape analysis is performed with a lower prompt energy threshold of 0.9 MeV in order to check the stability of the above result and study the sensitivity to $`\overline{\nu }_{geo}`$. With this threshold, the total background is estimated to be $`2.91\pm 1.12`$ events, most of which come from accidental and spallation events. The systematic error is 6.0%, which is smaller than that for the final event sample due to the absence of an energy threshold effect. When the maximum likelihood is calculated, the $`\overline{\nu }_{geo}`$ fluxes from 238U and 232Th are treated as free parameters. The best fit in this analysis yields $`\mathrm{sin}^22\theta =0.91`$ and $`\mathrm{\Delta }m^2=6.9\times 10^5\mathrm{eV}^2`$. These results and the allowed region of the oscillation parameters are in good agreement with the results obtained above. The numbers of $`\overline{\nu }_{geo}`$ events for the best fit are 4 for 238U and 5 for 232Th, which corresponds to $``$40 TW radiogenic heat generation according to model Ia in geonue . However, for the same model, $`\overline{\nu }_{geo}`$ production powers from 0 to 110 TW are still allowed at 95% C.L. with the same oscillation parameters.

If three neutrino generations are considered, the $`\overline{\nu }_e`$ survival probability depends on two mixing angles $`\theta _{12}`$ and $`\theta _{13}`$. In the region close to the best fit KamLAND solution the survival probability is, to a very good approximation, given by

$$P(\overline{\nu }_e\overline{\nu }_e)\mathrm{cos}^4\theta _{13}\left[1\mathrm{sin}^22\theta _{12}\mathrm{sin}^2\frac{\mathrm{\Delta }m_{12}^2L}{4E_\nu }\right],$$

with $`\mathrm{\Delta }m_{12}\mathrm{\Delta }m^2`$ from the 2-flavor analysis above. The CHOOZ experiment CHOOZ  established an upper limit of $`\mathrm{sin}^22\theta _{13}<0.15`$, or $`\mathrm{cos}^4\theta _{13}0.92`$. The best fit KamLAND result would correspond approximately to $`0.86<\mathrm{sin}^22\theta _{12}<1.0`$.

In summary KamLAND has used measurements at large distances ($``$180 km) to demonstrate, for the first time, reactor $`\overline{\nu }_e`$ disappearance at a high confidence level ($`99.95\%`$). Since one expects a negligible reduction of $`\overline{\nu }_e`$ flux from the SMA, LOW and VAC solar neutrino solutions, the LMA region is the only remaining oscillation solution consistent with the KamLAND result and CPT invariance. The allowed LMA region is further reduced by these results. Future measurements with greater statistical precision and reduced systematic errors will enable KamLAND to provide a high-precision measurement of the neutrino oscillation parameters.

The KamLAND experiment is supported by the COE program of Japanese Ministry of Education, Culture, Sports, Science and Technology, and funding from the United States Department of Energy. The reactor data are provided by courtesy of the following electric associations in Japan; Hokkaido, Tohoku, Tokyo, Hokuriku, Chubu, Kansai, Chugoku, Shikoku and Kyushu Electric Power Companies, Japan Atomic Power Co. and Japan Nuclear Cycle Development Institute. Kamioka Mining and Smelting Company has provided service for activities in the mine.


text:
# The Amplitude Modulation transform

## Introduction

This note has been motivated by the following question:

> Let $`0=x_0<x_1<\mathrm{}<x_n<\mathrm{}`$ be the sequence of local maxima of the *sinc* function $`\mathrm{sinc}(x)=\mathrm{sin}(x)/x.`$ Is the sequence $`1=\mathrm{sinc}(x_0),\mathrm{sinc}(x_1),\mathrm{},\mathrm{sinc}(x_n),\mathrm{}`$ decreasing?

This question is not difficult to answer. Indeed, at a critical point $`x_i,`$

$$\mathrm{sinc}^{}(x_i)=\frac{x_i\mathrm{cos}(x_i)\mathrm{sin}(x_i)}{x_i^2}=0,$$

which implies that

(0.1)
$$\mathrm{cos}(x_i)=\frac{\mathrm{sin}(x_i)}{x_i}=\mathrm{sinc}(x_i).$$

One can also write the above equation as:

(0.2)
$$x_i=\mathrm{tan}(x_i),$$

or, equivalently:

(0.3)
$$\mathrm{arctan}(x_i)=x_i.$$

Combining equations (0.1) and (0.3), we obtain:

(0.4)
$$\mathrm{sinc}(x_i)=\mathrm{cos}(\mathrm{arctan}(x_i))=\frac{1}{\sqrt{1+x_i^2}},$$

so the decrease of $`\mathrm{sinc}(x_i)`$ is immediate.

## 1. The *Amplitude Modulation* transform

The formula (0.4) suggest the following:

###### Definition 1.1.

The *Amplitude Modulation* transform $`𝒜(f)`$ of a function $`f:`$ is the set of functions whose values at the critical points of $`f(x)\mathrm{sin}(x)`$ agrees with those of $`f(x)\mathrm{sin}(x).`$

###### Remark 1.2.

In fact, if in the definition of the $`𝒜`$ transform we replace the multiplier $`\mathrm{sin}(x)`$ by $`\mathrm{sin}(x+k),`$ we obtain the same function $`𝒜f.`$ This is an observation of W. D. Smith, and it allows us to replace the definition above with the more pleasant definition below:

###### Definition 1.3.

The *Amplitude Modulation* transform $`𝒜(f)`$ of a function $`f:`$ is the function whose values at the critical points of $`f(x)\mathrm{sin}(x+k)`$ agrees with those of $`f(x)\mathrm{sin}(x+k)`$ for *all* values of the phase parameter $`k.`$

The discussion in the Introduction can thus be summarized as follows:

###### Theorem 1.4.

The function $`\frac{1}{\sqrt{1+x^2}}`$ is the $`𝒜`$ transform of $`\frac{1}{x}.`$

To get an analogous result for a general function $`f(x),`$ we perform the same sort of computation as in the Introduction:

(We will use the notation $`f_x`$ for the derivative of $`f`$ for typographical reasons.) The critical points of $`f(x)\mathrm{sin}(x)`$ are the points where:

$$\frac{df(x)\mathrm{sin}(x)}{dx}=0.$$

Expanding, we see that $`f(x)\mathrm{cos}(x)+f_x(x)\mathrm{sin}(x)=0,`$ and so

$$\mathrm{cot}(x)=\frac{f_x(x)}{f(x)},$$

so that

(1.1)
$$x=\mathrm{arccot}\left(\frac{f_x(x)}{f(x)}\right),$$

while

(1.2)
$$f(x)\mathrm{sin}(x)=\mathrm{tan}(x)\mathrm{sin}(x)f_x(x).$$

Combining Eq. (1.1) and Eq. (1.2) we see that at the critical points:

(1.3)
$$f(x)\mathrm{sin}(x)=\pm \frac{f^2(x)}{\sqrt{f^2(x)+f_x^2(x)}},$$

which we can summarize in

###### Theorem 1.5 (Theorem-Definition).

The function $`𝒜(f)`$ is defined by

$$𝒜(f)(x)=\frac{f^2(x)}{\sqrt{f^2(x)+f_x^2(x)}},$$

Here are some examples: As we have seen before, if $`f(x)=1/x,`$ then

$$𝒜(f)(x)=\frac{1}{\sqrt{1+x^2}}.$$

If $`f(x)=x^\alpha ,`$ then

$$𝒜(f)(x)=\frac{x^{\alpha +1}}{\sqrt{x^2+\alpha ^2}}.$$

If $`f(x)=\mathrm{exp}(x),`$ then

$$𝒜(f)(x)=\frac{\mathrm{exp}(x)}{\sqrt{2}},$$

while if $`f(x)=\mathrm{exp}(x),`$ then

$$𝒜(f)(x)=\frac{\mathrm{exp}(x)}{\sqrt{2}}.$$

If $`f(x)=\mathrm{exp}(g(x)),`$ then

(1.4)
$$𝒜(f)(x)=f(x)\frac{1}{\sqrt{1+g_x^2}}.$$

## 2. Some algebraic observations

We will need to recall a definition:

###### Definition 2.1.

A function $`y=f(x)`$ is called *algebraic* if there exists a two-variable polynomial $`P,`$ such that $`P(x,y)=0.`$

And some well-known results:

###### Theorem 2.2.

If $`z`$ is an algebraic function of $`y`$ and $`y`$ is an algebraic function of $`x,`$ then $`z`$ is an algebraic function of $`x.`$ Further, if $`y_1`$ and $`y_2`$ are algebraic functions of $`x,`$ then so are $`y_1y_2`$ and $`y_1+y_2.`$ Finally, if $`y`$ is an algebraic function, then so is $`y^{}.`$

###### Proof Sketch and lightning introduction to elimination theory.

The proofs of all the assertions follow from the following basic fact: two univariate polynomials $`P_1`$ and $`P_2`$ over a domain $`R`$ with unity have no common zeros if and only if their greatest common divisor is $`1,`$ or, equivalently, there exist polynomials $`Q_1`$ and $`Q_2,`$ such that

(2.1)
$$Q_1P_1+Q_2P_2=1.$$

Since Eq. (2.1) is a system of linear equations for the coefficients of $`Q_1`$ and $`Q_2,`$ the existence of $`Q_1`$ and $`Q_2`$ as above is easily seen to be equivalent to the non-vanishing of the determinant of the linear system. This determinant is the so-called *resultant* of the polynomials $`P_1`$ and $`P_2.`$ Now, if we have two polynomial equations $`P(x,y)=0`$ and $`Q(x,y)=0,`$ they can regarded as two polynomials in $`y`$ whose coefficients are polynomials in $`x,`$ and so the set of $`x`$-coordinates of the points in the common zero-set of $`P`$ and $`Q`$ all have the property that at those points, the resultant of the two equations vanishes. The resultant is a polynomial in $`x,`$ and so $`y`$ has been eliminated from consideration, hence the name “elimination theory.” Now, to proceed with the proof of the Theorem 2.2: If $`z`$ is an algebraic function of $`y`$ and $`y`$ is an algebraic function of $`x,`$ then there are equations $`P(x,y)=0`$ and $`Q(y,z)=0.`$ Eliminating $`y`$ from the two equations, we see that $`z`$ is algebraic. If $`y_1`$ and $`y_2`$ are algebraic, and $`y=y_1y_2,`$ then we have the the three equations (the ones satisfied by $`y_1`$ and $`y_2`$ and $`y=y_1y_2.`$ We can eliminate first $`y_1`$ and then $`y_2,`$ to show that $`y`$ is algebraic, similarly with $`y_1+y_2.`$

Finally, to show that the derivative is algebraic, we differentiate $`P(x,y)=0`$ implicitly, to obtain $`Q(x,y,y^{})=0.`$ Eliminating $`y`$ from the two equations we obtain the algebricity of $`y^{}.`$

###### Remark 2.3.

For considerably more detail on the subject of elimination, please see .

As corollaries of the above Theorem, we see that

###### Corollary 2.4.

If $`f(x)`$ is an algebraic function, then so is $`𝒮(f).`$ If $`f,g`$ are algebraic functions, then if $`h(x)=f(x)\mathrm{exp}(g(x)),`$ it follows that $`𝒜(h)(x)=k(x)h(x),`$ where $`k`$ is algebraic.

## 3. Inverse problems

The first obvious inverse problem is the following:

> Which functions $`g(x)`$ are $`𝒜`$ transforms?

Construction of the inverse transform is equivalent to the solution of the ODE

(3.1)
$$f^{}(x)=\pm f(x)\sqrt{\frac{f^2(x)}{g^2(x)}1}.$$

The choice of plus or minus is already troubling, as is the fact that the right hand side is frequently not Lipschitz, so the usual Picard existence theorem for ODE does not apply everywhere, and uniqueness fails spectacularly: the functions $`1/\mathrm{sin}x`$ and $`1`$ have the same $`𝒜`$ transform. This example also demonstrates that the initial value problem can develop singularities in finite time. Nevertheless, some things can be said. First:

###### Lemma 3.1.

The transformed function $`𝒜(f)`$ has a critical point whenever $`f`$ has a critical point. Furthermore, at such a critical point $`x,`$ $`𝒜(f)(x)=f(x),`$ and the last equality only holds at a critical point of $`f.`$

###### Proof.

A simple computation. ∎

We can thus simplify our life by attempting to solve 3.1 on an interval $`[a,b]`$ where $`g`$ is monotone, and in addition, $`0<g(x)<M.`$ We can pick between the two equations:

(3.2) $`f^{}(x)`$ $`=`$ $`f(x)\sqrt{{\displaystyle \frac{f^2(x)}{g^2(x)}}1}`$
(3.3) $`f^{}(x)`$ $`=`$ $`f(x)\sqrt{{\displaystyle \frac{f^2(x)}{g^2(x)}}1}`$
Now, by the Lemma 3.1 we know that we have local existence and uniqueness of solutions, and so the only thing we need check is that singularities do not develop in finite time. To do this we analyze two separate cases:

* Case 1. $`g(x)`$ is *decreasing* on $`[a,b].`$ In this case we take Eq. (3.3). Local existence and uniqueness is assured by the Picard theorem (see \[1, Chapter 1\]). We pick the initial value $`f(a)`$ at will (as long as it is bigger than $`g(a).`$) Since $`f^{}(x)`$ is always negative we know that $`f(x)<f(a),`$ and since we know that $`g(x)>0`$ on $`[a,b]`$ we know that $`f(x)>g(x)>0,`$ so it follows that we have a solution on $`[a,b].`$
* Case 2. $`g(x)`$ is *increasing*. In this case we start at the right endpoint $`b,`$ and use Eq. (3.2), and then construct the solution going right to left. The reasoning in Case 1 goes through verbatim.

What happens if $`g(x)`$ has critical points on $`[a,b]`$? In that case, it is fairly obvious that we can construct a “weak inverse,” but anything more seems to require much more work. In case the reader is dissatisfied with the nonexplicit nature of our construction, (s)he will perhaps be mollified by the observation that the following problem can be solved explicitly:

> Given a function $`r(x)>1`$ on $`[a,b],`$ construct a positive $`f(x),`$ such that $`f(x)/𝒜(f)(x)=r(x).`$

Using Eq. (1.4), it is easy to see that

(3.4)
$$f(x)=\mathrm{exp}\left(\pm C_a^x𝑑t\sqrt{r^2(t)1}\right)$$

is the desired solution.

## 4. Questions

The most natural question is:

###### Question 4.1.

Given a smooth $`g(x),`$ is there a natural way to construct an $`f(x)`$ such that $`g(x)=𝒜(f)`$

Changing categories:

###### Question 4.2.

Is $`𝒜`$ invertible on the set of algebraic functions?

or

###### Question 4.3.

Suppose $`f`$ satisfies a first order linear differential equation with algebraic coefficients. Is the same always true of $`𝒜(f)\mathrm{?}`$