Computation of eigenmodes
on a compact hyperbolic 3space
Abstract
Measurements of cosmic microwave background (CMB) anisotropy are ideal experiments for discovering the nontrivial global topology of the universe. To evaluate the CMB anisotropy in multiplyconnected compact cosmological models, one needs to compute the eigenmodes of the LaplaceBeltrami operator. Using the direct boundary element method, we numerically obtain the lowlying eigenmodes on a compact hyperbolic 3space called the Thurston manifold which is the second smallest in the known compact hyperbolic 3manifolds. The computed eigenmodes are expanded in terms of eigenmodes on the unit threedimensional pseudosphere. We numerically find that the expansion coefficients behave as Gaussian pseudorandom numbers for lowlying eigenmodes. The observed gaussianity in the CMB fluctuations can partially be attributed to the Gaussian pseudorandomness of the expansion coefficients assuming that the Gaussian pseudorandomness is the universal property of the compact hyperbolic spaces.
Yukawa Institute for Theoretical Physics
Kyoto University, Kyoto 6068502, Japan
1 Introduction
In recent years, there has been a great interest in properties of CMB
anisotropy in multiplyconnected cosmological
models [2, 3, 4, 5].
Most of these studies deal with flat models or noncompact
hyperbolic models for
which the eigenmodes are known explicitly. Since no closed
analytic expression of eigenmodes is known for compact
hyperbolic (CH) models, so far, analysis of the CMB anisotropy in CH models
has been considered to be quite difficult although they have
interesting properties which are strikingly different from that of
multiplyconnected flat models.
For instance, in low
adiabatic models, the largeangular fluctuations can be produced at
periods after the last scattering as the
curvature perturbations decay in the curvature dominant era. Therefore, the
argument of the suppression of the largeangular fluctuations due
to the ”modecutoff” in the multiplyconnected flat models cannot
simply be applicable to the multiplyconnected hyperbolic models.
Because the effect of the multiplyconnectedness becomes
significant as the volume of the space becomes small, it is very
important to study whether the”small” universe scenario
is plausible. For instance, the Weeks manifold and the
Thurston manifold have volume where
denotes the curvature radius, and they are the smallest and the second
smallest compact hyperbolic manifolds, respectively. For a technical
reason, we study the properties of eigenmodes on the Thurston
manifold. Since each type(scalar, vector and tensor) of perturbations
evolves independently on
the locally homogeneous and isotropic FRW background space, in
space,
the perturbations are given by the eigenmodes and the
time evolution of the perturbations on the locally
homogeneous and isotropic FRW space. The periodic
boundary conditions on the eigenmodes drastically change
the nature of the CMB fluctuations on the topological identification
scale while on smaller scale they asymptotically converge to that
in the standard locally and globally homogeneous and isotropic
FRW background space. Thus computation of eigenmodes are very
important in understanding the properties of CMB fluctuations
^{1}^{1}1It should be noted that the computation of
eigenmodes is also essential in the
framework of spectral geometry [6]..
Computing the eigenmodes of the
LaplaceBeltrami operator in CH spaces (manifolds)
is equivalent to solving the Helmholtz
equation with appropriate periodic boundary conditions in the
universal covering space.
A number of numerical methods have been used for solving the
Helmholtz equation such as the finite element methods and the
finite difference methods [7, 8, 9].
A numerical method called the ”direct boundary element method”
(DBEM) has been used by Aurich and Steiner for finding out
the eigenmodes of the LaplaceBeltrami operator in a
twodimensional compact multiplyconnected
space for studying the statistical properties of the eigenmodes in
highlyexcited states or
equivalently the semiclassical wavefunctions [10].
We find that this pioneering
work for ”quantum chaology” , the study of
the imprints of classical chaos in the quantum
mechanical counterparts is very useful for the study
of the CMB anisotropy in CH cosmological models as well.
The advantage of the DBEM is that it reduces the dimensionality of the
problem by one which leads to economy in the numerical task.
Since one needs to discretize only the boundary, generation of
meshes is much easier than the other methods. Furthermore, as we shall see
later, the DBEM is suitable for expanding the eigenfunctions that are
continued onto the whole Poincar ball by the periodic
boundary conditions in terms of
eigenfunctions on the pseudosphere^{2}^{2}2A set of the continued
eigenfunctions is a subset of all eigenfunctions on the universal
covering space.. In order to compute a set of expansion
coefficients, one needs to compute the values of the corresponding
eigenfunction on a hyperbolic sphere with appropriate
radius. If one does not care about the normalization of the
eigenfunctions, unlike the FEM, the computation of the eigenfunctions
on the whole fundamental domain is not
necessary. In the DBEM, computation of the eigenfunction at an arbitrary point
either inside or outside of the fundamental domain
can be done by using the values of the eigenfunction and the
normal derivatives on the boundary.
As the classical dynamical systems in CH spaces are strongly
chaotic, one can naturally assume that the imprint of the
classical chaos is hidden in the
corresponding quantum systems in some way.
It has been found that the expansion coefficients with a
certain basis behave as if they
are random Gaussian numbers in some
classically chaotic systems ([10, 11]), which is consistent
with the prediction of randommatrix theory([12]).
Since the CMB temperature fluctuations in CH spaces are written
in terms of expansion
coefficients and the eigenfunctions on the universal covering space
plus initial fluctuations, if the random behavior of the expansion
coefficients is confirmed,
the origin of the random gaussianity
in the CMB temperature fluctuations can be partially explained
in terms of the geometric property of the universe.
In this paper, we introduce the DBEM for solving the Helmholtz equation.
Then we apply the DBEM for computing the lowlying eigenmodes of the
LaplaceBeltrami operator in the Thurston manifold.
The computed eigenfunctions are naturally continued onto the whole
Poincar
ball because of the periodic boundary conditions and are expanded in
terms of eigenmodes on the simplyconnected pseudosphere.
Statistical properties
of the expansion coefficients are examined, since they are key factors
in understanding of the CMB anisotropy in CH models.
2 The direct boundary element method (DBEM)
The boundary element methods (BEM) use free Green’s function as
the weighted function, and the Helmholtz equation is
rewritten as an integral equation defined on the
boundary using Green’s theorem. Discretization of the
boundary integral equation yields a
system of linear equations. Since one needs the discretiztion on only
the boundary, BEM reduces the dimensionality of the
problem by one which leads to economy in the numerical task.
To locate an eigenvalue, the DBEM ^{3}^{3}3The DBEM uses only
boundary points in evaluating the integrand in
Eq.(5). The indirect
methods use internal points in evaluating the integrand in
Eq.(5) as well as the boundary points. requires
one to compute many determinants of the corresponding boundary
matrices which are dependent on the wavenumber .
Firstly, let us consider the Helmholtz equation with certain boundary
conditions,
(1) 
which is defined on a bounded Mdimensional connected and simplyconnected domain which is a subspace of a Mdimensional Riemannian manifold and the boundary is piecewise smooth. , and is the covariant derivative operator defined on . A function in Sobolev space is the solution of the Helmholtz equation if and only if
(2) 
where is an arbitrary function in Sobolev space called weighted function and is defined as
(3) 
Next, we put into the form
(4) 
where ’s are linearly independent
squareintegrable functions. Numerical
methods such as the finite element methods try
to minimize the residue function for a fixed
weighted function by changing the coefficients .
In these methods, one must
resort to the variational principle to find the ’s which
minimize .
Now we formulate the DBEM
which is a version of BEMs. Here we search ’s for the space
.
First, we slightly
modify Eq.(2) using the Green’s theorem
(5) 
where and ; the surface element is given by
(6) 
where denotes the Mdimensional LeviCivita tensor. Then Eq.(2) becomes
(7) 
As the weighted function v, we choose the fundamental solution which satisfies
(8) 
where , and is Dirac’s delta function. is also known as the free Green’s function whose boundary condition is given as
(9) 
where is the geodesic distance between and . Let be an internal point of . Then we obtain from Eq.(7) and Eq.(8),
(10) 
Thus the values of eigenfunctions at internal points can be computed using only the boundary integral. If , we have to evaluate the limit of the boundary integral terms as becomes divergent at (see appendix A). The boundary integral equation is finally written as
(11) 
or in another form,
(12) 
where and . Note that we assumed that the boundary surface at is sufficiently smooth. If the boundary is not smooth, one must calculate the internal solid angle at (see appendix A). Another approach is to rewrite Eq.(10) in a regularized form [14]. We see from Eq.(11) or Eq.(12) that the approximated solutions can be obtained without resorting to the variational principle. Since it is virtually impossible to solve Eq.(12) analytically, we discretize it using boundary elements. Let the number of the elements be N. We approximate by some loworder polynomials (shape function) on each element as where and denote the coordinates on the corresponding standard element ^{4}^{4}4It can be proved that the approximated polynomial solutions converge to as the number of boundary elements becomes large [15, 16]..
Then we have the following equation:
(13) 
where {u} and {q} are Ndimensional vectors which consist of the boundary values of an eigenfunction and its normal derivatives, respectively. [H] and [G] are NN dimensional coefficient matrices which are obtained from integration of the fundamental solution and its normal derivatives multiplied by and , respectively. The explicit form of [H] and [G] for constant elements are given in section 4. Note that the elements in [H] and [G] include implicitly. Because Eq.(13) includes both and , the boundary element method can naturally incorporate the periodic boundary conditions:
(14) 
where ’s are the facetoface identification maps defined on the boundary(see appendix B). The boundary conditions constrain the number of unknown constants to N. Application of the boundary condition (14) to Eq.(13) and permutation of the columns of the components yields
(15) 
where NN dimensional matrix is constructed from and and N dimensional vector is constructed from ’s and ’s. For the presence of the nontrivial solution, the following relation must hold,
(16) 
Thus the eigenvalues of the LaplaceBeltrami operator acting on the space are obtained by searching for ’s which satisfy Eq.(16).
3 Computation of lowlying eigenmodes
In this section, we apply the DBEM for computing
the lowlying eigenmodes on the Thurston manifold .
We have chosen for a technical reason that the
fundamental domain of is much simpler than that of
the Weeks manifold (generation of meshes is much simpler). See appendix B
for understanding the basic aspects of threedimensional
hyperbolic geometry.
The Helmholtz equation in the Poincar coordinates
is written as
(17) 
where and are the Laplacian and the gradient on the corresponding threedimensional Euclidean space, respectively. Note that we have set the curvature radius without loss of generality. By using the DBEM, the Helmholtz equation (17) is converted to an integral representation on the boundary. Here Eq.(12) can be written in terms of Euclidean quantities as
(18) 
where . The fundamental solution is given as [17, 18]
(19) 
where and . Then Eq.(18) is discretized on the boundary elements as
(20) 
where N denotes the number of the boundary elements. An example of elements on the boundary of the fundamental domain in the Poincar coordinates is shown in figure 1. These elements are firstly generated in Klein coordinates in which the meshgeneration is convenient.
The maximum length of the edge in these elements is 0.14. The condition that the corresponding de Broglie wavelength is longer than the four times of the interval of the boundary elements yields a rough estimate of the validity condition of the calculation as . On each , u and q are approximated by low order polynomials. For simplicity, we use constant elements:
(21) 
Substituting Eq.(21) into Eq.(20), we obtain
where
(22) 
The singular integration must be carried out
for II components as the fundamental solution diverges at
(). This is not an intractable problem. Several numerical
techniques have already been proposed by some authors [19, 20].
We have applied Hayami’s method
to the evaluation of the singular integrals [20]. Introducing
coordinates similar to spherical coordinates
centered at , the singularity is
canceled out by the Jacobian which makes the integral regular.
Let be the generators of the discrete
group which identify a boundary face
with another boundary face :
(23) 
The boundary of the fundamental domain can be divided into two regions and and each of them consists of boundary elements,
(24) 
The periodic boundary conditions
(25) 
reduce the number of the independent variables to N, i.e. for all , there exist and such that
(26) 
Substituting the above relation into Eq.(22), we obtain
(27) 
where and and matrices and are written as
(28) 
Eq. (27) takes the form
(29) 
where NN dimensional matrix is constructed from and and N dimensional vector is constructed from and . For the presence of the nontrivial solution, the following relation must hold,
(30) 
Thus the eigenvalues of the LaplaceBeltrami operator in a CH space are obtained by searching for ’s which satisfy Eq.(30). In practice, Eq.(30) cannot be exactly satisfied as the test function which has a locally polynomial behavior is slightly deviated from the exact eigenfunction. Instead, one must search for the local minima of det[A(k)]. This process needs long computation time as depends on k implicitly. Our numerical result () is shown in table 1.
k  

5.41  1 
5.79  1 
6.81  1 
6.89  1 
7.12  1 
7.69  1 
8.30  1 
8.60  1 
8.73  1 
9.26  2 
9.76  1 
9.91  1 
9.99  1 
The first ”excited state” which corresponds to is important for the understanding of CMB anisotropy. Our numerical result is consistent with the value obtained from Weyl’s asymptotic formula
(31) 
assuming that no degeneracy occurs.
One can interpret the first excited state as the mode that has
the maximum de Broglie wavelength . Because of the
periodic boundary conditions, the de Broglie wavelength
can be approximated by the
”average diameter” of the fundamental domain defined as a sum of
the inradius and the outradius ^{5}^{5}5The inradius
is the radius of the largest
simplyconnected sphere in the fundamental domain, and the outradius
is the radius of the smallest sphere that can enclose the
fundamental domain. , for
the Thurston manifold.,
which yields just
less than the numerical value. From these estimates,
supercurvature modes in small CH spaces ()
are unlikely to be observed.
To compute the value of eigenfunctions inside the fundamental domain,
one needs to solve Eq.(29). The singular decomposition
method is the most suitable numerical method for solving any linear
equation with a singular matrix , which can be
decomposed as
(32) 
where and are unitary matrices and D is a diagonal matrix. If in is almost zero then the complex conjugate of the ith row in V is an approximated solution of Eq.(29). The number of the ”almost zero” diagonal elements in D is equal to the multiplicity number.
Substituting the values of the eigenfunctions and their normal
derivatives on the
boundary into Eq.(10), the values of the
eigenfunctions inside the fundamental domain can be computed.
Eigenfunctions and
in Poincar
coordinates plotted as (), where
are shown in
figure 2. The eigenfuctions we computed are all
realvalued. Note that the
nondegenerated eigenfunctions must be realvalued.
The numerical accuracy of the obtained eigenvalues is roughly estimated as
follows. First, let us write the obtained numerical solution in terms
of the exact solution as
and ,
where and are the exact eigenvalue and
eigenfunction, respectively. The singular decomposition method
enables us to find the best approximated solution which satisfies
(33) 
where is a Ndimensional vector and denotes the Euclidean norm. It is expected that the better approximation gives the smaller . Then Eq. (33) can be written as,
(34) 
Ignoring the terms in second order, Eq.(34) is reduced to
(35) 
Since it is not unlikely that is anticorrelated to , we obtain the following relation by averaging over ,
(36) 
where denotes the averaging over
Thus one can estimate the expected deviation of
the calculated eigenvalue from and .
We numerically find that for
and for
. The other deviation values lie in between and
.
By computing the second derivatives, one
can also estimate the accuracy of the computed eigenfunctions.
The accuracy parameter is defined as
(37) 
where is normalized (). We see from figure
3 that the accuracy becomes worse
as the evaluation point approaches the boundary. However, for points with
hyperbolic distance between the evaluating point and the
nearest boundary, the errors are very small
indeed: . This result is considered to be
natural because the characteristic
scale of the boundary elements is for our 1168 elements.
If , the integrands in Eq.(10) become appreciable on
the neighborhood of the nearest boundary
point because the free Green’s
function approximately diverges on the point.
In this case, the effect of the deviation from the exact
eigenfunction is significant. If , the integrand on
all the boundary points contributes almost equally to the
integration so that the local deviations are cancelled out.
As we shall see in the next section, expansion coefficients
are calculated using the values of eigenfunctions on a sphere.
Since the number of evaluating points which
are very close to the boundary is negligible on the sphere,
expansion coefficients can be computed with
relatively high accuracy.
4 Statistical properties of eigenmodes
Properties of eigenmodes of the LaplaceBeltrami operator
are determined by the Helmholtz
equation. Therefore, at first glance it does not seem to make a sense
to study the statistical properties of the eigenmodes.
However, if one recognizes the LaplaceBeltrami operator in
a CH space as the Hamiltonian in a quantum mechanical system,
each eigenmode can be interpreted as a wavefunction in an eigenstate.
Since the corresponding classical
system is known to be a typical chaotic system (Ksystem), it is
natural to assume that the imprints of classical chaos is hidden in
the corresponding quantum system.
Recent studies have demonstrated that some of the statistical
properties of energy spectrum are in accordance with the universal
prediction of randommatrix theory(RMT) [21, 22].
RMT also predicts that the
squared expansion coefficients of an eigenstate with respect to
a generic basis are Gaussian distributed [10, 11, 12].
In the limit , obeys the statistics
given for three universality classes of the orthogonal (GOE, ), unitary
(GUE, ) and symplectic (GSE, ) ensembles, each
distribution function is given by
(38) 
In our case, as the timereversal symmetry of the Hamiltonian implies,
one expects that obeys the GOE prediction.
In order to apply the GOE prediction to the statistical properties of
eigenstates on CH spaces, one needs to find a set of orthonormal
basis but no closed analytic expression is known for any CH spaces.
To avoid the problem, Aurich and Steiner noticed that the
wavefunctions on the hyperbolic octagons can be continued onto the
universal covering space ,
and eigenstates can be expanded in terms of circularwaves
[10].
They numerically found that the squared expansion coefficients
obeys the GOE prediction in highly excited states of
a hyperbolic asymmetrical octagon model.
We extend their method to threedimensional CH models where we consider only
lowlying modes. First, we normalize the obtained 14 eigenfunctions
on the Thurston manifold. The eigenfunctions are naturally continued
onto the whole
unit Poincar ball by the periodic boundary condition.
As a ”generic basis”, we consider a set of orthonormal
eigenfunctions (Tcomplete functions) on the unit
pseudosphere which is isometric to the Poincar ball,
(39)  
where , and denote the associated Legendre function, the spherical harmonics and gamma function, respectively. can be written in terms of the hypergeometric function [23],
(40) 
Eigenfunctions can be expanded in terms of ’s as
(41) 
Note that each has no components with
because ’s are
complete and linearly independent.
At first glance the computation of in
Eq. (41) seems cumbersome as the domain of the integration
extends over the whole pseudosphere.
Fortunately, one can obtain by evaluating
twodimensional integrals. can be written as
(42) 
which is satisfied for the arbitrary value of . In practice the numerical instability occurs in the region where the absolute value of is too small. In our computation, the values of are chosen as shown in table 2.
0 l 8  8 l 13  13 l 18  

k8  0.53  1.3  1.6 
k8  0.53  1.1  1.3 
Thus can be computed if one obtains the values of eigenfunctions on the sphere
(43) 
with radius .
In order to compute the values of eigenfunctions on the sphere
with radius longer than the inradius ,
the points outside of the fundamental domain must be pulled back to the
inside, since Eq.(10) is valid only if is a set of
coordinates of an internal point.
The plots of eigenfunctions
on a sphere are shown in figure 4 and
figure 5. Apparent structure of the eigenfunctions on the sphere
seems complicated. However, some regular patterns are hidden in the
structure due to the periodic boundary conditions. Actually, there
are pairs of highly
correlated points on the sphere, since
any partial surface of the sphere that is enclosed by copies
of the boundary of the fundamental domain pulled back to
inside the fundamental
domain by the corresponding element of the discrete isometry group
intersects another partial surface that is pulled
back to inside the fundamental domain.
To evaluate the correlation pattern, let us
estimate how often a sphere with radius
intersects the copies of the fundamental domain.
The approximate number of the copies of the
fundamental domain inside
the sphere with radius (in proper length) is given by
(44) 
From this formula, in the case of the Thurston manifold, if . Because the sphere intersects the fundamental domain at random, the copies of the fundamental domain on the sphere stick out their half portions on average. Therefore, the approximate number of the copies that intersect the sphere is given by
(45) 
where . This estimate gives if . Approximating each eigenmode by de Broglie waves, we obtain the corresponding fluctuation scale in steradian on the sphere,
(46)  
Thus correlation patterns are observed in pairs of patches with typical size . When , angular fluctuation scales are given as , respectively. , for
Next, we extract a set of independent variables from ’s. In general, any is related to as
(47) 
where
(48) 
If is real, from Eq.(47),
(49)  
therefore,
(50) 
Thus can be written in terms of . To extract a set of independent variables from ’s, we rewrite Eq.(49) as follows
(51)  
where
(52) 
and
(53)  
Thus the real eigenfuctions can be written in terms of real independent coefficients and realvalued ,
(54) 
where
(55) 
and
(56) 
Now we turn to the statistical properties of the coefficients . As in [10], we consider the cumulative distribution of following quantities,
(57) 
where is the mean of ’s and is the variance. The cumulative distribution is compared to the cumulative RMT distribution functions which are directly derived from Eq.(38),
(58)  
where is the incomplete gamma function. To test the goodness of fit between the computed cumulative distribution function and that predicted by RMT, we use KolmogorovSmirnov statistic which is the least upper bound of all pointwise differences [24],
(59) 
where is the empirical cumulative distribution function defined by
(60) 
where are the computed values of a random sample which consists of elements. If is ”close” to , the observed must be so small that it falls within the range of possible fluctuations which depend on the size of the random sample. For the random variable for any , it can be shown that the probability of is given by [25]
(61) 
where
(62) 
From observed maximum difference , we obtain
the significant level which is equal to the probability
of . If is found to be large enough,
the hypothesis is verified. The
computed cumulative distributions of and the
GOE() prediction for four examples are plotted
in figure 6,
and the maximum difference and the significant levels for
and are
shown in table 3. Note that the last digit in
is not guaranteed, since Eq.(61) is an asymptotic formula.
We see from figure 6 and
table 3 that the agreement with GOE prediction is
remarkably good. The Gaussian behavior for highly excited states
is naturally expected as the semiclassical wavefunctions must reflect
the chaotic nature in the corresponding classical systems. However,
the Gaussian behavior for lowlying modes is not apparent at all.
It is possible that the
nongaussianity behavior is rather prominent as these modes have
fluctuations on large scale and that reflects the
pure quantum mechanical
behavior^{6}^{6}6The typical scale of angular fluctuation of mode
is approximately and the typical scale of radial fluctuation of
mode is approximately given by
.. Nevertheless, our numerical results serve to strengthen
the hypothesis that the expansion coefficients behave as Gaussian
pseudorandom numbers even for lowlying modes.
k  

5.41  9.42  23.4%  4.29  51.8% 
5.79  3.88  99.3%  4.29  52.1% 
6.81  5.95  78.6%  3.47  77.8% 
6.89  9.34  24.1%  2.36  98.8% 
7.12  7.12  57.2%  2.59  96.8% 
7.69  10.23  15.9%  4.84  36.7% 
8.30  6.38  70.8%  3.88  65.0% 
8.60  5.63  83.8%  2.09  99.7% 
8.73  9.46  22.9%  2.78  94.3% 
9.26a  7.21  55.6%  3.46  78.1% 
9.26b  5.99  77.9%  6.11  13.5% 
9.76  7.41  52.0%  4.57  43.8% 
9.91  8.90  29.3%  3.23  84.7% 
9.99  4.27  98.0%  3.72  70.0% 
ave.  7.23  56.3%  3.69  68.8% 
Next we examine the randomness of ’s. Because ’s are actually determined by the Helmholtz equation, it is appropriate
to describe ’s as pseudorandom numbers. We apply the
run test for testing randomness (see [24]).
Suppose that we have observations
of the random variable and m observations of the random variable
. The combination of those variables into observations
placed in ascending order of magnitude yields
xxx yy xx yyy x y xx yy,
where denotes an observation of and denotes an observation
of . Each underlined group which consists of successive
values of or is called run. The statistics of
number of runs are used for testing whether and have
the same the distribution function. Regardless of the type of the distribution
function, the run number is known to behave as Gaussian random
numbers in the limit .
The run test is also used as a test
for randomness. Let be the observed values of a
random variable . For simplicity, assume that is even. The
median divides the observed values into a lower and an upper half.
It is represented as if it falls below the median, and it is
represented as if it falls above the median. For instance, a
sequence
UUU L UU LL U U L UU,
has 8 numbers of runs (). The critical region for testing
the hypothesis of randomness is of the form or where
and is
readily given by the Gaussian distribution function. The significant
level is the probability of or . As the
KolmogolovSmirnov test, is given by the observed .
The run numbers and the significant levels are shown
in table 4. High significant levels are again observed except
for the one at for . As the
corresponding is larger than the averaged value, this may be due
to the cyclic effect. On the whole, it is
concluded that ’s behave as if they are random variables.
k  

5.41  62  85.5%  185  67.3% 
5.79  58  46.5%  174  39.9% 
6.81  69  14.4%  196  11.4% 
6.89  60  71.5%  168  14.0% 
7.12  69  14.4%  184  75.2% 
7.69  57  36.1%  191  29.2% 
8.30  63  71.5%  177  59.8% 
8.60  59  58.3%  184  75.2% 
8.73  70  10.0%  201  3.5% 
9.26a  56  27.3%  177  59.8% 
9.26b  55  20.1%  182  91.6% 
9.76  59  58.3%  182  91.6% 
9.91  70  10.0%  196  11.4% 
9.99  58  46.5%  179  75.2% 
ave.  61.8  40.7%  184  50.4% 
5 Summary
In this paper, we have demonstrated that the DBEM is eminently suitable for
computing eigenmodes on CH spaces and we obtain some lowlying eigenmodes on
a CH space called Thurston manifold which is the second smallest
in the known CH manifolds and we have studied the statistical properties of
these eigenmodes.
The lowlying eigenmodes are expanded in terms of eigenmodes on the
pseudosphere, and we find that the expansion coefficients behave as
if they are Gaussian random numbers. Why are they so random even for
lowlying modes? It should be pointed out that the
randomness of the expansion coefficients for lowlying eigenmodes is
not the property of the eigenmodes themselves but rather the property of the
images of eigenmodes on the whole universal covering space, since the
fluctuation scales for lowlying eigenmodes are comparable to the
the size of the fundamental domain. We conjecture that the origin of
the random behavior
of eigenmodes comes from the almost randomly
distributed images of a set of points in the universal covering
space.
Computation of eigenmodes is essential in simulating the CMB in
CH cosmological models. As the DBEM needs only a set of facetoface
identification maps and the discretization of the corresponding
fundamental domain, it can be applied to other CH spaces
straightforwardly. However, the computation of the modes with small
fluctuation scale is still a difficult task as the
number of modes increases as .
Nevertheless, the contribution of the modes
with small fluctuation to the temperature correlation of
CMB can be estimated by assuming that
the expansion coefficients for excited states ()
also behave as Gaussian
pseudorandom numbers as well as that for lowlying modes.
The assumption is numerically
confirmed in a twodimensional CH model [10].
If the observed Gaussian
pseudorandomness is found to be the universal behavior in CH spaces
for lowlying modes as well as excited modes,
the origin of the gaussianity of the CMB fluctuations can be partially
explained. This is because the amplitude of the CMB fluctuation is written in
terms of:
1. expansion coefficients of the initial fluctuation
in terms of eigenmodes on the CH space
2. expansion coefficients of eigenmodes on the CH space
that are extended onto the whole pseudosphere
in terms of eigenmodes on the pseudosphere.
Acknowledgments
I would like to thank my advisor Professor Kenji Tomita for his many helpful discussions and continuous encouragement. I would also like to thank Professor Toshiro Matsumoto for his extensive advice on the boundary element methods and Professor Jeff Weeks and the Geometry Center in University of Minnesota for providing me the data of CH spaces. Numerical computation in this work was carried out by HP Exemplar at the Yukawa Institute Computer Facility. I am supported by JSPS Research Fellowships for Young Scientists, and this work is supported partially by GrantinAid for Scientific Research Fund (No.9809834).
References
 [1]
 [2] Stevens D, Scott D and Silk J 1993 Phys. Rev. Lett. 71 20
 [3] OliveiraCosta A de, Smoot G F and Starobinsky A A 1996 Astrophys. J. 468 457
 [4] Levin J L, Barrow, J D, Bunn E F and Silk J 1997 Phys. Rev. Lett. 79 974
 [5] Levin J L, Scannapieco E and Silk J astroph/9802021
 [6] Seriu M 1996 Phys. Rev. D 53 6902
 [7] Aurich R and Steiner F 1989 Physica D 39 169
 [8] Grunewald F and Huntebrinker F 1996 Experimental Mathematics 5 57
 [9] Cornish N J and Turok N G 1998 Class.Quant.Grav. 15 2699
 [10] Aurich R and Steiner F 1993 Physica D 64 185
 [11] Haake F and Zyczkowski K 1990 Phys. Rev. A 42 1013
 [12] Brody T A, Flores J, French J B, Mello P A, Pandy A and Wong S S M 1981 Rev. Mod. Phys. 53 385
 [13] Brebbia C A and Dominguez J 1992 Boundary Elements  An introductory Course Second Edition (London: Computational Mechanics Publications)
 [14] Tanaka M, Sladek V and Sladek J 1994 ”Regularization techniques applied to boundary element methods” Applied Mechanics Reviews Vol.47 457
 [15] Tabata M 1994 Bibunhouteishikino Suchikaihou II(Numerical Solution of Differential Equations II) (Tokyo: Iwanami Shoten Publishers)
 [16] Johnson C 1987 Numerical Solution of Partial Differential Equations by the Finite Element Methods (Cambridge: Cambridge University Press)
 [17] Elstrodt E, Grunewald F and Mennicke J 1983 Uspekhi Mat.Nauk 38 No.1 119
 [18] Tomaschitz R 1989 Physica D 34 42
 [19] Telles J C F 1989 International Journal for Numerical Methods in Engineering 24 959
 [20] Hayami K and Brebbia C A 1987 Boundary Elements IX Vol.1 375
 [21] Mehta M L 1991 Random Matrices (New York:Academic Press)
 [22] Bohigas O 1991 ”Random Matrix Theories and Chaotic Dynamics” in Proceedings of the 1989 Les Houches School on Chaos and Quantum Physics ed Giannoni A et al (Amsterdam: Elsevier)
 [23] Harrison E R 1967 Rev. Mod. Phys. 39 No.4 862
 [24] Hogg R V and Tanis E A 1977 Probability and Statistical Inference (New York: Macmillan Publishing Co., Inc.)
 [25] Birnbaum Z W 1962 Introduction to Probability and Mathematical Statistics (New York: Haper & Brothers)

[26]
Thurston W P 1979 The geometry and topology of three manifolds
Princeton Lecture Notes (Now available via internet:
http://www.msri.org/gt3m/)  [27] Thurston W P 1997 ThreeDimensional Geometry and Topology (Princeton University Press)
 [28] Beardon A E 1983 The Geometry of Discrete Groups Graduate Texts in Mathematics 91 (New York: SpringerVerlag)
 [29] Fomenko A T and Kunii T L 1997 Topological Modeling for Visualization (Tokyo: SpringerVerlag)

[30]
Gabai D, Meyerhoff R and Thurston N 1996 Homotopy Hyperbolic
3Manifolds Are Hyperbolic MSRI preprint 1996058 available at
http://www.msri.org/ 
[31]
Weeks J R 1998 SnapPea: A computer program for creating and studying
hyperbolic 3manifolds available at
http://www.geom.umn.edu/software/download/snappea.html  [32] Weeks J R 1985 PhD thesis (Princeton: Princeton University)
 [33] Matveev S V and Fomenko A T 1988 Russian Math Surveys 43 No.1 3
 [34] Hodgson C and Weeks J 1994 Experimental Math. 261
Appendix A Boundary integral equation
Here, we derive the boundary integral equation (11) in section 1. For simplicity,