Numerical study of quantum resonances in chaotic scattering.
This paper presents numerical evidence that in quantum systems with chaotic classical dynamics, the number of scattering resonances near an energy E scales like h-(D(KE)+1)/2 as h  0. Here, KE denotes the subset of the energy surface bounded for all time under the flow generated by the classical Hamiltonian H and D (KE) denotes its fractal dimension. Since the number of bound states in a quantum system with n degrees of freedom scales like h-n, this suggests that the quantity (D (KE) represents the effective number of degrees of freedom in chaotic scattering problems. The calculations were performed using a recursive refinement technique for estimating the dimension of fractal repellors in classical Hamiltonian scattering, in conjunction with tools from modern quantum chemistry and numerical linear algebra.
Introduction
Quantum mechanics identifies the energies of stationary
states in an isolated physical system with
the eigenvalues of its Hamiltonian operator. Because
of this, eigenvalues play a central role in
the study of bound states, such as those describing
the electronic structures of atoms and molecules. 1
When the corresponding classical system allows
escape to infinity, resonances replace eigenvalues
as fundamental quantities: The presence of a resonance
at
, with E real and
> 0, gives
rise to a dissipative metastable state with energy E
and decay rate
, as described in [37]. Such states
are essential in scattering theory. 2
Department of Mathematics, University of California,
Berkeley, CA 94720. E-mail: kkylin@math.berkeley.edu.
The author is supported by a fellowship from the Fannie and
John Hertz Foundation.
1 For examples, see [5].
2 Systems which are not effectively isolated but interact only
weakly with their environmentcan also exhibit resonant behav-
ior. For example, electronic states of an "isolated" hydrogen
atom are eigenfunctions of a self-adjoint operator, but coupling
the electron to the radiation field turns those eigenstates into
metastable states with finite lifetimes. This paper does not deal
with dissipative systems and is only concerned with scattering.
An important property of eigenvalues is that one
can count them using only the classical Hamiltonian
function
Planck's constant
the number N eig of eigenvalues in
where n denotes the number of degrees of freedom
and vol () phase space volume. This re-
sult, known as the Weyl law, expresses the density
of quantum states using the classical Hamiltonian
generalization to resonances
is currently known.
In this paper, numerical evidence for a Weyl-
like power law is presented for resonances in a
two-dimensional model with three symmetrically-
placed gaussian potentials. A conjecture, based on
the work of Sjostrand [27] and Zworski [35], states
that the number of resonances
with
asymptotically
lies between C 1 h
flow generated by H:
(2)
If D (KE ) depends continuously on E and
sufficiently small, then D (KE1 )
and the number of resonances in such a
region is comparable to
for any E 2 [E
The sets K and KE are trapped sets and consist
of initial conditions which generate trajectories
that stay bounded forever. In systems where
3 For a beautiful exposition of early work on this and related
themes, see [16]. For recent work in the semiclassical context,
see [7].
fH  Eg is bounded for all E, the conjecture reduces
to the Weyl asymptotic h n . Note that the
conjecture implies that
effective number of degrees of
freedom for metastable states =
for both quantum and classical chaotic scattering.
The notion of dimension requires some com-
ment: The "triple gaussian" model considered
here has very few trapped trajectories, and K and
KE (for any energy E) have vanishing Lebesgue
measures. Thus, D(K) is strictly less than
and D (KE ) < In fact, the sets K and
KE are fractal, as are trapped sets in many other
chaotic scattering problems. Also, in this paper,
the term "chaotic" always means hyperbolic; see
Sjostrand [27] or Gaspard [12] for definitions.
This paper is organized as follows: First, the
model system is defined. This is followed by
mathematical background information, as well as
a heuristic argument for the conjecture. Then,
numerical methods for computing resonances and
fractal dimensions are developed, and numerical
results are presented and compared with known
theoretical predictions.
Notation. In this paper, H denotes the Hamiltonian
function
H the corresponding
Hamiltonian operator h 2
is the usual Laplacian and
acts by multiplication.
Triple Gaussian Model
The model system has degrees of freedom;
its phase space is R 4 , whose points are denoted by
First, it is convenient to define
G
Similarly, put
G
G
in two dimensions.
define H by
y
-22
Figure

1: Triple gaussian potential
where the potential Vm is given by
(R cos 2k
sin 2k
That is, it consists of m gaussian "bumps" placed
at the vertices of a regular m-gon centered at the
origin, at a distance R > 0 from the origin. This
paper focuses on the case because it is
the simplest case that exhibits nontrivial dynamics
in two dimensions. However, the case
also relevant because it is well-understood: See
Miller [21] for early heuristic results and Gerard
and Sjostrand [13] for a rigorous treatment. Thus,
double gaussian scattering serves as a useful test
case for the techniques described here.
The quantized Hamiltonian b
H is similarly defined


Figure

1.
3 Background
This section provides a general discussion of resonances
and motivates the conjecture in the context
of the triple gaussian model. However, the notation
reflects the fact that most of the definitions
and arguments here carry over to more general
systems with n degrees of freedom. The reader
should keep in mind that for the triple gaussian
model.
There exists an extensive literature on resonances
and semiclassical asymptotics in other set-
tings. For example, see [9, 10, 11, 34] for detailed
studies of the classical and quantum mechanics of
hard disc scattering.
3.1 Resonances
Resonances can be defined mathematically as fol-
lows:
real z, where I
is the identity operator. This one-parameter family
of operators R(z) is the resolvent and is meromorphic
with suitable modifications of its domain
and range. The poles of its continuation into the
complex plane are, by definition, the resonances
of b
H . 4
Less abstractly, resonances are generalized
eigenvalues of b
H. Thus, we should solve the time-independent
Schrodinger equation
to obtain the resonance  and its generalized
eigenfunction . In bound state computations,
one approximates as a finite linear combination
of basis functions and solves a finite-dimensional
version of the equation above. To carry out similar
calculations for resonances, it is necessary that
lie in a function space which facilitates such ap-
proximations, for example L 2 .
Let and  solve (10). Then e i
the time-dependent Schrodinger equation
ih
@
@t
It follows that Im () must be negative because
metastable states decay in time. Now suppose,
for simplicity, that solutions of
(10) with energy E behave like e i
2Ex for
large x > 0. Substituting
yields e i
2x , which grows exponentially because
Im
< 0. Thus, finite rank approximations
of b
H cannot capture such generalized
eigenfunctions. However, if we make the formal
substitution x 7! xe i , then the wave function
becomes exp
. Choosing
=E) forces to decay exponentially

This procedure, called complex scaling, transforms
the Hamiltonian operator b
H into the scaled
operator b
H . It also maps metastable states
with decay rate
tan (2) to genuine L 2
4 For more details and some references, see [37].
5 The analysis in higher dimensions requires some care, but
the essential result is the same.

Figure

2: Illustration of complex scaling: The
three lines indicate the location of the rotated continuous
spectrum for different values of , while
the box at the top of the figure is the region in
which resonances are counted. Eigenvalues which
belong to different values of  are marked with
different styles of points. As explained later, only
eigenvalues near the region of interest are com-
puted. This results in a seemingly empty plot.
eigenfunctions  of b
H . The corresponding resonance
becomes a genuine eigenvalue: b
. Furthermore, resonances of b
H will be invariant
under small perturbations in , whereas
other eigenvalues of b
H will not. The condition
=E) implies that, for small
and fixed E, the method will capture a resonance
if and only if
. We
can perform complex scaling in higher dimensions
by substituting r 7! re i in polar coordinates.
In algorithmic terms, this means we can compute
eigenvalues of b
H for a few different values
of  and look for invariant values, as demonstrated
in Figure 2. In addition to its accuracy and
flexibility, this is one of the advantages of complex
scaling: The invariance of resonances under
perturbations in  provides an easy way to check
the accuracy of calculations, mitigating some of
the uncertainties inherent in computational work. 6
Note that the scaled operator b
H is no longer self-
adjoint, which results in non-hermitian finite-rank
approximations and complex eigenvalues.
This method, first introduced for theoretical
purposes by Aguilar and Combes [1] and Balslev
6 For a different approach to computing resonances, see [33]
and the references there.
and Combes [3], was further developed by B. Simon
in [26]. It has since become one of the main
tools for computing resonances in physical chemistry
[22, 31, 32, 24]. For recent mathematical
progress, see [18, 27, 28] and references therein.
For reference, the scaled triple-gaussian operator
H is
where
Note that these expressions only make sense because
G
x0 (x) is analytic in x, x 0 , and .
3.2 Fractal Dimension
Recall that the Minkowski dimension of a given set
U
where U fy g. A simple
calculation yields
log (vol (U
log (1=)
if the limit exists.
Texts on the theory of dimensions typically
begin with the Hausdorff dimension because it
has many desirable properties. In contrast, the
Minkowski dimension can be somewhat awk-
ward: For example, a countable union of
zero-dimensional sets (points) can have positive
Minkowski dimension. But, the Minkowski dimension
is sometimes easier to manipulate and almost
always easier to compute. It also arises in the
heuristic argument given below.
For a detailed treatment of different definitions
of dimension and their applications in the study of
dynamical systems, see [8, 23].
3.3 Generalizing the Weyl Law
The formula
makes no sense in scattering problems because the
volume on the right hand side is infinite for most
choices of E 0 and E 1 , and this seems to mean that
there is no generalization of the Weyl law in the
setting of scattering theory. However, the following
heuristic argument suggests otherwise:
As mentioned before, a metastable state corresponding
to a resonance
has
a time-dependent factor of the form e i
wave packet whose dynamics
is dominated by  (and other resonances near it)
would therefore exhibit temporal oscillations of
frequency O(E=h) and lifetime O(h=
tically, then, the number of times the particle
"bounces" in the "trapping region" 7 before escaping
should be comparable to E
.
In the semiclassical limit, the dynamics of the
wave packet should be well-approximated by a
classical trajectory. Let T (x; the
time for the particle to escape the system starting
at position (x; y) with momentum (p x ; p y ). The
diameter of the trapping region is O(R), and typical
velocities in the energy surface are
O(
E) (mass set to unity), so the number of times
a classical particle bounces before escaping should
be O(T
E=R). This suggests that, in the limit
E=R  E=
and consequently
T  R
Fix
for fixed energies E 0 and Equation (17) implies
that T  R
analogy with the
Weyl law,
vol
follows as an approximation for the number of
quantum states with the specified energies and decay
rates.
Now, the function 1=T is nonnegative for all
and vanishes on K [E0 ;E1
g. Assuming that 1=T is suffi-
7 For our triple gaussian system, that would be the triangular
region bounded by the gaussian bumps.
ciently regular, 8 this suggests
1=T (x;
where
denotes distance to K [E0 ;E1 ] . It
follows that N res should scale like
vol
2o
For small
0 , this becomes
for some constant C, by (14). Choosing
and assuming that D (KE ) decreases monotonically
with increasing E (as is the case in Figure
22), we obtain
res  C 0 h
If sufficiently small, then
D (KE
and
res  h
In [27], Sjostrand proved the following rigorous
upper bound: For
1=C,

holds for all  > 0. When the trapped set is of pure
dimension, that is when the infimum in Equation
is achieved, one can take Setting
h gives an upper bound of the form (24).
In his proof, Sjostrand used the semiclassical argument
above with escape functions and the Weyl
inequality for singular values. Zworski continued
this work in [35], where he proved a similar result
for scattering on convex co-compact hyperbolic
surfaces with no cusps. His work was motivated
by the availability of a large class of examples
with hyperbolic flows, easily computable di-
mensions, and the hope that the Selberg trace formula
could help obtain lower bounds. But, these
hopes remain unfulfilled so far [14], and that partly
motivates this work.
8 In fact, this is numerically self-consistent: Assume that
1=T vanishes to order  (with  not necessarily equal to 2) on
K, and assume the conjecture. Then the number of resonances
would scale like h (2n D(K))= , from which one can solve
for . With the numerical data we have, this indeed turns out
to be 2 (but with significant fluctuations).
Also, if 1=T does not vanish quadratically everywhere on
K, variations in its regularity may affect the correspondence
between classical trapping and the distribution of resonances.
Computing Resonances
Complex scaling reduces the problem of calculating
resonances to one of computing eigenvalues.
What remains is to approximate the operator b
by a rank N operator b
HN; and to develop appropriate
numerical methods. For comparison, see
[22, 31, 32, 24] for applications of complex scaling
to problems in physical chemistry.
4.1 Choice of Scaling Angle.
One important consideration in resonance computation
is the choice of the scaling angle : Since
we are interested in counting resonances in a box
it is necessary to choose
tan 1
so that the continuous spectrum of
H is shifted out of the box [E

Figure

2).
In fact, the resonance calculation uses
This choice of  helps avoid the pseudospectrum
[30, 36]:
Let A be an N  N matrix, and let R(z) be
the resolvent It is well known that
when A is normal, that is when A commutes with
its adjoint A  , the spectral theorem applies and the
inequality
holds ((A) denotes the spectrum of A). When
A is not normal, no such inequality holds and
jjR(z)jj can become very large for z far from
(A). This leads one to define -pseudospectrum:
Using the fact that A is a matrix, one can show that
is equal to the set
That is, the -pseudospectrum of A consists of
those complex numbers z which are eigenvalues
of an -perturbation of A.
The idea of pseudospectrum can be extended to
general linear operators. In [30], it is emphasized
that for non-normal operators, the pseudospectrum
can create "false eigenvalues" which make the accurate
numerical computation of eigenvalues diffi-
cult. In [36], this phenomenon is explained using
semiclassical asymptotics. Roughly speaking, the
pseudospectrum of the scaled operator b
H is given
by the closure of
of its symbol H , which is the scaled Hamiltonian
function
in this case. Choosing  to be comparable to h ensures
that the imaginary part of H is also comparable
to h, which keeps the pseudospectrum away
from the counting box [E larger
would contribute a larger  2 term to the imaginary
part of H and enlarge the pseudospectrum.
As one can see in Figures 31 - 34, the invariance
of resonances under perturbations in  also helps
filter out pseudospectral effects.
This consideration also points out the necessity
of the choice
To avoid pseudospectral
effects,  must be O (h). On the other
hand, if approximations
may fail to capture resonances in the region
of interest.
4.2 Eigenvalue Computation
Suppose that we have constructed b
HN; . In
the case of eigenvalues, the Weyl law states that
as h ! 0, since our system has
degrees of freedom. Thus, in order to capture
a sufficient number of eigenvalues, the rank N
of the matrix approximation must scale like h 2 .
In the absence of more detailed information on the
density of resonances, the resonance computation
requires a similar assumption to ensure sufficient
numerical resolution.
Thus, for moderately small h, the matrix has
rapidly becomes prohibitive
on most computers available today. Fur-
thermore, even if one does not store the entire ma-
numerical packages like LAPACK [2] require
auxiliary storage, again making practical
calculations impossible.
Instead of solving the eigenvalue problem
directly, one solves the equivalent
eigenvalue problem
Efficient implementations of the Arnoldi algorithm
[19] can solve for the largest few eigenvalues
0 of
. But
this method allows one to compute a subset of the
spectrum of b
HN; near a given  0 .
Such algorithms require a method for applying
the matrix
to a given vector v
at each iteration step. In the resonance computa-
tion, this is done by solving ( b
for w by applying conjugate gradient to the normal
equations (see [4]). 9 The resonance program,
therefore, consists of two nested iterative methods:
An outer Arnoldi loop and an inner iterative linear
solver for ( b
v. This computation
uses which provides a flexible and efficient
implementation of the Arnoldi method.
To compute resonances near a given energy E,
the program uses  instead of
This helps control the condition number
of b
HN;  0 and gives better error estimates and
convergence criteria. 11
4.3 Matrix Representations
4.3.1 Choice of Basis
While one can discretize the differential operator
H via finite differences, in practice it is better to
represent the operator using a basis for a subspace
of the Hilbert space should better represent
the properties of wave functions near infinity
and obtain smaller (but more dense) matrices.
Common basis choices in the chemical literature
include so-called "phase space gaussian" [6]
and "distributed gaussian" bases [15]. These bases
are not orthogonal with respect to the usual L 2 inner
product, so one must explicitly orthonormal-
ize the basis before computing the matrix repre-
9 That is, instead of solving
A  v. This is necessary because b
HN; is non-hermitian, and
conjugate gradient only works for positive definite matrices.
This is not the best numerical method for non-hermitian prob-
lems, but it is easy to implement and suffices in this case.
for details on the package, as well as an overview
of Krylov subspace methods.
11 Most of the error in solving the matrix equation ( b
concentrates on eigenspaces of ( b
with large eigenvalues. These are precisely the desired eigen-
values, so in principle one can tolerate inaccurate solutions.
However, the calculation requires convergence criteria and error
estimates for the linear solver, and using a > 0, say a
turns out to ensure a relative error of about 10 6 after about
17-20 iterations of the conjugate gradient solver. Since we only
wanted to count eigenvalues, a more accurate (and expensive)
computation of resonances was not necessary.

Figure

3: A sinc function with
sentation of b
H . In addition to the computational
cost, this also requires storing the entire matrix
and severely limits the size of the problem one can
solve. Instead, this computation uses a discrete-
variable representation (DVR) basis [20]:
Consider, for the moment, the one dimensional
problem of finding a basis for a "good" subspace
of L 2 (R). Fix a constant x > 0, and for each
integer m, define
x
sin
x
(This is known as a "sinc" function in engineering
literature [25]. See Figure 3.) The Fourier transform
of  m;x is
One can easily verify that f m;x g forms an orthonormal
basis for the closed subspace of L 2
functions whose Fourier transforms are supported
in [ =x;=x].
To find a basis for corresponding space of band-limited
functions in L 2 (R 2 ), simply form the tensor
products
mn
The basis has a natural one-to-one correspondence
with points (m 1 regular
lattice of grid points in a box [X
covering the spatial region of interest. (See Figure
4.) Using this basis, it is easy to compute matrix
elements
H .
y
grid
points

Figure

4: Illustration of resonance program parameters
in configuration space: The lower-left
corner of the mesh is (X while the upper
right corner is (X The mesh contains
x  N y grid points, and a basis function mn is
placed at each grid point. Stars mark the centers
of the potentials, the circles have radius
Parameters for the classical computation
are depicted in Figure 5.
4.3.2 Tensor Product Structure
An additional improvement comes from the separability
of the Hamiltonian: Each term in the
scaled Hamiltonian b
H splits into a tensor product

dxI y (36)
I
x
dy 2
G
G
where I x and I y denote identity operators on
copies of L 2 (R). Since the basis fmn g consists
of tensor products of one dimensional bases, b
is also a short sum of tensor products. Thus, if
we let N x denote the number of grid points in
the x direction and let N y denote the number of
grid points in the y direction, then
and b
HN; is a sum of five matrices of the form
A
x
A y , where A x is N x N x and A y is N y N y .
Such tensor products of matrices can be applied
to arbitrary vectors efficiently using the outer
product representation. 12 Since the rank of b
12 The tensor product of two column vectors v and w can be
is y and N x  N y in these computa-
tions, we can store the tensor factors of the matrix
using O(N ) storage instead of O(N 2 ), and
apply b
HN; to a vector in O N 3=2
of O N 3
. The resulting matrix is not sparse,
as one can see from the matrix elements for the
Laplacian below.
Note that this basis fails to take advantage of the
discrete rotational symmetry of the triple gaussian
Hamiltonian. Nevertheless, the tensor decomposition
provides sufficient compression of information
to facilitate efficient computation.
4.3.3 Matrix Elements
It is straightforward to calculate matrix elements
for the Laplacian on
There is no closed form expression for the matrix
elements of the potential, but it is easy to perform
numerical quadrature with these functions. For ex-
ample, to compute
Z
G(x)m (x)n (x) dx (40)
for
Vmn
where the stepsize  should satisfy   x=2. It
is easy to show that the error is bounded by the
sum of

which controls the aliasing error, and
x
exp
which controls the truncation error.
4.3.4 Other Program Parameters
The grid spacing x implies a limit on the maximum
possible momentum in a wave packet formed
representedas vw T . We then have
B)(v
(Bw) T , which extendsby linearity to
by this basis. In order to obtain a finite-rank opera-
tor, it is also necessary to limit the number of basis
functions.
The resonance computation used the following
parameters:
1. are chosen to cover the
region of the configuration space for which
2. Let L
the dimension of the computational domain.
The resonance calculation uses
basis functions, with N
2h and
2h .
3. This gives
which limits the maximum momentum in a
wave packet to jp x j
and jp y j
2E.
5 Trapped Set Structure
5.1 Poincar e Section
Because the phase space for the triple gaussian
model is R 4 and its flow is chaotic, a direct computation
of the trapped set dimension is difficult.
Instead, we try to compute its intersection with a
Poincare section:
Let E be a fixed energy, and recall that R is the
distance from each gaussian bump to the origin.
Choose R 0 < R so that the circles C k of radius R 0
centered at each potential, for do not
intersect. The angular momentum p  with respect
to the kth potential center is defined by p
be the submanifold
(see

Figure

5), where the coordinates (; p  ) in the
submanifold P k are related to ambient phase space
coordinates (x;
and the radial momentum p r is
s
Figure

5: A typical trajectory: Stars mark the potential
centers. In this case, R = 1:4 and
The circles drawn in the figure have radius 1, and
the disjoint union of their cotangent bundles form
the Poincare section. Trajectories start on the circle
centered at bump 0 (the bumps are, counter-
clockwise, 0, 1, and 2) with some given angle
and angular momentum p  . This trajectory generates
the finite sequence ( _
bolic sequences are discussed later in the paper.)
An illustration of resonance computation is depicted
in Figure 4. The dashed line is the time-reversed
trajectory with the same initial condi-
tions, generating the sequence (1; 2; 0; 2; _
0).
Note that this implicitly embeds P into the energy
surface and the radial momentum p r is
always positive: The vector (p x
from the center of C k .
The trapped set is naturally partitioned into two
subsets: The first consists of trajectories which
visit all three bumps, the second of trajectories
which bounce between two bumps. The second set
forms a one-dimensional subspace of KE , so the
finite stability of the Minkowski dimension 13 implies
that the second set does not contribute to the
dimension of the trapped set. More importantly,
most trajectories which visit all three bumps will
also cut through P .
One can thus reduce the dimension of the problem
by restricting the flow to KE \ P , as follows:
Take any point (; p  ) in P k , and form the corresponding
point Equation
(45). Follow along the trajectory  t (x;
If the trajectory does not escape, eventually it
must encounter one of the other circles, say C
Generically, trajectories cross C k 0 twice at each
encounter, and we denote the coordinates (
(in of the outgoing intersection by
e
If a trajectory escapes from the trapping region,
we can symbolically assign 1 to e
. The map e
then generates stroboscopic recordings of the flow
t on the submanifold P , and the corresponding
discrete dynamical system has trapped set KE \
P . So, instead of computing  t on R 4 , one only
needs to compute e
on P . By symmetry, it will
suffice to compute the dimension of e
Pushing e
KE forward along the flow  t adds
one dimension, so D (KE
KE
+1. Being
a subset of the two-dimensional space P 0 , e
KE is
easier to work with.
Readers interested in a more detailed discussion
of Poincare sections and their use in dynamics are
referred to [29]. For an application to the similar
but simpler setting of hard disc scattering, see [9,
12]. Also, Knauf has applied some of these ideas
in a theoretical investigation of classical scattering
by Coulombic potentials [17].
5.2 Self-Similarity
Much is known about the self-similar structure of
the trapped set for hard disc scattering [9, 12];
13 That is,
see [8].
theta
ptheta

Figure

Points in P 0 which do not go to 1 after
one iteration of e
. The horizontal axis is  and the
vertical axis is p  .
less is known about "soft scatterers" like the triple
gaussian system. However, computational results
and analogy with hard disc scattering give strong
support to the idea that K (and hence e
KE ) is self-
similar. 14 Consider Figures 6 - 12: They show
clearly that e
KE is self-similar. (In these images,
However, it is also clear
that, unlike objects such as the Cantor set or the
Sierpinski gasket, e
KE is not exactly self-similar.
5.3 Symbolic Dynamics
The computation of D
KE
uses symbolic se-
quences, which requires a brief explanation: For
any point (; p  ), let s i denote the third component
of e
integer i. Thus,
s i is the index k of the circle C k that the trajectory
intersects at the ith iteration of e  (or the
jijth iteration of e
occuring only at
the ends. Let us call sequences satisfying these
conditions valid.
For example, the trajectory in Figure 4 generates
the valid sequence ( _
1),where the
dot over 0 indicates that the initial point (; p  ) of
the trajectory belongs to P 0 . Thus, we can label
collections of trajectories using valid sequences,
and label points in P with "dotted" sequences.
14 More precisely, self-affine.
theta
ptheta

Figure

7: Points in P 0 which do not go to 1 after
one iteration of e  1 . The horizontal axis is  and
the vertical axis is p  .
2.5 3 3.5 4
theta
ptheta

Figure

8: The intersection of the sets in Figures
6 and 7. These points correspond to symmetric
sequences of length 3.
-0.050.05theta
ptheta

Figure

9: The lower-right "island" in Figure 8,
magnified. The white cut-out in the middle is the
subset corresponding to symmetric sequences of
length 5.
3.55 3.6 3.65
theta
ptheta

Figure

10: The cut-out part of Figure 9, magni-
fied. Recall that these correspond to symmetric
sequences of length 5; compare with Figure 8.
3.2 3.4 3.6 3.8 40.10.30.5
theta
ptheta

Figure

11: The upper-right island in Figure 8. The
white cut-out in the middle is, again, the subset
corresponding to symmetric sequences of length
5.
ptheta

Figure

12: The cut-out part of Figure 11, magni-
fied. Recall that these correspond to symmetric
sequences of length 5. Compare with Figures 8
and 10.
Clearly, trapped trajectories generate bi-infinite
sequences. 15
The islands in Figures 8 - 11 correspond to
symmetric sequences centered at 0, of the form
track of the symbolic sequences generated by
each trajectory, one can easily label and isolate
each island. This is a useful property from the
computational point of view.
5.4 Dimension Estimates
To compute the Minkowski dimension using
Equation (15), we need to determine when a given
point is within  of e
KE . This is generally impossi-
ble: The best one can do is to generate longer and
longer trajectories which stay trapped for increasing
(but finite) amounts of time.
Instead, one can estimate a closely related quan-
tity, the information dimension, in the following
way: Let e
denote the set of all points in P 0
corresponding to symmetric sequences of length
centered at 0. That is, e
consists of
all points in P 0 which generate trajectories (both
forwards and backwards in time) that bounce at
least k times before escaping. The sets e
decrease
monotonically to e
e
KE .
One can then estimate the information dimension
using the following algorithm:
1. Initialization: Cover P 0 with a mesh L 0 with
points and mesh size  0 .
2. Recursion: Begin with e
K (1)
E , which consists
of four islands corresponding to symmetric
sequences of length 2
8). Magnify each of these islands and compute
the sub-islands corresponding to symmetric
sequences of length 5 (see Figures 9
and 11). Repeat this procedure to recursively
compute the islands of e
E from those of
e
sufficiently
large that each island of e
diameter smaller than the mesh size  0 of L 0 .
3. Estimation: Using the islands of e
In hard disc scattering, the converse holds for sufficiently
large R: To each bi-infinite valid sequence there exists a
trapped trajectory generating that sequence. This may not hold
in the triple gaussian model, and in any case it is not necessary
for the computation.
theta
ptheta

Figure

13: This figure illustrates the recursive
step in the dimension estimation algorithm: The
dashed lines represent while the solid lines
represent a smaller mesh centered on one of the
islands. The N 0  N 0 mesh L 0 remains fixed
throughout the computation, but the smaller N 1
constructed for each island of e
to the value of k specified by the algorithm.
mate the probability
vol
vol
for the (ij)th cell of L 0 . We can then compute
the dimension via
KE
log
which reduces to (15) when the distribution
is uniform because  0  1=N 0 .
Under suitable conditions (as is assumed to be
the case here), the information dimension agrees
with both the Hausdorff and the Minkowski di-
mensions.
The algorithm begins with the lattice L 0 with
which one wishes to compute the dimension. It
then recursively computes e
E for for increasing
values of k, until it closely approximates e
KE relative
to the mesh size of L 0 . It is easy to keep track
of points belonging to each island in this com-
putation, since each island corresponds uniquely
See [23] for a discussion of the relationship between these
dimensions, as well as their use in multifractal theory.

Figure

14: These are the eigenvalues of b
f0:0624; 0:0799; 0:0973g. This calculation used
an 102  108 grid, and 90 out of
eigenvalues were computed.
to a finite symmetric sequence. Note that while
the large mesh L 0 remains fixed throughout the
computation, the recursive steps require smaller
around each island of e
to the value of k specified by the algorithm. See

Figure

13.
6 Numerical Results
6.1 Resonance Counting
As an illustration of complex scaling, Figures 14
0:025]. Eigenvalues of b
HN; for different
values of  are marked by different styles of
points, and the box has depth h and width 0:2, with
These plots may seem
somewhat empty because only those eigenvalues
of b
HN; in regions of interest were computed.
Notice the cluster of eigenvalues near the bottom
edge of the plots: These are not resonances because
they vary under perturbations in . Instead,
they belong to an approximation of the (scaled)
continuous spectrum.
It is more interesting to see log (N res ) as a function
of log (h) and R. This is shown in Figures
19 and 20. Using least-squares regression, we can
extract approximate slopes for the curves in Figure
19; these are shown in Table 2 and plotted in

Figure

15: Eigenvalues for
using 112  119 grid
and 98 out of eigenvalues.

Figure

using 123  131 grid
and 107 out of eigenvalues.

Figure

17: Eigenvalues for
using 135  144 grid
and 116 out of

Figure

18: Eigenvalues for
using
out of
eigenvalues.
3.7 3.8 3.9 43.64
log(#res)

Figure

19: log (N res ) as a function of log (h),
for h varying from 0:017 to 0:025 and
0:05  k, with 0  k  6. (The lowest curve corresponds
to while the highest curve corresponds
to
R
log(#res)

Figure

20: log (N res ) as a function of R, for different
values of log (h): The highest curve corresponds
to while the lowest curve corresponds
to
R
scaling
exponents

Figure

21: The slopes extracted from Figure 19, as
a function of R. The dotted curve is a least-squares
regression of the "noisy" curve.
1.6 1.2986 1.2725 1.2511
1.7 1.2893 1.2636 1.2524

Table

1: Estimates of D(KE)+1as a function of R.

Figure

21.
6.2 Trapped Set Dimension
For comparison, D(KE)+1is plotted as a function
of R in

Figure

22. The figure contains curves
corresponding to different energies E: The top
curve corresponds to 0:4, the middle curve
0:5, and the bottom curve
also contains curves corresponding to different
program parameters, to test the numerical convergence
of dimension estimates. These curves were
computed with  2
recursion depth k (corresponding to symmetric
sequences of length 2  6 the
caption contains the values of N 0 and N 1 for each
curve. For reference, Table 1 contains the dimension
estimates shown in the graph. It is important
to note that, while the dimension does depend on
R
(D+1)/2

Figure

22: This figure shows D(KE)+1as a function
of R: The top group of curves have
the middle 0:5, and the bottom
Solid curves marked with circles represent computations
3 , and 1
2 . Dashed curves marked
with X's represent computations where N
14142, whereas dashed curves marked with triangles
represent computations where N
and 71. The recursion depth k 0 in all these
figues is 6. The curve does not appear
to have completely converged but suffices for our
purpose here.
R slope D(KE)+1relative error
1.6 1.3055 1.2725 0.025256

Table

2: This table shows the slopes extracted
from

Figure

19, as well as the scaling exponents
one would expect if the conjecture were true
computed at
errors are also shown.
E and R, it only does so weakly: Relative to its
value, D(KE)+1is very roughly constant across
the range of R and E computed here.
6.3 Discussion

Table

contains a comparison of D(KE)+1
(for
as a function of R, versus the scaling exponents
from Figure 21. Figure 23 is a graphical
representation of similar information. This figure
shows that even though the scaling curve in Figure
21 is noisy, its trend nevertheless agrees with
the conjecture. Furthermore, the relative size of
the fluctuations is small. At the present time, the
source of the fluctuation is not known, but it is possibly
due to the fact that the range of h explored
here is simply too large to exhibit the asymptotic
behavior clearly. 17

Figures

plots of log (N res ) versus
log (h), for various values of R. Along with
the numerical data, the least-squares linear fit and
the scaling law predicted by the conjecture are also
plotted. In contrast with Figure 23, these show
clear agrement between the asymptotic distribution
of resonances and the scaling exponent predicted
by the conjecture.
6.4 Double Gaussian Scattering
Finally, we compute resonances for the double
gaussian model (setting in (8). This case is
17 But see Footnote 8.
The conjecture only supplies the exponents for power
laws, not the constant factors. In the context of these logarithmic
plots, this means the conjecture gives us only the slopes,
not the vertical shifts. It was thus necessary to compute an y-intercept
for each "prediction" curve (for the scaling law predicted
by the conjecture) using least squares.
R
slopes
scaling
exponents)

Figure

23: Lines with circles represent D(KE )+1as functions of R, for 0:6g. The
dashed curve with triangles is the scaling exponent
curve from Figure 21, while the solid curve
with stars is the linear regression curve from that
figure. Relative to the value of the dimension, the
fluctuations are actually fairly small: See Table 2
for a quantitative comparison.
3.7 3.8 3.9 43.43.6-log(hbar)
log(#res)
R=1.4

Figure

24: For
data, circles least squares regression, and
stars the slope predicted by the conjecture. h
ranges from 0:025 down to 0:017.
log(#res)

Figure

25: Same for
from 0:025 to 0:017.
3.7 3.8 3.9 43.53.7-log(hbar)
log(#res)

Figure

3.7 3.8 3.9 43.53.73.9
log(#res)

Figure

3.7 3.8 3.9 43.63.8-log(hbar)
log(#res)
R=1.6

Figure

log(#res)

Figure

3.7 3.8 3.9 43.73.9-log(hbar)
log(#res)

Figure


Figure

31: Resonances for two-bump scattering
with
interesting for two reasons: First, there exist rigorous
results [13, 21] against which we can check
the correctness of our results. Second, it helps determine
the validity of semiclassical arguments for
the values of h used in computing resonances for
the triple gaussian model.
The resonances are shown in Figures 31 - 37:
In these plots, ranges from 0:035
to 0:015. One can observe apparent pseudospectral
effects in the first few figures [30, 36]; this is
most likely because the scaling angle  used here
is twice as large as suggested in Section 4.1, to exhibit
the structure of resonances farther away from
the real axis.
To compare this information with known results
[13, 21], we need some definitions: For a given
Z x1 (E)
where the limits of integration are
Let (E) denote the larger (in absolute value)
eigenvalue of D e
(0; 0); log () is the Lyapunov
exponent of e
, and is easy to compute numerically
in this case. Note that for two-bump scat-
tering, each energy E determines a unique periodic
trapped trajectory, and C(E) is the classical
action computed along that trajectory.

Figure

32: Resonances for two-bump scattering
with
-0.06
-0.020.020.06

Figure

33: Resonances for two-bump scattering
with
-0.06

Figure

34: Resonances for two-bump scattering
with

Figure

35: Resonances for two-bump scattering
with

Figure

Resonances for two-bump scattering
with

Figure

37: Resonances for two-bump scattering
with
-0.4
F
Figure

38: Lattice points for
Since these expressions are analytic, they have
continuations to a neighborhood of the real line -
C(E) becomes a contour integral. In [13], it was
shown that any resonance
must satisfy
O
where m and n are nonnegative integers. (The 1in m+ 1comes from the Maslov index associated
with the classical turning points.) This suggests
that we define the map F
where
Re (C())
and
Im (C())
h log ( (Re ()))
F should map resonances to points on the square
integer lattice, and this is indeed the case: Figures
images of resonances under
F , with circles marking the nearest lattice points.
The agreement is quite good, in view of the fact
that we neglected terms of order h 2 in Equation
(52).
Conclusions
Using standard numerical techniqes, one can compute
a sufficiently large number of resonances for
the triple gaussian system to verify their asymptotic
distribution in the semiclassical limit h ! 0.
F
Figure

39: Lattice points for
-0.4
F
Figure

40: Lattice points for
-0.4
F
Figure

41: Lattice points for
-0.4
F
Figure

42: Lattice points for
F
Figure

43: Lattice points for
-0.4
F
Figure

44: Lattice points for
This, combined with effective estimates of the
fractal dimension of the classical trapped set, gives
strong evidence that the number of resonances
N res in a box [E sufficiently
small h,
res  h
as one can see in Figure 23 and Table 2. Fur-
thermore, the same techniques, when applied to
double gaussian scattering, produce results which
agree with rigorous semiclassical results. This
supports the correctness of our algorithms and the
validity of semiclassical arguments for the range
of h explored in the triple gaussian model. The
computation also hints at more detailed structures
in the distribution of resonances: In Figures 14 -
18, one can clearly see gaps and strips in the distribution
of resonances. A complete understanding
of this structure requires further investigation.
While we do not have rigorous error bounds
for the dimension estimates, the numerical results
are convincing. It seems, then, that the primary
cause for our failure to observe the conjecture in a
"clean" way is partly due to the size of
could study resonances at much smaller values of
h, the asymptotics may become more clear.

Acknowledgments

Many thanks to J. Demmel and B. Parlett for crucial
help with matrix computations, and to X. S.
Li and C. Yang for ARPACK help. Thanks are
also due to R. Littlejohn and M. Cargo for their
help with bases and matrix elements, and to F.
Bonetto and C. Pugh for suggesting a practical
method for computing fractal dimensions. Z. Bai,
W. H. Miller, and J. Harrison also provided many
helpful comments. Finally, the author owes much
to M. Zworski for inspiring most of this work.
KL was supported by the Fannie and John Hertz
Foundation. This work was supported in part by
the Applied Mathematical Sciences subprogram
of the Office of Energy Research of the U. S. Department
of Energy under Contract DE-AC03-76-
SF00098. Computational resources were provided
by the National Energy Research Scientific Computing
Center (NERSC), the Mathematics Department
at Lawrence Berkeley National Laboratory,
and the Millennium Project at U. C. Berkeley.



--R

"A class of analytic perturbations for one body Schrodinger Hamiltonians,"
LAPACK User's Guide
"Spectral properties of many-body Schrodingeroperators with dilation analytic interactions,"
Templates for the Solution of Linear Sys- tems: Building Blocks for Iterative Meth- ods
Quantum Mechanics of One- and Two-Electron Atoms
"Semiclassical Gaussian basis set method for molecular vibrational wave functions,"
Spectral Asymptotics in the Semi-Classical Limit
Fractal Geometry: Mathematical Foundations and Applications.
"Scattering from a classically chaotic repellor,"
"Semiclassical quantization of the scattering from a classically chaotic repellor,"
"Exact quantization of the scattering from a classically chaotic repellor,"

"Semiclassical resonances generated by a closed trajectory of hyperbolic type,"
"Wave trace for Riemann surfaces,"
"On distributed Gaussian bases for simple model multidimensional vibrational problems,"
"Can one hear the shape of a drum?"
"The n-Centre Problem of Celestial Mechanics"

ARPACK User's Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods
"Gen- eralized discrete variable approximation in quantum mechanics,"
"Classical-limit Green's function (fixed-energy propagator) and classical quantization of nonseparable systems,"
"Tunneling and state specificity in unimolecular reactions,"
Dimension Theory in Dynamical Systems: Contemporary Views and Ap- plications
"Collisional breakup in a quantum system of three charged particles,"

"The definition of molecular resonance curves by the method of exterior complex scaling,"
"Geometric bounds on the density of resonances for semi-classical problems,"
"Complex scaling and the distribution of scattering poles,"
Structure and Interpretation of Classical Mechanics.
"Pseudospectra of Linear Operators,"
"Model studies of mode specificity in unimolecular reaction dynamics,"
"Mode specificity in unimolecular reaction dynamics: the Henon-Heiles potential energy surface,"
"Numerical Computation of the Scattering Frequencies for Acoustic Wave Equations,"
"Quantum mechanics and semi- classics of hyperbolic n-disk scattering sys- tems,"
"Dimension of the limit set and the density of resonances for convex co-compact hyperbolic surfaces,"
"Numerical linear algebra and solvability of partial differential equations"
"Resonances in Geometry and Physics,"
--TR
Circuits,signals,and systems
Numerical computational of the scattering frequencies for acoustic wave equations
Pseudospectra of Linear Operators
Structure and interpretation of classical mechanics
scattering resonances;semiclassical asymptotics;chaotic trapping;fractal dimension