text:
# Introduction

## Introduction

One of the most ancient and difficult problems of physics is that of the origin of the arrow of time, i.e., of the nature of the difference between the past and the future. It is conventional to search for a solution to this problem in dynamics, i.e., time evolution of a state of a physical system. The solution may be given by a dynamics which is asymmetric, or orientable in the sense of the sequence of states.

A deterministic dynamics does not involve such an orientability. Therefore it seems advisable to search for the solution in an indeterministic dynamics. But, according to an established opinion, in the standard indeterministic quantum dynamics there is no arrow of time as well .

A crucial feature of an orientable dynamics is the prevalence of retrodiction over prediction. An essential mathematical concept inherent in indeterminism is that of randomness \[4-6\]. It is the latter that seems to be the main impediment to the construction of a dynamics with the above-mentioned feature.

The simplest way of constructing an orientable dynamics is to construct one in which a past-directed time evolution would have been deterministic. In doing this, the pivotal point is the condition that a future-directed time evolution represent the standard indeterministic, i.e., probabilistic quantum dynamics.

In the present paper, that idea is realized. We introduce a jump-deciding mechanism for a two-state tangency vertex. The mechanism consists of the Planck clock (one with the Planck time period) and rules which decide a quantum jump by the reading of the clock. On the basis of this mechanism, a deterministic past-directed dynamics is constructed. The corresponding future-directed dynamics turns out to be indeterministic. The rules are selected in such a way that probabilities be standard quantum ones.

The construction outlined above exhibits a retrodictive universe. Its salient features are the following. A complete retrodiction and a partial prediction are involved. It is impossible to introduce initial conditions and to extend dynamics forward in time; thus the retrodictive universe is constructed at once as a whole—for all times from a maximal future to a maximal past. In this connection we quote Weyl : “The objective world simply is, it does not happen. Only to the gaze of my consciousness, crawling upward along the life line of my body, does a section of this world come to life as a fleeting image in space which continuously changes in time.”

## 1 The problems of the arrow of time and randomness

### 1.1 The problem of the time arrow

We quote Roger Newton : “One of the greatest puzzles of physics is the manifest discord between two facts: on one hand, all the fundamental equations and laws of physics are (essentially) invariant under time reversal;…on the other hand, we are all aware that at the macroscopic level, many physical processes flow in one time direction only. The unidirectional flow of time …is one of the most obvious features both of our consciousness and of the physical world, …which near the end of the nineteenth century presented physics with one of the most profound challenges, and about which there are, to this day, strong disagreements among physicists.”

There are five arrows of time :

1. The delay between cause and effect.

2. The biological, or cognitive arrow.

3. The second law of thermodynamics.

4. The cosmological arrow (the expansion of the universe).

5. The direction of the time parameter used in physics.

The fifth arrow is present in any dynamics. The second and third arrows are related by an informational aspect: information on the past is greater than on the future, so that (informational) entropy increases in time. It is this aspect that forms the basis for our attacking the time arrow problem.

### 1.2 The problem of randomness

The mathematical concept of randomness is inherent in an indeterministic, i.e., probabilistic dynamics. A physical problem related to the concept is the impossibility of empirically verifying or falsifying the randomness of results of a generic quantum dynamical process. We quote Beltrami : “…there can be no formal proof that a sufficiently long string is random …In response to the persistent question ‘Is it random?’ the answer must now be ‘probably, but I’m not sure’; ‘probably’ because most numbers are, in fact, random …and ‘not sure’ because the randomness of long strings is essentially undecidable.” Thus the question of the randomness of quantum dynamics cannot be, strictly speaking, decided empirically. So there remains the possibility of choosing a decision which would provide a resolution of some theoretical problems.

### 1.3 Time nonorientability of the standard indeterministic quantum dynamics

There exists no arrow of time in the standard indeterministic quantum dynamics. This is exhibited by means of a fully time-symmetric construction from which conventional quantum mechanics may be derived  (see a detailed treatment in ).

Since the conventional indeterministic (probabilistic) dynamics is based on the concept of randomness, it seems reasonable to try to abandon the latter in constructing an orientable dynamics: we do not see what else may be done.

## 2 A deterministic past-directed quantum jump
dynamics

### 2.1 Tangency vertex

Let us consider a tangency vertex , i.e., a confluence or/and branch point of levels without crossing, in more detail. Assume that without quantum jumps the Hamiltonian $`H(t)`$ is a $`C^{\mathrm{}}`$ operator-valued function of time, so that for a vertex

$$\frac{d^nH_{\mathrm{ver}}^{}}{dt^n}=\frac{d^nH_{\mathrm{ver}}^+}{dt^n}=\frac{d^nH_{\mathrm{ver}}}{dt^n},n=0,1,2,\mathrm{}$$
(2.1.1)
holds. In an infinitesimal neighborhood of the vertex,

$$H(t_{\mathrm{ver}}+\mathrm{\Delta }t)=H[g(t_{\mathrm{ver}}+\mathrm{\Delta }t)]=H[g(t_{\mathrm{ver}})+\dot{g}(t_{\mathrm{ver}})\mathrm{\Delta }t]$$
(2.1.2)
and

$$H_{\mathrm{ver}}(t_{\mathrm{ver}}\mathrm{\Delta }t)=\underset{l}{\overset{1,n^{}}{}}\epsilon _l(t_{\mathrm{ver}}\mathrm{\Delta }t)P_l(t_{\mathrm{ver}}\mathrm{\Delta }t),\mathrm{\Delta }t0$$
(2.1.3)
For any projector $`E(t)`$

$$E\frac{dE}{dt}E=0$$
(2.1.4)
is fulfilled, so that

$$\frac{d^n}{dt^n}\left(E\frac{dE}{dt}E\right)=0,n=0,1,2,\mathrm{}$$
(2.1.5)
Hence it is easily seen that from

$$\frac{d^mE^{}}{dt^m}=\frac{d^mE^+}{dt^m}\mathrm{for}m=0,1,\mathrm{},n$$
(2.1.6)
follows

$$\left(E\frac{d^{n+1}E}{dt^{n+1}}E\right)^{}=\left(E\frac{d^{n+1}E}{dt^{n+1}}E\right)^+$$
(2.1.7)
Now for a tangency vertex, we obtain by the method of

$$\frac{d^nP_{\mathrm{tan}}^{}}{dt^n}=\frac{d^nP_{\mathrm{tan}}^+}{dt^n}=\frac{d^nP_{\mathrm{tan}}}{dt^n},n=0,1,2,\mathrm{}$$
(2.1.8)
where

$$P_{\mathrm{tan}}^{}=\underset{l}{}P_l^{}=P_{\mathrm{tan}}$$
(2.1.9)
and

$$\frac{d^n\epsilon _l^{}}{dt^n}=\frac{d^n\epsilon _l^{}^+}{dt^n}=\frac{d^n\epsilon _{\mathrm{tan}}}{dt^n},n=0,1,2,\mathrm{}$$
(2.1.10)
Thus, at a tangency vertex, there is contact of order $`n=\mathrm{}`$.

### 2.2 A two-state jump-deciding mechanism in a past-directed
dynamics

We assume that the degeneracy of levels occurs only at vertices, so that there are only two-state vertices. Therefore we consider a two-state tangency vertex. For the latter, we introduce a jump-deciding mechanism. The mechanism consists of a Planck clock and rules deciding a jump by the reading of the clock. The Planck clock is one with the Planck time $`t_\mathrm{P}`$ period. (There is no other natural time interval.) We will measure time in units of $`t_\mathrm{P}`$ so that the period of the clock is 1.

There are two states $`i=1,2`$ for $`t<t_{\mathrm{tan}}`$ and two states $`f=1,2`$ for $`t>t_{\mathrm{tan}}`$. We consider transitions, i.e., quantum jumps $`if`$ in a future-directed dynamics and transitions $`fi`$ in a past-directed dynamics. It is the latter that has to be deterministic, so that let $`f=a`$ be an actual state before a quantum jump $`ai`$. Quantum probabilities are

$$p(f|i)=p_{fi}=p_{if}=|f|i|^2$$
(2.2.1)
with

$$\underset{f}{}p(f|i)=\underset{i}{}p(f|i)=1$$
(2.2.2)
Let

$$ϵ_{i=1}<\epsilon _{i=2}$$
(2.2.3)
Introduce some $`\overline{t}`$,

$$0\overline{t}1$$
(2.2.4)
A deterministic past-directed dynamics is defined by the following jump-deciding rule for $`ai`$:

$$\mathrm{if}0t_{\mathrm{tan}}<\overline{t}\mathrm{then}a1,\mathrm{if}\overline{t}t_{\mathrm{tan}}<1\mathrm{then}a2$$
(2.2.5)
### 2.3 A future-directed dynamics

In order to provide a correct probabilistic future-directed dynamics, we introduce

$$p_{\mathrm{min}/\mathrm{max}}=(\mathrm{min}/\mathrm{max})_i\{p(a|i)\},p_{\mathrm{min}}1/2p_{\mathrm{max}}$$
(2.3.1)
$$t_1=\frac{p_{\mathrm{min}}}{1/2+p_{\mathrm{min}}},t_2=\frac{1/2}{1/2+p_{\mathrm{min}}},t_1t_2,t_1+t_2=1$$
(2.3.2)
and put

$$\overline{t}=t_i\mathrm{where}p(a|i)=p_{\mathrm{min}}$$
(2.3.3)
The rule (2.2.5) is equivalent to these:

$$\mathrm{if}0t_{\mathrm{tan}}<t_1\mathrm{then}a1,\mathrm{if}t_2t_{\mathrm{tan}}<1\mathrm{then}a2$$
(2.3.4)
$$\mathrm{if}t_1t_{\mathrm{tan}}<t_2\mathrm{then}ai\mathrm{where}p(a|i)=p_{\mathrm{max}}>1/2(\mathrm{here}t_1<t_2)$$
(2.3.5)
Now if an initial state is $`i`$, then, in view of (2.3.4),

$$t_{i1}t_{\mathrm{tan}}<t_{i+1}$$
(2.3.6)
where

$$0=t_0t_1t_2t_3=1$$
(2.3.7)
and

$$t_{i+1}t_{i1}=t_2$$
(2.3.8)
Therefore we should have for probabilities

$$p_{\mathrm{min}}=\frac{t_1/2}{t_2},p_{\mathrm{max}}=\frac{(t_2t_1)+t_1/2}{t_2}=\frac{t_2t_1/2}{t_2}$$
(2.3.9)
since $`t_1t_{\mathrm{tan}}<t_2`$ implies $`p(a|i)=p_{\mathrm{max}}`$, whereas if $`t_{\mathrm{tan}}<t_1`$ or $`t_{\mathrm{tan}}t_2`$ then $`p(a|i)=p_{\mathrm{min}}`$ and $`p(a|i)=p_{\mathrm{max}}`$ are equiprobable. The choice (2.3.3) provides (2.3.9).

### 2.4 The Planck clock and two fundamental theories

For the sake of generality, we might choose any time period $`t_{\mathrm{per}}`$ rather than the Planck time $`t_\mathrm{P}`$. We have put

$$t_{\mathrm{per}}=t_\mathrm{P}$$
(2.4.1)
for lack of any other natural time interval. There is an added reason for incorporating the Planck clock into quantum jump dynamics, i.e., putting (2.4.1).

The Planck time

$$t_\mathrm{P}=\left(\frac{\mathrm{}G}{c^5}\right)^{1/2}=5.3906\times 10^{44}\mathrm{s}$$
(2.4.2)
may be introduced as a natural unit for all physical quantities. Put

$$c=1$$
(2.4.3)
then

$$\mathrm{}G=t_\mathrm{P}^2$$
(2.4.4)
Next, there are two systems of units:

$$\mathrm{}=1,G=t_\mathrm{P}^2$$
(2.4.5)
(natural units) and

$$G=1,\mathrm{}=t_\mathrm{P}^2$$
(2.4.6)
(“geometrized units” ). In the natural units, general relativity involves $`t_\mathrm{P}`$: it appears in the Einstein equation, but the standard quantum theory does not involve $`t_\mathrm{P}`$: it does not appear in the Schr$`\ddot{\mathrm{o}}`$dinger equation and quantum probabilities. By contrast, in the geometrized units, it is quantum theory rather than general relativity that involves $`t_\mathrm{P}`$.

The Planck clock incorporates $`t_\mathrm{P}`$ into quantum theory via quantum jump dynamics. Now, in the natural units, general relativity and quantum theory enjoy equal rights with respect to $`t_\mathrm{P}`$, which links those fundamental theories in addition to the construction introduced in .

## 3 The retrodictive universe and its salient features

### 3.1 The retrodictive universe

The construction accomplished above exhibits a universe which is naturally called retrodictive. Let us consider its salient features.

### 3.2 A complete retrodiction and a partial prediction

Since the past-directed dynamics is deterministic, in the retrodictive universe there exists a complete retrodiction, to which the universe owes its name.

An informational aspect of retrodiction is this: A complete information on the universe’s states in the past reduces to and may be obtained from the information contained in the description of a present state and in the laws of the past-directed dynamics.

The future-directed dynamics is indeterministic and reproduces the dynamics of the standard probabilistic quantum theory. Notwithstanding this fact, there exists a partial prediction provided by (2.3.5): if

$$t_1t_{\mathrm{tan}}<t_2$$
(3.2.1)
then in the future-directed dynamics the jump $`if`$ happens to the $`f`$ for which

$$p(f|i)=p_{\mathrm{max}}>1/2$$
(3.2.2)
But this partial prediction does not violate the probabilistic relations (2.3.9).

### 3.3 The impossibility of introducing initial conditions

The past-directed dynamics constructed backwards in time starting from some state at $`t=t_{\mathrm{initial}}`$, $`\mathrm{}<t_{\mathrm{initial}}<\mathrm{}`$, cannot be extended forward in time: we may run into a situation where

$$i=1,t_2t_{\mathrm{tan}}<1\mathrm{or}i=2,0t_{\mathrm{tan}}<t_1$$
(3.3.1)
which is inconsistent with the rule (2.3.4). Thus it is impossible to construct a future-directed dynamics starting from some initial conditions. It is only a past-directed dynamics that may be constructed in the case of the retrodictive universe.

### 3.4 The entirety of the retrodictive universe

The impossibility of constructing a future-directed dynamics starting from some initial conditions given at some $`t_{\mathrm{initial}}`$, i.e., solving the Cauchy problem, implies the entirety of the retrodictive universe: the latter is determined at once for all times from a maximal future to a maximal past. The universe exists but does not evolve in time . A seeming evolution is the result of the indeterministic character of the future-directed dynamics.

This conclusion cracks the problem of initial conditions for the universe: there are none.

### 3.5 The arrow of time

The arrow of time is inherent in the retrodictive universe as well as a future-directed increase of informational entropy.

## Acknowledgments

I would like to thank Alex A. Lisyansky for support and Stefan V. Mashkevich for helpful discussions.


format: text

text:
# Multiple-photon resolving fiber-loop detector

## Abstract

We show first reconstructions of the photon-number distribution obtained with a multi-channel fiber-loop detector. Apart from analyzing the statistics of light pulses this device can serve as a sophisticated post-selection device for experiments in quantum optics and quantum information. We quantify its efficiency by means of the Fisher information and compare it to the efficiency of the ideal photodetector.

A major drawback of common detectors of weak light fields is their lack of photon-number resolution. Due to the nonlinear character of amplifying process the response of such detectors is not sensitive to the strength of the input signal. The only two detection events in such a case are “click” and “no click” that correspond to the presence or absence of the signal. A device capable of photon-number resolution would contribute both to fundamental research in quantum optics and to implementation of quantum communication and information protocols. While such devices were recently indeed constructed Kim et al. (1999); Kurtsiefer et al. (2000); Brattke et al. (2001); Santori et al. (2001); Pelton et al. (2002); Stevenson et al. (2002), they require operation under extreme conditions at present and therefore did not become a common laboratory tool yet. Also, their photon-number resolution is still limited to only few photons.

Another way of circumventing this problem is splitting the input pulse using a multiport device followed by an array of conventional binary detectors as was proposed in Paul et al. (1996); Kok and Braunstein (2001). In the ideal case of many output ports the input pulse gets perfectly split and each photon is detected separately. However reasonable performance of such device would require a very large number of beamsplitters and detectors which would results in a bulky and costly detection device. It has been suggested Banaszek and Walmsley (2003) to replace the complicated multiport device by a fiber loop and a single photodetector, see Fig. 1.

After each round-trip part of the incoming pulse gets transmitted to a conventional binary detector. This results in a time resolved series of detections, each of them corresponding to a different output port of the multiport device. Such a multiple photon resolving device has recently been built in our laboratory Haderka et al. . The variable ratio coupler inserted at the entrance to the fiber-loop delay line, see Fig. 1, is used to adjust the transmission probabilities of the output channels to a certain extent and thus tweak the overall resolving power of our instrument.

The purpose of this communication is twofold. First, we will analyze the performance of the fiber-loop detector and compare it to the efficiency of the ideal photodetection device. Second, we will show how to reconstruct the photon statistics of the input pulse from the data measured at the fiber-loop detector via the maximum-likelihood principle and apply this technique to experimental data.

The fiber-loop detector is used to count photons contained in the input pulse in an indirect way. Therefore it is natural to relate its performance to that of the “textbook photon counter”—the ideal photodetector of quantum efficiency $`\overline{\eta }`$. The concept of the equivalent efficiency is extremely useful in this context. The equivalent efficiency of the fiber-loop detector is defined to be equal to the quantum efficiency of the ideal photodetector that gives the same measurement error. This concept was used by Braunstein and Nemoto Braunstein and Nemoto (2002) to analyze another indirect photon-counting device—the homodyne detection scheme. Since their motivation was to use the homodyne detection for signaling they chose the mutual information as their measure of the measurement accuracy. Our motivation is different. We are rather interested in the analysis of the input light statistics. Let $`\rho _n`$ be the true photon-number distribution (i.e. the diagonal elements of the density matrix) of the input pulse. After $`N`$ identically prepared pulses have been detected we can make an estimate of the photon number distribution $`\overline{\rho }_n`$. Due to the finite measurement time and the presence of noise it will usually differ from the true distribution. We will take the mean quadratic distance $`d=_n(\rho _n\overline{\rho }_n)^2`$ between $`\rho _n`$ and $`\overline{\rho }_n`$ as our measure of the estimation error. Let us note that for density matrices diagonal in the Fock basis this distance equals the Hilbert-Schmidt distance between the true and estimated states. Its reciprocal divided by $`N`$, $`I=1/(dN)`$, quantifies the information about the input pulse gained from the measurement per pulse. Naturally, this distance depends not only on the measurement device but also on the chosen estimation strategy. Let us calculate the error of the optimal estimation: It is well known that no matter how clever estimation procedure is adopted, its error cannot be smaller than the Cramér-Rao bound Rao (1945); Cramér (1946). This bound applied to $`I`$ implies that

$$I=\left[N\underset{n}{}(\mathrm{\Delta }\rho _n)^2\right]^1(N\text{Tr}\text{F}^1)^1,$$
(1)
where $`F`$ is the Fisher information matrix,

$$\text{F}_{kl}=\frac{}{\rho _k}\mathrm{log}P(𝐧|\rho )\frac{}{\rho _l}\mathrm{log}P(𝐧|\rho )_𝐧,$$
(2)
and $`P(𝐧|\rho )`$ is the conditional probability of detecting the outcome n provided the true statistics of the input pulse is $`\rho `$. The averaging is done over the all possible outcomes n of the measurement. The probability $`P(𝐧|\rho )`$ completely describes the measurement apparatus. Provided the individual detection events are independent, $`P(𝐧|\rho )`$ is a multinomial distribution,

$$P(𝐧|\rho )p_n^{Nf_n}.$$
(3)
Here $`p_n`$ is the probability that the input pulse gives rise to a detection in a particular channel $`n`$. The Fisher information matrix then simplifies to

$$\text{F}_{kl}=\frac{1}{N}\underset{n}{}\frac{1}{p_n}\frac{p_n}{\rho _k}\frac{p_n}{\rho _l},$$
(4)
In the case of the ideal detector of efficiency $`\stackrel{~}{\eta }`$ the independent outcomes are just the different counted numbers of photons and the corresponding probabilities are given by the well-known Bernouli distribution,

$$p_n=\underset{mn}{}\left(\genfrac{}{}{0pt}{}{m}{n}\right)\stackrel{~}{\eta }^n(1\stackrel{~}{\eta })^{mn}\rho _m.$$
(5)
Now we will analyze the fiber-loop detector which can be thought of as a multiport device. Such a multiport, for simplicity limited to only three output channels, is shown in Fig. 2.

The simplest situation occurs when the coincidences between detections at different output channels are not registered. Then an $`s`$-channel device reduces to $`s`$ independent binary detectors of generally different quantum efficiencies. As was shown in Mogilevtsev (1998), such a detection scheme already yields enough information for the reconstruction. In this case the independent events are zero-detection events in the individual channels having probabilities

$$p_j^{}=\underset{m}{}[1\eta _j]^m\rho _m.$$
(6)
A significant gain in accuracy of the reconstruction can be expected when all possible coincidences are taken into account. A detection event of such a complex observation is best recorded in the binary notation where “1” or “0” appearing at a certain position corresponds to a “click” or “no click,” respectively, in the given output channel. The most simple nontrivial case of such a measurement is the fiber-loop detector with only two output channels (the pulse goes round the loop only once). This is equivalent to letting the transmissivity $`T_2`$ of the second beamsplitter in Fig. 2 go to unity. In this case we have four detection events. Their probabilities read

$$\begin{array}{cc}\hfill p_{00}& =\underset{m}{}\rho _m[1\eta _1T\eta _2(1T)]^m,\hfill \\ \hfill p_{10}& =\underset{m}{}\rho _m[1\eta _2(1T)]^mp_{00},\hfill \\ \hfill p_{01}& =\underset{m}{}\rho _m(1\eta _1T)^mp_{00},\hfill \\ \hfill p_{11}& =1p_{00}p_{10}p_{01},\hfill \end{array}$$
(7)
and $`T=T_1`$. It is clear that only weak input signals can fully be analyzed by the simple measurement (7). In the following we will assume that the input density matrix can be truncated and thus fully specified by giving its first few elements. The most simple case is a pulse containing at most two photons,

$$\rho =(1\rho _1\rho _2)|00|+\rho _1|11|+\rho _2|22|.$$
(8)
A successful estimation consist in identifying the single photon and two photon contributions. Using Eqs. (4), (5), and (7) the information (1) about the state (8) yielded by the fiber-loop detector can easily be calculated. As the resulting expression is rather complicated and not suitable for discussion we will consider only the leading term of its expansions in $`\eta `$s. Even best binary detectors have $`\eta `$ significantly smaller than unity guaranteeing the rapid convergence of the series. Setting $`\eta _1=\eta _2=\eta `$ the information reads,

$`I_{\text{loop}}`$ $``$ $`{\displaystyle \frac{2\tau \eta ^2}{5\rho _2}},\rho _20,`$ (9)
$`I_{\text{loop}}`$ $``$ $`{\displaystyle \frac{2\eta }{\rho _1}}+2(1\tau )\eta ^2,\rho _2=0,`$ (10)
and $`\tau =T(1T)`$. As can be seen the information is determined by the last non-vanishing $`\rho _n`$. Naturally, states with large multi-photon contribution are more difficult to estimate. If $`\rho _2`$ significantly differs from zero it is optimal to split the input pulse between the output channels equally, $`T=1/2`$, see also Fig. (3). However for states with negligible two-photon content $`\rho _20`$ the information exhibits a local minimum at this point, see Fig. (3). This minimum is caused by the term in Eq. (10) which is quadratic in $`\eta `$ and since $`\eta `$ is always significantly smaller than unity the minimum s a shallow one. In an experiment where the true state is not known, $`T=1/2`$ would be the best choice. The corresponding information gain of the ideal detector is

$`I_{\text{ideal}}`$ $``$ $`{\displaystyle \frac{\overline{\eta }^2}{5\rho _2}},\rho _20,`$ (11)
$`I_{\text{ideal}}`$ $``$ $`{\displaystyle \frac{2\overline{\eta }}{\rho _1}},\rho _2=0.`$ (12)
As can be seen, for states with negligible two-photon content the ideal and fiber-loop detectors show the same performance. This is not surprising because for single photon states there is no difference between the ideal and binary detector. However for states with significant two-photon content the information the fiber loop detector becomes inferior to the ideal one. The equivalent efficiency of the fiber-loop detector yields in the worst case can be calculated by comparing Eqs. (9) and (11),

$$\eta _{\text{eq}}\eta /\sqrt{2}.$$
(13)
The dependence of the equivalent efficiency calculated by the exact inversion of the Fisher matrix on the quantum efficiency of the binary detector used in the fiber-loop detector is shown in Fig. 4.

As one can see Eq. (13) makes an excellent approximation. Let us close this part noting that for weak input fields the information gain of the fiber-loop and ideal detectors differ only by the factor of $`\sqrt{2}`$. Further improvement can be achieved by increasing the number of output channels. Provided the noise of the detectors can be neglected, splitting output channels into two or more new channels leads to a refinement of the probability operator measure describing the detector. Any such refinement can only increase the information gain. Such many-channel fiber-loop detectors will be discussed elsewhere.

In real experiment experimenter’s control of the parameters of the fiber-loop detector is severely limited for many reasons, so it is not possible to choose the optimal configuration of the device.

Prior to its utilization the device has to be calibrated. We use coherent pulses whose statistics is known to be Poissonian for this purpose. From the probabilities of zero counts at the individual output channels the transmissions of beamsplitters in the substituting scheme are readily calculated. In the next step we calculate probabilities $`p_n`$ of all possible outcomes of the experiment. We start from the probabilities of “zero-arbitrary” detection events. Such events consist of detecting nothing in given output channels while we do not care about the rest; they read

$$\stackrel{~}{p}_j=\underset{k=0}{\overset{\text{dim}}{}}\left(1\underset{l=0}{\overset{s}{}}ϵ_{jl}\eta _l\right)^k\rho _k,$$
(14)
where dim is the chosen cutoff of the Hilbert space and $`ϵ_{jl}`$ is zero if there was an arbitrary event (denoted $``$) detected at the $`l`$-th channel and unity if there was no detection, and $`\eta _l=n_l/N`$ is the detection efficiency of the $`l`$-th channel. Now the probability of having a coincidence at the first and third channels of a three-channel detector can be written in terms of the single-count probabilities as follows: $`p_{101}=\stackrel{~}{p}_0p_{100}p_{001}p_{000}`$. The latter can, in turn, be expressed as $`p_{100}=\stackrel{~}{p}_{00}p_{000}`$ and $`p_{001}=\stackrel{~}{p}_{00}p_{000}`$, respectively. The probabilities of all $`2^s`$ detection events possible with an $`s`$-channel loop detector can be easily calculated using this simple recursive procedure. Having $`p_n`$ determined one can proceed with the reconstruction. Notice that $`p_n`$ are linear combinations of density matrix elements $`\rho _n`$,

$$p_j=\underset{i}{}c_{ji}\rho _i,j=1,\mathrm{},2^s,$$
(15)
$`c_{ji}=p_j/\rho _i`$ being defined by the parameters of the fiber-loop detector. The problem (15) is a typical linear and positive problem. The latter property follows from the non-negativity of $`\rho _n`$. Of course, the experimentally observed frequencies will generally differ from the theoretical probabilities, $`f_jp_j`$. Because of this the problem (15) has usually no exact solution and has to be solved in the statistical sense, for example via minimization of some statistical measure of distance between $`f_j`$ and $`p_j`$. We use the approach based on minimizing the Kullback-Leibler divergence (also called relative entropy) $`d=_jf_j\mathrm{log}p_j_jf_j\mathrm{log}f_j`$ between data and theory. The state $`\rho _n`$ minimizing $`d`$ is at the same time the maximum-likelihood estimator of the input state and therefore attains the Cramér-Rao bound asymptotically for large $`N`$. The minimization can be carried on e.g. via the iterative Expectation-Maximization algorithm Dempster et al. (1977); Vardi and Lee (1993):

$$\rho _k^{(n+1)}=\rho _k^{(n)}\underset{j}{}\frac{f_jc_{jk}}{_kc_{ji}\rho _i^{(n)}}.$$
(16)
Left panel in Fig. 5 shows the reconstruction of the photon statistics of a laser pulse containing $`0.6`$ photons on the average. The agreement between the measured and theoretical statistics is very good. In this case only first three output channels were sufficient for the reconstruction. The equivalent efficiency of the three-channel detector calculated for Poissonian light is shown on the right in Fig. 5. We would like to emphasize that this is the efficiency of the real experimental apparatus where a significant fraction of light is either lost inside the apparatus or not detected due to small quantum efficiency of detectors Haderka et al. . Nonetheless, in a wide range of intensities our detector is the equal of the ideal detector of quantum efficiency of about $`25\%`$.

More intense pulses require the use of a larger number of output channels. We have tested our fiber-loop detector on different kinds of light pulses. The statistics has been artificially changed by changing the intensity of subsequent Poissonian pulses using a random number generator. The statistics $`P(I)`$ of intensity fluctuations determines the true photon statistics of light pulses through the composition rule

$$\rho _n=\frac{I^n}{n!}e^IP(I)𝑑I.$$
(17)
Experimental results obtained with different light sources are summarized in Fig. 6.

Notice that the composite nature of $`\rho `$ nicely shows in the upper right panel. Also the super-Poissonian character of the input pulse is clearly seen in the lower left panel. Several tens of thousands of pulses were used for each reconstruction.

In conclusion we have demonstrated the photon-counting capability of the fiber-loop detector and quantified its efficiency by means of the Fisher information. The possibility to reconstruct the statistics of the input light pulses of moderate intensities have also been shown.

This work was partially supported by grants LN00A015, and RN19972003012 of the Czech Ministry of Education, and by Czech-Italian project No. 29, “Decoherence and Quantum Measurement.”


text:
# Bounds for Entanglement of formation of two mode squeezed thermal States

## Abstract

The upper and lower bounds of entanglement of formation are given for two mode squeezed thermal state. The bounds are compared with other entanglement measure or bounds. The entanglement distillation and the relative entropy of entanglement of infinitive squeezed state are obtained at the postulation of hashing inequality.

Keywords: bound of entanglement measure, two mode squeezed thermal state,

PACS: 03.67.-a,03.65.Bz

Quantum entanglement plays an essential role in all branches of the emerging field of quantum information. Recently a great deal of attention has been devoted to the quantum information processing with canonical continuous variables and various protocols for quantum communication and computation are developed based on continuous variables   . With continuous variable optical system these protocols can be implemented. The efficiency of entanglement manipulation protocols critically depends on the quality of the entanglement that one can generate. For practical purpose  one needs to determine the entanglement degradation in transmission of two mode squeezed vacuum state light through absorbing fibers. It is therefore essential to be able to quantify the amount of entanglement in systems with continuous variables.

Of all the continuous variable quantum states, quantum Gaussian states are of great practical importance, since these comprise essentially all the experimentally realizable continuous variable states, theoretically they might also be simple enough to be analytically tractable for the problems involved. Gaussian state is completely specified by its mean $`\eta `$ and its correlation matrix $`\alpha `$, where the mean can be dropped by local unitary operation so that is irrelevant for entanglement problems. The density operator of Gaussian state can be characterized  by its quantum characteristic function $`\chi \left(z\right)=\mathrm{exp}[i\eta ^Tz\frac{1}{2}z^T\alpha z]`$ , and any Gaussian state of two modes can be transformed into what we called the standard form, using local unitary operation only  . The corresponding characteristic function has mean $`\eta =0`$ and the correlation matrix has the simple form

$$\alpha =\left[\begin{array}{cccc}b_1,\hfill & 0,\hfill & c_1,\hfill & 0\hfill \\ 0,\hfill & b_1,\hfill & 0,\hfill & c_2\hfill \\ c_1,\hfill & 0,\hfill & b_2,\hfill & 0\hfill \\ 0,\hfill & c_2,\hfill & 0,\hfill & b_2\hfill \end{array}\right],$$
(1)
A particular simple subset of two mode Gaussian states is the two mode squeezed thermal states (TMST), with $`b_1=b_2=b,`$ $`c_1=c_2=c`$. The density operator of TMST is $`\rho _{ST}=S_2^+\left(r\right)\rho _0S_2\left(r\right),`$ where $`S_2\left(r\right)=\mathrm{exp}[r(a_1^+a_2^+a_1a_2)]`$ is the two mode squeezed operator with squeezing parameter $`r0`$, and $`\rho _0`$ is the product of two thermal states A and B whose average photon numbers are $`N_A=N_B=N`$. The relationship between the two expressions of TMST is $`b=(N+1/2)\mathrm{cosh}2r,`$ $`c=(N+1/2)\mathrm{sinh}2r`$. One can use the unit scale parameters $`\lambda =\mathrm{tanh}r`$ and $`v=N/(N+1)`$, then the inseparability criterion has a very simple form of $`\lambda >v`$ .

Even for the two parameter TMST, however, until now the amount of the entanglement in several main measures remains an open question. In this paper, I will study the upper and lower bounds for the entanglement of formation of TMST. First I will establish the lower bound of entanglement of formation, then the upper bound, then compare them with other entanglement bounds and measure, at last draw the conclusions.

Now let us consider the lower bound of entanglement of formation. It was proved that entanglement of formation is nonincreasing under local operation and classical communication . This presents us a nature way to obtain the lower bound of the entanglement of formation of a mixed state. By performing a special kind of local operation and (or) classical communication to the state, converting it to a state whose entanglement of formation is well known, the entanglement of formation of the resultant state will provide a lower bound to that of the original state.

For TMST, the simplest way of conversion is to map it on a two-qubit system by means of local operation. One of the operation for such purpose is with the Hermitian spin one-half operators $`S_j^\mu ,`$ $`\mu =A,B`$ for A, B parts respectively, $`S_1^\mu +iS_2^\mu =2_{m=0}^{\mathrm{}}|2m_{\mu \mu }2m+1|,S_3^\mu =_{m=0}^{\mathrm{}}\left(1\right)^m|m_{\mu \mu }m|,`$which satisfy the Pauli matrix algebra: $`[S_i^\mu ,S_j^\nu ]=2i\epsilon _{ijk}\delta _{\mu \nu }S_k^\mu `$.

Let us denote the mapped state with $`\rho ^{}`$ and assign the following qubit density matrix $`\rho _A^{}`$ to the original reduced density matrix $`\rho _A`$, $`\rho _A^{}=\frac{1}{2}\left(I_A+𝐒^A\sigma \right),`$where $`𝐒^A\sigma =_{i=1}^3Tr\left(\rho _AS_i^A\right)\sigma _i`$, $`\sigma _i`$ are Pauli matrices and $`I_A`$ is the identity operator on the Hilbert space of qubit A.

The map of bipartite density matrix is

$$\rho ^{}=\frac{1}{4}\left(I_AI_B+𝐒^A\sigma I_B+I_A𝐒^B\sigma +\underset{i,j}{}t_{ij}\sigma _i\sigma _j\right)$$

, with $`t_{ij}=Tr\left(\rho S_i^AS_j^B\right)`$ .

For TMST density matrix $`\rho _{ST}`$, its reduced state $`\rho _A`$ $`=Tr_B\rho _{ST}`$ turns out to be a one mode thermal state of the form $`\left(1v^{}\right)_{m=0}^{\mathrm{}}v^m|mm|`$ with $`v^{}=\frac{v+\lambda ^2}{1+v\lambda ^2}`$ . It is clear that $`Tr\left(\rho _AS_{1,2}^A\right)=0`$, and a direct calculation shows that $`Tr\left(\rho _AS_3^A\right)=\frac{1v^{}}{1+v^{}}`$. Meanwhile the two mode squeezing operator $`S_2\left(r\right)`$ will generate or annihilate the same number of bosons in mode A and mode B, it will not change the difference of the occupation numbers of the two modes, one readily gets $`t_{ij}=0`$, for $`ij`$. As for the diagonal elements $`t_{ii}`$, from the definition of $`t_{ij}`$ , it follows that $`t_{11}=t_{22}=2Tr_{m,n=0}^{\mathrm{}}\rho _{ST}|2m,2n2m+1,2n+1|.`$Making use of the coherent state representation of $`\rho _{ST}`$ ,after some algebra, one obtains $`t_{33}=\left(\frac{1v}{1+v}\right)^2`$ and

$$t_{11}=2C_0\underset{m,n=0}{\overset{\mathrm{}}{}}\underset{l=0}{\overset{2\mathrm{min}(m,n)}{}}\sqrt{\frac{2n+1}{2m+1}}\left(\begin{array}{c}2m+1\\ l+1\end{array}\right)\left(\begin{array}{c}2n\\ l\end{array}\right)\tau ^{2l+1}\omega ^{2m+2n2l},$$
(2)
with $`C_0=\frac{\left(1v\right)^2\left(1\lambda ^2\right)}{1v^2\lambda ^2}`$, $`\tau =\frac{\lambda \left(1v^2\right)}{1v^2\lambda ^2}`$ , $`\omega =\frac{v\left(1\lambda ^2\right)}{1v^2\lambda ^2}`$.

The qubits map of TMST is

$$\rho ^{}=\left[\begin{array}{cccc}\frac{1+v^2\lambda ^2}{\left(1+v\right)^2\left(1+\lambda ^2\right)},\hfill & 0,\hfill & 0,\hfill & \frac{t_{11}}{2}\hfill \\ 0,\hfill & \frac{v}{\left(1+v\right)^2},\hfill & 0,\hfill & 0\hfill \\ 0,\hfill & 0,\hfill & \frac{v}{\left(1+v\right)^2},\hfill & 0\hfill \\ \frac{t_{11}}{2},\hfill & 0,\hfill & 0,\hfill & \frac{v^2+\lambda ^2}{\left(1+v\right)^2\left(1+\lambda ^2\right)}\hfill \end{array}\right]$$
(3)
The concurrence  of $`\rho ^{}`$ will be $`C=\mathrm{max}\{0,t_{11}2v/\left(1+v\right)^2\}`$. Denote $`x=\frac{1}{2}\left(1+\sqrt{1C^2}\right)`$, the entanglement of formation of $`\rho ^{}`$ reads

$$E\left(\rho ^{}\right)=h\left(x\right)=x\mathrm{log}x(1x)\mathrm{log}(1x).$$
(4)
Throughout this paper the base of $`\mathrm{log}`$ is $`2`$. One then obtains the lower bound of entanglement of formation : $`E_{lf}\left(\rho _{ST}\right)=E\left(\rho ^{}\right)`$.

We now turn to the upper bound for entanglement of formation of quantum Gaussian state and specifically TMST. The entanglement of formation is defined as follows. Given a density matrix $`\rho `$ of bipartite quantum system, consider all possible pure state decompositions of $`\rho `$, that is, all ensembles of states $`|\psi _i`$ with probabilities $`p_i`$, such that $`\rho =_ip_i|\psi _i\psi _i|`$. For each pure state $`|\psi `$, the entanglement $`E\left(\psi \right)`$ is defined as the entropy of the reduced state of $`|\psi `$. The entanglement of formation of the mixed state $`\rho `$ is then defined as the average entanglement of the pure states of the decomposition, minimized over all decompositions of $`\rho :`$ $`E\left(\rho \right)=\mathrm{min}_ip_iE\left(\psi _i\right)`$.

From the definition, it is known that any ensemble as a decomposition of $`\rho `$ has an average entanglement no less than the entanglement of formation of $`\rho `$, and so that is a candidate in providing an upper bound for the entanglement of formation of $`\rho `$ . The aim of us is to construct a series of decompositions, and within the series find out the decomposition of minimum average entanglement. Because the series is just a subset of the set of all possible decompositions, the minimum within the series is a local minimum. It provides a relatively tight upper bound for global minimum, the entanglement of formation of $`\rho `$.

Among all the decompositions, Gaussian state decomposition is mathematically easy to deal with. Rather than decompose a Gaussian state to its pure Gaussian state ensemble, let us generate mixed Gaussian state from a pure Gaussian seed. In the viewpoint of quantum channel, the way of generating mixed Gaussian state from a pure Gaussian state is just the way of transmitting the pure Gaussian state through quantum Gaussian channel . So that the problem of how to generate mixed Gaussian state is physical, at least in principle. Suppose a pure Gaussian state is given with correlation matrix $`\alpha _s`$ and zero mean. Applying local unitary displacement operation on it yields another pure Gaussian state, with the same correlation matrix but a nonzero mean value or displacement $`\eta `$. The two states have the same amount of entanglement because local unitary operation does not change entanglement of pure state share by two parts A and B. If the displacement $`\eta `$ has a Gaussian probability distribution, then the mixture of all these Gaussian states will be a Gaussian state. Let $`P`$ be a Gaussian probability distribution with zero mean (for convenience) and covariance matrix $`M`$ for variables $`(\eta _A^q,\eta _A^p,\eta _B^q,\eta _B^p)`$ =$`\eta ^T`$. The quantum characteristic function of the mixed state $`\rho _P`$ is

$$\mathrm{exp}(i\eta ^Tz\frac{1}{2}z^T\alpha _sz)P\left(d\eta \right)=\mathrm{exp}[\frac{1}{2}z^T\left(M+\alpha _s\right)z].$$
(5)
The correlation matrix of $`\rho _P`$ is $`\alpha _P=M+\alpha _s`$. Normalizable of probability measure requires that $`M`$ should be positive semi-definite, $`M=0`$ is included to correspond to $`\delta `$ function type of probability distribution. So that arbitrary Gaussian mixed state $`\rho _P`$ can be decomposed into an ensemble of pure Gaussian states as far as the difference of the correlation matrices is positive semi-definite. That is

$$M=\alpha _P\alpha _s0.$$
(6)
So if the separable pure Gaussian state correlation matrix $`\alpha _s`$ is found under the restriction of Eq. (6), state $`\rho _P`$ should be separable. If it is impossible to find such an $`\alpha _s`$, then one has to turn to find an entangled pure Gaussian state correlation matrix $`\alpha _s`$ to fulfill Eq.(6). The minimum of the entanglement of those pure entangled Gaussian states will give out the upper bound of entanglement of formation of state $`\rho _P`$.

Consider the entangled TMST $`\rho _{ST}`$ with $`\lambda >v`$. Let

$$\alpha _s=\frac{1}{2}\left[\begin{array}{cccc}\mathrm{cosh}2x,\hfill & 0,\hfill & \mathrm{sinh}2x,\hfill & 0\hfill \\ 0,\hfill & \mathrm{cosh}2x,\hfill & 0,\hfill & \mathrm{sinh}2x\hfill \\ \mathrm{sinh}2x,\hfill & 0,\hfill & \mathrm{cosh}2x,\hfill & 0\hfill \\ 0\hfill & \mathrm{sinh}2x,\hfill & 0,\hfill & \mathrm{cosh}2x\hfill \end{array}\right],$$
(7)
be exactly the squeezed vacuum state with squeezing parameter $`x`$. $`M0`$ requires that $`b\frac{1}{2}\mathrm{cosh}2x0,(b\frac{1}{2}\mathrm{cosh}2x)^2(c\frac{1}{2}\mathrm{sinh}2x)^20.`$This can be written as $`\mathrm{cosh}2\left(rx\right)\frac{1+v^2}{1v^2}.`$Denote $`v=\mathrm{tanh}w`$, one gets $`rxw,`$ or $`xrw`$. The minimum is $`x_0=rw`$. Let $`y=\mathrm{sinh}^2x_0,`$ in explicit form of $`\lambda ,v`$, one has $`y=\frac{\left(\lambda v\right)^2}{\left(1\lambda ^2\right)\left(1v^2\right)}`$. The upper bound of entanglement of formation then will be

$$E_{uf}\left(\rho _{ST}\right)=g\left(y\right)=\left(y+1\right)\mathrm{log}\left(y+1\right)y\mathrm{log}y$$
(8)
Let us compare the entanglement bounds obtained above with other bounds and entanglement measure. First of them is the upper bound of relative entropy of entanglement. The entanglement measure is defined by the minimization of the distance of the entangled Gaussian state to the set of separable states measured by the relative entropy, its bound was obtained by restricting the minimization set to that of separable Gaussian states. We can further restrict it to the set of separable TMST $`\stackrel{~}{\rho }_{ST}=S_2\left(\stackrel{~}{r}\right)\rho _0\left(\stackrel{~}{v}\right)S_2^+\left(\stackrel{~}{r}\right)`$, with $`\stackrel{~}{\lambda }=\mathrm{tanh}\stackrel{~}{r}\stackrel{~}{v}`$. Then the minimization can be carried out. The upper bound will be

$$E_{ur}\left(\rho _{ST}\right)=2g\left(N\right)2\mathrm{log}\left(1\stackrel{~}{v}^{}\right)[\frac{1+v}{1v}\mathrm{cosh}\left(r\stackrel{~}{w}^{}\right)1]\mathrm{log}\stackrel{~}{v}^{}.$$
(9)
The minimum reaches at $`\stackrel{~}{\lambda }=\stackrel{~}{v}=\stackrel{~}{v}^{}=\mathrm{tanh}\stackrel{~}{w}^{}`$ for some $`\stackrel{~}{v}^{}`$.

The second one we should compare with is entanglement measure of logarithmic negativity , for TMST, it reads

$$E_{LN}(\rho _{ST})=\mathrm{log}[(1+\lambda )(1v)]\mathrm{log}[(1\lambda )(1+v)].$$
(10)
The last one used in comparing is coherent information. The coherent information of bipartite state $`\sigma `$ with reductions $`\sigma _A`$ and $`\sigma _B`$ is defined as $`I^\mu \left(\sigma \right)=S\left(\sigma _\mu \right)S\left(\sigma \right),`$for $`S\left(\sigma _\mu \right)S\left(\sigma \right)0`$ and $`I^\mu \left(\sigma \right)=0`$ otherwise, where $`S`$ denotes the von Neumann entropy. And if for any bipartite state the one-way distillable entanglement is no less than coherent information which is called hashing inequality, then one obtains Shannon-like formulas for the capacities . The hypothetical hashing inequality reads

$$D_{}\left(\sigma \right)I^B\left(\sigma \right)$$
(11)
where $`D_{}\left(\sigma \right)`$ is forward classical communication aided distillable entanglement of the state. I do not mean to prove this inequality, but rather to check if it is violated for some cases in continuous variable system. It is known that two way distillation of entanglement is no less than that of one way, and entanglement of formation is no less than relative entropy of entanglement, the later is no less than distillation of entanglement. In short one has $`E_f\left(\sigma \right)E_r\left(\sigma \right)D_{}\left(\sigma \right)D_{}\left(\sigma \right)`$. Combining with hashing inequality one then obtains

$$E_{uf}\left(\sigma \right)E_f\left(\sigma \right)I^B\left(\sigma \right),E_{ur}\left(\sigma \right)E_r\left(\sigma \right)I^B\left(\sigma \right).$$
(12)
For TMST $`\rho _{ST}`$, the entropy of its reduced state is $`g\left(N^{}\right)`$, with $`N^{}=\frac{v^{}}{1v^{}}=\frac{v+\lambda ^2}{\left(1v\right)\left(1\lambda ^2\right)}`$. The coherent information is

$$I^B\left(\rho _{ST}\right)=\mathrm{max}\{g\left(N^{}\right)2g\left(N\right),0\}$$
(13)
By comparison of all bounds and measure of TMST I find that: (1) the upper bound of entanglement of formation and of relative entropy of entanglement are no less than coherent information, so it can be anticipated that hashing inequality is also true for TMST or more generally for bipartite Gaussian state; (2) As the squeezing parameter tends towards zero, the lower and upper bounds of entanglement of formation tend to coincide with each other. (3) As the squeezing parameter tends towards infinitive, the upper bound of relative entropy of entanglement tends to coincide with the coherent information, so if the hashing inequality is right, the relative entropy of entanglement, the distillation entanglement and coherent information are all equal for infinitive squeezing state. (4) The entanglement measure of logarithmic negativity is at the topmost of all the bounds.

The main points of this paper are as follow: (1) The lower bound of entanglement of formation of TMST is drawn from the fact that local operation can only decreases if not preserves the entanglement, at lower squeezing side, it almost coincides with the upper bound. This result leads to the possible determination of the entanglement of formation at lower squeezing level. (2)As for the upper bound of entanglement of formation of quantum Gaussian state, an inequality is derived, for the case of TMST, an analytical expression is given out. (3)By carrying out the minimization of the upper bound of relative entropy of entanglement for TMST I determine the entanglement distillation and the relative entropy of entanglement of TMST with infinitive squeezing at the postulation of hashing inequality.


text:
# REFERENCES

Comment on ”Phase Diagram of La2-xSrxCuO4 Probed in the Infrared: Imprints of Charge Stripe Excitations”

Recently Lucarelli et al. have reported temperature-dependence of the in-plane optical reflectivity of La2-xSrxCuO4 over a wide doping range, focusing on the infrared peaks at 30 cm-1 (for $`x`$=0.12), 250 cm-1 and 510 cm-1. They interpreted the first peak (30 cm-1) as a signature of charge stripe ordering, while the latter two (250 cm-1 and 510 cm-1) are attributed to the polaronic charge excitations. However, careful readers would notice that the reported spectra are largely different from those so far measured on the same system. As we illustrate below, all these peaks are caused by an uncontrolled leakage of the $`c`$-axis reflectivity into the measured spectra.

First, we show that the absorption peaks at 250 cm-1 and 510 cm-1 are nothing but the $`c`$-axis phonon modes (A2u). The reported TO-phonon frequencies for two of the three A2u-modes coincide well with the above two frequencies. In Fig.1 we compare the inset spectra of Fig.1c in ref. with the purely $`c`$-polarized spectrum of La1.96Sr0.04CuO4. This clearly demonstrates that Lucarelli et al. observed the spectra mixed with the $`c`$-axis component. A similar mixing is more or less observed in most of their samples except for x=0 and 0.26.

Second, the $`c`$-component mixing seriously affects the reflectivity values below $``$200 cm-1. At low tempreatures, the reflectivity is close to unity for $`Eab`$, while it decreases with lowering $`\omega `$ to less than 0.5 for $`Ec`$, reflecting the incoherent charge dynamics in the $`c`$-axis direction. If a $`c`$-component is mixed into the measured in-plane spectrum, then the measured reflectivity tends to decrease with decreasing $`\omega `$, creating an artificial new absorption peak in the Kramers-Kronig transformed conductivity spectrum. Note that the spectra reported by the other groups do not show such a pronounced peak below 100 cm-1. Although Lucarelli et al. mentioned that their spectrum for $`x`$=0.12 is consistent with the result by Dumm et al., the gigantic peak at 30 cm-1 is not seen in the latter.

Finally, we point out a non-systematic doping dependence of reflectivity spectrum seen in Fig.1 of ref.. For example, the 510 cm-1 peak that is substantially weakened at $`x`$=0.15 develops again for x=0.19, which strongly suggests an accidental mixing of $`c`$-component. As to the origin of $`c`$-component mixing, there are several possibilities such as polarizer leakage, mis-cutting and multi-domain structure of crystals grown by traveling-solvent-floating-zone method. From the non-systematic spectral change with $`x`$, it is speculated that the source of $`c`$-component and the mixing rate may be different from sample to sample in ref.. It should be noted that even a small amount of admixture of the $`c`$-component seriously affects the optical spectrum, whereas neutron scattering and/or transport measurements are more robust against a few percent mixture of different crystal angles.

In summary, the three infrared peaks observed in ref. in addition to the in-plane phonon peaks do most certainly originate from the $`c`$-axis component mixed into the in-plane spectra. Therefore, neither the $`d`$-band scenario nor the charge stripe dynamics can be deduced from these experimental results.
S. Tajima1), S. Uchida2), D. van der Marel3) and D. N. Basov4)
1)Superconductivity Research Laboratory, ISTEC, Tokyo 135-0062, Japan
2)Dept. of Physics, The University of Tokyo, Tokyo 113-8656, Japan
3)Material Science Center, University of Groningen, Nijenborgh 4, 9747AG, Groningen, The Netherland
4)Dept. of Physics, University of California at San Diego, La Jolla, California 92093-0319


text:
# Level spacing statistics of classically integrable systems -Investigation along the line of the Berry-Robnik approach-

## I Introduction

For a bounded quantum system in the semiclassical limit, statistical properties of energy levels at a given energy have been intensively studied. Universal behaviors are found in statistics of unfolded energy levels, which are sequence of numbers uniquely determined by the energy levels using the mean level density obtained from the Thomas-Fermi rule(see Appendix A). For quantum systems whose classical counterparts are strongly chaotic, the unfolded level statistics are well characterized by the random matrix theory which gives level spacing distribution obeying the Wigner distribution.

For quantum systems whose classical counterparts are integrable(such systems will be referred to as classically integrable systems), level statistics were theoretically studied by Berry-Tabor, Sinai, Molchanov, Bleher, Connors and Keating, and have been the subject of many numerical investigations. Still its mechanism is not well understood, the level spacing distributions are believed to be the Poisson (exponential) distribution for generic classically integrable systems.

As suggested e.g., by Hannay (see the discussion of ), one possible explanation would be as follows: For an integrable system of $`f`$ degrees-of-freedom, almost every orbit is generically confined in each inherent torus, and whole region in the phase space is densely covered by invariant tori as suggested by the Liouville-Arnold theorem. In other words, the phase space of the integrable system consists of infinitely many tori which have \] infinitesimal volumes in Liouville measure. Then, the energy level sequence of the whole system is a superposition of sub-sequences which are contributed from those regions. Therefore, if the mean level spacing of each independent subset is large, one would expect the Poisson distribution as a result of the law of small numbers.

As shown by the pioneering work of Berry and Robnik, this would be the case for the nearly integrable systems consisting of large number of independent components. Quantum energy eigenstates are considered to be ’a superposition’ of various classical trajectories, which are connected by the tunneling effect. In the semiclassical limit, because of the suppression of the tunneling, energy eigenstates are expected to be localized on individual region in phase space explored by a typical trajectory of classical system, such as tori and chaotic regions, and to form independent components. The level spacing distribution is, then, regarded as a product of superposition of energy levels contributed respectively from those independent components. Based on this view, Berry and Robnik observed that the level spacing distribution of the system consisting of a single regular component, described by a Poisson distribution, and $`N`$ chaotic components, described by the Wigner distributions of equal strength, approaches the Poisson distribution when $`N+\mathrm{}`$.

The localization of energy eigenfunctions onto individual phase space structures was studied by Berry as a semiclassical behavior of the Wigner function. For the classically integrable system, it is shown that the Wigner function tends to delta function on tori in the semi-classical limit. Such localization phenomenon of eigenfunctions has been observed numerically for several systems. For example in Ref., it was shown that the Husimi distribution function of each energy eigenstate is well localized either on the classical chaotic regions or on the regular regions filled with classical tori and that the energy eigenvalues corresponding to eigenfunctions supported by distinct phase-space regions obey distinct statistics.

Those works seem to imply that the existence of infinitely many independent components is an essential ingredient of the appearance of Poissonian level spacing distribution and that the semiclassical limit is one of the mechanisms providing infinitely many components. However, in some classically integrable/nearly-integrable systems, where one might expect infinitely many components, deviations from the Poissonian distribution have been observed. Therefore, it is interesting to explore the consequences only of the existence of infinitely many independent components.

In this paper, along the line of thoughts of Berry and Robnik, we investigate the distribution of unfolded energy levels when energy levels consist of infinitely many independent components, and show the possibility of deviation from the Poisson distribution. Hereafter, for the sake of simplicity, unfolded energy levels will be referred to as energy levels or levels.

We briefly review the Berry-Robnik theory. It relates the statistics of the energy level distribution to the phase-space geometry by assuming that the sequence of the energy spectrum is given by the superposition of statistically independent subspectra, which are contributed respectively from eigenfunctions localized onto the invariant regions in phase space. Formation of such independent subspectra is a consequence of the condensation of energy eigenfunctions on disjoint regions in the classical phase space and of the lack of mutual overlap between their eigenfunctions, and, thus, can only be expected in the semi-classical limit where the Planck constant tends to zero, $`\mathrm{}0`$. This mechanism is sometimes referred to as the principle of uniform semi-classical condensation of eigenstates, which is based on an implicit state by Berry. The validity of this assumption is checked numerically in the semi-classical(high energy) region where the Planck volume is much smaller than the phase volume of each invariant region.

In the Berry-Robnik approach, the overall level spacing distribution is derived along a line of mathematical framework by Mehta, as follows: Consider a system whose classical phase space is decomposed into $`N`$ disjoint regions. The Liouville measures of these regions are denoted by $`\rho _i(i=1,2,3,\mathrm{},N)`$ which satisfy $`_{i=1}^N\rho _i=1`$. Let $`E(S)`$ be the gap distribution function which stands for the probability that an interval $`(0,S)`$ contains no level. $`E(S)`$ is expressed by the level spacing distribution $`P(S)`$ as follows;

$$E(S)=_S^{\mathrm{}}𝑑\sigma _\sigma ^{\mathrm{}}P(x)𝑑x.$$
(1)
When the entire sequence of energy levels is a product of statistically independent superposition of $`N`$ sub-sequences, $`E(S;N)`$ is decomposed into those of sub-sequences, $`E_i(S;\rho _i)`$,

$$E(S;N)=\underset{i=1}{\overset{N}{}}E_i(S;\rho _i).$$
(2)
In terms of the normalized level spacing distribution $`p_i(S;\rho _i)`$ of sub-sequence, $`E_i(S;\rho _i)`$ is given by

$$E_i(S;\rho _i)=\rho _i_S^{\mathrm{}}𝑑\sigma _\sigma ^{\mathrm{}}p_i(x;\rho _i)𝑑x.$$
(3)
and $`p_i(S;\rho _i)`$ is assumed to satisfy

$$_0^{\mathrm{}}Sp_i(S;\rho _i)𝑑S=\frac{1}{\rho _i}.$$
(4)
Note that the spectral components are not always unfolded automatically in general even when the total spectrum is unfolded. However, in the sufficient small interval $`[ϵ,ϵ+\mathrm{\Delta }ϵ]`$, each spectral component obeys a same scaling law (see Appendix A) and thus is unfolded automatically by an overall unfolding procedure. Equations (2) and (4) relate the level statistics in the semiclassical limit with the phase-space geometry.

In most general cases, the level spacing distribution might be singular. In such a case, it is convenient to use its cumulative distribution function $`\mu _i`$;

$$\mu _i(S)=_0^Sp_i(x;\rho _i)𝑑x.$$
(5)
The corresponding quantity of the overall level spacing distribution is

$$M(S;N)=_0^SP(x;N)𝑑x,$$
(6)
where $`P(S;N)`$ is the level spacing distribution function corresponding to $`E(S;N)`$.

In addition to equations (2) and (4), we assume two conditions for the statistical weights:

* Assumption (i): The statistical weights of independent regions uniformly vanishes in the limit of infinitely many regions;

$$\underset{i}{\mathrm{max}}\rho _i0\text{as}N+\mathrm{}.$$
(7)
* Assumption (ii): The weighted mean of the cumulative distribution of energy spacing,

$$\mu (\rho ;N)=\underset{i=1}{\overset{N}{}}\rho _i\mu _i(\rho ),$$
(8)
converges in $`N+\mathrm{}`$ to $`\overline{\mu }(\rho )`$

$$\underset{N+\mathrm{}}{lim}\mu (\rho ;N)=\overline{\mu }(\rho ).$$
(9)
The limit is uniform on each closed interval: $`0\rho S`$.

In the Berry-Robnik theory, the statistical weights of individual components are related to the phase volumes of the corresponding invariant regions. This relation is satisfactory if the Thomas-Fermi rule for the phase space fractions still holds in general. Here we do not specify their physical meaning and deal with them as parameters.

Under the assumptions (i) and (ii), eqs.(2) and (4) lead to the overall level spacing distribution whose cumulative distribution function is given by the following formula in the limit of $`N+\mathrm{}`$,

$$M_{\overline{\mu }}(S)=1\left(1\overline{\mu }(S)\right)\mathrm{exp}\left[_0^S\left(1\overline{\mu }(\sigma )\right)𝑑\sigma \right],$$
(10)
where the convergence is in the sense of the weak limit. When the level spacing distributions of individual components are sparse enough, one may expect $`\overline{\mu }=0`$ and the level spacing distribution of the whole energy sequence reduces to the Poisson distribution,

$$M_{\overline{\mu }=0}(S)=1\mathrm{exp}\left(S\right).$$
(11)
In general, one may expect $`\overline{\mu }0`$ which corresponds to a certain accumulation of the levels of individual components.

In the following sections, the above statement is proved and the limiting level spacing distributions are classified into three classes. One of them is the Poisson distribution as discussed in the original work by Berry and Robnik. The others are not Poissonian. We give examples of non-Poissonian limiting level spacing distributions in section III. In the concluding section, we discuss some relations between our results and other related works.

## II Limiting Level Spacing Distribution

### A Derivation of the limiting level spacing distribution

In this section, starting from eqs.(2) and (4), and the assumptions (i) and (ii) introduced in the previous section, we show that, in the limit of infinitely many components $`N+\mathrm{}`$, the level spacing distribution converges weakly to the distribution with the cumulative distribution function:

$$M_{\overline{\mu }}(S)=1\left(1\overline{\mu }(S)\right)\mathrm{exp}\left[_0^S\left(1\overline{\mu }(\sigma )\right)𝑑\sigma \right].$$
(12)
According to Helly’s theorem, this is equivalent to show that the cumulative distribution function $`M(S;N)`$ converges to $`M_{\overline{\mu }}(S)`$. The convergence is shown as follows.

Following the procedure by Mehta(see appendix A.2 of Ref.), we rewrite the gap distribution function $`E(S;N)`$ in terms of the cumulative level spacing distribution functions $`\mu _i(S)`$ of independent components:

$$E(S;N)=\underset{i=1}{\overset{N}{}}\left[\rho _i_S^+\mathrm{}𝑑\sigma \left(1\mu _i(\sigma )\right)\right]=\underset{i=1}{\overset{N}{}}\left[1\rho _i_0^S𝑑\sigma \left(1\mu _i(\sigma )\right)\right].$$
(13)
The second equality follows from Eq.(4), integration by parts and $`lim_{\sigma +\mathrm{}}\sigma \left(1\mu _i(\sigma )\right)=0`$, which results from the existence of the average. Then, the overall cumulative level spacing distribution function $`M(S;N)`$ is given by

$`M(S;N)`$ $`=`$ $`1+{\displaystyle \frac{d}{dS}}E(S;N)`$ (14)
$`=`$ $`1E(S;N){\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\rho _i\rho _i\mu _i(S)}{1\rho _i_0^S𝑑\sigma \left(1\mu _i(\sigma )\right)}}.`$ (15)
First we consider the behavior of $`E(S;N)`$. Since the convergence of $`_{i=1}^N\rho _i\mu _i(\sigma )\overline{\mu }(\sigma )`$ for $`N+\mathrm{}`$ is uniform on each interval $`\sigma [0,S]`$ by Assumption (ii),

$`\mathrm{log}E(S;N)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{log}\left[1\rho _i{\displaystyle _0^S}𝑑\sigma (1\mu _i(\sigma ))\right]`$ (16)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left[\rho _i{\displaystyle _0^S}𝑑\sigma (1\mu _i(\sigma ))+O(\rho _i^2)\right]`$ (17)
$`=`$ $`{\displaystyle _0^S}𝑑\sigma \left[1\mu (\sigma ;N)\right]+{\displaystyle \underset{i}{\overset{N}{}}}O(\rho _i^2)`$ (19)
$`{\displaystyle _0^S}𝑑\sigma \left[1\overline{\mu }(\sigma )\right]\text{as}N+\mathrm{},`$
where we have used $`|\mu _i(\sigma )|1`$, $`\mathrm{log}(1+ϵ)=ϵ+O(ϵ^2)`$ in $`ϵ1`$, and the following property obtained from Assumption (i)

$$|\underset{i=1}{\overset{N}{}}O(\rho _i^2)|C\underset{i}{\mathrm{max}}\rho _i\underset{i=1}{\overset{N}{}}\rho _i=C\underset{i}{\mathrm{max}}\rho _i0\text{as}N+\mathrm{}.$$
(20)
with $`C`$ a positive constant. The $`N+\mathrm{}`$ limit of the sum in the right-hand side of equation (15) can be calculated in a similar way. Indeed, as $`1/(1ϵ)=1+O(ϵ)`$ in $`ϵ1`$, one has

$`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\rho _i\rho _i\mu _i(S)}{1\rho _i_0^S𝑑\sigma (1\mu _i(\sigma ))}}`$ $`=`$ $`1{\displaystyle \underset{i=1}{\overset{N}{}}}\rho _i\mu _i(S)+{\displaystyle \underset{i=1}{\overset{N}{}}}O(\rho _i^2)`$ (22)
$`1\overline{\mu }(S)\text{as}N+\mathrm{}.`$
Therefore, we have the desired result:

$`\underset{N\mathrm{}}{lim}M(S;N)`$ $`=`$ $`M_{\overline{\mu }}(S)=1\left(1\overline{\mu }(S)\right)\mathrm{exp}\left[{\displaystyle _0^S}\left(1\overline{\mu }(\sigma )\right)𝑑\sigma \right]`$ (23)
### B Properties of the limiting level spacing distribution

Since $`\mu _i(S)`$ is monotonically increasing and $`0\mu _i(S)1`$, $`\overline{\mu }(S)`$ has the same properties. Then, $`1\overline{\mu }(S)0`$ for any $`S0`$ and one has

$$\frac{1}{S}_0^S𝑑\sigma (1\overline{\mu }(\sigma ))1\overline{\mu }(+\mathrm{})\text{as}S+\mathrm{}.$$
(24)
The limit classifies the following three cases:

* Case 1, $`\overline{\mu }(+\mathrm{})=0`$: The limiting level spacing distribution is the Poisson distribution. Note that this condition is equivalent to $`\overline{\mu }(S)=0`$ for $`{}_{}{}^{}S`$ because $`\overline{\mu }(S)`$ is monotonically increasing.
* Case 2, $`0<\overline{\mu }(+\mathrm{})<1`$: For large value of $`S`$, the limiting level spacing distribution is well approximated by the Poisson distribution, while, for small value of $`S`$, it may deviates from the Poisson distribution.
* Case 3, $`\overline{\mu }(+\mathrm{})=1`$: The limiting level spacing distribution deviates from the Poisson distribution for $`{}_{}{}^{}S`$, in such a way that the cumulative distribution function approaches 1 as $`S+\mathrm{}`$ more slowly than does the Poisson distribution. This case will be referred to as a sub-Poisson distribution.

One has Case 1 if the individual level spacing distributions are derived from scaled distribution functions $`f_i`$ as

$$\mu _i(S)=\rho _i_0^Sf_i\left(\rho _ix\right)𝑑x,$$
(25)
where $`f_i`$ satisfy

$$_0^+\mathrm{}f_i(x)𝑑x=1,_0^+\mathrm{}xf_i(x)𝑑x=1$$

and are uniformly bounded by a positive constant $`D`$: $`|f_i(S)|D`$ ($`1iN`$ and $`S0`$). Indeed, one then has

$`|\mu (S;N)|`$ $``$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\rho _i^2{\displaystyle _0^S}\left|f_i\left(\rho _ix\right)\right|𝑑x`$ (26)
$``$ $`DS{\displaystyle \underset{i=1}{\overset{N}{}}}\rho _i^2DS\underset{i}{\mathrm{max}}\rho _i{\displaystyle \underset{i=1}{\overset{N}{}}}\rho _i0\overline{\mu }(S).`$ (27)
This includes the case studied by Berry and Robnik, where the level spacing distribution is a superposition of a single regular component and $`N`$ equivalent chaotic components, and the latter is expressed by the product of the scaled distributions as in Eq.(2). Indeed, one has

$$E^{BR}(S;N)=\mathrm{exp}\left(\rho _0S\right)\underset{i=1}{\overset{N}{}}E_i^{\text{WIGNER}}(S;\rho _i),$$
(28)
where the statistical weights are $`\rho _i=\frac{1\rho _0}{N}`$ and the individual level spacing distributions $`f_i`$ corresponding to the gap distributions $`E_i^{\text{WIGNER}}(S;\rho _i)`$ are given by

$$f_i(x)=\frac{\pi x}{2}\mathrm{exp}\left[\frac{\pi }{4}x^2\right].$$
(29)
In addition, this would be the case when the system consists of identical $`N`$ components where the level spacing distribution is described by a scaled form as Eq.(25). Such a case is expected when there is a symmetry such as the regular polygonal billiards.

We remark that, when the limiting function $`\overline{\mu }(S)`$ is differentiable, the asymptotic level spacing distribution admits the following density:

$$P_{\overline{\mu }}(S)=\left[(1\overline{\mu }(S))^2+\overline{\mu }^{}(S)\right]\mathrm{exp}\left[_0^S\left(1\overline{\mu }(\sigma )\right)𝑑\sigma \right].$$
(30)
## III Example

As an example of the deviation from the Poisson distribution, we study the quantum systems whose energy levels are described by using positive integer numbers, $`m`$ and $`i`$ as follows,

$$ϵ_{m,i}=m^2+\alpha i^2,$$
(31)
where $`\alpha `$ is the system parameter. Such energy levels are given, for instance, by the rectangular billiard system whose aspect ratio of two sides is characterized by $`\alpha `$. In this paper, the unfolding transformation of the energy levels $`\{ϵ_{m,i}\}\{\overline{ϵ}_{m,i}\}`$ (see Appendix A) is done by using the leading Weyl term of the cumulative mean number of energy levels,

$$\overline{ϵ}_{m,i}\mathrm{\#}(ϵ_{m,i})=\frac{\pi }{4\sqrt{\alpha }}ϵ_{m,i}.$$
(32)
For a given energy interval $`[\overline{ϵ},\overline{ϵ}+\mathrm{\Delta }\overline{ϵ}]`$, $`i`$ or $`m`$ can be regarded as an index which classifies energy levels into spectral components. In this paper, $`i=1,2,3,\mathrm{},N`$,

$$N=\left[\sqrt{\frac{4\sqrt{\alpha }(1+\gamma )\overline{ϵ}\pi }{\alpha \pi }\text{ }}\right],$$
(33)
is adopted for classification where $`\gamma \mathrm{\Delta }\overline{ϵ}/\overline{ϵ}`$, and $`[x]`$ stands for the maximum integer which does not exceed $`x`$. The relative weight of each component, $`\rho _i(i=1,2,3,\mathrm{},N)`$, is given by

$$\rho _i=\{\begin{array}{c}\frac{4(1+\gamma )}{N\pi \gamma (1+\frac{1}{\alpha N^2})}\left(\sqrt{1+\frac{1\alpha i^2}{\alpha N^2}}\sqrt{\frac{1}{1+\gamma }+\frac{1}{\alpha N^2}\left(\frac{1}{1+\gamma }\alpha i^2\right)}\right)+O\left(\frac{1}{N^2}\right)\text{if}i<\sqrt{\frac{N^2+\frac{1}{\alpha }}{1+\gamma }},\hfill \\ \frac{4(1+\gamma )}{N\pi \gamma (1+\frac{1}{\alpha N^2})}\sqrt{1+\frac{1\alpha i^2}{\alpha N^2}}+O\left(\frac{1}{N^2}\right)\text{if}\sqrt{\frac{N^2+\frac{1}{\alpha }}{1+\gamma }}iN.\hfill \end{array}$$
(34)
As easily seen, $`\rho _i`$ satisfies the assumption (i);

$$\underset{i}{\mathrm{max}}\rho _i\frac{4}{\pi }\sqrt{1+\frac{1}{\gamma }}\times \frac{1}{N}\left(1+\frac{1}{\alpha N^2}\right)^{1/2}+O\left(\frac{1}{N^2}\right)0\text{as}N+\mathrm{}.$$
(35)
Note that the limit of infinitely many components, $`N+\mathrm{}`$, corresponds to the high energy limit, $`\overline{ϵ}+\mathrm{}`$ (see Eq.(33)), which is equivalent to the semiclassical limit. In this limit, the statistical weight of each sub-spectrum becomes sparse, since each element of $`\mu (S;N)`$, $`\rho _i\mu _i(S)`$, tends to zero: $`\rho _i\mu _i(S)\mathrm{max}_j\rho _j0`$.

In the billiard system, each spectral component obeys a same scaling law (see Appendix A), and thus is unfolded automatically by an overall unfolding procedure.

Figures 1(a)–1(c) show numerical results of the level spacing distribution $`P(S)`$ for three values of $`\alpha `$. In case that $`\alpha `$ is far from rational, $`P(S)`$ is well approximated by the Poisson distribution (Fig.1(a)). While in case that $`\alpha `$ is close to a rational, $`P(S)`$ shows large deviation from the Poisson distribution (Figs.1(b) and 1(c)). When $`\alpha `$ is a rational expressed as $`\alpha =p/q`$, where $`p`$ and $`q`$ are coprime positive integers, $`P(S)`$ is not smooth and becomes a sum of delta functions in the limit of $`ϵ+\mathrm{}`$, which are separated by a same step $`X`$,

$$X=\frac{\pi }{4\sqrt{pq}}.$$

Figures 2(a)–2(c) show the function $`|\mathrm{log}(1M(S;N))|`$ for three values of $`\alpha `$ corresponding to figures 1(a)–1(c), respectively. The dotted line in each figure corresponds to the Poisson distribution, $`M(S)=1exp(S)`$. As shown in figures 2(b) and 2(c), the distribution function for small value of $`S`$ clearly deviates from the Poisson distribution. However, for large value of $`S`$, it approaches a line whose slope is $`1`$(see the dashed line), and thus, the level spacings for large value of $`S`$ obey the Poisson distribution.

In order to compare the non-Poissonian distributions and the classification given in the previous section, we consider,

$$D(S;N)=\frac{1}{S}\mathrm{log}\left[1_0^S\left(1M(\sigma ;N)\right)𝑑\sigma \right].$$
(36)
When $`N+\mathrm{}`$, $`D(S;N)`$ approaches $`\frac{1}{S}_0^S(1\overline{\mu }(\sigma ))𝑑\sigma `$ and this function distinguishes the three cases as follows: In Case 1 i.e., where the level spacing obeys the Poisson distribution, $`lim_{N+\mathrm{}}D(S;N)=1`$. In Case2, $`lim_{N+\mathrm{}}D(S;N)c`$ as $`S+\mathrm{}`$ $`(c1)`$, and in Case 3, where the sub-Poisson distribution is expected, $`lim_{N+\mathrm{}}D(S;N)0`$ as $`S+\mathrm{}`$.

Figure 3 shows $`D(S;N)`$ for different values of $`N`$. From this, one can think that $`D(S;N)`$ for $`N=61905,S10`$ well approximates $`lim_{N+\mathrm{}}D(S;N)`$.

Figure 4 shows $`D(S;N)`$ for the three values of $`\alpha `$ corresponding to figures 1(a)–1(c), respectively. In case that $`P(S)`$ is well characterized by the Poisson distribution (figure 1(a)), the corresponding function $`D(S;N)`$ agrees with $`1`$. In case that $`P(S)`$ deviates from the Poisson distribution ( figures 1(b) and 1(c) ), $`D(S;N)`$ approaches a number less than $`1`$ for $`S+\mathrm{}`$. Therefore these results correspond to the Case 2. In this model, we have not yet observed the clear evidence of Case 3. Such a case is expected when there is stronger accumulation of the energy levels of individual components.

## IV Conclusion and Discussion

In this paper, along the line of thoughts of Berry and Robnik, we investigated the level spacing distribution of systems with infinitely many independent components and discussed its deviations from the Poisson distribution. In the semiclassical limit, reflecting infinitely fine classical phase space structure, individual energy eigenfunctions are expected to be well localized in the phase space and give independent contribution to the level statistics. Keeping this expectation in mind, we considered a situation where the system consists of infinitely many components and each of them gives an infinitesimal contribution. And by applying the arguments of Mehta, and Berry and Robnik, the limiting level spacing distribution was obtained which is described by a single monotonically increasing function $`\overline{\mu }(S)`$ of the level spacing $`S`$. The limiting distribution is classified into three cases; Case 1: Poissonian if $`\overline{\mu }(+\mathrm{})=0`$, Case 2: Poissonian for large $`S`$, but possibly not for small $`S`$ if $`0<\overline{\mu }(+\mathrm{})<1`$, and Case 3: sub-Poissonian if $`\overline{\mu }(+\mathrm{})=1`$. Thus, even when energy levels of individual components are statistically independent, non-Poisson level spacing distribution is possible.

Note that the singular level spacing distribution can be taken into account in terms of non-smooth cumulative distribution $`M_{\overline{\mu }}`$. Such a singularity is expected when there is strong accumulation of the energy levels of individual components. For certain class of systems, such accumulation is observable. One example is shown in section III where the results show clear evidence of Case 2. Another example is studied by Shnirelman, Chirikov and Shepelyansky, and Frahm and Shepelyansky for a certain type of systems which contain quasi-degeneracy result from inherent symmetry(time reversibility). As is well known, the existence of quasi-degeneracy leads to the sharp Shnirelman peak at small level spacings.

One of the interesting features of the level statistics is the level clustering which is described by a nonvanishing value of the level spacing distribution function at $`S=0`$. Level clustering is expected for integrable systems or mixed systems, but not for strongly chaotic systems due to the level repulsion. For certain class of systems, rigorous results are available; Molchanov analyzed the energy levels of a one-dimensional Schrödinger operator with random potential and showed that the level clustering arose from the localization of the eigenfunction in the semiclassical limit. Minami also analyzed a one-dimensional Schrödinger operator with $`\delta `$-potentials and reported the same result. The limiting level spacing distribution obtained in this paper possesses level clustering property. Indeed, when $`\overline{\mu }(S)`$ is differentiable, the level spacing distribution function (30) has nonvanishing value at $`S=0`$: $`P_{\overline{\mu }}(0)=1+\overline{\mu }^{}(0)>0`$.

It is also interesting to extend the Berry-Robnik distribution (28) for the level statistics of the nearly integrable system with two degree-of-freedom. The classical phase space of this system consists of regular and chaotic regions. Since the regular regions corresponding to the system consist of infinitely many independent regions, original proposal for the gap distribution by Berry and Robnik would be replaced by

$$E_{\overline{\mu }}(S;N)=\mathrm{exp}\left[\rho _0_0^S(1\overline{\mu }(\sigma ))𝑑\sigma \right]\underset{i=1}{\overset{N}{}}E_i^{\text{RMT}}(S;\rho _i),$$
(37)
where $`E_i^{\text{RMT}}(S;\rho _i)`$ is the gap distribution function of the $`i`$th chaotic component obtained from the random matrix theory. The above distribution formula is classified into cases; Case 1’, $`\overline{\mu }(+\mathrm{})=0`$: Berry-Robnik distribution, Case 2’, $`0<\overline{\mu }(+\mathrm{})<1`$: Berry-Robnik distribution for large $`S`$, but possibly not for small $`S`$, and Case 3’, $`\overline{\mu }(+\mathrm{})=1`$: A distribution function obtained by the superposition of spectral components obeying the sub-Poisson statistics and the Random matrix theory. From this classification, one can see that the new formula (37) admits deviations from the Berry-Robnik distribution when $`\overline{\mu }(+\mathrm{})0`$.

For nearly integrable systems with two degree-of-freedom, Prosen and Robnik showed numerically that the Berry-Robnik formula (28) well approximates the level spacing distributions in the high energy region that is called the Berry-Robnik regime, while it deviates in the low energy region and approximates the Brody distribution. They studied this behavior in terms of a fractional power dependence of the spacing distribution near the origin at $`S=0`$, which could be attributed to the localization properties of eigenstates on chaotic components. From the above classification, one can see that the condition $`\overline{\mu }(+\mathrm{})=0`$ should be satisfied in the Berry-Robnik regime. While Case 2’ and Case 3’ in the above classification might propose another possibilities. When the spectral components corresponding to regular regions show strong accumulation, the level spacing statistics obeys the distribution formula (37) with $`0<\overline{\mu }(+\mathrm{})1`$, and shows deviations from the Berry-Robnik distribution. Such possibilities will be studied elsewhere.
Acknowledgments
The authors would like to thank Professor M. Robnik, Professor Y.Aizawa, Professor A.Shudo, and Dr G. Veble for their helpful advice. The authors also thank the Yukawa Institute for Theoretical Physics at Kyoto University. Discussions during the YITP workshop YITP-W-02-13 on ”Quantum chaos: Present status of theory and experiment” were useful to complete this work. This work is partly supported by a Grant-in-Aid for Scientific Research (C) from the Japan Society for the Promotion of Science.

## A Unfolding of spectrum

The unfolding transformation of each energy level $`\{ϵ_n\}\{\overline{ϵ}_n\}`$ is done by using the cumulative mean number of levels up to the energy $`ϵ`$,

$$\mathrm{\#}([0,ϵ])=_0^ϵd(x)𝑑x.$$
(A1)
In the above equation, $`d(ϵ)`$ exhibits the density of energy levels obtained by the Thomas-Fermi rule:

$$d(ϵ)=\frac{V(ϵ)}{(2\pi \mathrm{})^f},V(ϵ)=\delta \left(ϵH(𝐪,𝐩)\right)d^f𝐪d^f𝐩,$$
(A2)
where $`(2\pi \mathrm{})^f`$ is the Planck volume of the system with $`f`$ degrees-of-freedom, $`V(ϵ)`$ is the phase volume on the energy surface, $`\delta (ϵH)`$ is the delta function, $`H(𝐪,𝐩)`$ is the classical Hamiltonian function, and $`(𝐪,𝐩)`$ are the coordinates and momenta in the phase space. The unfolding transformation of spectrum $`\{ϵ_n\}\{\overline{ϵ}_n\}(n=1,2,3\mathrm{})`$ is defined by

$$\overline{ϵ}_n=\mathrm{\#}([0,ϵ_n]).$$
(A3)
Here, we consider the unfolding procedure of the $`i`$-th sub-spectrum. Since the phase volume of the $`i`$-th component is $`\rho _i(ϵ)V(ϵ)`$, the density of each sub-spectrum is then described by

$$d_i(ϵ)=\frac{\rho _i(ϵ)V(ϵ)}{(2\pi \mathrm{})^f}=\rho _i(ϵ)d(ϵ)$$
(A4)
When the energy interval $`[ϵ,ϵ+\mathrm{\Delta }ϵ]`$ is sufficient small, the phase space geometry on the energy surface does not change in general. In other words, $`\rho _i(ϵ)`$ is approximated by a constant value in this ivterval, and the cumulative mean number of the $`i`$-th sub-spectrum is thus described as

$`\mathrm{\#}_i([ϵ,ϵ+\mathrm{\Delta }ϵ])`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi \mathrm{})^f}}{\displaystyle _ϵ^{ϵ+\mathrm{\Delta }ϵ}}\rho _i(e)V(e)𝑑e`$ (A5)
$``$ $`{\displaystyle \frac{\rho _i}{(2\pi \mathrm{})^f}}{\displaystyle _ϵ^{ϵ+\mathrm{\Delta }ϵ}}V(e)𝑑e=\rho _i\mathrm{\#}([ϵ,ϵ+\mathrm{\Delta }ϵ]).`$ (A6)
Therefore, in the asymptotic limit $`\mathrm{\Delta }ϵ0`$, one can see that each spectral component obeys a same scaling law and is unfolded automatically by an overall unfolding transformation:

$$\overline{ϵ}_n=\mathrm{\#}(ϵ_n)=\frac{1}{\rho _i}\mathrm{\#}_i(ϵ_n).$$
(A7)
The billiard system is very convenient since the phase space geometry does not change for variety of $`ϵ`$, and $`\mathrm{\#}_i(ϵ)=\rho _i\mathrm{\#}(ϵ)`$ is satisfactory even when the size of energy interval does not small.