%Paper: hep-lat/9209005
%From: uwolff%surya11.cern.ch@AVM.CC.FSU.EDU (Ulli Wolff)
%Date: Wed, 2 Sep 92 11:43:12 +0200

% content: LATEX-file, 1 figure as separate .ps file appended;
% their beginning can be located by searching for   %Fig.1:
% the .ps file has to be removed before invoking LATEX
%
%
\documentstyle[12pt]{article}
\begin{document}

\newcommand{\be}{ \begin{equation}}
\newcommand{\ee}{ \end{equation}  }
\newcommand{\bi}{\bibitem}
\newcommand{\s}{ \sigma }
\newcommand{\rag}{ \rangle }
\newcommand{\lag}{ \langle }
\newcommand{\rd}{ \mbox{\rm d} }
\newcommand{\rD}{ \mbox{\rm D} }
\newcommand{\re}{ \mbox{\rm e} }

\renewcommand{\thefootnote}{\fnsymbol{footnote}}

%\renewcommand{\topfraction}{0.8}
%\renewcommand{\bottomfraction}{0.8}
\renewcommand{\floatpagefraction}{0.8}
\begin{titlepage}
\begin{flushright} CERN-TH 6636/92 \end{flushright}
\vspace{1ex}
\vspace{0.75cm}
\begin{center}
{ \LARGE High Precision Simulation Techniques \\
for Lattice Field Theory
\footnote[2]{Talk given at the 4th International Conference on Computational
Physics PC92,
Prague, August 1992}}
\vspace{0.75cm}

{\large Ulli Wolff \\
CERN, Theory Division,\\
CH-1211 Gen\`eve 23, Switzerland}
\date{}
%\maketitle
\end{center}
\vspace{1cm}
\thispagestyle{empty}
\begin{abstract}\normalsize
An overview is given over the recently developed and now widely used
Monte Carlo algorithms with reduced or eliminated critical slowing down.
The basic techniques are overrelaxation, cluster algorithms and
multigrid methods. With these tools one is able to probe much closer than
before the universal continuum behavior of field theories on the lattice.
\end{abstract}
\vspace{2.cm}
\begin{flushleft} CERN-TH 6636/92 \\ September 1992 \end{flushleft}
\vspace{1ex}
\end{titlepage}
%
%
\section{Introduction}
\noindent
The interest in statistical systems close to criticality is
shared by a large community including  condensed
matter physicists and particle theorists.  An important tool in the
study of such systems are numerical experiments in the form of Monte
Carlo simulations which complement analytical results that are
available for special systems and limiting cases. One of the
central problems in such simulations is the degradation
of most known simulation algorithms as criticality --- the cutoff or
continuum limit in quantum field theories describing particles ---
is approached.  The interdisciplinary search of improved techniques
is an on-going effort, but it already has yielded some very positive
results in recent years some of which we want to briefly review here.
Our focus will be on spin systems with variables belonging
to continuous manifolds ($\s$-models) and on pure lattice gauge theory.
This selection thus leaves out the vast field of discrete Ising-like
systems$^1$
as well as the enormous problems (and potential gains from algorithmic
improvement)
faced in QCD simulations with fermions.$^2$

The problem of critical slowing down (CSD)
near criticality is schematically described by the dynamical
scaling law
\be \label{z}
\tau = c \; \xi^z.
\ee
Here $\xi$ is some physical scale or correlation length in lattice units
and $\tau$ is a time scale in number of iterations.
A $\tau$ and corresponding $z$ can characterize
the time scale for equilibration or
the rate at which
statistically independent estimates for observables of interest are
produced. Then for the expectation value of $A$,
\be \label{avA}
\lag A[s] \rag = \frac{1}{Z}\int \rD s \; \re^{-\beta H[s(x)]} \; A[s],
\ee
the accuracy in estimating $A$ improves as
\be
\delta_A =  \sqrt{\mbox{variance}_A \times 2\tau/\#\mbox{ of iterations}}\;.
\ee
In (\ref{avA}) $\rD s$ = $\prod_x \rd\mu(s(x))$ means integration
over all fields or spins $s(x)$ at lattice sites $x$
with some measure, $Z$ is the partition function and
$\beta H[s]$ is the inverse bare coupling or temperature times
some action or Hamiltonian.

The algorithms to be described here have been designed to lower $z$
from its ``traditional'' value of about two for standard local Metropolis
methods to $z \simeq 1$ or even $z \simeq 0$ in some cases.
They thus improve the efficiency
of simulations by {\it one to two powers} of $\xi$.

Beyond the introduction this article will continue with a section on
the technique of embedded variables (common to many algorithms),
and short descriptions of cluster algorithms, multigrid methods and
overrelaxation followed by some conclusions.


\section{Updating Embedded Variables}
\noindent
Let us imagine some group of transformations $T\in G$ that
act on the configurations $s \hat{=} \{s(x)\}$,
\be
s \rightarrow T s,
\ee
such that the measure is left invariant, $\rD s = \rD T s$.
Such transformations can be local, $T s = \{ t(x) s(x) \}$,
and then $t(x)$ is a field similar to $s$ itself with values
in local group factors that make up $G$. One may however also
consider more global changes of $s$, where $T$ acts for instance on some
cluster of spins from $s$ with one and the same rotation.
For a given configuration $s$ we can now think of a statistical system
with configurations $T\in G$ and an effective or induced Boltzmann
factor $\exp(-\beta H[T s])$ whose simulation can also be considered.
Assume now that we have a Monte Carlo algorithm for this system
which is characterized by transition probabilities $w(s;T_1 \to T_2)$
that preserve $\exp(-\beta H[Ts])$ and
depend on s only through the effective Boltzmann factor such that
 $w(Ts;T_1\to T_2)$ = $w(s;T_1 T \to T_2 T)$ holds.
Then, if ervalid algorithm for the original field $s(x)$:
\begin{itemize}
\item for momentarily fixed $s=s_1$ put $T_1=Id$ as initial configuration,
\item update $T_1\to T_2$ with the $w$--algorithm,
\item take $s_2=T_2s_1$ as a new $s$.
\end{itemize}

To prove this assertion we first note that the overall transition
probability is given by
\footnote[2]{For simplicity we take a discrete group $G$ here;
in the continuous case the sum over all elements has to be replaced by
an invariant group integration.}
\be
W(s_1 \to s_2) = \sum_T w(s;Id \to T) \; \delta(s_2-Ts_1).
\ee
We have to show that $W$ preserves the Boltzmann weight $\exp(-\beta H[s])$.
To this end we transform
\begin{eqnarray}
\int \rD s_1 \; \re^{-\beta H[s_1]} \; W(s_1 \to s_2) &=& \nonumber\\
\int \rD s_1 \sum_{T'} \re^{-\beta H[s_1]} \;
 w(s_1;Id \to T') \; \delta(s_2-T's_1) \; &=& \nonumber\\
\frac{1}{|G|}\int \rD s_1 \sum_{T,T'} \re^{-\beta H[Ts_1]} \;
 w(s_1;T \to T') \; \delta(s_2-T's_1) \; &=& \nonumber\\
\frac{1}{|G|}\int \rD s_1 \sum_{T'} \re^{-\beta H[T's_1]} \;
 \delta(s_2-T's_1) \; &=&
\re^{-\beta H[s_2]}.
\end{eqnarray}
To arrive at the third line  changes of variables $s_1\to Ts_1$
and $T' \to T'T^{-1}$ are carried out and $T$ is averaged over.
In the last line  we absorbed $T'$ into $s_1$ and the $\delta$-function
is used.

In many cases,
in particular if $G$ is a lower dimensional manifold than the original
configuration space, the moves described are not ergodic. This can often
be improved by using families of different embeddings between which
one switches
in a random or deterministic order. In cases where this is not sufficient
one can always blend in conventional update steps to achieve ergodicity.



\section{Cluster Algorithms}
\noindent
Cluster algorithms$^3$ for continuous spins
are an example of the successful use ovariables. Their drawback has been up to
 now that
the $O(n)$ invariant
$n$-vector models are the only continuous variable $\s$-models
where they are powerful.\footnote[2]{There are principal
reasons for this limitation which
come close to a no-go theorem.$^4$}
Here however, according to accumulated
numerical evidence$^5$, they really achieve $z\simeq 0$ with even
a very small coefficient $c$ in (\ref{z}). With the additional advantage
of variance reduced estimators for Green functions, these models have
become an ideal testing ground for nonperturbative physics$^6$,
in particular
in two dimensions, where they are asymptotically free for $n\ge 3$.
Also the xy-model ($n=2$) is of great interest$^7$ as it allows for
checks of the Kosterlitz Thouless scenario.
We therefore now specialize our general setup, Eq.(\ref{avA}), to
the $O(n)$ models, where
\be
s(x)\in R^n, \quad \rD s= \prod_x \rd^n s(x) \; \delta(s(x)^2-1), \quad
-H[s] = \sum_{<xy>} s(x) \cdot s(y).
\ee

The key to efficient cluster algorithms is the embedding of Ising
spins and the use of global update techniques for them.
A family of embeddings is labelled by an $n$-component unit vector
$r$ of the same type as the local spins. The group involved
is $G_r \simeq Z_2^{\# \mbox{sites}}$ corresponding to one
Ising spin $\s(x)=\pm 1$ on each site, $T \leftrightarrow \{\s(x)\}$.
They act as reflections,
\be \label{r}
s'=Ts, \quad s'(x)=s_{\perp}(x) + \s(x) s_{\|}(x),
\ee
where $\perp$ and  $\|$ refer to the $r$-direction.
The induced Hamiltonian
\be
-H[Ts] = \sum_{<xy>} r\cdot s(x) \; r\cdot s(y) \; \s(x) \s(y)
 +\mbox{terms independent of } \s
\ee
is recognized to describe a {\it ferromagnetic} Ising model with random
bond strengths $J_{<xy>} = r\cdot s(x) \; r\cdot s(y)$. When they are
multiplied along closed loops the result is always positive or zero
due to their factorized form, which shows the absence of frustration.
We also note that there is no magnetic field coupling to $\s$.
This can be regarded as being due to the reflections in (\ref{r})
being part of the $O(n)$ global symmetry which leads to a global
$Z_2$ invariance for $\s$.
For these embedded random systems the Swendsen-Wang algorithm$^8$
or its single cluster variant$^3$ work very well, actually better
than in the standard Ising model where some remaining CSD is still
detectable.

The algorithm just described is ergodic if $r$ is chosen at random
such that all directions can appear. In practical realizations it is
usually convenient to always take $r=(1,0,\ldots,0)$. Then
the embedding requires fewer operations, and ergodicity is restored
by globally $O(n)$-transforming the whole configuration with a random
rotation or reflection after a certain number of updates.

\section{Multigrid Techniques}
\noindent
Multigrid Monte Carlo (MGMC) techniques have been proposed$^9$ to
eliminate CSD by allowing for efficient moves of a critical
system on all scales. They work well on nearly gaussian systems$^{10}$ as
they do in the related problem of linear difference equations
where the method has made its first appearance.
Here we will not review the truly recursive  MGMC$^{11}$
approach but a simpler unigrid variant that was proposed recently.$^{12}$
It has become the optimal method for $\s$-models other than the
$O(n)$ family.

We now present the method for the $O(3)$-model.$^{12}$
The actual updates are performed here (and in realizations for
other $\s$-models) on embedded $U(1)$ spins of the xy-model type.
To an $O(3)$ generator corresponding to some rotation, as for
example
\be
i\lambda = \left( \begin{array}{ccc} 0 & 1 & 0 \\ -1~ & 0 & 0 \\
           0 & 0 & 0 \end{array} \right) \; ,
\ee
we couple local angles $\alpha(x)$ by
\be
Ts \hat{=} \left\{ \re^{-i\alpha(x) \lambda} s(x) \right\}.
\ee
This induces a Hamiltonian for $\alpha(x)$,
\be
-H[Ts]= \sum_{<xy>} \mbox{Re} \left( J_{<xy>} \re^{i(\alpha(x)-\alpha(y)}
\right) +\mbox{terms independent of } \alpha
\ee
with complex bond strengths
\be
J_{<xy>} = (s^1(x)+i s^2(x)) (s^1(y)-i s^2(y))=J_{<yx>}^{\ast},
\ee
where $s^a, a=1,2,3,$ denotes the components of $s$.
The random bonds generated for the embedded xy-model are again
ferromagnetic due to their factorized form. As for $U(1)$ gauge fields
their orientation has to be properly taken into account, of course.

We now need an algorithm for $\alpha(x)$. For the version of MGMC
of Ref.~12 elementary moves are performed on
$B\times B \times \ldots \times B$
subblocks of the lattice. For such blocks one has a priory fixed
profiles of kernels $K(x)$ that vanish outside the block and are smooth.
Possible choices
in arbitrary dimensions are pyramids or the lowest mode
sine waves with the block as a cavity. The kernels  appear in the
nonlocal Metropolis proposals
\be
\{\alpha(x)\} \to \{ \alpha(x) + \Psi K(x) \}.
\ee
These are accepted or rejected in the usual way, and $\Psi$ is drawn
from a symmetric distribution with a width adjusted for reasonable
acceptance. For a lattice of length $L$, which should be a power of two,
one has to hierarchically cover the lattice with blocks of sizes
$B = L/2, L/4, \ldots , 1$ with the last choice producing just local
updates. It has turned out to be important that either the superimposed
block lattice or equivalently
the field configuration is randomly translated between
updates, such that the block boundaries are not at fixed positions.
Furthermore the generator $\lambda$ is randomized to achieve ergodicity.

For $O(n)$ models this MGMC algorithm is presumably
inferior to cluster methods, although detailed comparisons would
be somewhat hardware dependent.
The importance of the technique derives however from the fact,
that it can also be used for $CP^n$-models and for $SU(n)$ valued
spins. The main change is that appropriate generators
(from $U(n)$ and $SU(n)$ respectively) have to be substituted for $\lambda$.
The resulting embedded $U(1)$ system now {\it can have frustration}
depending on the configuration of the ``background'' spins $s$.
Practical tests for the $SU(3) \times SU(3)$ chiral model and for
the $CP^3$ system have shown that these frustrations do not
seriously slow down the evolution.$^{13}$ In all three cases $z = 0.2(1)$
has been found.

\section{Overrelaxation}
\noindent
In contrast to the two previous algorithms overrelaxation (OR) achieves
an improvement (down to $z \approx 1$) with still local updates
only, and hence it is as fast or even somewhat faster per sweep
than standard algorithms.
It also immediately carries over to gauge fields.
We now present OR in its ``hybrid'' form that found many
applications recently rather than the original ``tunable'' version$^{14}$.
We consider the local update problem at site $x$ with local Boltzmann
weight
\be \label{loBo}
\re^{\beta s(x) \cdot M(x)} \quad
\mbox{where } \quad M(x) = \sum_{y=
\mbox{\scriptsize n.n. of } x} s(y),
\ee
again for an $O(n)$ model spin for illustration.
A local heatbath procedure amounts to choosing a new $s(x)$ independently
of the old one with the weight (\ref{loBo}). For OR we need additional
microcanonical steps producing transitions $s_1(x) \to s_2(x)$ such that
\begin{itemize}
\item $s_1(x), s_2(x)$ have the same local weight (\ref{loBo}) and thus
the energy is unchanged,
\item $s_1(x), s_2(x)$ are as far from each other as possible.
\end{itemize}
For our example this principle leads to the change
\be
s_2(x) = -s_1(x) + 2 \frac{M \; M\cdot s_1(x)}{M^2}.
\ee
Experience has shown that both for vectorization and for the nonuniversal
it is best to group the local updates such as to do a maximum number
of independent ones in parallel, which usually amounts to checkerboard
ordering.
Now an OR iteration consists of $N$ microcanonical sweeps followed
by one heatbath or other standard ergodic step.
{}From exact gaussian analysis$^{15}$ and from numerical experiments$^{16}$
it is known that to achieve $z\approx 1$ it is necessary do let $N$ grow
proportional to $\xi$ as criticality is approached. Often $N=\xi$
is a reasonably good trial value. The goal is to achieve a roughly constant
autocorrelation time when measured in iterations which
implies $z\approx 1$ when referring to sweeps.

In particular, this kind of OR is the present method of choice
for SU(2) lattice gauge fields
(either fundamental or embedded to move SU(3) fields).
The local problem for a link variable
$U_{\mu}(x) \in SU(2)$ coincides with the $O(4)$ case when $SU(2)$
matrices are expanded in terms of the unit- and the Pauli-matrices.

\begin{figure}
\vspace{9.5cm}
\caption{Autocorrelation time in sweeps for
four dimensional $SU(2)$ lattice gauge theory
in a finite size scaling limit.}
\end{figure}
We close with the example of
 a recent simulation of the $SU(2)$ gauge theory$^{17}$
 where it was possible to
determine the relation between $\tau$ and a scale in lattice units
for a whole range of scales as shown in Fig.~1.
In this study a renormalized coupling
constant was held fixed which is equivalent to a finite size scaling limit
at fixed $L\sqrt{K}$ with the
string tension $K$ assuming the r\^ole
of a correlation length.
The time $\tau$ refers to independently estimating the
renormalized coupling.
The line in the plot represents a fit with the
form (\ref{z}) giving $z=1.0(1)$ and $c=0.5(1)$.
For further details on algorithmic and physical aspects of this work
we have to refer to Ref.~17.
\section{Conclusions}
\noindent
We have presented some of the accelerated algorithms for the Monte
Carlo simulation of spin
and gauge theories that have been discovered and tested in recent years.
As a result, critical continuum behavior can be studied now much more
accurately, especially in two and three dimensions. In four dimensions,
which is the most interesting case for particle physics, the situation
will become similar as larger systems will be studied on future computers.
In the presence of Goldstones modes, the new techniques are crucial
already now.

When algorithms  of the multigrid type with $z < 1$ will be found for
gauge theories, it will strongly depend on the overhead inflicted,
at which system size they really become profitable.
In simulations of the $CP^3$ model, for instance, it has been found
that on vector machines the crossover between OR and MGMC occurs only
for correlation lengths $\xi = 20 \ldots 30$.$^{13,18}$

\vspace{1ex}


\noindent
{\large\bf Acknowledgements} \newline
I would like to thank Martin Hasenbusch and Steffen Meyer for
correspondence and discussions on their multigrid technique.
%
\newpage
\begin{thebibliography}{99}
\bi{1.} P. G. Lauwers, ``Multiscale Monte Carlo algorithms in statistical
mechanics and quantum field theory'', Lecture notes of the 3rd European
Conference on Multigrid Methods, Bonn, Germany, October 1990 and
preprint Bonn-HE-90-16.
\bi{2.} R. C. Brower, R. G. Edwards, C. Rebbi and E. Vicari,
{ Nucl. Phys.} { B366} (1991) 689;
T. Kalkreuter, G. Mack and M. Speh,
``Blockspin and multigrid for staggered fermions in nonAbelian gauge fields'',
preprint DESY-91-070, July 1991;
 P. G. Lauwers, R. Ben-Av and S. Solomon,
{ Nucl. Phys.} { B374} (1992) 249.
\bi{3.} U. Wolff,
{ Phys. Rev. Lett.} { 62} (1989) 361 and
{ Nucl.Phys. B} (Proc. Suppl.) { 17} (1990) 93.
\bi{4.} S. Caracciolo, R. G. Edwards, A. Pelissetto and A. D. Sokal,
`` Wolff type embedding algorithms for general nonlinear sigma models'',
Tallahassee preprint FSU-SCRI-92-65;
A. D. Sokal, { Nucl.Phys. B} (Proc. Suppl.) { 20} (1991) 55.
\bi{5.} R. Edwards and A. Sokal,
{ Phys. Rev.} { D40} (1989) 1374;
U. Wolff, { Nucl. Phys.} { B334} (1990) 581.
\bi{6.}  M. L\"uscher, P. Weisz and U. Wolff,
{ Nucl. Phys.} { B359} (1991) 221.
\bi{7.} U. Wolff,
{ Nucl. Phys.} { B322} (1989) 759;
M. Hasenbusch, M. Marcu and K. Pinn, ``High precision verification of the
Kosterlitz Thouless scenario'', preprint DESY-92-010, January 1992.
\bi{8.} R. H. Swendsen and J.-S. Wang,
{ Phys. Rev. Lett.} { 58} (1987) 86.
\bi{9.}  G. Parisi, in { Progress in Gauge Field Theory },
(Carg\`ese 1983), ed. G.~t'Hooft et. al. (Plenum, New York, 1984).
\bi{10.} J. Goodman and A. D. Sokal, { Phys. Rev. Lett.}
{  56} (1986) 1015.
\bi{11.} J. Goodman and A. D. Sokal, { Phys. Rev.}
 { D40} (1989) 2035.
\bi{12.} M. Hasenbusch,  S. Meyer and G. Mack,
{ Nucl.Phys. B} (Proc. Suppl.) { 20} (1991) 110.
\bi{13.} M. Hasenbusch and  S. Meyer, { Phys. Rev. Lett.}
 { 68}
 (1992) 435 and
``Testing Accelerated Algorithms
 in the Lattice $CP^3$ Model'',
Kaiserslautern preprint Th 91-23, January 1992.
\bi{14.} S. L. Adler,
{ Phys. Rev.} { D23} (1981) 2901.
\bi{15.}  U. Wolff, ``Dynamics of Hybrid Overrelaxation in the
Gaussian Model'',
CERN preprint 6408/92, March 1992, to appear in
{ Phys. Lett. B.}
\bi{16.} J. Apostolakis, G. C. Fox and C. F. Baillie,
{ Nucl. Phys. B} (Proc. Suppl.) { 20} (1991) 678.
\bi{17.} M. L\"uscher, R. Sommer, P. Weisz and U. Wolff,
 ``Computation of the Running Coupling in the SU(2)
Yang-Mills Theory'' ,
preprint CERN-TH.6566/92, July 1992.
\bi{18.}
U. Wolff, { Phys. Lett.} { B284} (1992) 94.
\end{thebibliography}

\end{document}

%Fig.1:
%!PS-Adobe-2.0
%%Creator: gnuplot
%%DocumentFonts: Courier
%%BoundingBox: 50 50 403 554
%%Pages: (atend)
%%EndComments
/gnudict 40 dict def
gnudict begin
/Color false def
/gnulinewidth 5.000 def
/vshift -73 def
/dl {10 mul} def
/hpt 31.5 def
/vpt 31.5 def
/vpt2 vpt 2 mul def
/hpt2 hpt 2 mul def
/Lshow { currentpoint stroke moveto
  0 vshift rmoveto show } def
/Rshow { currentpoint stroke moveto
  dup stringwidth pop neg vshift rmoveto show } def
/Cshow { currentpoint stroke moveto
  dup stringwidth pop -2 div vshift rmoveto show } def
/DL { Color {setrgbcolor [] 0 setdash pop}
 {pop pop pop 0 setdash} ifelse } def
/BL { stroke gnulinewidth 2 mul setlinewidth } def
/AL { stroke gnulinewidth 2 div setlinewidth } def
/PL { stroke gnulinewidth setlinewidth } def
/LTb { BL [] 0 0 0 DL } def
/LTa { AL [1 dl 2 dl] 0 setdash 0 0 0 setrgbcolor } def
/LT0 { PL [] 0 1 0 DL } def
/LT1 { PL [4 dl 2 dl] 0 0 1 DL } def
/LT2 { PL [2 dl 3 dl] 1 0 0 DL } def
/LT3 { PL [1 dl 1.5 dl] 1 0 1 DL } def
/LT4 { PL [5 dl 2 dl 1 dl 2 dl] 0 1 1 DL } def
/LT5 { PL [4 dl 3 dl 1 dl 3 dl] 1 1 0 DL } def
/LT6 { PL [2 dl 2 dl 2 dl 4 dl] 0 0 0 DL } def
/LT7 { PL [2 dl 2 dl 2 dl 2 dl 2 dl 4 dl] 1 0.3 0 DL } def
/LT8 { PL [2 dl 2 dl 2 dl 2 dl 2 dl 2 dl 2 dl 4 dl] 0.5 0.5 0.5 DL } def
/M {moveto} def
/L {lineto} def
/P { stroke [] 0 setdash
  currentlinewidth 2 div sub moveto
  0 currentlinewidth rlineto  stroke } def
/D { stroke [] 0 setdash  2 copy  vpt add moveto
  hpt neg vpt neg rlineto  hpt vpt neg rlineto
  hpt vpt rlineto  hpt neg vpt rlineto  closepath  stroke
  P  } def
/A { stroke [] 0 setdash  vpt sub moveto  0 vpt2 rlineto
  currentpoint stroke moveto
  hpt neg vpt neg rmoveto  hpt2 0 rlineto stroke
  } def
/B { stroke [] 0 setdash  2 copy  exch hpt sub exch vpt add moveto
  0 vpt2 neg rlineto  hpt2 0 rlineto  0 vpt2 rlineto
  hpt2 neg 0 rlineto  closepath  stroke
  P  } def
/C { stroke [] 0 setdash  exch hpt sub exch vpt add moveto
  hpt2 vpt2 neg rlineto  currentpoint  stroke  moveto
  hpt2 neg 0 rmoveto  hpt2 vpt2 rlineto stroke  } def
/T { stroke [] 0 setdash  2 copy  vpt 1.12 mul add moveto
  hpt neg vpt -1.62 mul rlineto
  hpt 2 mul 0 rlineto
  hpt neg vpt 1.62 mul rlineto  closepath  stroke
  P  } def
/S { 2 copy A C} def
end
%%EndProlog
%%Page: 1 1
gnudict begin
gsave
50 50 translate
0.070 0.070 scale
90 rotate
0 -5040 translate
0 setgray
/Courier findfont 220 scalefont setfont
newpath
LTa
LTb
1584 1808 M
1647 1808 L
6873 1808 M
6810 1808 L
1452 1808 M
(2) Rshow
1584 2415 M
1647 2415 L
6873 2415 M
6810 2415 L
1452 2415 M
() Rshow
1584 2845 M
1647 2845 L
6873 2845 M
6810 2845 L
1452 2845 M
() Rshow
1584 3179 M
1647 3179 L
6873 3179 M
6810 3179 L
1452 3179 M
(5) Rshow
1584 3452 M
1647 3452 L
6873 3452 M
6810 3452 L
1452 3452 M
() Rshow
1584 3683 M
1647 3683 L
6873 3683 M
6810 3683 L
1452 3683 M
() Rshow
1584 3882 M
1647 3882 L
6873 3882 M
6810 3882 L
1452 3882 M
() Rshow
1584 4059 M
1647 4059 L
6873 4059 M
6810 4059 L
1452 4059 M
() Rshow
1584 4216 M
1647 4216 L
6873 4216 M
6810 4216 L
1452 4216 M
(10) Rshow
2183 771 M
2183 834 L
2183 4489 M
2183 4426 L
2183 551 M
(6) Cshow
3129 771 M
3129 834 L
3129 4489 M
3129 4426 L
3129 551 M
(8) Cshow
3862 771 M
3862 834 L
3862 4489 M
3862 4426 L
3862 551 M
(10) Cshow
4461 771 M
4461 834 L
4461 4489 M
4461 4426 L
4461 551 M
(12) Cshow
4968 771 M
4968 834 L
4968 4489 M
4968 4426 L
4968 551 M
(14) Cshow
5406 771 M
5406 834 L
5406 4489 M
5406 4426 L
5406 551 M
(16) Cshow
5793 771 M
5793 834 L
5793 4489 M
5793 4426 L
5793 551 M
(18) Cshow
6140 771 M
6140 834 L
6140 4489 M
6140 4426 L
6140 551 M
(20) Cshow
LTb
1584 771 M
6873 771 L
6873 4489 L
1584 4489 L
1584 771 L
220 2630 M
currentpoint gsave translate 90 rotate 0 0 moveto
(tau) Cshow
grestore
4228 331 M
(L/a) Cshow
4228 4709 M
(Fig.1) Cshow
LT0
LT0
3862 3035 D
4461 3434 D
4968 3613 D
5406 3873 D
6140 4180 D
2183 2240 D
3129 2692 D
3862 3109 D
3862 2984 M
3862 3083 L
3831 2984 M
3893 2984 L
3831 3083 M
3893 3083 L
4461 3393 M
4461 3474 L
4430 3393 M
4492 3393 L
4430 3474 M
4492 3474 L
4968 3574 M
4968 3650 L
4937 3574 M
4999 3574 L
4937 3650 M
4999 3650 L
5406 3833 M
5406 3912 L
5375 3833 M
5437 3833 L
5375 3912 M
5437 3912 L
6140 4133 M
6140 4225 L
6109 4133 M
6171 4133 L
6109 4225 M
6171 4225 L
2183 2218 M
2183 2263 L
2152 2218 M
2214 2218 L
2152 2263 M
2214 2263 L
3129 2667 M
3129 2716 L
3098 2667 M
3160 2667 L
3098 2716 M
3160 2716 L
3862 3077 M
3862 3140 L
3831 3077 M
3893 3077 L
3831 3140 M
3893 3140 L
LT1
1584 2102 M
1584 2102 L
1637 2127 L
1691 2151 L
1744 2175 L
1798 2200 L
1851 2224 L
1905 2248 L
1958 2273 L
2011 2297 L
2065 2321 L
2118 2346 L
2172 2370 L
2225 2394 L
2279 2418 L
2332 2443 L
2385 2467 L
2439 2491 L
2492 2516 L
2546 2540 L
2599 2564 L
2652 2589 L
2706 2613 L
2759 2637 L
2813 2662 L
2866 2686 L
2920 2710 L
2973 2735 L
3026 2759 L
3080 2783 L
3133 2808 L
3187 2832 L
3240 2856 L
3294 2881 L
3347 2905 L
3400 2929 L
3454 2954 L
3507 2978 L
3561 3002 L
3614 3027 L
3668 3051 L
3721 3075 L
3774 3100 L
3828 3124 L
3881 3148 L
3935 3173 L
3988 3197 L
4042 3221 L
4095 3245 L
4148 3270 L
4202 3294 L
4255 3318 L
4309 3343 L
4362 3367 L
4415 3391 L
4469 3416 L
4522 3440 L
4576 3464 L
4629 3489 L
4683 3513 L
4736 3537 L
4789 3562 L
4843 3586 L
4896 3610 L
4950 3635 L
5003 3659 L
5057 3683 L
5110 3708 L
5163 3732 L
5217 3756 L
5270 3781 L
5324 3805 L
5377 3829 L
5431 3854 L
5484 3878 L
5537 3902 L
5591 3927 L
5644 3951 L
5698 3975 L
5751 4000 L
5805 4024 L
5858 4048 L
5911 4073 L
5965 4097 L
6018 4121 L
6072 4145 L
6125 4170 L
6178 4194 L
6232 4218 L
6285 4243 L
6339 4267 L
6392 4291 L
6446 4316 L
6499 4340 L
6552 4364 L
6606 4389 L
6659 4413 L
6713 4437 L
6766 4462 L
6820 4486 L
6826 4489 L
stroke
grestore
end
showpage
%%Trailer
%%Pages: 1


