text:
# Discrete Wigner functions and the phase space representation of quantum teleportation

## I Introduction

Quantum teleportation  is a scheme by which an unknown quantum state is transported between two parties by only transmitting classical information between them. This remarkable task requires also that the two parties share the two halfs of a composite system prepared in a specific entangled state. In recent years this scheme has been studied in great detail. Interest in teleportation (whose experimental feasibility has also been recently demonstrated in simple systems ) is motivated by more than one reason. In fact, teleportation is one of the most remarkable processes requiring the use of entanglement as essential resource . Moreover, teleportation can also be conceived as a primitive for universal quantum computation . In the seminal work of Bennett et al  the teleportation of the state of a qbit (a two level system) was described. Shortly afterwards, Vaidman  showed how to generalize this method to include systems with a space of states of arbitrary dimensionality. In this way, the idea of teleporting states of systems with continuous variables was first proposed. A concrete analysis of the teleportation of the state of a continuous system (the electromagnetic field) was first discussed by Braunstein and Kimble . In their work, these authors proposed (and later performed ) interesting experiments to accomplish teleportation of continuous variables. For this case, it was natural to describe the whole procedure in terms of phase space distributions . However, this is not the case for systems with a finite dimensional Hilbert space, where the use of the phase space representation is not so common. The description of the usual teleportation protocol in phase space has been recently presented by Koniorczyk, Buzek and Jansky  using the discrete version of Wigner functions originally introduced by Wooters . This approach, as mentioned in , can only be used when the dimension of the Hilbert space of the system to be teleported is a prime number. In this paper we extend the results presented in  to the case where the space of states has arbitrary dimensionality. For this purpose we use a different definition for the discrete Wigner function that turns out to be very convenient to analyze situations where entanglement between subsystems of arbitrary dimensionality is an important issue.

Several methods exist to represent the quantum state of a system with an $`N`$ dimensional Hilbert space in phase space. As mentioned above, Wooters introduced a discrete version of the Wigner function that has all the desired properties only when $`N`$ is a prime number. His phase space is an $`N\times N`$ grid (if $`N`$ is prime) and a Cartesian product of such spaces corresponding to prime factors of $`N`$ in the most general case. On the other hand, a different approach to define a Wigner function for a system with an $`N`$ dimensional Hilbert space was introduced by Leonhardt  that rediscovered results previously used by Hannay and Berry  and others. This method, that has the property of being well defined for arbitrary values of $`N`$, was used in several contexts  and recently applied to analyze the phase space representation of quantum computers and algorithms . In this case, the phase space is constructed as a grid of $`2N\times 2N`$ points where the state is represented in a redundant manner (only $`N\times N`$ of them are truly independent). In this paper, we use a hybrid approach allowing us to capture the most useful features of both Wooters and Leonhardt methods. Thus, to represent a quantum state of a bipartite system we use, following Wooters, a phase space which is a Cartesian product of two grids. Each one of these grids has, as it does in Leonhardt approach, $`2N_i\times 2N_i`$ points (where $`N_i`$ is the dimensionality of the Hilbert space of the $`i`$–th subsystem).

The paper is organized as follows: In Section II we first review the usual approach to define Wigner functions for systems with an $`N`$ dimensional Hilbert space. Then we introduce a convenient generalization that can be adopted in order to analyze bipartite (or multi–partite) systems. We discuss some general properties of the Wigner function and analyze the phase space representation of a family of entangled states (generalized Bell–EPR states). In Section 3 we show how to describe teleportation of the quantum state of a system with an $`N`$ dimensional Hilbert space using the Wigner function. In Section IV we describe the procedure by which a direct measurement of the Wigner function of a composite system can be done. In Section V we summarize our conclusions.

## II Phase space representation of composite quantum systems

We describe here the formalism of Wigner functions to represent a composite quantum system in phase space. For simplicity, we consider the composite system to be formed by subsystems each one of which has an $`N`$ dimensional Hilbert space We first describe the properties of discrete Wigner functions for one of the subsystems and later discuss the phase space representation of the composite system. It is clear that one can always split a system into subsystems in many ways. For example, when $`N`$ is a composite number one can choose to separate each subsystem into even smaller subsystems each one of which has a space of states with the dimensionality of the prime factors of $`N`$. Adopting this description is simply a matter of physical convenience. Here, we assume that the $`N`$ dimensional subsystems are the relevant elementary components and that, in some physically interesting regime, the entanglement between them can be manipulated.

### A Discrete Wigner functions

To represent the quantum state in phase space we must first define the notions of position and momentum. To do this, we introduce a basis of the Hilbert space $`B_x=\{|n,n=0,\mathrm{},N1\}`$ which we arbitrarily interpret as the position basis (with periodic boundary conditions: $`|n+N=|n`$). Given the position basis $`B_x`$, we introduce the conjugate momentum basis $`B_p=\{|k,k=0,\mathrm{},N1\}`$ by means of the discrete Fourier transform. The states of $`B_p`$ can be obtained from those of $`B_x`$ as

$$|k=\frac{1}{\sqrt{N}}\underset{n}{}\mathrm{exp}(i2\pi nk/N)|n.$$
(1)
As in the continuous case, position and momentum are related by the discrete Fourier transform. The correct semi-classical limit corresponds to the large $`N`$ limit since the dimensionality of the Hilbert space is related to an effective Planck constant as $`N=1/2\pi \mathrm{}`$.

Displacement operators in position and momentum, denoted as $`\widehat{U}`$ and $`\widehat{V}`$, are defined as: :

$`\widehat{U}^m|n`$ $`=`$ $`|n+m,\widehat{U}^m|k=\mathrm{exp}(2\pi imk/N)|k,`$ (2)
$`\widehat{V}^m|k`$ $`=`$ $`|k+m,\widehat{V}^m|n=\mathrm{exp}(i2\pi mn/N)|n.`$ (3)
Commutation relations between $`\widehat{U}`$ and $`\widehat{V}`$ directly generalize the ones corresponding to finite translations in the continuous case:

$$\widehat{V}^p\widehat{U}^q=\widehat{U}^q\widehat{V}^p\mathrm{exp}(i2\pi pq/N).$$
(4)
A reflection operator $`\widehat{R}`$ can also be defined as the one acting in the position basis as $`\widehat{R}|n=|n`$ (again, this operation is to be understood mod $`N`$). $`\widehat{R}`$ is related to the Fourier transform $`U_{FT}`$ (where $`n^{}|U_{FT}|n=\mathrm{exp}(i2\pi nn^{}/N)`$) since $`\widehat{R}=U_{FT}^2`$.

To represent the state in phase space we use the Wigner function defined as the following expectation value :

$$W(\alpha )=\text{Tr}(\widehat{A}(\alpha )\widehat{\rho })$$
(5)
where $`\alpha `$ denotes a phase space point ($`\alpha =(q,p)`$) and $`\widehat{A}(\alpha )`$ are the so–called “phase space point operators” defined in terms of displacements and reflections as

$$\widehat{A}(\alpha )=\frac{1}{2N}\widehat{U}^q\widehat{R}\widehat{V}^p\mathrm{exp}(i\pi pq/N).$$
(6)
It is important to mention that, in order for $`W(\alpha )`$ to have all the desired properties the phase space should be a grid of $`2N\times 2N`$ points which for the rest of the paper will be denoted as $`G_{2N}`$ (i.e., $`G_{2N}`$ is the set of points $`\alpha =(q,p)`$ where $`q`$ and $`p`$ take values between $`0`$ and $`2N1`$). It will also be useful to denote the first $`N\times N`$ sub-grid as $`G_N`$ (i.e., $`G_N`$ is the set of points $`\alpha =(q,p)`$ where $`q`$ and $`p`$ take values between $`0`$ and $`N1`$).

The Wigner function (5) obeys three defining properties: (P1) It is real valued, which is a consequence of the fact that the operators $`\widehat{A}(\alpha )`$ are hermitian by construction. (P2) The Wigner function can be used to compute expectation values between states as

$$\text{Tr}[\rho _A\rho _B]=N\underset{\alpha G_{2N}}{}W_A(\alpha )W_B(\alpha ).$$
(7)
This is a consequence of the completeness of the set of operators $`\widehat{A}(\alpha )`$, that satisfy

$$\text{Tr}[\widehat{A}(\alpha )\widehat{A}(\alpha ^{})]=\frac{1}{4N}\delta _N(q^{}q)\delta _N(p^{}p)$$
(8)
where $`\alpha ,\alpha ^{}G_N`$ and $`\delta _N(q)\frac{1}{N}_{n=0}^{N1}\text{e}^{2\pi iq/N}`$ is the periodic delta function (which is zero unless $`q=0\text{mod}N`$).

As $`\widehat{A}(\alpha )`$ form a complete orthonormal basis of the space of operators, one can expand the density matrix $`\rho `$ in this basis and show that the Wigner function are nothing but the coefficients of such expansion. Thus:

$`\widehat{\rho }`$ $`=`$ $`4N{\displaystyle \underset{\alpha G_N}{}}W(\alpha )\widehat{A}(\alpha )`$ (9)
$`=`$ $`N{\displaystyle \underset{\stackrel{~}{\alpha }G_{2N}}{}}W(\stackrel{~}{\alpha })\widehat{A}(\stackrel{~}{\alpha }).`$ (10)
The last expression, where the sum is over $`\alpha G_{2N}`$, can be obtained from (9) by noticing that the contribution of each of the four $`N\times N`$ sub-grids are identical. This can be shown by using the fact that (for $`\sigma _q,\sigma _p=0,1`$):

$$\widehat{A}(q+\sigma _qN,p+\sigma _pN)=\widehat{A}(q,p)(1)^{\sigma _pq+\sigma _qp+\sigma _q\sigma _pN}.$$
(11)
Finally, $`W(\alpha )`$ satisfies a third crucial property. Consider a line $`L`$ in the grid $`G_{2N}`$ (a line $`L`$ is the set of all points $`\alpha =(q,p)`$ such that $`apbq=c`$ for given integers $`a`$, $`b`$ and $`c`$). Then the sum of $`W(\alpha )`$ over all points lying on any line $`L`$ is always positive. This property (P3) is a consequence of the following fact: Adding all phase space operators over a line $`L`$ defined by the equation $`apbq=c`$ one obtains a projection operator. Thus $`\widehat{A}_L=_{\alpha L}\widehat{A}(\alpha )`$ is a projector onto an eigenspace of the phase space translation operator

$$\widehat{T}(a,b)=U^aV^b\mathrm{exp}(i\pi ab/N)$$
(12)
with eigenvalue $`\mathrm{exp}(i\pi c/N)`$. (If $`T(a,b)`$ does not have $`\mathrm{exp}(i\pi c/N)`$ as one of its eigenvalues the projector $`\widehat{A}_L`$ vanishes). The simplest example of the use of this property is the following: consider the horizontal lines $`p=c`$. The sum of $`W(\alpha )`$ over these lines vanishes if $`c`$ is odd (because $`T(1,0)=U`$ has eigenvalues $`\mathrm{exp}(i\pi k/N)`$ when $`k`$ is an even integer). On the other hand, when $`c`$ is even, the sum of $`W(\alpha )`$ is equal to the probability of measuring a momentum equal to $`c/2`$. Thus, $`_qW(q,p)=p/2|\widehat{\rho }|p/2`$ if $`p`$ is even (and zero otherwise).

Let us summarize the results presented so far: The Wigner function for systems with an $`N`$ dimensional Hilbert space is defined in (5) as the expectation value of the phase space operator $`\widehat{A}(\alpha )`$ given in (6). This definition is such that $`W(\alpha )`$ is real, it can be used to compute inner products between states and it gives all the correct marginal distributions when added over any line in the phase space, which is a grid $`G_{2N}`$ with $`2N\times 2N`$ points. The size of the phase space grid is important to obtain a Wigner function with all the desired properties. The values of $`W(\alpha )`$ on the sub-grid $`G_N`$ are enough to reconstruct the rest of the phase space (since the set $`\widehat{A}(\alpha )`$ is complete when $`\alpha `$ belong to the grid $`G_N`$). However, the redundancy introduced by the doubling of the number of sites in $`q`$ and $`p`$ is essential when one imposes the condition that all the marginal distributions should be obtained from the Wigner function.

As an example, we show here the Wigner function of a position eigenstate $`\widehat{\rho }=|q_0q_0|`$ (see  for other interesting examples):

$$W(q,p)=\frac{1}{2N}\delta _N(q2q_0)(1)^{p(q2q_0)_N}.$$
(13)
Here, $`z_N`$ denotes $`z`$ modulo $`N`$. $`W(\alpha )`$ is non–zero only on two vertical strips located at $`q=2q_0`$ modulo $`N`$. When $`q=2q_0`$ we have $`W(2q_0,p)=1/2N`$ while for $`q=2q_0\pm N`$ the modulating factor $`(1)^p`$ produces oscillations. They can be interpreted as the interference between the $`q=2q_0`$ strip and its mirror image induced by the periodic boundary conditions. The fact that $`W(\alpha )`$ becomes negative in this interference strip is essential to recover the correct marginal distributions. Adding the values of $`W(q,p)`$ along a vertical line gives the probability of measuring $`q/2`$ which should be $`1`$ for $`q=2q_0`$ and zero otherwise.

### B Wigner function for composite systems

We now consider a composite system with Hilbert space $`^{(1,2)}=^{(1)}^{(2)}`$ (for simplicity we assume that the dimension of both spaces $`^{(i)}`$ is the same but the method can be generalized if this is not the case). To represent the states of this composite system in phase space we generalize the approach described in the previous subsection in an obvious way: We consider the phase space for the composite system as the Cartesian product of the ones for the subsystems (as in the classical case) and use the phase space point operators defined as

$$\widehat{A}(\alpha _1,\alpha _2)=\widehat{A}(\alpha _1)\widehat{A}(\alpha _2).$$
(14)
The set $`\{\widehat{A}(\alpha _1,\alpha _2)`$, with $`\alpha _iG_N\}`$, is a complete orthonormal basis of the space of operators on $`^{(1,2)}`$ since:

$`\text{Tr}_{1,2}(\widehat{A}(\alpha _1,\alpha _2)\widehat{A}(\alpha _1^{},\alpha _2^{}))={\displaystyle \frac{1}{(4N)^2}}`$ $`\delta _N(\alpha _1\alpha _1^{})`$ (15)
$`\times `$ $`\delta _N(\alpha _2\alpha _2^{}),`$ (16)
where $`\alpha _i,\alpha _i^{}G_N`$ and $`\delta _N(\alpha )=\delta _N(q)\delta _N(p)`$.

The Wigner function for the composite system is defined as the expectation value of the above operators:

$$W(\alpha _1,\alpha _2)=\text{Tr}(\widehat{A}(\alpha _1,\alpha _2)\rho )$$
(17)
This function has the same properties than the one for each subsystem. In fact, the three properties (P1–P3) generalize trivially to this case. Reality (P1) is once again an obvious consequence of the hermitian nature of the phase space point operators. Completeness of the operators $`\widehat{A}(\alpha _1,\alpha _2)`$ enable us to expand the total density matrix in this basis and write:

$`\rho `$ $`=`$ $`(4N)^2{\displaystyle \underset{\alpha _1,\alpha _2G_N}{}}W(\alpha _1,\alpha _2)\widehat{A}(\alpha _1,\alpha _2),`$ (18)
$`=`$ $`N^2{\displaystyle \underset{\alpha _1,\alpha _2G_{2N}}{}}W(\alpha _1,\alpha _2)\widehat{A}(\alpha _1,\alpha _2).`$ (19)
Once again, the first line (where both $`\alpha _1,\alpha _2G_N`$) can be transformed into the second line where the sum can be extended to the grid $`G_{2N}`$ by noticing that the contribution over sub-grids are identical (due to the relations obtained by using (11)). Using this, it is simple to show that inner products between two states of the composite system can be computed from the Wigner functions as (P2):

$$\text{Tr}(\rho _A\rho _B)=N^2\underset{\alpha _1,\alpha _2G_{2N}}{}W_A(\alpha _1,\alpha _2)W_B(\alpha _1,\alpha _2).$$
(20)
The third property (P3) is valid as well but it is worth discussing it with more detail. The Wigner function turns out to be positive when added over any “slice” of the total phase space. A “slice” in phase space can be defined (following Wooters ) as a set of lines $`\{L_1,L_2\}`$ (one line for each subsystem). These sets are called slices since in the continuous limit, the set of phase space points satisfying the equations $`a_1p_1b_1q_1=c_1`$ (that defines $`L_1`$) and $`a_2p_2b_2q_2=c_2`$ (defining $`L_2`$) form a two dimensional manifold immersed in the four dimensional phase space. The fact that the Wigner function is positive when added over all points belonging to the slice $`\{L_1,L_2\}`$ is obvious: Thus, adding $`\widehat{A}(\alpha _1,\alpha _2)`$ over all points where $`\alpha _1`$ and $`\alpha _2`$ respectively belong to the lines $`L_1`$ and $`L_2`$ one always obtains a tensor product of two projectors. These operators project onto the eigenspaces of the operators

$`T_1`$ $`=`$ $`T(a_1,b_1)\text{I},`$ (21)
$`T_2`$ $`=`$ $`\text{I}T(a_2,b_2)`$ (22)
with eigenvalue $`\mathrm{exp}(i\pi c_1/N)`$ and $`\mathrm{exp}(i\pi c_2/N)`$ respectively. We omit the proof of this property, that can be done by copying the one presented in  (see also below for a simple proof of a related property). Slices defined as above, by picking a line for each subsystem, will be denoted “separable slices”. As we have just seen, Wigner functions when added over separable slices can be used to compute probabilities for the outcomes of measurements of separable observables. These are properties that are measured by independent experiments on the two separate subsystems.

It is possible, however, to consider more general non–separable slices in the phase space. As we will show below, by adding phase space operators on these non–separable slices we will obtain projectors onto entangled states. For this we should use a general kind of manifold $`L_{1,2}`$ that can be defined as the set of points $`(\alpha _1,\alpha _2)`$ satisfying the condition $`a_1p_1b_1q_1+a_2p_2b_2q_2=c_{12}`$ (notice that calling $`L_{12}`$ a line can be misleading: in the continuous limit it is a three dimensional manifold). It is simple to show that by adding phase space operators over all points belonging to $`L_{1,2}`$ one also obtains a projection operator. It is rather instructive to see this proof: For this, we just have to use the fact that the Fourier transform of the subsystem’s phase space point operators $`\widehat{A}(\alpha _i)`$ is a translation :

$$T(a,b)=\underset{q,p=0}{\overset{2N1}{}}\widehat{A}(q,p)\mathrm{exp}(i\frac{2\pi }{2N}(apbq)).$$
(23)
Using this, we can compute the sum of $`\widehat{A}(\alpha _1,\alpha _2)`$ over $`L_{1,2}`$ as follows:

$`\widehat{A}_{L_{1,2}}`$ $``$ $`{\displaystyle \underset{(\alpha _1,\alpha _2)S_{1,2}}{}}\widehat{A}(\alpha _1,\alpha _2)`$ (24)
$`=`$ $`{\displaystyle \underset{\alpha _1,\alpha _2G_{2N}}{}}\widehat{A}(\alpha _1)\widehat{A}(\alpha _2)`$ (25)
$`\times \delta _{2N}(a_1p_1b_1q_1+a_2p_2b_2q_2c_{12})`$ (26)
$`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{\lambda =0}{\overset{2N1}{}}}T^\lambda (a_1,b_1)T^\lambda (a_2,b_2)\text{e}^{i\frac{\pi }{N}\lambda c_{12}},`$ (27)
where to obtain the last line we represented the delta–function as a sum of exponentials and used (23). As the translation operators are unitary and cyclic one can always express them in terms of their eigenstates and eigenvalues as

$$T(a_i,b_i)=\underset{n=0}{\overset{N1}{}}|\varphi _{i,n}\varphi _{i,n}|\mathrm{exp}(i2\pi n/N).$$
(28)
Using this equation one finally obtains:

$`\widehat{A}_{L_{1,2}}={\displaystyle \underset{n,m}{}}|\varphi _{1,n}\varphi _{1,n}|`$ $``$ $`|\varphi _{2,m}\varphi _{2,m}|`$ (29)
$`\times `$ $`\delta _{2N}(n+mc_{1,2}/2).`$ (30)
This explicitly shows that $`A_{L_{1,2}}`$ is a projector onto the eigenspaces of a collective operator of the bipartite system. Thus, (30) projects onto eigenspaces of

$$\widehat{T}_{1,2}=\widehat{T}(a_1,b_1)\widehat{T}(a_2,b_2)$$
(31)
with eigenvalue $`\mathrm{exp}(i\pi c_{1,2}/N)`$. It is clear that, when $`N`$ is even the projector is non-vanishing only if $`c_{1,2}`$ is even. The most important conclusion is that, generically, the sum of phase space point operators over non–separable manifold $`L_{12}`$ will correspond to a projector over an entangled state. In the previous case $`\widehat{T}_{1,2}`$ is a collective operator generating simultaneous phase space translations of both systems. In the continuous limit it is clear that this operator is generated by linear combinations of the momenta and coordinates of the two subsystems.

In general, a manifold $`L_{12}`$ will not be associated with a one dimensional projector. For example, the manifold $`L_{+,p_\beta }`$ defined by the equation $`p_1+p_2=2p_\beta `$ corresponds to an $`N`$ dimensional subspace. The same is true for the manifold $`L_{,q_\beta }`$ defined as the set of points satisfying $`q_1q_2=2q_\beta `$. The intersection between these two sets will be denoted as $`L_\beta `$ (defined as the set of points belonging to both $`L_{+,p_\beta }`$ and $`L_{,q_\beta }`$) and corresponds to a one dimensional projector over an entangled state (see below).

Wigner functions of separable states have very different features than those of entangled states: In fact, if the density matrix is a tensor product $`\rho =\rho _1\rho _2`$ then $`W(\alpha _1,\alpha _2)`$ is a product of the form $`W(\alpha _1,\alpha _2)=W_1(\alpha _1)W_2(\alpha _2)`$. More generally, if the state is separable (i.e., the density matrix is a convex sum of tensor products) then the Wigner function is a convex sum of products of the above form. For entangled states (states which are not separable) this is not the case as will be explicitly seen below. It is also useful to notice that reduced Wigner functions can be computed for one subsystem by adding the total Wigner function over the complementary half of the phase space. Thus, summing $`W(\alpha _1,\alpha _2)`$ over $`\alpha _2`$ is equivalent to taking the partial trace over the second subsystem since

$`W_1(\alpha _1)`$ $`=`$ $`{\displaystyle \underset{\alpha _2G_{2N}}{}}W(\alpha _1,\alpha _2)`$ (32)
$`=`$ $`\text{Tr}_1(\widehat{A}(\alpha _1)\rho _1),`$ (33)
where $`\rho _1`$ is the reduced density matrix of the first system obtained from the total density matrix as $`\rho _1=\text{Tr}_2\rho `$. Finally, it is also useful to notice other properties of the composite Wigner function when we sum it over half of the phase space. For example, if $`\rho _A`$ and $`\rho _B`$ are two states of the composite system, then

$`F(\alpha _1,\alpha _1^{})`$ $``$ $`{\displaystyle \underset{\alpha _2G_N}{}}W_A(\alpha _1,\alpha _2)W_B(\alpha _1^{},\alpha _2)`$ (34)
$`=`$ $`{\displaystyle \frac{1}{4N}}\text{Tr}_2(\text{Tr}_1(\rho _A\widehat{A}(\alpha _1))\text{Tr}_1(\rho _B\widehat{A}(\alpha _1^{})),`$ (35)
an equation that will useful later.

### C Wigner functions for Bell states

Let us first introduce a complete basis of entangled states (generalized Bell states ). First we define the state $`|\mathrm{\Theta }_0`$ as

$$|\mathrm{\Theta }_0=\frac{1}{\sqrt{N}}\underset{n=0}{\overset{N1}{}}|n_1|n_2.$$
(36)
This pure state for the composite system is maximally entangled since the reduced density matrix of each subsystem is proportional to the identity. A complete basis of entangled states can be defined from $`|\mathrm{\Theta }_0`$ as follows:

$$|\mathrm{\Theta }_\beta =\frac{1}{\sqrt{N}}\underset{n=0}{\overset{N1}{}}\text{e}^{i2\pi p_\beta n/N}|n_1|nq_\beta _2.$$
(37)
These states can all be obtained from $`|\mathrm{\Theta }_0`$ by using one of the following equivalent expression:

$`|\mathrm{\Theta }_\beta `$ $`=`$ $`T_1(q_\beta ,p_\beta )\text{e}^{i\pi q_\beta p_\beta /N}\text{I}_2|\mathrm{\Theta }_0`$ (38)
$`=`$ $`V_1^{p_\beta }U_2^{q_\beta }|\mathrm{\Theta }_0`$ (39)
Above, the notation $`\beta =(q_\beta ,p_\beta )`$ was used. When $`\beta G_N`$ these states form a complete orthonormal basis of the Hilbert space $`^{(1,2)}`$. In fact, these $`N^2`$ vectors satisfy

$$\mathrm{\Theta }_\beta |\mathrm{\Theta }_\beta ^{}=\delta _N(\beta \beta ^{}).$$
(41)
It is interesting to notice that these are the common eigenstates of the following complete set of commuting operators:

$`U_+`$ $`=`$ $`U_1U_2,`$ (42)
$`V_{}`$ $`=`$ $`V_1V_2^{}.`$ (43)
In fact, $`V_{}`$ and $`U_+`$ commute and that Bell states (37) obey the following identities:

$`U_+|\mathrm{\Theta }_\beta `$ $`=`$ $`\mathrm{exp}(i2\pi p_\beta /N)|\mathrm{\Theta }_\beta ,`$ (44)
$`V_{}|\mathrm{\Theta }_\beta `$ $`=`$ $`\mathrm{exp}(i2\pi q_\beta /N)|\mathrm{\Theta }_\beta .`$ (45)
These expressions allow us to better understand the nature of Bell states: $`U_+`$ displaces both systems in position by the same amount while $`V_{}`$ displaces them in momentum in opposite direction. As Bell states are common eigenstates of these operators, such states can be interpreted as corresponding to the eigenstates of the total momentum and relative position (note that in the continuum limit $`U_+=\mathrm{exp}(i(P_1+P_2)\delta x/\mathrm{})`$ and $`V_{}=\mathrm{exp}(i\delta p(Q_1Q_2)/\mathrm{})`$. In this sense, Bell states (37) are precisely the ones used by Einstein, Podolsky and Rosen  to present their argument against completeness of quantum mechanics as a description of nature.

Having this in mind, one expects the phase space representation of Bell states to be rather simple. This is indeed the case:

$`W_\beta (\alpha _1,\alpha _2)`$ $`=`$ $`\text{Tr}(|\mathrm{\Theta }_\beta \mathrm{\Theta }_\beta |\widehat{A}(\alpha _1)\widehat{A}(\alpha _2))`$ (46)
$`=`$ $`W_0(\alpha _12\beta ,\alpha _2)`$ (47)
where

$`W_0(\alpha )`$ $`=`$ $`\text{Tr}(|\mathrm{\Theta }_0\mathrm{\Theta }_0|\widehat{A}(\alpha _1)\widehat{A}(\alpha _2))`$ (48)
$`=`$ $`{\displaystyle \frac{1}{(2N)^2}}\delta _N(q_{\alpha _1}q_{\alpha _2})\delta _N(p_{\alpha _1}+p_{\alpha _2})`$ (49)
$`\times (1)^{(q_{\alpha _1}p_{\alpha _1}+q_{\alpha _2}p_{\alpha _2})/N}`$ (50)
Thus, the Wigner function of $`|\mathrm{\Theta }_\beta `$ is nonzero only when $`q_{\alpha _1}q_{\alpha _2}=2q_\beta `$ and $`p_{\alpha _1}+p_{\alpha _2}=2p_\beta `$ (modulo $`N`$). The oscillations, whose origin we described above for a simpler case (13), can also be interpreted as due to the interference with the mirror images created by the boundary conditions . Notice that these are precisely the equations defining the manifold $`L_\beta `$. Thus, the projectors onto Bell state $`|\mathrm{\Theta }_\beta `$ is simply the sum of phase space point operators over all points belonging to $`L_\beta `$.

Before explicitly discussing teleportation it is useful to present some further results related to Bell states and their connection to phase space point operators. A complete basis of the space of operators on $`^{(1,2)}`$ can be constructed using Bell states: Thus, the operators $`\widehat{B}(\beta _1,\beta _2)=|\mathrm{\Theta }_{\beta _1}\mathrm{\Theta }_{\beta _2}|`$ form a complete orthogonal set (with $`\beta _1,\beta _2G_N`$). The change of basis between this set and phase space point operators is

$$\widehat{B}(\beta _1,\beta _2)=(4N)^2\underset{\alpha _1,\alpha _2G_N}{}K(\beta _1,\beta _2|\alpha _1,\alpha _2)\widehat{A}(\alpha _1,\alpha _2).$$
(51)
where the coefficients $`K(\beta _1,\beta _2|\alpha _1,\alpha _2)`$ are, in general, complex numbers (when $`\beta _1=\beta _2`$ we have $`K(\beta ,\beta |\alpha _1,\alpha _2)=W_\beta (\alpha _1,\alpha _2)`$). The precise form of this coefficients can be easily obtained but will not be needed here. The inverse of (51) is also useful and reads

$$\widehat{A}(\alpha _1,\alpha _2)=\underset{\beta _1,\beta _2}{}\stackrel{~}{K}(\alpha _1,\alpha _2|\beta _1,\beta _2)\widehat{B}(\beta _1,\beta _2)$$
(52)
Simple relations between the coefficients of (51) and (52) exist. In particular, one can show that

$`\stackrel{~}{K}(\alpha _1,\alpha _2|\beta _1,\beta _2)`$ $`=`$ $`K(\beta _2,\beta _1|\alpha _1,\alpha _2)`$ (53)
$`=`$ $`K^{}(\beta _1,\beta _2|\alpha _1,\alpha _2)`$ (54)
Finally, we mention that Bell states satisfy the following identity:

$$\text{Tr}_1(\widehat{A}(\alpha )\text{I}_2|\mathrm{\Theta }_0\mathrm{\Theta }_0|)=\frac{1}{N}\text{I}_1\widehat{A}^T(\alpha ),$$
(55)
where the transpose of the phase space point operator appears in the right hand side.

## III Teleportation in phase space

We show here how the usual teleportation protocol can be described in phase space. We consider three identical subsystems labeled by the integers $`j=1,2,3`$. The aim is to teleport the initial state of system $`1`$, which is characterized by an arbitrary Wigner function $`W(\alpha _1)`$. For this we initially prepare systems $`2`$ and $`3`$ in one of the Bell states (for simplicity we use $`|\mathrm{\Theta }_0_{2,3}`$ as the initial state). Thus, the initial density matrix of the combined three-partite system is

$`\rho _{1,2,3}`$ $`=`$ $`\rho _1|\mathrm{\Theta }_0\mathrm{\Theta }_0|_{2,3}`$ (56)
$`=`$ $`N^3{\displaystyle \underset{\alpha _1,\alpha _2,\alpha _3G_{2N}}{}}W(\alpha _1)W_0(\alpha _2,\alpha _3)`$ (57)
$`\times \widehat{A}(\alpha _1)\widehat{A}(\alpha _2)\widehat{A}(\alpha _3)`$ (58)
where the Wigner function for the Bell state $`|\mathrm{\Theta }_\beta `$ is given in (50). After preparing this initial state the teleportation protocol proceeds as follows: First we perform a measurement of system $`1`$ and $`2`$ that projects them into the Bell basis. Physically, as discussed above, this corresponds to a collective measurement that determines the total momentum $`p_1+p_2`$ and the relative coordinate $`q_1q_2`$ for these two subsystems. After this measurement the state of the system is projected into the resulting state $`|\mathrm{\Theta }_\beta _{1,2}`$, where $`\beta =(q_\beta ,p_\beta )`$ ($`p_\beta `$ and $`q_\beta `$ are the measured values of the total momentum and distance). To explicitly write down the resulting state, it is better to rewrite equation (58) expressing the phase space point operator $`\widehat{A}(\alpha _1)\widehat{A}(\alpha _2)`$ in terms of the Bell operators $`B(\beta _1,\beta _2)`$ as in (51). Thus,

$`\rho _{1,2,3}`$ $`=`$ $`(4N)^3{\displaystyle \underset{\genfrac{}{}{0pt}{}{\alpha _jG_N}{\beta _kG_N}}{}}W(\alpha _1)W_0(\alpha _2,\alpha _3)\stackrel{~}{K}(\alpha _1,\alpha _2|\beta _1,\beta _2)`$ (59)
$`\times \widehat{B}(\beta _1,\beta _2)\widehat{A}(\alpha _3)`$ (60)
From this equation it is obvious that, after the Bell measurement of the first two subsystems, only the terms with $`\beta _1=\beta _2=\beta `$ in the above expression survive. Therefore, the state of the third subsystem is

$$\rho _3^{}=4N\underset{\alpha _3G_N}{}W^{}(\alpha _3)\widehat{A}(\alpha _3)$$
(61)
where the term multiplying the phase space point operator in (60) can be identified as the new Wigner function of the third system (up to a normalization constant). In turn, this Wigner function can be seen to be defined by an expression involving a sum that contains the initial Wigner function of the first system. Thus, this can be written as

$$W^{}(\alpha _3)=\underset{\alpha _1G_N}{}Z(\alpha _3,\alpha _1)W(\alpha _1).$$
(62)
The matrix $`Z(\alpha _3,\alpha _1)`$ simply tells us how to construct the final Wigner function for the third system from the initial Wigner function of the first one. Again, the explicit form of this matrix is read from (60):

$$Z(\alpha _3,\alpha _1)=(2N)^4\underset{\alpha _2G_N}{}W_0(\alpha _3,\alpha _2)W_\beta (\alpha _1,\alpha _2).$$
(63)
The above expression is easy to evaluate since it contains a sum over half of the total phase space. Thus, we can use the relation (35) to simplify it. Doing this (and using (55)) one discovers that $`Z(\alpha _3,\alpha _1)`$ is just the trace of a product of two phase space point operators acting on the second subsystem which are evaluated at points $`\alpha _1`$ and $`\alpha _32\beta `$. Therefore, taking into account the orthogonality conditions for phase space point operators, one obtains

$$Z(\alpha _3,\alpha _1)=\delta _N(\alpha _3\alpha _12\beta ).$$
(64)
This means that the state of the third subsystem has a Wigner function that is displaced in phase space by an amount $`\beta `$ with respect to the initial state of the first system. Therefore, to recover the initial state one must displace the third system in phase space by applying the evolution operator $`\widehat{T}(\beta )=\widehat{T}(q_\beta ,p_\beta )`$ (in fact, one can show that when the operator $`\widehat{T}(\beta )`$ acts on the state of the system, the Wigner function is simply shifted in phase space by the amounts $`(2a,2b)`$ the factors of $`2`$ being originated in the fact that the grid has $`2N\times 2N`$ points ). Obviously, the recovery operation depends on $`\beta `$, which is revealed by the Bell measurement performed on the first two subsystems. In this way the final state of the third system is always identical to the unknown initial state of the first subsystem, which is the goal of the teleportation protocol.

## IV Measuring the Wigner Function

The Wigner function, as mentioned above, is in one to one correspondence with the quantum state of the system. Thus, experimentally determining the value of $`W(\alpha _1,\alpha _2)`$ in every phase space point is equivalent to completely determining the state of the system. Experimentally reconstructing the Wigner function has been the goal of a series of experiments, all dealing with continuous systems . In general, these experiments are aimed at determining first marginal distribution of some observables and later reconstructing the Wigner function by means of a Radon–like transform. However, it is interesting (and useful) to realize that there is no need to fully determine the quantum state to evaluate the Wigner function at a given phase space point. Indeed, this was proposed originally by Davidovich and Lutterbach for the Wigner function of the state of the electromagnetic field stored in a superconducting cavity . Using a minor variation of this method the direct measurement of the Wigner function at the origin of phase space has been recently performed in a cavity QED experiment. More recently, this method was generalized and shown to be applicable to the measurement of discrete (or continuous) Wigner functions of generic systems in . Here, we will describe how this tomographic scheme can be generalized to directly measure the Wigner function we presented in this paper. The fact that this Wigner function can be measured in an efficient way (i.e., in a number of steps which scales polynomially with the dimensionality of the system) may be important if, for example, one is able to relate interesting physical quantities (like some entanglement measures) with phse space observables (work is in progress in this direction).

The efficient strategy to measure the Wigner function of a composite system at (any) given phase space point is, as mentioned above, a direct generalization of the idea originally proposed in  to measure the Wigner function of an $`N`$ dimensional system. The basic ingredient can be described in terms of the following quantum algorithm. Consider a system initially prepared in a quantum state $`\widehat{\rho }`$. We put this system in contact with an ancillary qbit prepared in the state $`|0`$. This ancillary qbit plays the role of a “probe particle” in a scattering–like experiment. The algorithm is: i) Apply an Hadamard transform to the ancillary qbit (where $`H|0=(|0+|1)/\sqrt{2}`$, $`H|1=(|0|1)/\sqrt{2}`$), ii) Apply a “controlled–$`\widehat{M}`$” operator (if the ancilla is in state $`|0`$ this operator acts as the identity for the system but if the state of the ancilla is $`|1`$ it acts as the unitary operator $`\widehat{M}`$ on the system), iii) Apply another Hadamard gate to the ancilla and finally perform a weak measurement on this qbit detecting its polarization (i.e., measuring the expectation values of Pauli operators $`\sigma _z`$ and $`\sigma _y`$). It is easy to show that the above algorithm has the following remarkable property:

$$\sigma _z=Re(Tr(\widehat{M}\widehat{\rho })),\sigma _y=Im(Tr(\widehat{M}\widehat{\rho })).$$
(65)
Thus, the final polarization measurement of the ancillary qbit reveals a property determined both by the initial state $`\widehat{\rho }`$ and the unitary operator $`\widehat{M}`$.

In  we discussed how to view this simple algorithm as the basic tool to construct a rather general tomographer (and also a rather general spectrometer). In particular, we showed how to use it to measure the Wigner function of a simple system. Here, we show how to adapt it to measure the Wigner function of the composite system we have been discussing so far. This can be done by applying the algorithm shown in Figure 1.

From the previous discussion it is clear that the above algorithm is such that by measuring the polarization of the ancillary qbit we determine the Wigner function. Indeed, this follows from the identity

$$\sigma _z=4N^2\text{Tr}(\rho \widehat{A}(\alpha _1)\widehat{A}(\alpha _2))=4N^2W(\alpha _1,\alpha _2).$$
(66)
As phase space point operators (6) are simply a product of displacement operators (which implement addition of one, modulo $`N`$) and reflections (which are the square of the Fourier transform) the network of Figure 1 can be implemented efficiently (i.e., it involves a number of elementary gates that grows polynomially with $`\mathrm{log}(N)`$).

## V Conclusion

In this paper we used a hybrid approach to construct a Wigner function to represent quantum states of a composite system in phase space. The function we defined has interesting features enabling us to study situations where entanglement between subsystems plays an important role. Thus, the hybrid method captures some of the most useful properties of the Wigner functions defined by Wooters  and Leonhardt (and others). For a bipartite system this function depends upon two phase space coordinates $`W(\alpha _1,\alpha _2)`$. The phase space is a Cartesian product of the phase spaces of the subsystems, as it is the case in Wotters proposal. However, each phase space grid has $`2N\times 2N`$ points, as suggested by Leonhardt and others. For separable states the Wigner function is, in general, a convex sum of products of independent functions for each subsystem. Thus, this Wigner function is a natural tool to study entanglement between subsystems. In this paper we showed that basis of entangled states can be identified with non–separable slices in the phase space (the basis formed by Bell states is one such example). We also showed that $`W(\alpha ,\alpha ^{})`$ is measurable by a simple scattering–like experiment where an ancillary particle successively interacts with the two subsystems.

###### Acknowledgements.

JPP thanks L. Davidovich and Marcos Saraceno for useful discussions. He also thanks Cecilia Lopez for carefully reading the manuscript. This work was partially supported with grants from Ubacyt, Anpcyt and Fundación Antorchas. JPP is a fellow of CONICET.


format: text

text:
# Periodic orbit quantization of the closed three-disk billiard as an example of a chaotic system with strong pruning

## I Introduction

The relation between quantum spectra, on the one hand, and the dynamics of the corresponding classical systems, on the other, is a problem of fundamental importance in physics, which has attracted attention ever since the early days of quantum mechanics. A key to understanding this relation in chaotic systems was Gutzwiller’s discovery (cf. Gut90 ; Gut01 ) that the semiclassical density of states in these systems can be written as a sum over all (isolated) periodic orbits of the classical system. However, practical applications of Gutzwiller’s trace formula, like other periodic orbit sums, are greatly impeded by the fact that the sums usually diverge in the domain where the physical eigenvalues are located, mainly as a consequence of the rapid proliferation of periodic orbits with growing period. Various techniques were developed over the years to overcome the convergence problem of the periodic orbit sums Cvi89 ; Aur92 ; Ber90 , but most of them depend on special properties of the systems, such as ergodicity or the existence of complete symbolic dynamics.

As an alternative, the technique of harmonic inversion of semiclassical trace formulae has been introduced as a very general method for extracting eigenvalues from periodic orbit sums for both chaotic and regular systems Mai97b ; Mai98 ; Mai99a ; Mai99b ; Wei00 ; Mai00 . Harmonic inversion circumvents the convergence problems, and allows one to obtain the semiclassical eigenvalues from a relatively small number of periodic orbits. Moreover, the method does not require any special properties of the systems, and therefore is applicable also to the particularly challenging case of systems with strong pruning, i.e., systems with a highly incomplete symbolic coding of the periodic orbits. In this paper we demonstrate the power of the harmonic inversion technique in the semiclassical quantization of the closed three-disk billiard system, for which all other semiclassical quantization methods to date have failed in determining more than the few lowest eigenvalues, chiefly because of the extremely rapid proliferation of periodic orbits with increasing action and the strong pruning typical of this system. From a technical point of view our study will show that harmonic inversion is numerically stable even when handling huge (on the order of several million) periodic orbit sets.

The three-disk billiard system consists of three equally spaced hard disks of unit radius. The existence or nonexistence of periodic orbits and the behavior of the periodic orbit parameters in the three-disk system sensitively depend on the distance $`d`$ between the centers of the disks. Large disk separations, especially $`d=6`$, have served as a test case for periodic orbit quantization of chaotic systems in many investigations in recent years. In particular, the system has been employed as prototype example for the usefulness of periodic orbit quantization by cycle expansion techniques Cvi89 ; Eck93 ; Eck95 ; Wir99 . For large disk separations, the assumptions of the cycle expansion, namely that the contributions from long orbits are shadowed by those of short orbits, are almost exactly fulfilled. This is no longer true for small disk separations. As the disks approach each other, the convergence of the cycle expansion becomes slower and slower, until it finally breaks down.

Periodic orbit quantization by harmonic inversion does not depend on specific properties, such as shadowing of orbits etc., and therefore is expected to work well also for small disk separations. In fact, the harmonic inversion technique has successfully been applied to the “standard” literature disk separation $`d=6`$ as well as to the small distance $`d=2.5`$, where the convergence of the cycle expansion is already slow Mai98 ; Mai99a . In this paper we concentrate on the limiting case of mutually touching disks, $`d=2`$, where the system turns into a bound system, and the requirements of the cycle expansion and other semiclassical methods are no longer fulfilled at all.

For the closed three-disk system, the lowest exact quantum eigenvalues were first calculated by Tanner et al. Tan91  and by Scherer Sche91 . A remarkable step towards the semiclassical quantization of the closed three-disk system was achieved by combining the conventional cycle expansion with a functional equation Tan91 . Based on periodic orbit corrections to the mean density of states, approximations to the very lowest eigenvalues could be obtained from a small set of periodic orbits up to cycle length $`l=3`$, where the symbolic dynamics is still complete (including formally a zero length orbit). For the very lowest eigenvalues, the results were found to be in good agreement with the exact quantum values. However, the method failed for higher eigenvalues, where only the mean density of states could be reproduced. Since the problems with the method arose from the strong pruning of orbits in this system, they are of a fundamental nature, and could not be overcome by including more orbits. No other semiclassical method has been successful in obtaining higher eigenvalues in the closed three-disk billiard system yet. In this paper, we will demonstrate that harmonic inversion passes even this acid test of semiclassical quantization.

In Sec. II we discuss the construction of the periodic orbit signal for the closed three-disk system, which is a nontrivial task because of, firstly, the extremely rapid proliferation of periodic orbits and, secondly, the strong pruning of the symbolic dynamics. Once the signal has been constructed, we refer in Sec. III to the harmonic inversion technique for the extraction of the semiclassical eigenvalues. The numerical recipe of the harmonic inversion method has already been established in the literature Mai97b ; Mai98 ; Mai99a ; Mai99b ; Wei00 ; Mai00 , however, to make the article self-consistent we recapitulate the important steps in Sec. III.1. The semiclassical results are discussed in Secs. III.2 and III.3.

## II Periodic orbits of the closed three-disk system

We start by summarizing the basic properties of the three-disk system relevant for periodic orbit quantization. Furthermore, the pruning of orbits and the distribution of the orbit parameters are investigated, and the resulting strategies for the numerical search for periodic orbits are discussed.

### II.1 Symbolic code and symmetry reduction

The periodic orbits of the three-disk system can be labeled by a ternary symbolic code, which is complete for sufficiently large disk separations. The classical dynamics of the three-disk system has been studied, among others, by Gaspard and Rice Gas89 . If the disks are labeled by the numbers 1, 2, 3, each periodic orbit is characterized by a sequence of these numbers, indicating the disks the particle collides with during one period of the orbit. For sufficiently separated disks, there is a one-to-one correspondence between the symbolic code and the periodic orbits of the system: For every sequence there exists one unique periodic orbit (with the restriction that consecutive repetitions of the same symbol are forbidden and circular shifts of a sequence describe the same orbit). As was first pointed out by Hansen Han93 , at disk separation $`d=2.04821419`$ “pruning” sets in, i.e., part of the periodic orbits become unphysical as they begin to run through one of the disks, and thus the symbolic code is no longer complete. The number of pruned orbits rapidly becomes larger and larger if the disks continue to approach each other. In the limiting case of touching disks, $`d=2`$, the system exhibits strong pruning. Examples of pruned orbits will be discussed in Section II.3.

In our calculations, we make use of the symmetry reduction introduced by Cvitanović and Eckhardt in Ref. Cvi89 : The three-disk system is invariant under the symmetry operations of the group $`C_{3v}`$, i.e., reflections at three symmetry lines and rotations by $`2\pi /3`$ and $`4\pi /3`$. The periodic orbits fall into three classes of distinct symmetry: orbits invariant under reflections at one of the symmetry lines (multiplicity 3), orbits invariant under rotations by $`2\pi /3`$ and $`4\pi /3`$ (multiplicity 2), and orbits with no symmetry (multiplicity 6). The quantum states are grouped in the three irreducible subspaces $`A_1`$, $`A_2`$, and $`E`$, where the states of the $`A_1`$ ($`A_2`$) subspace are symmetric (antisymmetric) under reflection at the symmetry lines, respectively, and the states of the $`E`$ subspace are invariant under rotations by $`2\pi /3`$ and $`4\pi /3`$.

Following Refs. Cvi89 ; Eck95 , one can map the system onto a fundamental domain, which consists of a one-sixth slice of the full system, with the symmetry axes acting as straight mirror walls (see Fig. 1).

The periodic orbits of the full system can be described completely in terms of the periodic orbits in the fundamental domain. The symmetry reduced periodic orbits are labeled by a binary code, where the symbol ‘0’ represents backscattering (or change between clockwise and anti-clockwise scattering) and the symbol ‘1’ stands for scattering to the third disk in the original full domain picture. For example, the shortest full domain orbit 12 maps onto the 0 orbit, and the 123 orbit maps onto the 1 orbit.

### II.2 Semiclassical density of states

The quantum resonances of the three different subspaces $`A_1`$, $`A_2`$ and $`E`$ can be obtained separately from the periodic orbits in the fundamental domain by introducing appropriate weight factors for the orbits in Gutzwiller’s trace formula, as was shown by Cvitanović and Eckhardt Cvi93 . We concentrate on the $`A_1`$ subspace, for which each orbit has a weight factor equal to 1.

As for all billiard systems, the shape of the periodic orbits is independent of the wave number $`k=\sqrt{2mE}/\mathrm{}`$, and the action scales as

$$S/\mathrm{}=ks,$$
(1)
where the scaled action $`s`$ is equal to the physical length of the orbit. We consider the density of states as a function of the wave number

$$\rho (k)=\frac{1}{\pi }\mathrm{Im}g(k),$$
(2)
with a scaled response function $`g(k)`$. Applying the Gutzwiller trace formula to the closed three-disk system yields for the $`A_1`$ subspace

$`g(k)=g_0(k)\mathrm{i}{\displaystyle \underset{\mathrm{po}}{}}{\displaystyle \frac{s_{\mathrm{po}}\mathrm{e}^{\mathrm{i}\frac{\pi }{2}\mu _{\mathrm{po}}}}{r|det(M_{\mathrm{po}}1)|^{1/2}}}\mathrm{e}^{\mathrm{i}ks_{\mathrm{po}}}`$ (3)
$`=g_0(k)\mathrm{i}{\displaystyle \underset{\mathrm{po}}{}}(1)^{l_s}{\displaystyle \frac{s_{\mathrm{po}}}{r|(\lambda _{\mathrm{po}}1)(\frac{1}{\lambda _{\mathrm{po}}}1)|^{1/2}}}\mathrm{e}^{\mathrm{i}ks_{\mathrm{po}}},`$
where $`M_{\mathrm{po}}`$ and $`\mu _{\mathrm{po}}`$ are the monodromy matrix and the Maslov index of the orbit, respectively, $`l_s`$ is the symbol length, $`s_{\mathrm{po}}`$ is the scaled action, and $`\lambda _{\mathrm{po}}`$ denotes the expanding stability eigenvalue of the orbit (i.e., the eigenvalue with an absolute value larger than one). The sum runs over all symmetry reduced periodic orbits including multiple traversals. Here, $`r`$ denotes the repetition number with respect to the corresponding primitive orbit.

In practice, only the primitive periodic orbits have to be determined. The parameters of the $`r^{\mathrm{th}}`$ repetition of the primitive orbit (here characterized by the index 0) are then given by $`l_s=rl_{s0}`$, $`s=rs_0`$ and $`\lambda =\lambda _0^r`$.

### II.3 Numerical search for periodic orbits

For extracting the quantum resonances of the system from Eq. (II.2) by harmonic inversion, all periodic orbits up to a maximum scaled action have to be included. The parameters of the periodic orbits – scaled action and stability eigenvalues – have to be determined numerically. We calculate the primitive periodic orbits using the symbolic code as input. For simplicity, the calculations are carried out in the full domain, and the results are then translated back into the symmetry reduced system. The disks are “connected” according to the code, starting with arbitrary reflection points on the disks as initial condition. The reflection points are then varied in such a way that the total length (i.e., the action) of the orbit reaches a minimum. All orbit parameters can then be calculated from the reflection points. The scaled action is given by the length of the (symmetry reduced) orbit, and the stability eigenvalue $`\lambda `$ can be determined by an algorithm proposed by Bogomolny Bog88 . In addition to action and stability, one can also determine averages of different classical quantities (distance from the center of the system, angular momentum, etc.), which are needed for the cross-correlation technique (see Section III).

For disk separations smaller than the pruning limit $`d=2.04821419`$, it has to be checked whether or not the orbits are physical, i.e., whether or not they stay completely outside the disks. In our numerical calculations, we found two different classes of pruned orbits. An example of the first class is given in Figure 2: As the disks approach each other, a section of the orbit connecting two of the disks gets inside the third disk.

The pruned orbit still corresponds to a unique minimum of the total length when the reflection points on the disks are varied for given symbolic code. This is not the case for the second class, an example of which is shown in Figure 3: In this case, as the disks come closer, the reflection angle at one specific reflection point approaches $`\pi `$. As the value $`\pi `$ is reached, the minimum of the action with respect to this specific reflection point splits into two equal minima and one local maximum. Accordingly, the orbit splits into three unphysical orbits with the same symbolic code.

In one case, the orbit is now reflected on the inside of the disk (see Fig. 3c), which corresponds to a reflection angle larger than $`\pi `$ (i.e., the orbit has penetrated into the disk in the vicinity of the reflection point). This orbit does not correspond to a minimum of the total length but to a saddle point. In the other two cases, Figs. 3d and 3e, the reflection point moves along the disk in such a way that the reflection law is no longer fulfilled but the orbit just passes the “reflection point” in a straight line and penetrates into the disk. These two cases correspond to two equally deep minima of the total length. The two orbits have the same shape but only differ as to which penetration point is considered as “reflection point”, indicated by the squares in Figs. 3d and 3e. In fact, the shape of the orbit is equal to that of a pruned orbit with a (symmetry reduced) symbol length that is shorter by 1 and which belongs to the first class of pruned orbits as defined above. (In this orbit the reflection point in question is simply missing, i.e., none of the penetration points is considered as reflection point.)

Besides the existence or nonexistence of periodic orbits, the distance between the disks also strongly influences the distribution of the periodic orbit parameters. This is illustrated in Figure 4, which compares the orbit parameters of the shortest primitive periodic orbits found for the open three-disk system with disk separation $`d=6`$ and for the case of touching disks, $`d=2`$. For $`d=6`$, the action and stability of the orbits are mainly determined by the symbol length, the orbit parameters are nicely bunched in easily distinguishable intervals, and the number of orbits up to a given action is relatively small. However, the picture changes completely when the disks approach each other: As all orbits become shorter, the total number of orbits up to a given action increases rapidly. The parameters of the orbits are no longer simply determined by the symbol length, and their overall behavior becomes much more complicated, as is evident from Figure 4.

For $`d=2`$, the orbits can be grouped in “channels” with the same “tail” (end figures) but growing number of leading ‘0’s in the code. (A sequence of $`n`$ leading ‘0’s in the code will be denoted by $`0^n`$ in the following). These orbits have the same basic shape but run deeper and deeper into the corner formed by two touching disks, bouncing back and forth between the two disks (see Fig. 5).

In each channel, the action of the orbits grows very slowly with increasing symbol length, while the expanding stability eigenvalue increases relatively fast (in fact exponentially) with every additional collision. Since adding a leading ‘0’ to the symbolic code does not change the action considerably, there is a huge number of orbits with very long symbol lengths but relatively small, and comparable, actions. In fact, most channels break off because of pruning. E.g., the $`0^n1`$ series, labeled (A) in Fig. 4b, breaks off after six orbits, and the $`0^n11`$ series, marked (B), already breaks off after three orbits. Our analysis strongly suggests that there exist only two infinite channels: The $`0^n101`$ and the $`0^n111`$ series, dubbed (C) in Fig. 4b (these channels are only shown up to 20 leading ‘0’s in the code). The shapes and actions of these orbits converge very slowly towards the shape and action of the limiting orbit which starts exactly in the touching point of two disks and is directly reflected back by the opposite disk (cf. Fig. 5). Although all other channels break off because of pruning, the number of orbits in a channel up to a given action often becomes very large (or, in the two cases mentioned, even infinite). Note that the 0 orbit, bouncing back and forth between two disks, does not exist in the closed three-disk system.

While for large disk separations finding all orbits up to a given action simply corresponds to calculating all orbits up to a certain maximum symbol length, the complicated distribution of periodic orbit parameters in the closed three-disk system renders the search for the relevant periodic orbits up to a given action a nontrivial task. Due to the extremely rapid proliferation of periodic orbits with increasing length $`s`$, we were not able to determine all orbits up to a given length $`s_{\mathrm{max}}`$, but we had to introduce restrictions for the orbits to be included. As the amplitude in the Gutzwiller formula depends on the inverse of the stability parameter, a reasonable cut-off criterion is the size of the stability of the orbit. On the other hand, since in each channel the absolute value of the expanding stability eigenvalue grows exponentially fast, the main contributions of each channel to the periodic orbit sum come from the orbits with small symbol lengths. Furthermore, the Maslov indices of adjacent orbits in each channel differ by two, and therefore the periodic orbit amplitudes alternate in sign, which leads to approximate cancellations of terms in the semiclassical trace formula. As a second criterion, we therefore restricted the number of consecutive symbols ‘0’ in the symbolic code, thus selecting the most relevant orbits of each channel.

The orbits were calculated channel by channel, starting from the shortest orbit in the channel and adding more and more leading zeros to the code. The calculation of the channel was broken off if one of the following conditions was fulfilled:

* The action exceeds the maximum action $`s_{\mathrm{max}}`$.
* The absolute value of the stability eigenvalue exceeds the given value $`\lambda _{\mathrm{max}}`$.
* The maximum number of consecutive symbols ‘0’ is reached.
* The orbits become pruned.

For the set of periodic orbits shown in Figure 4b, we searched for orbits with physical length $`s<s_{\mathrm{max}}=5.0`$. The maximum for the absolute value of the stability eigenvalue was chosen to be $`\lambda _{\mathrm{max}}=10^9`$. The maximum number of successive zeros in the code was restricted to 20 for orbits with length $`s<3.1`$ and to 12 for orbits with $`3.1<s<5.0`$.

## III Periodic orbit quantization of the closed three-disk system by harmonic inversion

The fundamental problem, and challenge, now is to extract semiclassical eigenvalues in a numerically stable way from the huge periodic orbit set presented in Fig. 4b. To this end, we resort to harmonic inversion of periodic orbit signals Mai97b ; Mai98 ; Mai99a ; Mai99b ; Wei00 ; Mai00 . Here, we review the basic ideas and present an extension of the method introduced in Mai00  to harmonic inversion of cross-correlated periodic orbit sums.

### III.1 The method

The starting point is to introduce a weighted density of states in terms of $`k`$

$$\varrho _{\alpha \alpha ^{}}(k)=\frac{1}{\pi }\mathrm{Im}g_{\alpha \alpha ^{}}(k),$$
(4)
with

$$g_{\alpha \alpha ^{}}^{\mathrm{qm}}(k)=\underset{m}{}\frac{b_{\alpha n}b_{\alpha ^{}n}}{kk_m+\mathrm{i}ϵ},$$
(5)
where $`k_m`$ is the eigenvalue of the wave number of eigenstate $`|m`$ and

$$b_{\alpha m}=m|\widehat{A}_\alpha |m$$
(6)
are the diagonal matrix elements of a chosen set of $`D`$ linearly independent operators $`\widehat{A}_\alpha `$, $`\alpha =1,2,\mathrm{},D`$. The Fourier transform of (5) yields a $`D\times D`$ cross-correlated signal

$`C_{\alpha \alpha ^{}}^{\mathrm{qm}}(s)`$ $`=`$ $`{\displaystyle \frac{\mathrm{i}}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}g_{\alpha \alpha ^{}}(k)\mathrm{e}^{\mathrm{i}sk}dk`$ (7)
$`=`$ $`{\displaystyle \underset{m}{}}b_{\alpha m}b_{\alpha ^{}m}\mathrm{e}^{\mathrm{i}k_ms}.`$
A semiclassical approximation to the weighted quantum response function (5) is given by the semiclassical response function (3) weighted with the classical averages

$$a_{\alpha ,\mathrm{po}}=\frac{1}{s_{\mathrm{po}}}_0^{s_{\mathrm{po}}}A_\alpha (𝐪(s),𝐩(s))ds,$$
(8)
with $`A_\alpha (𝐪,𝐩)`$ the Wigner transform of the operator $`\widehat{A}_\alpha `$ Mai99c ; Hor00 . By Fourier transformation we obtain the semiclassical approximation to the quantum signal (7), which reads

$`C_{\alpha \alpha ^{}}^{\mathrm{sc}}(s)`$ $`=`$ $`{\displaystyle \underset{\mathrm{po}}{}}{\displaystyle \frac{a_{\alpha ,\mathrm{po}}a_{\alpha ^{},\mathrm{po}}s_{\mathrm{po}}\mathrm{e}^{\mathrm{i}\frac{\pi }{2}\mu _{\mathrm{po}}}}{r|det(M_{\mathrm{po}}\mathrm{𝟏})|^{1/2}}}\delta \left(ss_{\mathrm{po}}\right)`$ (9)
$``$ $`{\displaystyle \underset{\mathrm{po}}{}}𝒜_{\alpha \alpha ^{}}^{\mathrm{po}}\delta \left(ss_{\mathrm{po}}\right),`$
where $`r`$ is the repetition number counting the traversals of the primitive orbit, and $`M_{\mathrm{po}}`$ and $`\mu _{\mathrm{po}}`$ are the monodromy matrix and Maslov index of the orbit, respectively. Semiclassical approximations to the eigenvalues $`k_m`$, and eventually also to the diagonal matrix elements $`m|\widehat{A}_\alpha |m`$, are obtained by adjusting the semiclassical cross-correlated periodic orbit signal (9) to the functional form of the quantum signal (7). This can be done with the filter-diagonalization method for the signal processing of cross-correlation functions Nar97 ; Man98 . Here, we resort to an alternative method introduced in Ref. Mai00 , which will now be extended for the harmonic inversion of cross-correlated periodic orbit sums.

The special form of the periodic orbit signal (9) as a sum of $`\delta `$ functions allows for an even simpler procedure than filter-diagonalization, viz. analytical filtering. In the following we will apply a rectangular filter, i.e., $`f(k)=1`$ for frequencies $`k[k_0\mathrm{\Delta }k,k_0+\mathrm{\Delta }k]`$, and $`f(k)=0`$ outside the window. Applying this filter to the semiclassical signal $`C_{\alpha \alpha ^{}}^{\mathrm{sc}}(s)`$ in (9) we obtain the band-limited (bl) cross-correlated periodic orbit signal,

$`C_{\alpha \alpha ^{}}^{\mathrm{sc},\mathrm{bl}}(s)`$ $`=`$ $`{\displaystyle \frac{\mathrm{i}}{2\pi }}{\displaystyle _{k_0\mathrm{\Delta }k}^{k_0+\mathrm{\Delta }k}}g_{\alpha \alpha ^{}}^{\mathrm{sc}}(k)\mathrm{e}^{\mathrm{i}s(kk_0)}dk`$ (10)
$`=`$ $`{\displaystyle \frac{\mathrm{i}}{2\pi }}{\displaystyle \underset{\mathrm{po}}{}}𝒜_{\alpha \alpha ^{}}^{\mathrm{po}}{\displaystyle _{k_0\mathrm{\Delta }k}^{k_0+\mathrm{\Delta }k}}\mathrm{e}^{\mathrm{i}sk_0\mathrm{i}(ss_{\mathrm{po}})k}dk`$
$`=`$ $`\mathrm{i}{\displaystyle \underset{\mathrm{po}}{}}𝒜_{\alpha \alpha ^{}}^{\mathrm{po}}{\displaystyle \frac{\mathrm{sin}[(ss_{\mathrm{po}})\mathrm{\Delta }k]}{\pi (ss_{\mathrm{po}})}}\mathrm{e}^{\mathrm{i}s_{\mathrm{po}}k_0}.`$
The introduction of $`k_0`$ into the arguments of the exponential functions in (10) causes a shift of frequencies by $`k_0`$ in the frequency domain. Note that $`C_{\alpha \alpha ^{}}^{\mathrm{sc},\mathrm{bl}}(s)`$ is a smooth function and can be easily evaluated on an arbitrary grid of points $`s_m<s_{\mathrm{max}}`$ provided the periodic orbit data are known for the set of orbits with classical action $`s_{\mathrm{po}}<s_{\mathrm{max}}`$.

Applying now the same filter as used for the semiclassical periodic orbit signal to the quantum one, we obtain the band-limited quantum signal

$`C_{\alpha \alpha ^{}}^{\mathrm{qm},\mathrm{bl}}(s)={\displaystyle \frac{\mathrm{i}}{2\pi }}{\displaystyle _{k_0\mathrm{\Delta }k}^{k_0+\mathrm{\Delta }k}}g_{\alpha \alpha ^{}}^{\mathrm{qm}}(k)\mathrm{e}^{\mathrm{i}s(ikk_0)}dk`$ (11)
$`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}b_{\alpha m}b_{\alpha ^{}m}\mathrm{e}^{\mathrm{i}(k_mk_0)s},|k_mk_0|<\mathrm{\Delta }k.`$
In contrast to the signal $`C_{\alpha \alpha ^{}}^{\mathrm{qm}}(s)`$ in Eq. (7), the band-limited quantum signal consists of a finite number of frequencies $`k_m`$, $`m=1,\mathrm{},M`$, where in practical applications $`M`$ can be of the order of $``$ (50-200) for an appropriately chosen frequency window, $`\mathrm{\Delta }k`$. The problem of adjusting the band-limited semiclassical signal in Eq. (10) to its quantum mechanical analogue in Eq. (11) can now be written as a set of nonlinear equations

$`C_{\alpha \alpha ^{}}^{\mathrm{sc},\mathrm{bl}}(n\tau )c_{\alpha \alpha ^{},n}={\displaystyle \underset{m=1}{\overset{M}{}}}b_{\alpha m}b_{\alpha ^{}m}\mathrm{e}^{\mathrm{i}k_m^{}n\tau },`$ (12)
$`n=0,1,\mathrm{},2N1,`$
for the unknown variables, viz. the shifted frequencies, $`k_m^{}k_mk_0`$, and the parameters $`b_{\alpha m}`$, $`m=1,\mathrm{},M`$. Note that $`M`$ is related to the signal length $`2N`$ and the dimension $`D`$ of the cross-correlation matrix by $`M=ND`$. The signal (12) can be evaluated on an equidistant grid, $`s=n\tau `$, with relatively large step size $`\tau \pi /\mathrm{\Delta }w`$, and the discrete signal points $`c_{\alpha \alpha ^{},n}C_{\alpha \alpha ^{}}^{\mathrm{sc},\mathrm{bl}}(n\tau )`$ are called the “band-limited” cross-correlated periodic orbit signal.

We now wish to solve the nonlinear system, Eq. (12), which can be written as

$$c_{\alpha \alpha ^{},n}=\underset{m=1}{\overset{M}{}}b_{\alpha m}b_{\alpha ^{}m}z_m^n,$$
(13)
where $`z_m\mathrm{exp}(\mathrm{i}k_m^{}\tau )`$ and $`b_{\alpha m}`$ are, generally complex, variational parameters. We assume that the number of frequencies in the signal is relatively small ($`M50`$ to $`200`$). Although the system of nonlinear equations is, in general, still ill-conditioned, frequency filtering reduces the number of signal points, and hence the number of equations. In Ref. Mai00  three different methods, viz. linear predictor (LP), Padé approximant (PA), and signal diagonalization (SD) have been employed to solve Eq. (13) for a one-dimensional signal, i.e., $`D=1`$. The SD method can be generalized in a straightforward manner to the harmonic inversion of cross-correlated periodic orbit sums. The problem of solving the nonlinear set of equations (13) can be recast in the form of the generalized eigenvalue problem Wal95 ; Man97 ; Man98

$$𝐔𝑩_m=z_m𝐒𝑩_m,$$
(14)
where the elements of the $`M\times M`$ operator matrix $`𝐔`$ and overlap matrix $`𝐒`$ depend trivially upon the signal points:

$`U_{\alpha i,\alpha ^{}j}`$ $`=`$ $`c_{\alpha \alpha ^{},i+j+1};`$
$`S_{\alpha i,\alpha ^{}j}`$ $`=`$ $`c_{\alpha \alpha ^{},i+j};i,j=0,\mathrm{},N1.`$ (15)
The matrices $`𝐔`$ and $`𝐒`$ in Eq. (14) are complex symmetric (i.e., non-Hermitian), and the eigenvectors $`𝑩_m`$ are orthogonal with respect to the overlap matrix $`𝐒`$,

$$\left(𝑩_m|𝐒|𝑩_m^{}\right)=N_m\delta _{mm^{}},$$
(16)
where the brackets define a complex symmetric inner product $`(a|b)=(b|a)`$, i.e., no complex conjugation of either $`a`$ or $`b.`$ The overlap matrix $`𝐒`$ is not usually positive definite and therefore the $`N_m`$’s are, in general complex, normalization parameters. An eigenvector $`𝑩_m`$ cannot be normalized for $`N_m=0`$. The parameters $`b_{\alpha m}`$ in Eq. (13) are obtained from the eigenvectors $`𝑩_m`$ via

$$b_{\alpha m}=\frac{1}{\sqrt{N_m}}\underset{\alpha ^{}=1}{\overset{D}{}}\underset{n=0}{\overset{N1}{}}c_{\alpha \alpha ^{},n}𝑩_{m,\alpha ^{}n}.$$
(17)
The parameters $`z_m`$ in Eq. (13) are given as the eigenvalues of the generalized eigenvalue problem (14), and are simply related to the frequencies $`k_m^{}`$ in Eq. (12) via $`z_m=\mathrm{exp}(\mathrm{i}k_m^{}\tau )`$.

The resolution of the results depends on the signal length. For a one-dimensional signal, the method requires a signal length of $`s_{\mathrm{max}}4\pi \overline{\rho }(k)`$, with $`\overline{\rho }(k)`$ the mean density of states, to resolve the frequencies Mai97b . This means that all periodic orbits up to the scaled action $`s_{\mathrm{max}}`$ have to be included. The advantage of using the cross-correlation approach is based on the insight that the total amount of independent information contained in the $`D\times D`$ signal is $`D(D+1)`$ multiplied by the length of the signal, while the total number of unknowns (here $`b_{\alpha m}`$ and $`k_m`$) is $`(D+1)`$ times the total number of poles $`k_m`$. Therefore the informational content of the $`D\times D`$ signal per unknown parameter is increased (compared to the one-dimensional signal) by roughly a factor of $`D`$, and the cross-correlation approach should lead to a significant improvement of the resolution. The power of this method for periodic orbit quantization has already been proven for the example of an integrable system Mai99b ; Wei00 . We will now demonstrate how the method works in the case of three touching disks, as an example of a chaotic system with extremely rapid proliferation of orbits and strong pruning.

### III.2 Semiclassical eigenvalues of the closed three-disk system by harmonic inversion of a one-dimensional signal

We first investigate a one-dimensional signal ($`N=1`$) by simply choosing $`\widehat{A}_1=\mathrm{𝟏}`$, i.e., the unity operator. The signal is constructed from the set of periodic orbits shown in Fig. 4b, which contains about 5 million primitive orbits.

We analyzed the signal up to length $`s_{\mathrm{max}}=4.9`$ using the decimated signal diagonalization method Mai00 . The frequencies obtained by harmonic inversion of the signal are the semiclassical approximations to the eigenvalues $`k_n`$. Fig. 6 shows the results from harmonic inversion (solid lines), compared with the exact quantum results (dashed lines), both presented in terms of the energy $`E=(\mathrm{Re}k)^2/2`$. The exact quantum eigenvalues were taken from Scherer Sche91  and Wirzba Wir . The harmonic inversion results clearly reproduce the quantum eigenvalues up to $`E4500`$ with only small deviations. However, the signal length was not sufficient to resolve the eigenvalues in the region $`E>4500`$, where the density of states with respect to the wave number $`k`$ (which grows $`\sqrt{E}`$) becomes too large.

Apart from the values shown in Fig. 6, the harmonic inversion of the signal also yielded a number of unconverged frequencies. Converged frequencies were identified by having an imaginary part close to zero (since the eigenvalues of the wave number must be real for bound systems) and an amplitude close to the theoretical value $`m_k=1`$. The whole set of frequencies obtained is shown in Figure 7 together with the real part of their amplitudes. The solid vertical lines mark the positions of the exact quantum eigenvalues. The semiclassical values included in Fig. 6 are represented by filled circles in Fig. 7. The amplitudes are not as well converged as the frequencies (i.e. the eigenvalues) but partly show larger deviations from the theoretical value $`m_k=1`$. \[This can be explained by the observation that when the decimated signal diagonalization method is used for harmonic inversion the frequencies usually converge faster than the amplitudes.\] However, up to $`E4500`$, the eigenvalues can be clearly identified.

The small deviations of the eigenvalues obtained by harmonic inversion from the exact quantum values may have different reasons. Especially for the lowest eigenvalues, the major part is probably due to the semiclassical error. On the other hand, as explained above, we did not include all orbits in the signal but left out a large number of relatively unstable orbits. Although each of these orbits only gives a negligible small contribution to the periodic orbit sum, the number of excluded orbits may be so large that their contributions sum up in a way so as to have a visible effect on the density of states. Finally, even in the low-lying part of the spectrum, deviations may arise from the fact that the signal length was very short. The influence of the missing orbits and the relatively short signal length are reflected in the relatively poor convergence of the amplitudes.

In order to test the stability of the results with respect to the signal parameters, and to obtain an estimate of random errors, we performed the same calculation with various sets of different signal parameters. Figure 8 shows the harmonic inversion results for the eigenvalues as a function of the signal length $`s_{\mathrm{max}}`$. From all frequencies obtained the values with $`|\mathrm{Im}k|<1.0`$ and with amplitudes $`\mathrm{Re}m_k>0.2`$ were singled out, respectively. The filled circles mark amplitudes with $`\mathrm{Re}m_k>0.5`$. It can clearly be seen how the maximum eigenvalue up to which the spectrum can be resolved depends on the signal length. Fig. 8 shows that the results for the low-lying eigenvalues are quite stable with respect to variation of the signal length.

Compared with the results obtained by Tanner et al. Tan91  using the extended cycle expansion method, the analysis of a single signal by harmonic inversion has yielded about the same number of eigenvalues for the closed three-disk system. For the resolution of higher eigenvalues, the signal length was not sufficient. Extending the signal to significantly larger scaled actions is in practice not possible because of the extremely rapid increase of the number of orbits. It should be stressed, however, that, in contrast to the cycle expansion method, which runs into severe problems because of the pruning of orbits in this system, it is at least in principle possible to improve the harmonic inversion results by including more orbits. On the other hand, as discussed before, the harmonic inversion method offers the possibility to significantly reduce the signal length required for the resolution of eigenvalues by the construction and analysis of cross-correlated periodic orbit sums.

### III.3 Improvement of the resolution by harmonic inversion of cross-correlated periodic orbit sums

To increase the accuracy of the semiclassical eigenvalues, we now apply the cross-correlation technique to the closed three-disk system. Instead of the one-dimensional signal investigated in the previous section, we choose a set of three operators and construct a $`3\times 3`$ cross-correlated signal according to Eq. (9). As input, we again use the set of periodic orbits presented in Fig. 4 together with the averages of different classical quantities over the orbits.

Fig. 9 shows the results from the harmonic inversion of a $`3\times 3`$ signal of length $`s_{\mathrm{max}}=4.8`$ (solid lines), which was constructed using the operators $`\mathrm{𝟏}`$ (unity), $`r^4`$ and $`L^4`$, where $`r`$ and $`L`$ denote the distance from the center of the system and the absolute value of the angular momentum, respectively.

The eigenvalues are again presented in terms of the energy $`E=(\mathrm{Re}k)^2/2`$. For comparison, the dashed lines indicate the positions of the exact quantum eigenvalues. The converged eigenvalues for the wave number $`k`$ have been singled out from the whole set of frequencies obtained by the condition that they should have an imaginary part close to zero and an amplitude close to the theoretical value $`m_k=1`$. The amplitudes and the imaginary parts of all frequencies obtained are shown in Figure 10 by circles. In particular, the filled circles denote the frequencies included in Fig. 9. In Fig. 10a, the positions of the exact quantum eigenvalues are given by the solid vertical lines.

With the cross-correlated signal, one can now clearly identify eigenvalues up to $`E6500`$. In addition, the convergence of the lowest eigenvalues is improved as compared to the results from the single signal obtained in Section III.2, as can be seen from the amplitudes. It is not possible to determine to what extent the accuracy of the semiclassical eigenvalues has improved since the results can only be compared with the exact quantum eigenvalues and the size of the semiclassical error is unknown.

Again, we have performed the same calculation for various sets of parameters in the harmonic inversion scheme. Figure 11 shows the results from a $`3\times 3`$ and from a $`4\times 4`$ signal as functions of the signal length. The $`3\times 3`$ signal was constructed from the operators $`\mathrm{𝟏}`$ (unity), $`r^4`$ and $`L^4`$, and the $`4\times 4`$ signal contained the operators $`\mathrm{𝟏}`$, $`r^2`$, $`r^4`$ and $`L^4`$. In both diagrams, only the frequencies with an imaginary part $`|\mathrm{Im}k|<1.0`$ and an amplitude $`\mathrm{Re}m_k>0.5`$ were included. The filled circles mark values with $`|\mathrm{Im}k|<0.5`$. The positions of the exact quantum eigenvalues are again marked by the solid vertical lines.

As with the one-dimensional signal analyzed in the previous section, it can be observed that the results for the lowest eigenvalues are very stable with respect to the variation of the signal length as long as the signal is not too short. However, the very lowest frequencies show a tendency to split into two or more, which was a frequent observation in our calculations when the signal length and the matrix dimensions were chosen too large. The results for the higher eigenvalues depend more sensitively on the signal length as the density of states approaches the limiting resolution that can be achieved with the signal. In this region, the $`4\times 4`$ signal shows a better resolution than the $`3\times 3`$ signal.

With the $`3\times 3`$ signal and the two largest values of the signal length considered, the lowest eigenvalue was not obtained. The reason for this probably lies in inaccuracies at the end of the signal, as the signal length approaches the value $`s_{\mathrm{max}}=5.0`$: In the decimated signal diagonalization scheme, which we used for the harmonic inversion of the signal, the $`\delta `$ functions in the original signal (9) are replaced with sinc functions Mai00  \[$`\mathrm{sinc}x=(\mathrm{sin}x)/x`$\]. Therefore, every periodic orbit contributes to the signal not only exactly at its scaled action $`s=s_{\mathrm{po}}`$, but also in a range $`\mathrm{\Delta }s`$ around this point. If the signal length is close to $`s_{\mathrm{max}}=5.0`$, also orbits with actions slightly larger than this limit may give non-negligible contributions to the very end of the signal. Since we did only include orbits with $`s5.0`$, these contributions are missing in our signal.

Some eigenvalues between $`E=3500`$ and $`E=5000`$ seem to be particularly hard to obtain from the set of orbits used. Possibly, these eigenvalues are related to the orbits running very deep into the corners between the disks, which were not included in the signal. We conjecture that the missing orbits as well as the short signal length are again responsible for the deviations especially of the results for higher eigenvalues from the exact quantum values.

## IV Conclusion

We have successfully applied the technique of harmonic inversion to the semiclassical quantization of a system with an extremely rapid proliferation of periodic orbits and strong pruning. With the harmonic inversion of cross-correlated periodic orbit sums we were able to calculate accurate semiclassical energy eigenvalues of the closed three-disk system up to the region $`E6500`$ (a total of 29 levels). This is the first time that such high-lying eigenvalues have been correctly calculated by periodic orbit theory in this system. Our results demonstrate that the harmonic inversion of cross-correlated periodic orbit sums is indeed a powerful method of periodic orbit quantization and allows the stable handling even of huge periodic orbit sets. In contrast to other methods, which depend on the existence of a complete symbolic code, the harmonic inversion technique is not affected by the strong pruning of orbits. In the closed three-disk system, the only restriction for the practical application of the harmonic inversion method lies in the extremely large number of periodic orbits, which is a special feature of this system. It should also be noted that the semiclassical eigenvalues were obtained without resorting to the mean density of states, which, including all necessary correction terms, is nontrivial to calculate even in the case of the closed three-disk system Sche91 . The present results may stimulate future work on other challenging systems without special structural information on the quantum or classical level. One example with strong pruning of the symbolic dynamics is the three-dimensional generalization of the closed three-disk billiard, i.e., a system consisting of four touching spheres at the corners of a regular tetrahedron. By contrast to the closed three-disk system this is an open system where no functional equation can be applied. Experiments on chaotic light-scattering from the four spheres have recently attracted much attention Swe99 .

###### Acknowledgements.

We thank A. Wirzba for communicating to us quantum mechanical data for the closed three-disk system. This work was supported by Deutsche Forschungsgemeinschaft and Deutscher Akademischer Austauschdienst.