\section{Introduction - Motivation}
     One  of the most interesting open questions in non-linear
dynamics is the nature and evolution of transport in the
chaotic phase space regions of perturbed integrable Hamiltonian systems. Most
of the work in this field has been done for 2-D systems or
for the standard map, following two different
approaches. In the first approach, which may be designated as
{\it microscopic}, one studies the homoclinic and heteroclinic tangle
of stable and unstable
manifolds inside the chaotic region. This gives a ``complete''
qualitative picture of transport in phase space, since a chaotic
trajectory follows the unstable manifolds of the unstable
periodic orbits, jumping from one to another near the
heteroclinic points. On the other hand, invariant tori surrounding
stable periodic orbits act as ``barriers'' in the transport
process, modifying at the same time the topological structure of
phase space. However, it is not presently well understood how the
results from this kind of study can be used in a quantitative
description of the evolution of an ensemble of trajectories.

     In the second approach, which may be designated as {\it
macroscopic}, one assumes {\it a priori} that transport is a pure
random phenomenon (a Markovian process). Then the transport of an
ensemble of trajectories may be considered as ``classical
diffusion'' in action space.
Unfortunately, it is well known that transport in Hamiltonian
systems cannot be considered as a pure Markovian process,
due to the problem of
{\it stickiness}. Moreover the process may follow L\'evy rather
than classical statistics (Shlesinger et
al., 1993). In this case the macroscopic approach still works,
provided that the process is considered as {\it fractal} diffusion
in {\it normal} space (Zaslavsky, 1994). In this way it becomes evident that
L\'evy-like statistics are closely related to the self-similarity of
the phase space. This formalism results in a differential equation
with fractional partial derivatives, whose solution is formidably
difficult, even in the simple case of 2-D systems.

     There are rather few attempts to study transport in
Hamiltonian systems with more than two degrees of freedom, mainly
because the above mentioned approaches for 2-D systems cannot be
directly generalized. It is therefore important to note
that there exists yet a third approach for the description of
transport, which may be implemented in a straightforward way in
systems with more than two degrees of freedom, and this is {\it
normal diffusion} in {\it fractal space} (e.g. see West and
Deering, 1994, and references therein). The fractality of phase space
has been already demonstrated for 2-D systems and has been
associated with the self similarity of the hierarchical
structure of island families on a surface of section (e.g. see
Zaslavsky, 1994, BenKadda et al., 1997).

In the case of more than
two degrees of freedom the situation is more complicated, since
the topological structure of phase space in any region depends on
the number of local integrals of motion. Recently
Varvoglis et al. (1997) have found indications
of phase space fractality in the same model 3-D Hamiltonian
system studied in the present work, but they were unable, due to the
particular method of study they had selected,
to identify actual fractal or multifractal sets in phase space.
They have conjectured that the fractality is due to the
presence of {\it quasi-barriers} (Varvoglis and Anastasiadis,
1996), i.e. geometrical objects of lower than full dimensionality
(depending on the  number of local integrals of motion), such as
periodic orbits, invariant tori around them and cantori (Aubry,
1978; Percival, 1979).

From the above discussion it is obvious that the importance of the
fractality of phase space in the
diffusion process lies in the nature and in the rate of the diffusion:
In a simply connected space,
the diffusion is classical and follows Fick's law (transport
proportional to $\sqrt (t)$). In a fractal
space the diffusion is non-classical (presumably a L\'evy process)
and, in the case of a
Hamiltonian system, usually it has a lower than the Fickian rate.

On the other hand, the distribution of
periodic orbits perpendicular to the $x-y$ plane
in the present 3-D model potential
has been studied by Barbanis and Contopoulos (1995) (henceforth referred
to as BC) and Barbanis (1996). It was found that it is not
random: the perpendicular crossings (p.c.'s)
of orbits with multiplicities forming Farey sequences with the $x-y$ plane
are arranged along lines,
named in the present paper basic-lines (BL) and Farey tree lines (FTL).
These lines are not continuous but
they possess gaps; two of them have already been mentioned in BC.
Since the presence of Farey sequences
implies a form of self-similarity, it is possible that the BL's and FTL's
may be related to the above-mentioned quasi-barriers and that the gaps
may be related to any fractal properties of phase space.

In the present paper we use the model potential studied by
BC in order to test the conjecture put forward
by Varvoglis et al. (1997) about the origin of the fractality of phase space.
Furthermore, through this, we attempt to assess
the relation between
the topological structure of the phase space and transport, on the one hand,
and the system of BL's and FTL's, on the other.
Although this potential can model a true galaxy only locally,
its study is expected, nevertheless, to contribute to the understanding
of galactic evolution problems, based on the fact that evolution is intimately
connected to transport in phase and configuration space.
Transport in a chaotic region of a perturbed
integrable dynamical system, in turn, may be viewed as
(stochastic) diffusion, which depends on the value of the
Lyapunov Characteristic Number as well as on the topological structure of the
phase space. We believe that similar
behaviour would characterise any 3-D perturbed integrable dynamical system,
as are most of the realistic model galactic potentials.

The
purpose of this paper is threefold:

(a) to study in finer detail the
distribution of the periodic orbits along the BL's and FTL's.

(b) to examine whether these lines are connected, in any way,
to transport in phase space and are, thus, related to the
``quasi-barriers'' and

  to confirm the existence of true
irregular families of periodic orbits in a 3-D Hamiltonian system.

     This paper is organized as follows. In the next section we
present, for reasons of completeness, some remarks on
previous results. In Section 3 we show that the LCN's may ``map''
the fractality of phase space. In Section 4 we present our
results on the interconnection and discontinuities of BL's and FTL's.
In Section 5 we show numerically the existence of irregular families
of periodic orbits.
Finally in Section 6 we summarise and discuss our results.

\section{Remarks on previous results}

\subsection{Basic lines}
     The systematic exploration of the phase space structure of
3-D Hamiltonian systems was initiated by Magnenat (1982) and later
by Contopoulos \& Barbanis (1989), BC and Barbanis (1996),
through the study of periodic orbits of unit mass test
particles in a 3-D model potential, corresponding to the Hamiltonian
\begin{eqnarray}
    H = & {1 \over 2} \left( p_x^2 + p_y^2 + p_z^2 \right) +
    {1 \over 2} \left( A x^2 + B y^2 + C z^2 \right) \\
    & {-\epsilon x z^2 - \eta y z^2 = h} \nonumber
\end{eqnarray}
with a particular set of parameters (A = 0.9, B = 0.4, C =
0.225, $\epsilon = 0.560$, $\eta = 0.20$ and $h = 0.00765$). The
specific values for the parameters A, B, C were
selected by Magnenat (1982) because the 2-D system $x-z$ for
$\sqrt{A/C}$ = 2 is topologically
equivalent to the Inner Lindblad Resonance and the 2-D system $y-z$ for
$\sqrt{B/C}$ = 4/3
is a well studied 2-D dynamical system with extended chaotic regions.
For this set of parameters it was found in
BC that the distribution of the p.c.'s
of the periodic orbits is not accidental: the crossings are arranged along
particular lines on the $\xbar-\ybar$ plane (where $\xbar =
A^{1/2}x$, $\ybar = B^{1/2}y$), as it is evident in Fig.~\ref{figM1}
(taken from Fig.~5 of BC). Each orbit has either one p.c., marked
by a (x), or two p.c.'s, marked by (+). In this paper, for clarity
reasons, we use a point ($\cdot$) instead of (x). Each p.c. is
designated by a number (or a number and a letter), representing
the multiplicity of the orbit; different orbits with the same
multiplicity are differentiated by a  prime. The lines are formed
by the p.c.'s of the orbits, whose multiplicities form arithmetic
progressions with increment the multiplicity of the limit orbit.
They are termed in this paper {\it basic} if they connect orbits
of the unperturbed system (e.g. 1a, 1b in Fig. 1) or orbits with low
multiplicity, $m$ (e.g. $m \leq 5$).  For example the basic
sequences
\begin{itemize}
\item[] \hspace{0.5cm} 3(+), 5($\cdot$), 7(+), 9($\cdot$), ...2c(+)
\item[] and 3(+), 5$'$($\cdot$), 7(+), 9$'$($\cdot$), ...2c(+),
\end{itemize}
which form the lines C and C$'$, start with
an orbit of multiplicity 3 and increment 2, which is the
multiplicity of the limit orbit 2c.

\subsection{Farey tree lines}
     In BC it was recognized that the distribution of p.c.'s on the
$\xbar-\ybar$ plane presents some self-similar features. To begin with,  to
each basic sequence of orbits corresponds a number of new sequences of higher
order. Each member of a new sequence is a periodic orbit having
as multiplicity the sum of the multiplicities of two consecutive
orbits of the generating sequence. For example between the first
two orbits of the basic sequence C, i.e. 3(+) and 5($\cdot$), there is
the orbit 8(+). Two second order sequences are formed between 3
and 8 as well as between 8 and 5, i.e.
\begin{eqnarray}
   & 8(+), 11(\cdot), 14(+), 17(\cdot), 20(+), 23(\cdot), ...3(\cdot)
\nonumber\\
\mathrm{and} & \\
 & 8(+), 13(+), 18(+), 23(+), 28(+), ...5(\cdot)\nonumber
\end{eqnarray}
Both new sequences have as first member the orbit 8, but they
have different limit orbits and increments (i.e. 3 for the first
and 5 for the second, respectively). In the same way we can form
sequences of even higher order.
We call the orbits of these higher order
sequences {\it Farey tree orbits} and the lines on which
they lie {\it Farey tree lines}, because such sequences are
similar to the Farey sequences discussed by Niven and Zuckerman (1960).
For a good review on Farey trees see Efthymiopoulos et al. (1997).

The periodic orbits are used, in the present paper, as an additional tool
for probing the
topological structure of the phase space. In particular, periodic
 orbits of high multiplicities are
needed in order to search for self-similar structures in the BL's
and FTL's for a wide interval of
multiplicities (see Avnir et al., 1998). The importance of the
non-basic (Farey) lines lies exactly
on the above line of reasoning, as well as in the fact that they
might be necessary in order to
calculate correctly an appropriate diffusion coefficient. It should
be noted that we have restricted
our study to the periodic orbits crossing perpendicular the x-y plane
for simplicity reasons.

\begin{figure}
  \resizebox{\hsize}{!}{\includegraphics*{7919.f1}}
  \caption{The majority of the perpendicular crossings
  of the periodic orbits of the Hamiltonian (1) lie on particular lines on
  the $\xbar-\ybar$ plane inside  the boundary circle
  $\xbar^2+\ybar^2=2h$. Most of the lines are concentrated around
  the basic orbits
  1a and 1b. There is also a noticeable  spiral structure with focal
  point the orbit 1a. Four spirals  illustrate the area of this
  structure. Note that in this figure (reproduced from BC, Fig. 5) a different
  notation is used than
  in the rest of the present paper. Each
  symbol, $(\times)$ or $(+)$, shows a periodic orbit with one or two
  crossing points, respectively.}
  \label{figM1}
\end{figure}

\begin{figure}
\resizebox{\hsize}{!}{\includegraphics*{7919.f2}}
\caption{A segment of the BL between the p.c.'s of the periodic
orbits 3 and 5 (indicated by a bold line) and the network of various FTL's
(corresponding to the upper central part of Fig.~\ref{figM1}, line C).
The nonexistent part of C between 8 and 14
is indicated by a dashed line. The contour lines
of constant LCN's  $\lambda = 10^{-4}$ are superimposed and the
regions of $\lambda  > 10^{-4}$ are gray shaded}
\label{figM2}
\end{figure}

\subsection{Irregular orbits}
All periodic orbits of the same multiplicity that can be found through the
continuous variation of a parameter in the equations of motion belong
to the same
family. The graph giving one of the initial conditions of the periodic
orbits as a function of the parameter is the {\it characteristic}
of this family. Families, which do not exist for values of the
perturbing parameter below a critical value, appear only as pairs, are well
known in 2-D dynamical systems and are called {\it irregular
families}  (Contopoulos 1970; Barbanis 1986).
However, their existence has not been proven for 3-D systems.
Pairs of families that have been found in such systems without
a direct connection to a basic orbit, they were all of small
multiplicities and
it turns out that they have an indirect connection to a basic orbit
through other families bifurcating from it
(Barbanis, 1996).
It is desirable to know whether
families of periodic orbits in 3-D systems (presumably of high multiplicities)
do exist without {\it any} connection to a basic orbit.
Their presence might play a certain role in the fractality of phase
space and, therefore,
in trajectory transport.

\subsection{Lyapunov Numbers}
Varvoglis et al. (1997), while
studying numerically transport phenomena in the trajectories of
the Hamiltonian (1), have found, in an indirect way, that the phase
space shows a
fractal structure. In particular they presented numerical
evidence that the function $dV/dt$ has multifractal properties,
where $V(t)$ denotes essentially the coarse-grained volume of phase
space visited by a trajectory up to
time $t$. They found, also, that the evolution of the function $dV/dt$
is closely
related to the evolution of the function
\begin{equation}
\chi (t) = { 1\over t} \ln {d(t)\over d(0)}
\end{equation}
whose limit, as $t \rightarrow \infty$, is defined as the LCN. The
function $\chi (t)$ shows plateaus in the time intervals where $dV/dt$
is close to
zero and steeply increasing transient
segments in the time intervals where $dV/dt$ takes large positive
values. This fact shows that a trajectory is confined successively
in regions of phase space where the LCN converges to a limit and, therefore,
the fractal properties of $dV/dt$ should arise from the
(apparently) self-similar distribution of periodic orbits and
other quasi-barriers in phase space (Zaslavsky, 1994).

It should be emphasised that, in the present paper, LCN's are not
calculated in order to estimate
``rates of diffusion'', since LCN's alone
cannot describe completely this phenomenon in the case where there exist
invariant tori of considerable measure. We do,
however, calculate LCN's in
order to ``probe'' the phase space structure, i.e. to delineate, in an
independent from other methods
way, the ``topology'' of the phase space. This new method is based on
the property that
trajectories starting in a certain region, bounded by
``quasi-barriers'', are characterised by a certain
``Local Lyapunov Number''. By drawing the
contour lines of LCN's, one
draws, essentially, the geometric boundary of the region bounded by
the ``quasi-barriers''.

\section{LCN's and fractal properties}

The ``classical'' calculation of LCN's
requires a continuous integration of the corresponding
trajectories, until the function $\chi (t)$ has
reached a plateau. It is, therefore, obvious that this method cannot
be used in a ``mass production''
procedure. For this reason we decided to use a different approach, by
calculating the value of the function $\chi (t)$
at various integration time intervals, $\Delta t$, up to $\Delta t = 3\,10^5$.
Comparing the results for different $\Delta t's$ we have found that they
are qualitatively the same (i.e. they form level lines with increasing
complexity as one goes to finer details),
provided that $\Delta t > 10^5$.
So, even if $\chi (t)$ does not reach a plateau, its value at the end of the
corresponding time interval $\Delta t > 10^5$
gives a correct estimate, at least to order of magnitude, of the
``degree of stochasticity'' of the trajectory in the phase space region
restricted by the quasi-barriers. We designate this value
as the {\it Fixed Time LCN}. In what follows we use, for
convenience, the notation LCN. Through this kind of LCN's we attempt to
probe the topological structure on the $\xbar-\ybar$ plane as
follows:

We draw a dense mesh on the $\xbar-\ybar$ plane (for $0 < \xbar <
0.020$ and $0.010 < \ybar < 0.045$ at intervals $\Delta \xbar$ =
$\Delta \ybar = 0.0005$) and we use the nodes of this mesh as
initial conditions for trajectories starting perpendicular to this plane.
We calculate  the LCN of each trajectory at the end of various
time intervals and we draw the contour lines of constant LCN
on the $\xbar - \ybar$ plane at the LCN value $\lambda = 10^{-4}$. In
Fig.~\ref{figM2}
we show our results for the LCN's at
the end of a time interval $\Delta t = 3\,10^5$.
Fig.~\ref{figM2}  contains only a small part of the studied region. The
structure of the
contour lines in the whole region as well as their variation with $\Delta t$
is still under investigation and
will be the topic of another publication.

The emergent
picture is very interesting. One can immediately see that the
studied area may be divided into two regions: Most of it is
characterized by small values of LCN's (white regions);
inside the white regions one
can find small elongated regions characterized by high values
of LCN's (gray regions).
It should be noted that the position and elongation of these
regions is closely related to the position and the direction of BL's and
FTL's (see next Section). This fact may be interpreted in the
following way. Since a considerable fraction of the periodic orbits
along these lines are unstable, the exponential
divergence of
trajectories in the area around them is governed mainly by the
unstable manifolds of the orbits-members of the lines. The situation
is reminiscent of a {\it fractal repeller} as defined in Gaspard and
Baras (1995), i.e. a set of countably
infinite unstable periodic orbits, which has zero Lebesgue measure
and a finite Hausdorff
dimension. In our case the orbits-members of BL's and FTL's may be
considered as forming a sort of an incomplete
(since not all periodic orbits are unstable) {\it fractal repeller}.
In other words
the BL's and FTL's are closely associated with the ``quasi-barriers'',
which are geometrical objects that separate ``loosely'' different
phase space regions with
different Local LCN, as it has been discussed
by Varvoglis et al. (1997) and Tsiganis et al. (1998).

Since (I) there is evidence that BL's and FTL's have
fractal properties and (ii) there is a relation
between the BL's and FTL's, on one hand, and the
LCN's on the other, it should be interesting to examine
whether the contour lines have fractal properties as well. To do so we
calculate the LCN's in the region $-0.02< \xbar < 0.01, 0.05 < \ybar < 0.10$
using a coarse ($\Delta \xbar = \Delta \ybar = 0.005$)
and a fine mesh ($\Delta \xbar = \Delta \ybar = 0.001$). This
region was selected because it was studied, although with a much coarser mesh
( $\Delta \xbar = \Delta \ybar = 0.01$)
by Contopoulos and Barbanis (1989, Fig. 9).
We then plot
both contour lines at $\lambda = 10^{-4}$, on the same graph, as
shown in Fig.~\ref{figM3}.
It is easy to see that, while the dashed line, corresponding to the
coarser mesh, seems rather smooth, the continuous one, corresponding to the
finer mesh, shows ``tongues'' that oscillate
about the dashed curve, a picture reminiscent of the classical
examples of fractal curves. Of course this cannot be taken as a {\it
proof}, even numerical, that contour lines are fractal curves, since
the available data are not exhaustive (see also Avnir et al.,
1998). However, we believe that Fig.~\ref{figM3}, if considered
together with the
fact that LCN's reflect the distribution of the unstable manifolds of
unstable periodic orbits, is a strong indication that the contour
lines have, indeed, fractal properties.

\section{BL and FTL are not simple}

As it is evident from Fig.~\ref{figM1} and the relevant discussion in BC
and Barbanis (1996), the BL's formed by the
lowest order sequences look continuous and smooth.
Because of this it had been assumed that the
lines appearing in Fig.~\ref{figM1} contain all the higher order
sequences described in Section 2 as well, so that any fractal properties,
arising from the Farey tree character of the p.c. sequences, are
restricted {\it along} the lines and not {\it across} them. However,
the existence of two small
gaps in a BL, noticed by BC, gave the first evidence that, somehow,
this picture is
not correct. In order to examine, therefore, in depth the
situation, we proceed in the calculation of
periodic orbits belonging to several higher order sequences than those
appearing in
Fig.~\ref{figM1}.

We focus our interest on a small region, between the orbits 3(+) and
5($\cdot$) of the
line C (BC, Fig. 8), since it is in this area that the first
noticeable gap was observed.
Between these two orbits lies
the orbit 8(+). Therefore two second order sequences are formed with orbit
8 as first member and increments 3 and 5 respectively, i.e.
the sequences (2). In Fig.~\ref{figM2} we can see that the
first sequence
presents a large gap between orbits 8 and 14, where higher
order FTL's, branching away from the BL, may be observed.
Furthermore, as we discuss below, orbit
11 (which lies between orbits 8 and 14, instead of being part of the BL)
belongs to a higher order FTL. In contrast the second
sequence does not show any noticeable gaps.
The reason for this difference becomes obvious, if one
considers the BL's and FTL's in association with the degree of stochasticity
of the phase space where these orbits lie, estimated through
the calculation of LCN's.
In Fig.~\ref{figM2} one can immediately notice
that the first sequence crosses a region of high LCN's,
while the second is confined in a region of low LCN's.

Let us now examine what happens to the higher order sequences
between 11 and 8, on one hand, and 11 and 14 on the other, i.e.
11($\cdot$), 19(+), 27($\cdot$), 35(+), ..., 8(+) and
14(+), 25(+), 36(+), ..., 11($\cdot$).
As it is evident from Fig.~\ref{figM2},
the FTL corresponding to the first sequence is torn and split into
two branches, one beginning from
orbit 11 going upwards and the other beginning from orbit 8
going downwards.
Similar is the situation with the second sequence. The FTL is torn and split
into two branches, one going upwards from orbit 14 and the other
going downwards from orbit 11. The branch going downwards from orbit 11
is connected to the FTL
of the first sequence going upwards from orbit 11, resulting in a
composite S-shaped line (the three leftmost continuous
thin lines in Fig.~\ref{figM2}).
In this way we see that
the gap between orbits 8 and 14 is, in a way, ``filled''
with FTL's corresponding to higher order sequences. This is an example
of a case where
higher order FTL's do not lie on the ``parent'' BL.

It is interesting to note that the branch of the second FTL
going upwards from orbit 11 is developed inside a
narrow strip of ordered motion.

In Fig.~\ref{figM2} we notice also two ``vertical'' higher order FTL's
that emanate from the orbits 57($\cdot$) and 17($\cdot$), corresponding
to higher order sequences with increment 40. The first of them,
the one emanating from
orbit 57($\cdot$), is vertical to the local gradient of the contour lines
and lies in a region of low LCN's. The second one, the one emanating
from orbit 17($\cdot$), is vertical to the local gradient of the contour lines
as well but it seems to span a {\it high}
LCN's region. The most probable explanation is that the line
lies, in fact, in a narrow strip of low LCN's, with a width less
than the size of the mesh used. This is something that should
be expected for a geometrical object with fractal properties
and has been indeed encountered in some other
regions of the studied area.

From the above discussion emerges an intuitively appealing picture:
BL's and FTL's are, in a sense, lines representing some ordered features of the
dynamical system, which are deformed in the neighbourhood of the
{\it chaotic seas}
and tend to run perpendicular to the local gradient of the
contour lines of constant
LCN's. Therefore we understand that the structure of
the BL's and FTL's on the  $\xbar - \ybar$ plane is considerably more
complicated than it was assumed in BC, a fact which, as we showed in Section 3,
seems to play an important role in the nature of transport in the phase
space of Hamiltonian (1).

\section {Existence of irregular families}

In two previous papers (BC; Barbanis 1996) the authors investigated
the bifurcation and the evolution of known periodic orbits belonging to some
BL's. It was found that a small number of these families
form pairs that do not have any direct or indirect connection with the
periodic orbits of the unperturbed system. However, for a given multiplicity,
there are many other bifurcating families, which may play the role of
connecting bridges between those pairs.

\begin{figure}
\resizebox{\hsize}{!}{\includegraphics*{7919.f3}}
\caption{The contour lines of constant LCN at the value $\lambda = 10^{-4}$
for a coarse
mesh ($\Delta \xbar = \Delta \ybar = 0.005$, dashed line) and a fine mesh
($\Delta \xbar = \Delta \ybar = 0.001$, solid line). Regions of
 $\lambda > 10^{-4}$
and $\lambda < 10^{-4}$ are marked with a H and an L, respectively}
\label{figM3}
\end{figure}

Following these considerations, we studied the evolution of
the bifurcations of some {\it specific} multiplicities from the
basic orbits 1a and 1b. We have selected the multiplicities 11
and 29 for three reasons: (I) 11 and 29 are prime numbers, so that
there are no families of these multiplicities resulting from lower
multiplicity families, except from those with multiplicity one.
(ii) The number of bifurcating families is neither small nor
great. (iii) There are a lot of orbits of these multiplicities on
the various BL's and FTL's for $\epsilon = 0.560$ and $\eta = 0.20$.

The detailed discussion in Appendix A and the examples in Appendix B lead to
the conclusion that pairs of families
without any direct or indirect connection to a basic orbit of the
unperturbed system (e.g. the six pairs
appearing in Fig.~\ref{figA1}) exist in a 3-D system. This happens
mainly for high multiplicities, exactly
as one would expect considering the analogous case of irregular families
in 2-D systems. It should be noted also that, besides the above mentioned
{\it classical} case of irregular families, in this work we have found a new
kind of irregular families, which we call a {\it double pair}, i.e
a pair consisting of two single pairs (see Fig.~\ref{figB6} and its discussion
in Appendix B).

\section{Summary and conclusions}
Barbanis and Contopoulos (1995), in studying the
distribution of the p.c.'s of the periodic orbits in a well known 3-D model
potential with the
$\xbar-\ybar$ plane, have discovered a noteworthy order. Except for a
couple of gaps, the p.c.'s
of periodic orbits, whose multiplicities form Farey-tree sequences,
are arranged
along continuous and smooth-looking lines. The p.c.'s along these
lines have some fractal
properties, arising from the fact that they belong to Farey-tree sequences
of various orders. In the present work we tried to study the properties
of these lines in association with the values of the LCN's of the area
span by them.

We have found numerical evidence that the contour lines of
LCN's
show fractal properties. We have also found numerical evidence that the
splitting of initially continuous and smooth-looking
BL's and FTL's, as well as the ensuing formation of gaps, is not an exception,
as was implied in BC. As a rule, it is observed in regions
of the $\xbar - \ybar$ plane characterised by high values of the LCN's and
appears not only in BL's, but in higher order
FTL's as well.

Since our results are only numerical and involve only a small
number of
periodic orbits, the above results cannot be considered as firm proof
that the distribution of BL's and FTL's as well as the contour lines
of LCN's have, beyond any doubt, fractal properties.
However we feel that there is enough evidence that the above two
geometrical objects
(BL's and FTL's on the one hand and contour lines on the other)
are closely related and that, furthermore, they show strong indications
of fractal behaviour. Therefore we think that the study of some
selected regions of the $\xbar - \ybar$ should be done in more detail,
by calculating LCN's in even finer mesh
and by comparing the contour lines to higher order FTL's,
in order to establish to a higher level of confidence the fractal
nature of both geometrical objects.

As far as irregular periodic orbits are concerned, we have found
that the majority of the periodic orbits forming BL's and FTL's either
bifurcate directly from 1a or 1b or they are connected indirectly with
them. However,
by considering families of periodic orbits with large multiplicities,
we find that some
of them do not have any connection with 1a or 1b but they form a pair with
another family. In the present work we have encountered also the more complex
situation of a {\it double pair}, i.e. a pair whose members are pairs too.

\begin{acknowledgements}
The authors would like to thank Prof. L. Martinet for a critical
reading of the first version, which improved the final form of the
paper. H.V. would also like to thank Prof. J.H. Seiradakis and
Mr. K. Tsiganis for
several helpful discussions.
\end{acknowledgements}


\section{Introduction}

The dominant characteristics of the optical spectrum of starburst galaxies are
its nebular emission lines. These lines are formed in the surrounding
interstellar medium of the starburst that is photoionized by photons emitted by
stars with mass higher than 10 M$\odot$. These massive stars dominate the
ultraviolet and contribute significantly to the optical continuum, that comes
mainly from A, B and O stars. The most conspicuous features of the spectra of
these early stars are their strong hydrogen Balmer and neutral helium absorption
lines. In the spectrum of a starburst, the nebular emission lines are
superimposed on the stellar absorption lines, and usually emission dominates
over absorption. However, the contribution of the underlying absorption becomes
increasingly important for higher terms of the Balmer series (H$\gamma$,
H$\delta$,...) and some of the HeI lines ($\lambda$ 4471, 4387, 4026, ...). Very
often, the effect of the stellar population is also seen as absorption wings in
the H$\beta$, H$\gamma$ and H$\delta$ lines. The observation of these lines with
high spectral resolution makes it possible to estimate the contribution of the
underlying absorption. On the other hand, the analysis of the profiles
high-order of the Balmer series, which are dominated by absorption instead of
emission, can yield information on the properties of the stellar content of the
starbursts and on their evolutionary stage. Therefore, high-resolution
absorption profiles of the high-term Balmer and He I lines covering a wide range
in effective temperature and gravity are needed to predict the composite stellar
spectra of starburst galaxies.

Observational and synthetic stellar spectra are available in the literature.
However, the observational atlases have only intermediate or low spectral
resolution (Burstein et al. 1984; Jacoby, Hunter \& Christian 1984; Walborn \& Fitzpatrick 1990;  Cananzi, Augarde, \& Lequeux 1993) and the synthetic stellar
libraries do not cover the high-order Balmer series and HeI lines (Auer \&
Mihalas 1972; Kurucz 1979, 1993). For these reasons, a grid of stellar
atmospheres and synthetic spectra in the wavelength range from 3700 to 5000 \AA\
with a sampling of 0.3 \AA\ and covering 50000 to 4000 K has been computed using
the code of Hubeny (1988). In this paper, we present the grid of synthetic
profiles, which are compared with observations and synthetic profiles computed
by Auer \& Mihalas (1972) and Kurucz (1993). In a companion paper (Gonz\'alez
Delgado, Leitherer \& Heckman 1999; hereafter paper II), evolutionary synthesis
profiles of H and He absorption lines for star forming regions are presented.

\section{Stellar library}

\subsection{The Grid}

The grid includes the synthetic stellar spectra of the most relevant H and He I
lines from 3700 to 5000 \AA, in five different wavelength ranges (Table 1).
H$\epsilon$ is not synthesized because it is coincident with the CaII H line,
and this line shows a very strong interstellar component in individual stars
and in the integrated spectrum of a galaxy. HeII lines, He I $\lambda$5876, HeI $\lambda$6678 and H$\alpha$ are not synthesized because in hot stars the
profiles of these lines are affected by stellar winds (Gabler et al. 1989;
Herrero et al. 1992; Bianchi et al. 1994). However, HeI at wavelengths shorter
than 5000 \AA, H$\beta$ and the higher order-terms of the Balmer series are
only partially (mainly the core of the line but not the wings) or not at all
affected by sphericity and winds (Gabler et al. 1989). The spectra span a range
of effective temperature from 4000 to 50000 K, with a variable step from 500 K
to 5000 K, and a surface gravity log$g$=0.0 to 5.0 with a step of 0.5 (Table
2). The metallicity is solar.

The spectra are generated in three different stages with a set of computer
programs developed by Hubeny and collaborators (Hubeny 1988; Hubeny \& Lanz
1995a; Hubeny, Lanz, \& Jeffery 1995). First, the stellar atmosphere is
calculated; then, the stellar spectrum is synthesized, and finally the
instrumental and rotational convolution are performed.

\subsection{The model atmospheres of the grid}

 To generate a synthetic stellar spectrum, a model atmosphere is needed. For
$T_{eff}\geq$25000 K, the atmosphere is produced using version 193 of the
program TLUSTY (Hubeny 1988; Hubeny \& Lanz 1995a,b); for $T_{eff}\leq$25000 K
we use a Kurucz (1993) LTE atmosphere. TLUSTY calculates a plane-parallel,
horizontally homogeneous model stellar atmosphere in radiative and hydrostatic
equilibrium. The program allows departures from local thermodynamic equilibrium
and metal line blanketing, using the hybrid complete linearization and
accelerated lambda iteration (CL/ALI) method.

To reduce computational time, non-blanketed non-LTE (NLTE) models are computed
with TLUSTY. H and He are considered explicitly. The population of their levels
(9 atomic energy levels of HI, 14 levels of HeI and 14 levels of HeII are
considered) are determined by solving the corresponding statistical equilibrium
equations. 25 additional atoms and ions contribute to the total number of
particles and to the total charge, but not to the opacity.

The NLTE models are computed in three stages. First, an LTE-gray atmosphere is
generated. Here, the opacity is independent of wavelength, and the populations
of the energy levels are calculated assuming the local value of the temperature
and electron density. This model is used as a starting approximation for the
LTE model. Finally, the NLTE model is computed, where the gas and the radiation
are coupled. Here, departures from LTE are allowed for 39 energy levels.
Convection is suppressed in all the models. A depth-independent turbulence
velocity of 2 km s$^{-1}$ is assumed. Doppler broadening is assumed for all the line
transitions. The properties of the atmosphere are calculated at 54 depth
points.

We use Kurucz (1993) LTE atmospheres for $T_{eff}\leq$25000 K because for stars
cooler than B1, NLTE effects are not very important, These models are line
blanketed, and we take the models with a turbulence velocity of 2  km s$^{-1}$
and solar metallicity.

\subsection{The synthetic profiles}

The synthetic spectra are computed with the program SYNSPEC (Hubeny, Lanz, \&
Jeffery 1995). This program reads the input model atmosphere, either calculated
by TLUSTY or from the Kurucz models, and solves the radiative transfer
equation, wavelength by wavelength in a specified spectral range. The program
also uses an input line list that contains the transitions in the six specified
wavelength ranges that we synthesize here. The line list has a format similar
to the Kurucz \& Peytremann (1975) tables.

The continuum opacity is calculated exactly the same way as it was done in the
atmosphere model. For $T_{eff}\geq$ 25000 K, only continuum opacities from the
atomic energy levels of H and He are considered. The opacity sources are: 1)
photoionizations from all the explicit levels (39 levels in TLUSTY models); 2)
free-free opacity for all the explicit ions (HI, HeI and HeII in TLUSTY models
and H, Mg, Al, Si and Fe in the Kurucz models); 3) electron scattering. The
line opacity is calculated from the line list. Although the structure of the
atmosphere (for $T_{eff}\geq$ 25000 K) is calculated with only H and He, the
synthetic spectra are calculated assuming line blanketing from elements with
atomic number Z$\leq$28. The atoms and transitions that were treated in NLTE in
the atmosphere models are also treated in NLTE by SYNSPEC. For those lines that
originate between levels for which the population is calculated in LTE, SYNSPEC
uses an approximate NLTE, based on the second-order escape probability theory
(Hubeny et al. 1986).

The line profiles have the form of a Voigt function and take into account the
effect of natural, Stark, Van der Waals and thermal Doppler broadening. The
broadening parameters for H and He are those tabulated by Vidal, Cooper, \&
Smith (1973) for the first four members of the Balmer series, and up to H10 by
Butler (private communication). For HeI, the line broadening tables for
$\lambda$4471 are from Barnard, Cooper, \& Smith (1974) and those for
$\lambda$4026, 4387, 4922 are from Shamey (1969). The maximum distance of two
neighboring frequency  points for evaluating the spectrum is 0.01 \AA. A
turbulence velocity of 2 km s$^{-1}$ is assumed.

Finally, the program ROTIN performs the rotational and instrumental
convolution. A FWHM instrumental profile of 0.01 \AA\ and a rotational velocity
of 100 km s$^{-1}$ are assumed which is the typical rotational velocities  in
massive stars (Conti \& Ebbets, 1977). The final sampling of the spectra is 0.3
\AA.

\section{Results}

The synthetic spectra of the grid are available for retrieval at
http://www.iaa.es/ae/e2.html and http://www.stsci.edu/science/starburst/.
Figure 1 compares the spectra for $T_{eff}$=25000 K and log$ g$= 4.0 computed
with TLUSTY+SYNSPEC assuming NLTE with those computed with SYNSPEC and Kurucz
LTE atmospheres. The two spectra are essentially identical; this justifies the
use of LTE stellar atmosphere models for $T_{eff}\leq$25000 K. Figure 2 shows
the spectra of a typical O ($T_{eff}$=40000 K, log$ g$=4.0), B ($T_{eff}$=20000
K, log$ g$=4.0) and A ($T_{eff}$=10000 K, log$ g$=4.0) star. The most important
lines in the spectrum are labelled. The equivalent width of the most important
H and HeI lines have been measured in the synthetic spectra. Several
measurements were done with different spectral windows. This allows a
calibration of the contribution of weaker lines to the spectral index that
characterizes each of the Balmer lines. The spectral index was measured from
the continuum equal to 1 (this represents the real equivalent width of the
lines), but also from a pseudo-continuum that was determined by fitting a first
order polynomial (except for the spectral range 3700-3900 \AA, for which we
used a third order) to the windows defined in Table 3. This spectral index
simulates the equivalent width we can measure in observed stellar spectra.
Tables 4 to 7 show the equivalent widths of H$\delta$, H8, HeI $\lambda$4471
and HeI $\lambda$3819, respectively.

Figure 3 shows the equivalent widths of H$\beta$, H$\delta$ and two of the
high-order terms of the Balmer series (H8 and H9) as a function of $T_{eff}$
for main sequence stars (log$ g$= 4.0). The maximum equivalent width occurs at
$\sim$9000 K,  corresponding to an early-type A star. The plot indicates that
the high-order Balmer series lines show, like the lower terms, a strong
dependence on $T_{eff}$ and on gravity (Figure 4); thus, they are an efficient
tool for the determination of the fundamental stellar parameters. Figure 5
shows the equivalent width of the HeI lines as a function of $T_{eff}$ for main
sequence stars (log$ g$= 4.0) (see also Tables 6 and 7 for HeI $\lambda$4471
and  HeI $\lambda$3819, respectively). They also show a strong dependence on
$T_{eff}$ and gravity (Figure 6). The maximum is at $\sim$20000 K, corresponding
to an early-type B star. The increase of the equivalent width at
$T_{eff}\leq$10000 K is not due to HeI absorption, but to the contribution of
some metal lines (FeII $\lambda$4385, FeI $\lambda$4920, FeII $\lambda$4924)
which fall in the windows where the equivalent width is measured.

\section{Test of the profiles}

In this section the synthetic profiles of some of the H and HeI lines are
compared to the  Auer \& Mihalas (1972) and  Kurucz (1993) profiles. They are
also compared with observations. The goal is to test the advantages and
limitations of our grid with respect to previous work.

\subsection{Comparison to Auer-Mihalas models}

Auer \& Mihalas (1972) computed non-blanketed NLTE profiles of H (P$\alpha$,
H$\alpha$, H$\beta$, H$\gamma$ and L$\alpha$) and HeI ($\lambda$4026, 4387,
4471 and 4922) for  $T_{eff}$ from 30000 to 50000 K, and log$ g$=3.3 to 4.5.
Figure 7 compares the profiles of H$\gamma$ for $T_{eff}$=30000 and 40000 and
log$ g$=4.0 with the spectra of the grid. The profiles are very similar, the
small discrepancy comes from the effect of rotation (we assume a  rotation of
100 km s$^{-1}$ and none in the Auer-Mihalas models) and the inclusion of
atomic transitions in  the spectral range we have synthesized. The contribution
of these lines makes the profile of H$\gamma$ asymmetric in our spectra but not
in the Auer-Mihalas profiles. This represents  an improvement of our synthetic
profiles over those of Auer \& Mihalas because our grid is more similar to the
observed spectra of stars (see for example Walborn \& Fitzpatrick 1990).

Figure 8 shows the profile of HeI $\lambda$4471 for $T_{eff}$=40000 K and
log$ g$=4.0.  The apparent disagreement comes from the effect of the rotational
convolution performed in the grid.  In fact, the difference between the
profiles disappears when the comparison is done between  non-convolved
profiles. The effect of the rotational convolution is more drastic in the
profile of  the HeI lines than in the Balmer lines because HeI lines are weaker
and narrower than Balmer lines.

\subsection{Comparison to Kurucz models.}

For $T_{eff}\leq$ 25000 K, the grid has been generated using the stellar
atmosphere structure calculated by Kurucz (1993). Kurucz also synthesizes the
profiles of H$\alpha$, H$\beta$, H$\gamma$  and H$\delta$. However, the
synthesis does not include the atomic transitions  that fall within  $\pm$100
\AA\ of the center of the Balmer lines. Figure 9 compares the profiles  of our
grid with those of Kurucz (1993). They are very similar, the apparent
disagreement at $T_{eff}\leq$ 7000 K is due to the effect of the rotational
convolution and the inclusion of metallic lines in the spectral ranges
synthetized here. However, note that the profile of the Balmer lines before
performing the rotation is equal to the profile synthesized by Kurucz.

Our synthetic spectra assume fewer sources of continuum and line opacity than
the Kurucz stellar atmosphere. This inconsistency leads to differences between
the continuum  calculated by TLUSTY+SYNSPEC and the low resolution spectra
generated by Kurucz. However,  the difference is less than 8$\%$ in  the whole
spectral range synthesized here if $T_{eff}\geq$ 7000 K. For lower
effective temperature, the deviation is more important since the shape of the
continuum also changes.  However, the normalized profiles of the Balmer lines
of our grid and Kurucz (1993)  are the same for all effective temperatures
($T_{eff}\leq$ 25000 K). Thus, we can confidently use the normalized
profiles of the grid in our evolutionary synthesis code, where these profiles
are  calibrated in absolute flux.

\subsection{Comparison to observations}

Figure 10 compares the synthetic profiles from 3700 to 3900 \AA\ with the
spectra of observed stars from the stellar library of Jacoby, Hunter, \&
Christian (1984). The observed spectra are normalized assuming a
pseudo-continuum that was defined by fitting a third order polynomial to the
windows defined in Table 3. The synthetic spectra are binned to the resolution
of the observations ($\sim$4 \AA). There is good agreement between the
synthetic and observed profiles. The discrepancy between the observations of
the O5V star and synthetic profiles at  $\lambda\leq$3760 \AA\ is probably due
to uncertainties in the wavelength and flux calibration of the data.

We have also compared the synthetic spectra with two stars observed by one of
us for a different project. The stars HD24760 and HD31295 are classified as
B0.5V and A0V, respectively. They were observed with the 2.5m Isaac Newton
Telescope at the Roque de los Muchachos Observatory using the 500 mm camera of
the Intermediate Dispersion Spectrograph and a TEK CCD detector. The dispersion
is 0.4 \AA/pix. Figures 11 and 12 compare the observed normalized spectra with
the synthetic profiles with effective temperature of 27500 K and 9500 K,
respectively. In both cases, the gravity is log$ g$=4.0. The agreement between
observations and synthetic spectra is very good, even if an optimization of the
fit was not attempted; we have only taken the spectra of our grid with values
of the $T_{eff}$ and gravity closer to the characteristic values of B0.5V and
A0V stars.

\subsection{Effect of the metallicity}

Opacities play an important role in determing the properties of the structure
of the atmosphere. Balmer lines, however, are not very much affected by the
abundance of elements heavier than H and He if the temperature is higher than
7000 K. The comparison of the Balmer profiles synthesized by Kurucz (1993) for
metallicity between Z$\odot$ and Z$\odot$/10 shows that the profiles are
essentially the same in this range of metallicity if $T_{eff}\geq$7000 K. Thus,
we can use our normalized spectra to predict the synthetic spectra of a stellar
population younger than 1 Gyr if the metallicity is higher than Z$\odot$/10.

\section{Summary}

We have computed a grid of stellar atmosphere models and synthetic spectra
covering five spectral ranges between 3700 and 5000 \AA\ that include the
profiles of the Balmer (H13, H12, H11, H10, H9, H8, H$\delta$, H$\gamma$ and
H$\beta$) and HeI ($\lambda$3819, $\lambda$4009, $\lambda$4026, $\lambda$4120,
$\lambda$4144, $\lambda$4387, $\lambda$4471, and $\lambda$4922) lines with a
sampling of 0.3 \AA. The grid spans a range of effective temperature  50000
K$\geq T_{eff} \geq$4000 K, and gravity 0.0$\leq$log$ g\leq$5.0 at solar
metallicity.

The spectra are generated using a set of computer programs developed by Hubeny
et al. (1995a,b). The profiles are generated in three stages. First, for
$T_{eff}\geq$25000 K, we use the TLUSTY code (Hubeny 1988) to compute NLTE
stellar atmosphere models. We assume 9 energy levels of HI, 14 levels of HeI
and 14 levels of HeII explicitly in NLTE. For $T_{eff}\leq$25000 K, we take the
Kurucz (1993) LTE atmospheres. In the second stage, the synthetic spectra are
produced with the program SYNSPEC (Hubeny et al. 1995b) that solves the
radiative transfer equation using as input the model atmosphere and a line list
that contains information about atomic transitions in the relevant spectral
wavelength ranges. Although the NLTE models generated with TLUSTY are pure H
and He stellar atmospheres, SYNSPEC assumes line blanketing for elements
heavier than H and He.

Our grid of synthetic spectra has limitations due to the inconsistencies
between the continuum and line opacities assumed in the stellar atmosphere and
those assumed in calculating the transfer of the lines, which start to be
important for $T_{eff}\leq$7000 K. However, they reproduce very accurate the
Balmer lines generated by Kurucz (1993) even at $\leq$7000 K. The profiles of
the HeI lines are also very similar to those generated by Auer \& Mihalas (1972)
for $T_{eff}$ between 30000 and 50000 K.

This work was motivated by the importance of the HeI and the high-order Balmer
lines in the study of the stellar population of galaxies with active star
formation. These galaxies show Balmer and He lines in emission produced in
their HII regions that are super-imposed on the stellar absorption lines.
However, the equivalent width of the nebular lines decreases quickly with
decreasing wavelength, whereas the equivalent width of the stellar absorption
lines is constant with wavelength (e.g. the equivalent width of H8 ranges
between 2 and 12 \AA\ and H$\beta$ between 2.5 and 14 \AA). Thus, the
high-order Balmer series members are less contaminated by emission, and they
can be very useful to study the underlying stellar population in starburst
galaxies. The evolutionary synthesis profiles of H and He absorption lines for
starbursts and post-starbursts are presented in paper II. The full set of the
resulting models is available at our websites http://www.iaa.es/ae/e2.html and
http://www.stsci.edu/science/starburst/, or on request from the authors at
rosa@iaa.es.

{\bf Acknowledgments}

We thank Ivan Hubeny for his help during the initial phase of this project and
for making his code available to the community, and to Enrique P\'erez, Tim
Heckman and Angeles D\'\i az for their helpful suggestions and comments. The
spectra of some observed stars come from data of an ongoing project to build an
observational library of stellar spectra in collaboration with Angeles D\'\i az,
Pepe V\'\i lchez and Enrique P\'erez. This work was supported by the Spanish
DGICYT grant PB93-0139.

\clearpage


\section{Introduction}

   Pion production in nucleon-nucleon collisions has been given
increased attention due to the recent availability of precise
near-threshold data for $pp\rightarrow pp\pi^{0}$,
$\vec{p}{p}\rightarrow d\pi^{+}$ and $\vec{p}p\rightarrow pn\pi^+$.
These reactions involve a large momentum transfer,
but near threshold only a few partial waves can contribute.
Here the study of the short-range features
of the NN interaction is greatly facilitated because it is  possible
to experimentally separate the contributions from different
angular momentum states. In recent studies it was
found that within 10-20 MeV of threshold,  i.e. below 300 MeV, only Ss
for  $pp\rightarrow pp\pi^{0}$, \cite{meyer90,meyer92,bonder}
or at most two partial waves (Ss and Sp) for
$\vec{p}{p}\rightarrow d\pi^{+}$ and $\vec{p}p\rightarrow pn\pi^+$
\cite{korkmaz,heimberg,drochner,daehnick95,hardie,daehnick98,flammang}
 were important, and direct comparisons with
simple theoretical approaches can be made. The expected dominance
of Ss contributions ($l_{NN}=0$, $L_{\pi}=0$) for the final states was
seen, but at the same time it was found that for
$pp\rightarrow pp\pi^0$ traditional models under-predicted the cross
sections by factors of 3-5.

At higher energies the reaction $\vec{p}p\rightarrow pn\pi^+$
tends to be dominated by the delta resonance.
\cite{waltham85,shypit89,wicklund87}.
The present study in the near threshold region will provide information
on polarization observables in a region still subject to
relatively parameter-free microscopic calculations.
Its major motivation is the hope and expectation to reduce the
theoretical uncertainties for a process fundamental in nuclear physics.
Additional references to prior experimental and theoretical work can be
found in
\cite{flammang,hardie}.

   Much effort has gone into attempts to explain
 ${pp\rightarrow pp\pi^0}$, probably the theoretically most difficult
 of the three $\rm NN \rightarrow NN\pi$  branches.
Modern chiral perturbation work on this reaction has been discussed recently
in
\cite{cohen,sato}, but this approach has not yet explained the data.
It has become apparent that for $pp \rightarrow pn\pi^+$
the first order terms dominate and delta and  complicating heavy
meson exchange contribution are relatively small making this reaction
theoretically important.  At this time there is no refereed publication on
chiral
perturbation calculation; however, very recently a preprint has become
available that reports $\rm \chi pT$ calculations for $pp \rightarrow
pn\pi^+$
and $pp \rightarrow d\pi^+$.  These $\rm \chi pT $
calculations under-predict near threshold data by factors of 2 to 3,
but with less room to adjust poorly known parameters than
for ${pp\rightarrow pp\pi^0}$ \cite{darocha}. The varying
theoretical approaches proposed for the inclusion of off-shell
contributions,
``heavy meson''  exchange, and contributions involving the $\Delta$ isobar
are still under discussion \cite{Hanhart98a}.  Divergent results
suggest that conclusions about the adequacy of a
 theoretical pion production model cannot be based on agreement for
just a single pion production branch.  A satisfactory model must succeed in

 explaining all $\rm NN\rightarrow NN\pi$ branches.

Promising attempts to simultaneously describe data from more than one
reaction
have been made using the Hamiltonian approach . A  relativistic
phenomenological
model was offered by Shyam and Mosel \cite{shyam}.
Their approach naturally included higher angular momenta and
the exchange of heavier mesons and was able to fit $pp\rightarrow pn\pi^+$
and $pp\rightarrow pp\pi^0$ cross sections over a wide energy
range.   This model should also predict analyzing powers;
however, such predictions have not been published to date.
More complete microscopic calculations in the meson exchange
picture, involving
angular momentum states as large as $L_{\pi}=2$ are now being published by
the J\"ulich group \cite{haidenbauer,hanhart97,hanhart98b}.
Their approach
of \cite{hanhart97,hanhart98b} is largely  free of adjustable
parameters.  It gives a very good reproduction of the published
$pp\rightarrow pp\pi^0$ and $\vec{p} p \rightarrow pn\pi^+$ data.
However, it  has only qualitative
success in reproducing the new spin correlations observables
\cite{meyer98} for
$\vec{p} \vec{p} \rightarrow pp\pi^0$.  The spin correlation data
reported here for $\vec{p}\vec{p}\rightarrow pn\pi^+$
present a different test of the predictive power
of this approach.

Precision measurements of pion production near threshold became
possible with the advent of cooled beams and the
technology developed at IUCF and other facilities around the
world to use a stored cooled beam on a windowless internal target.
In the present experiment we used a polarized beam on a polarized
internal target to measure the $\vec{p}\vec{p}\rightarrow pn\pi^+$
reaction at 325, 350, 375 and 400~MeV.  Bombarding energies for the present
experiment were chosen to cover the range over which some higher
partial waves gradually become significant. This is the first investigation
of
spin correlation coefficients in this reaction near threshold.
We have measured $\rm A_{y}$, $\rm S \equiv  A_{xx}+A_{yy}$,
$\rm D\equiv  A_{xx}-A_{yy}$,
integrated over all kinematic variables and the observed range of the
 of the polar angles. We also determined
 the differential analyzing powers $\rm A_y(\theta)$.

\section{Experiment}

The present experiment was conducted at the Cooler ring of the Indiana
University Cyclotron Facility (IUCF). A beam of polarized protons
was accelerated to 197~MeV and stack-injected
into the ring.  There it was accumulated for 2-3 minutes,  yielding
orbiting beam currents of 100-300~$\rm \mu$A. The stored beam was then
accelerated in the ring to full energy. After about 5-10 minutes of data
taking the remaining beam was discarded and the cycle was repeated with
reversed
beam spin. The beam polarization was vertical, alternately up or down.

The target consisted of  $74\pm 4\% $ polarized atomic hydrogen
in a storage cell, where it was maintained by an atomic beam \cite{wise}.
The open-ended cell,
a movable cylinder of 12~mm diameter with a 25~$\rm \mu$m thick
aluminum wall, minimized the background induced by the beam halo
in the cell wall. Depolarization
due to wall collisions was minimized by coating the cell with Teflon.
The target gas density distribution along the storage cell was
of triangular shape with a spread of $\pm$12.5 cm about the center.
The thickness of the target was $\rm  \approx 10^{13}~atoms/cm^2$.
The target polarization was flipped every 2 sec, pointing in
sequence up, down, left, right. Each change of direction took less than
10~ms.

Elastic proton-proton scattering events were observed concurrently with
pion production in order to determine the luminosity and beam
and target polarization.
Four plastic scintillators at azimuthal angles of
$\rm \pm45^\circ$ and $\rm \pm135^\circ$  detected
coincident protons elastically scattered near $\rm \theta_{lab}=45^\circ$.
The
product PQ of the beam (P) and target (Q) polarization was
determined from the spin correlation
 $\rm A_{xx}-A_{yy}$, which is large and well known for pp
scattering.  Typically, we saw  $\rm PQ\approx 0.5$, but for the 350 MeV
run PQ dropped to  0.34.  The time-integrated luminosity was deduced
from the known pp elastic cross section.
\input{epsf}

\vspace{-4.cm}
\hspace{-3.cm}
\begin{figure}
\epsfysize=10.cm
\centerline{
\epsfbox{fig1.ps}
}
\vspace{-0.0cm}
\caption{Apparatus for the experiment: WC1 and WC2 are wire
chambers. E and K are segmented  plastic scintillators stacks that
determine the energy of the  charged reaction products.
V is the veto detector.
F is a timing detector. X is one of the four  detectors for the elastic pp
scattering
monitor.
\label{fig:Fig1}}
\end{figure}

The detector stack is shown in Fig.~1.
The trajectories of the reaction pions and protons were determined by  wire
chambers WC1 and WC2, spaced 23 cm apart, each consisting of two wire
planes.
The angular accuracy of the deduced trajectories was
limited by multiple scattering of the ejectiles and by the
resolution of the wire chambers.
The major contribution to the multiple scattering came from the 1.5~mm
thick plastic scintillator detector (F) which served as timing detector.
We estimate the angular uncertainty for the
 trajectories for protons and pions
of about $\sigma=\rm 0.5^\circ$ and $\sigma=1.0^\circ $, respectively.
The energy of the charged particles was measured by a stack of segmented
plastic scintillators of total thickness of 25.4 cm.
The energy resolution of the stack was about $\rm \sigma=3.5\%$.
Events are vetoed if a charged particle reaches the detector V behind the
K-detector.

Various types of events were analyzed for the present experiment.
The $\rm \vec{p}\vec{p}\rightarrow pn\pi^{+}$ candidates
had to have simultaneous counts in one or more F segments and two
or more E segments. They were
rejected if there was a count in the veto detector V. The veto
detector effectively eliminated (or else selected) very energetic
(elastic) projectiles.  For
the 400 MeV run it also rejected the most energetic
reaction protons and pions, complicating the interpretation of this run.
In addition to the pion trigger there were a ``monitor trigger,'' to study
pp
elastic scattering,  and a ``single prong'' trigger for diagnostic purposes.

More details of the experimental set-up are given in \cite{pia}.

\section{Analysis}

Coincident protons and pions were identified by their deposited
energy and their time of flight
between the E and  F detectors. Software cuts were used to
further eliminate background. We employed a particle identification
cut, a cut based on the trajectory origin, and a missing-mass cut.
Fig.~2 shows the distribution of the missing mass reconstructed
from the measured momenta of the proton and pion.
\begin{figure}
\epsfysize=11.cm
\centerline{
\epsfbox{fig2.ps}
}
\vspace{-1.5cm}
\caption{Distributions of the calculated mass $\rm m_x$ of the undetected
particle at 325~MeV bombarding energy, shown for the four
combinations of vertical beam and target
polarization. A sharp peak ($\approx \rm 3.5~MeV/c^2$ FWHM) is seen at the
neutron rest mass. The shaded region indicates the
assumed background distribution. }
\end{figure}

The events of interest are those
near the neutron mass. The  software cuts and the coincidence
requirements greatly reduce  background, but some of it remains under
the  neutron mass peak. This amounts to 3 to 8\% of the peak area
and primarily comes from pions produced in the walls of
the target cell.  To assess the background shape we took data
with pure $\rm N_2$ gas in the target cell. We deduced
the background under the missing mass peak by using a background
shape derived from the $\rm N_2$ missing mass spectrum. This
empirical shape was normalized to the visible background continuum
outside the neutron peak. To estimate the error from  background
corrections we varied the background subtraction by $\rm \pm$25\% .
The effect on the final results was found to be smaller than
the statistical errors.
At 325 MeV the resolution of the neutron missing mass peak
was  $\rm \sigma=1.6~MeV$.

A Monte Carlo simulation of the experiment was used to determine various
limiting effects of the apparatus and to derive corresponding corrections.
The code contained the detailed geometry
of the detector systems and the target density distribution. It included the
effects
of multiple scattering of the charged particles, pion decay, loss of energy
by the charged particles, and the energy resolutions of the scintillators.
The simulation provided a guideline for the expected angular and energy
distributions of pions and protons. This was important at 400 MeV,
because about 0.2\% of the coincident protons and 14.4\% of the pions,
intercepted by the detectors, did not stop
in the scintillators. In the final analysis we used only those
events where the proton and the pion energies were measured accurately
enough to produce a  missing mass value within $\rm \pm 3.5~MeV$ of the
neutron mass. The simulation was also needed
to study the acceptance of the detector system, which  did
not cover the entire phase space of the exit channel.
Very forward events were suppressed by the central hole in the detector
which was needed to accommodate the beam pipe for the circulating
beam.  Many large-angle events are lost because of the finite detector size.
Coincident events are detected only if the polar angle of the
proton is larger than $\rm 6^\circ$
and if the pion angle is $\rm 6^\circ \leq\theta_{lab}\leq 32^\circ$.
The event acceptance for $\rm p\pi^{+}$ coincidences
was 21.5\% at 325 MeV,  19.4\%  at 350 MeV, 18.2\%  at 375 MeV, and 15.0\%
at 400 MeV.  The resulting
loss in count rate is significant;
however, since target and projectile are identical particles
we sample a  phase space twice as large as the acceptance
suggests.

In the c.m. system low energy pions are seen for all
angles, but energetic pions are not sampled for polar angles near
$\theta_{c.m.}=90^{\circ}$.
Near threshold the analyzing power for pions in the $\rm p n \pi^{+}$
final state  varies about linearly with pion momentum, independent of
bombarding energy \cite{daehnick98}.
Therefore $\rm A_{y}$ data for c.m. angles between
$\rm65^{\circ} \ and \ 115^{\circ}$ not
only have low statistics, but they also require (model dependent)
corrections, which results in large error bars.
The integrated spin dependent cross section ratios S, D, and $A_{y}$
reported here lack data near $\theta_{c.m.}=90^{\circ}$,
and therefore only approximately represent the full integral over $\rm
\theta_{c.m.}$.
Nevertheless, accurate comparisons with theory are possible.  They
require a  corresponding integration range for the theoretical calculation,
i.e., the
restriction of pion lab angles to the range  $\rm 6^{\circ}$ to
$32^{\circ}$.

\section{Results and Discussion}

Fig.~2 compares four missing mass spectra with the same integrated
luminosity for the four spin-spin orientations with vertical beam and
target polarization.
 The spin-dependent  cross section $\Delta\sigma'_T$
 (defined as the observed cross sections with
the spins of the colliding protons opposite minus the
observed cross section with spins parallel) can be calculated
by using the cross ratio method\cite{meyer97b}:
\begin{equation}
\rm \frac{\Delta\sigma'_T}{\sigma_{unpol}} = \frac{2}{PQ} \
\frac{1-\sqrt{R}}{1+\sqrt{R}}
\end{equation}
where $\rm R = (Y_{uu}Y_{dd}/Y_{ud}Y_{du})$.  $\rm Y_{mn}$
denotes the yields
for the four spin combinations  and  $\rm \sigma_{unpol}$ is the
unpolarized cross section for the measured  range of the polar angle
$\rm \theta$.
Because of the limited $\rm \theta$ range of the detector the integrated
cross
section ratio $\rm \Delta\sigma'_T/\sigma_{unpol}$ will differ from
the standard ratio $\rm \Delta\sigma_T/\sigma_{total}$ if the spin
correlation coefficients have a strong $\rm \theta$  dependence.
The recent Hamiltonian calculation by Hanhart et al. \cite{hanhart98b}
suggests relatively small changes for our geometry (see Fig.~4).

For the present experiment spin-dependent differential cross section
as a function of the pion azimuthal angle $\rm \phi_{\pi}$
can be expressed as\cite{meyer97b}
\begin{eqnarray}
\sigma(\rm \phi_{\pi}) = \sigma_{tot}[1+(p_yA_y^B - q_yA_y^T)\cos\phi_{\pi}
\nonumber \\
  +p_yq_y(\frac{A_{xx}+A_{yy}}{2}-\frac{A_{xx}-A_{yy}}{2}\cos2\phi_{\pi})
\nonumber\\
  + terms\ involving\  p_x,p_z,q_x,q_z],
\end{eqnarray}
where $\rm A_y^B$ and $\rm A_y^T$ are respectively the beam and target
analyzing powers, and $\rm p_x,p_z,q_x,q_z$ are very small.
The values $A_{nn}$ are the spin correlation coefficients and
$\rm p_i(q_i)$ are the beam (target) polarization components. In this
equation
the observables $\rm A_y^B$, $\rm A_y^T$ and $A_{nn}$ are integrals
over all kinematic variables of the proton, the neutron and the pion
except over the pion azimuthal angle $\phi_{\pi}$. Spin correlation
observables
are deduced by sorting the
data as a function  $\rm \phi_{\pi}$ for a specific
combination of beam and target spin.
We show the $\phi$-dependence of the measured ratios in Fig. 3.

Except for the scales, the presentation
is identical to that in ref.\cite{meyer98} for $\rm
\vec{p}\vec{p}\rightarrow  pp\pi^{0}$.
The columns represent the four bombarding energies,
and the rows the different combinations $\rm R_i$ for the spin-dependent
yields. The definitions of the $\rm R_i$ are given in Fig. 3.
To obtain the plotted quantities, we calculate
$\rm W_i = (R_i-1)/(R_i+1)$ and divide by $\rm \sqrt{PQ}$ for the first
three rows and by PQ for the last two two rows.

\onecolumn
\begin{figure}
\epsfysize=20.0cm
\centerline{
\epsfbox{fig3.ps}
}
\vspace{-4.0cm}
\caption{Asymmetries for different spin combinations $\rm R_i$ (listed on
the
left) as a function of the azimuthal angle of the $\rm \pi^+$. The solid
curves represent a least-square fit using the theoretical
$\rm \phi_{\pi}$ dependence (listed on the right), varying $\rm a_B$,
$\rm a_T$, $\rm S\equiv (A_{xx}+A_{yy})$, and $\rm D\equiv
(A_{xx}-A_{yy})$.}
\end{figure}
\twocolumn

The first three rows are directly related to the beam and target
analyzing powers, whereas the
last two rows reflect the spin correlation coefficient combinations
$\rm S\equiv (A_{xx}+A_{yy})$, and $\rm D\equiv (A_{xx}-A_{yy})$. In this
figure $\rm a_B$ and $\rm a_T$ are related to the beam and target analyzing
powers $\rm A_y^B$ and $\rm A_y^T$ by $\rm a_B=A_y^BP/\sqrt{PQ}$, and
$\rm a_T=A_y^TQ/\sqrt{PQ}$.
The solid lines in Fig.~3 represent a least-square fit using the
$\rm \phi_{\pi}$ dependence from Eq.(2) (listed to the right of Fig.~3),
and varying the four parameters $\rm a_B$, $\rm a_T$, S, and D.

\begin{figure}
\epsfxsize=5.5cm
\centerline{
\epsfbox{fig4.eps}
}
\vspace{+0.5cm}
\caption{The integrated $\rm \vec{p}\vec{p}\rightarrow pn\pi^+$
 observables
$\rm \Delta\sigma'_T/\sigma_{unpol} = -(A_{xx}+A_{yy}$),
$\rm D\equiv A_{xx}-A_{yy}$, and $\rm A_y$ as listed in Table~1 are shown
as a function of  $\rm \eta$, the maximum
pion center-of-mass momentum in units of the pion mass. The
observables $\rm A_{xx}$, $\rm A_{yy}$ and $\rm A_{y}$ were
integrated over all kinematic variables except the
pion azimuthal angle. Solid curves are from Hanhart et al..
These calculations are for the full phase space and are not directly
comparable
with the data. The dashed curves approximate the theoretical calculations
for the actual phase space observed.}
\end{figure}

The integrated spin correlation coefficient S is related to the
spin-dependent  cross section by  $\rm \Delta\sigma'_T/\sigma_{tot} = -S$,
and we find that these two analysis methods agree.
The integrated analyzing power is
given by $\rm A_y=(-a_B\cdot a_T)^{1/2}$. The resulting polarization
observables are shown in Fig.~4 and listed in Table~1. These values,
but not those in Fig.~3, are corrected for background contributions.
The values for the integrated $\vec{p}\vec{p}\rightarrow pn\pi^+$
spin correlation coefficients in Fig.~4  are in remarkably
good agreement with recent
Hamiltonian predictions following ref. \cite{hanhart98b}.\\

In our experiment beam and target analyzing powers
are related by  $\rm A_y^B(\theta)=-A_y^T(\pi-\theta) $ and,
therefore, permit a study of the back angles for  $\rm A{_y}(\theta_{c.m.})$
with the present apparatus.
Fig.~5 shows the observed $\theta_{c.m.}$ dependence of the analyzing
powers $A_{y}(\theta_{c.m.})$.  As mentioned above,
the limited size of the detector reduced the usefulness of data
 near $\theta_{c.m.}= 90^{\circ}$.  However, a large enough
range in $\theta$ could be covered to demonstrate the change in the
angular distributions with energy and the need for
partial waves higher than Ss and Sp.
The dashed curves show calculations taken from \cite{hanhart98b}.
\begin{figure}
\epsfxsize=5.5cm
\centerline{
\epsfbox{fig5.eps}
}
\vspace{0.1cm}
\caption{The measured $\rm \vec{p}\vec{p}\rightarrow pn\pi^+$
analyzing powers $\rm {A_{y}}(\theta_{c.m.}, T)$ compared with fits (solid
curves) described in the text. The angular distributions suggest some
asymmetry
around $\rm \theta_{c.m.}=90^{\circ}$. Acceptable fits require the
inclusion of an associated Legendre polynomial with $L_{\pi}=2$. The
detector size did not permit reliable measurements for
$\rm -0.25 < cos(\theta_{c.m.}) < 0.25$. The dashed curves are calculations
by Hanhart et al.[18].}

\end{figure}

The solid curves in Fig. 5 represent fits with the relation
\begin{equation}
A_{y}(x)=(1/\sigma_{unpol}(x))*(a~P^1_{1}(x)
+b~P^1_{2}(x)+c~P^1_{3}(x)),
\end{equation}
where the expressions ${P^1}_{\nu}(x)$ are
associated Legendre polynomials and $x=cos(\theta)$.
$\rm P^1_{3}(x)$ was not needed or used for
the fits shown; however,  $\rm P^{1}_{2}(x)$,
i.e. the inclusion of an $L_{\pi}=2$ term, was essential to obtain a
reasonable representation of the data at the four energies.

We have found that the recent calculations of the J\"ulich group
are quite successful in predicting the integrated spin correlations
for $\vec{p}\vec{p}\rightarrow pn\pi^+$ in spite of their difficulties
with $\vec{p}\vec{p}\rightarrow pp\pi^\circ$.  This dramatic difference
in the otherwise similar reactions should help identify the reason for the
less
than satisfactory agreement of calculations and measurements
\cite{meyer98} for $\vec{p}\vec{p}\rightarrow pp\pi^\circ$.
A further indication of remaining shortcomings is found in the
difference between calculated and measured angular distributions
for  $\rm A_{y}(\theta_{c.m.})$ in Fig. 5. It appears that some
diagrams that are secondary in $\vec{p}\vec{p}\rightarrow pn\pi^+$
but important in $\vec{p}\vec{p}\rightarrow pp\pi^\circ$ have been
omitted or treated inadequately.

We are currently analyzing additional measurements of the
$\rm \vec{p}\vec{p}\rightarrow pn\pi^+$ reaction with
longitudinal beam and target polarization. Such  results
together with the present data may allow a deduction of
the dominant partial waves in a model independent way.\\

\begin{acknowledgements}
    We acknowledge the support of Drs.~M.~Dzemidzic, F.~Sperisen
and D.~Tedeschi in the early stages of the experiment.
We thank Dr. W. Haeberli for his advice and continuing interest and J.
Doskow for technical support.  We also
wish to thank the IUCF accelerator operations group for their
dedicated efforts. We are grateful to the authors of
ref. \cite{hanhart98b} for  making available to us calculations
for $\rm \vec{p}\vec{p}\rightarrow pn\pi^+$
obtained with their model.
This work was supported by the US National Science Foundation under
Grants PHY95-14566, PHY96-02872, PHY-97-22556, and by the department
of Energy under Grant DOE-FG02-88ER40438.
\end{acknowledgements}

\pagebreak


\section{Introduction}

The question whether accretion discs can be in a state of self-organised
criticality (SOC) (Bak et al. \cite{Bak:Tang}, Kadanoff et al.
\cite{Kadanoff:Nagel}) was raised explicitly by
Mineshige et al. (\cite{Mineshige:Takeuchi}), and strongly implicitly
by Bak et al. (\cite{Bak:Tang}).  Young \& Scargle
(\cite{Young:Scargle}) have raised the related question of transient chaos in
accretion systems.  A distinctive feature of SOC is flickering energy transport
with no characteristic lengthscale or time separation, displaying a $1/f$
power
spectrum; in this context, $1/f$ is
shorthand for inverse power law frequency dependence with unspecified
index.  In
SOC systems, for which mathematical sandpiles (Bak et al. \cite{Bak:Tang},
Kadanoff et al. \cite{Kadanoff:Nagel}) provide a paradigm, global transport
occurs as a result of self-organised avalanches which are triggered locally
when
the accretion of sand leads to a critical gradient being exceeded at a given
point.  This causes local redistribution of sand, which may lead to the
critical gradient being
exceeded at neighbouring points, resulting in further redistribution and,
cumulatively, to a global avalanche.  Following the avalanche, the system
returns
to a subcritical configuration; accretion then continues until it again
creates a
local excess gradient, triggering a further avalanche.  In the SOC-sandpile
paradigm, it is helpful to note that the word ``critical" is being used in two
senses: the sandpile has a critical gradient, redistribution being triggered
wherever this gradient is exceeded; and the self-organised global avalanches
which emerge from the integrated effects of local redistribution may be
scale-free, in which case they are linked to the generic field of critical
phenomena.

In mathematical models
and also some experimental realisations (Nagel \cite{Nagel:1992}, Feder
\cite{Feder:1995}, Frette et al. \cite{Frette:Christensen},
Christensen et al. \cite{Christensen:Corral}), avalanche statistics display
scale-free $1/f$ characteristics.  SOC sandpile algorithms are extremely
simple and possess the
attraction of any successful reduced system: it becomes unnecessary to
attempt to
model the detailed, and perhaps insuperably complex, microphysics of transport
in the real system.  Such an approach has been applied to the statistics of
solar X-ray
bursts (Lu \& Hamilton \cite{Lu:Hamilton}; for a recent treatment, see
MacKinnon \& Macpherson \cite{MacKinnon:Macpherson} and references therein
) and to transport phenomena in magnetically confined plasmas (Newman et
al. \cite{Newman:Carreras},
Carreras et al. \cite{Carreras:Newman},
Dendy \& Helander \cite{Dendy:Helander}, \cite{Dendy:Hel2}). It is clearly of interest to
establish whether,
in certain circumstances, the detailed modelling of astrophysical accretion
flows
could be substituted by a simple SOC sandpile paradigm.  This is a highly
cross-disciplinary question.  In the present paper, we aim to carry forward the
debate in two ways: by identifying in greater detail the classes of accretion
flow where SOC might play a role; and by examining theoretical arguments
that we
believe point clearly towards SOC in some accretion flows.

\section{Some sites of interest}

It was noted by Bak et al.  (\cite{Bak:Tang}) that flickering signals with $1/f$
power spectra have been observed for the X-ray variability of active galactic
nuclei (AGNs), specifically: 0.05-2keV X-rays from the Seyfert galaxy NGC4051
(Lawrence et al \cite{Lawrence:Watson}) and 2-7keV X-rays from the Seyfert
galaxy NGC5506 (McHardy \& Czerny \cite{McHardy:Czerny}), and for 10-140 keV
X-rays from the massive compact binary Cyg~X-1 (Nolan et al.
\cite{Nolan:Gruber}).  Similar $1/f$ spectra for X-ray variability from binary
accreting systems were noted by Mineshige et al.  (\cite{Mineshige:Takeuchi})
(neutron stars, Makashima \cite{Makashima:1988}) and by Geertsema \& Achterberg
(\cite{Geertsema:Achterberg}) (cataclysmic variables and dwarf novae, Wade \&
Ward \cite{Wade:Ward}).  While we believe this list can be extended, as we
discuss below, it is also important to be more specific about where, in these
diverse systems, SOC might be occurring.

Perhaps the simplest case is presented by a paradigmatic AGN. Let us take for
this the standard picture of a massive black hole fed by a cascade of
structures: the accretion disc, the molecular torus and on a larger scale the
galactic disc and its barred or spiral structure, which together generate a
clumpy and irregular transfer of gas and stars.  This irregular mass transfer
could be analogous to the random and discrete feeding of a sandpile with grains
of sand near its apex.  It is believed that, at least in thin discs, the
residence time $\tau_{res}$ of matter in the disc is much longer than
the free--fall time, $\tau_{ff}\approx (R^3/GM)^{1/2}$ at radius $R$:
\[\frac {\tau_{res}} {\tau_{ff}}\sim \frac {R^2} {\alpha H^2}\]
where $H\ll R$ is the thickness of the disc, and $\alpha\ll 1$ is the
well-known viscosity parameter. The condition $\tau_{res}\gg\tau_{ff}$,
which is implicit in most accretion disc models, appears to be
necessary, at least in principle, for the validity of a sandpile--type
approach to mass transfer within the disc.

Thus if the accretion disc, like the sandpile, is in a state
of SOC, there would be no clear link between the pattern of accretion from the
torus and the pattern of avalanches leading to mass transfer across the disc,
over its inner edge, and onto the black hole, giving rise to the X-ray signal.
The latter would automatically display $1/f$ flicker.  On larger space and
timescales relating to our own Galaxy, even steady gas inflow from the Galactic
bar could also result in avalanches to the inner regions, and these avalanches
could give rise to episodes of AGN activity during an otherwise quiescent phase,
with a duty cycle of a few percent, as inferred on statistical grounds for
other AGNs (Mezger et al., \cite{Mezger:Duschl}).  The sandpile-SOC paradigm thus provides a
candidate framework for summing up the consequences of the complex and intricate
physics that relates the largescale dynamics of the Galactic disc to the
activity of its central black hole.

Cataclysmic variables, and in particular dwarf novae, present further
opportunities for SOC. First, the flow of matter from the secondary
across the inner Lagrangian point could itself be an avalanching process
governed by SOC. In this case, if the radiation source was a hot spot
where the accreting mass flow reached the outer edge of an accretion
disc, it would flicker with $1/f$ statistics.  Steady mass flow from the
secondary would also be compatible with an SOC signal, however, provided
the latter originated from SOC mass transfers from the accretion disc
(or accretion column in the strongly magnetised regime) to the white
dwarf.  Either regime appears possible in principle, both for dwarf
novae and in the wider context of binary accreting systems.  In dwarf
novae, it is pointed out by van Amerongen et al.
(\cite{Amerongen:Kuulkers}) and Lasota et al.  (\cite{Lasota:Hameury})
that outbursts are due to suddenly increased accretion towards the white
dwarf, but that it is unclear whether the instability resides in mass
transfer from the secondary to the accretion disc, or across the
accretion disc and onto the white dwarf.  We note that, observationally,
flickering is more often strong during the quiescent phase of dwarf
novae than during major eruptions (Wade \& Ward \cite{Wade:Ward}).
Questions relating to the location of unstable flows in soft X-ray
transients (Lasota et al.  \cite{Lasota:Narayan}) are similar to those
in dwarf novae.  For the wind-fed X-ray binary pulsar GX301-2, a model
has been proposed (Orlandini \& Morfill \cite{Orlandini:Morfill})
involving ``noisy" accretion of blobs of matter formed by
magnetohydrodynamical (MHD) instability at the magnetospheric radius,
and not caused by inhomogeneities present in the stellar wind from the
optical companion.  This approach is somewhat reminiscent of that of
Baan (\cite{Baan:1977}), where accreted matter accumulates at the
magnetopause of a rotating neutron star until an interchange instability
is triggered, after which the released matter generates an X-ray burst.
Baan (\cite{Baan:1977}) suggested that, since most of the time between
bursts is a refilling time, an approximately linear relation should
exist between burst energy and the subsequent quiescent interval.
However, in an SOC model where randomly arising local instability is
sufficient to trigger a global avalanche, there would be no such
correlation.  This appears to be a key observational discriminant for
the possible presence of SOC in a given accreting system.  We also note
that, from the theoretical point of view, Frank {\it et al.}
\cite{Frank:King} have pointed out that a local instability in a given
annulus of the disc can only trigger large--scale instability accross
the disc if parameters in neighbouring annuli are such that the effect
of this instability in these annuli can in turn trigger local
instability there. This amounts to a prescription for a sandpile--type
approach, and hence for the possibility of SOC.

It seems clear from the foregoing that there is good observational and
interpretative motivation for testing for SOC in a broad range of accreting
astrophysical systems, encompassing flickering AGNs and certain distinct
locations within a variety of binary objects.  Firm identification of SOC would
yield information on the global consequences of the smaller-scale physics of
the accretion process, while short-circuiting the need for detailed modelling.
Let us now turn to theoretical arguments for expecting SOC.

\section{Theoretical considerations}

The identification of mechanisms responsible for the dissipation of shear flow
energy within an accretion disc is an active field of research, because of its
importance in determining transport across the disc and onto the compact
object.
If, as has been suggested (Bak et al \cite{Bak:Tang}, Mineshige et al.
\cite{Mineshige:Takeuchi}), the combined global effects of local transport
physics result in SOC, observational signatures of the type described above
will
emerge.  It is widely accepted (see, for example, Longair (\cite{Longair:1994})
and Narayan (\cite{Narayan:1997})) that anomalous viscosity caused by MHD
turbulence probably plays an important role in the flux of angular momentum
within accretion discs. MHD turbulence can arise naturally in accretion
discs, see for example the
instability mechanisms proposed by Tagger et al. (\cite{Tagger:Henriksen}),
Vishniac et al. (\cite{Vishniac:Jin}), and Balbus \& Hawley
(\cite{Balbus:Hawley}).  It has also been pointed out (Chen et al.
\cite{Chen:Abramowicz}) that there exist ranges of accretion rate $\dot{M}$ and
disc radius $R$ for which no stable steady state solution of the basic
equilibrium equations is possible.  In this case, it is suggested (Chen et al.
\cite{Chen:Abramowicz}, Narayan \cite{Narayan:1997}) that the flow is
forced into
a time-dependent variable mode, satisfying the required $\dot{M}$ in the mean.
This provides further motivation to test the applicability of sandpile-type
models.  In this section, we concentrate on the model for MHD turbulent energy
dissipation in accretion discs presented by Geertsema \& Achterberg
(\cite{Geertsema:Achterberg}), because it appears to
give rise to SOC.  While the link to cataclysmic variable and dwarf nova
observations (but not, explicitly, SOC) was made by  Geertsema \& Achterberg
(\cite{Geertsema:Achterberg}) themselves, we believe that this work has much
wider implications: for the general question of SOC in accretion flows; for the
connection between SOC and turbulence in general;
for the role of SOC in plasma physics, since it represents the first
instance where SOC
has been observed in a mathematical model derivable from the fundamental
equations of MHD; and for the role of SOC in terrestial experimental systems -
real sandpiles and ricepiles, as distinct from mathematical idealisations
thereof
- where uncertainties about its scope remain (Nagel \cite{Nagel:1992}, Feder
\cite{Feder:1995}, Frette et al. \cite{Frette:Christensen}, Christensen et al.
\cite{Christensen:Corral}).

Before considering the model of Geertsema \& Achterberg
(\cite{Geertsema:Achterberg}) in greater detail, let us turn to its
results.  Figure 12 of Geertsema \& Achterberg
(\cite{Geertsema:Achterberg}) shows the calculated times series of
energy dissipation events within the disc.  We note that this is
qualitatively very similar to the observed time series of energy
dissipation measured in an experimental ricepile displaying SOC, Fig 2c
of Frette et al.  (\cite{Frette:Christensen}), and in a related
mathematical model of Dendy \& Helander (\cite{Dendy:Helander}, fig.  3,
and \cite{Dendy:Hel2}, fig.5).  More quantitatively, Fig 13 of Geertsema
\& Achterberg (\cite{Geertsema:Achterberg}) shows the power spectrum of
energy dissipated by MHD disc turbulence, which displays the $1/f$
dependence characteristic of SOC; compare, for example, Fig 3 of Frette
et al.  (\cite{Frette:Christensen}) and Fig 3 of Christensen et al.
(\cite{Christensen:Corral}) which show measured spectra of energy
dissipation and particle transit times, respectively, in SOC ricepiles.

Given the clear indications of SOC emerging from the MHD turbulence model of
Geertsema \& Achterberg (\cite{Geertsema:Achterberg}), it will be of interest
to  establish  how  it  has arisen.  A full explanation must await diagnostic
analysis of the code runs generated by this model.  Pending this, we
conclude the
present section by seeking to identify some of the relevant salient features.
In outline, the model of Geertsema \& Achterberg (\cite{Geertsema:Achterberg})
is constructed  as  follows.

The accretion disc is regarded as a differentially rotating turbulent MHD
fluid,
and is modelled by a reduced system of equations reflecting the most important
features of three-dimensional MHD turbulence. The disc is assumed to be thin in
comparison with its diameter, and the flow is taken to be subsonic and hence
incompressible. The possible existence of a large-scale magnetic
field is neglected, and the only turbulent structures considered have a length
scale shorter than the height of the disc, making the turbulence essentially
three-dimensional. In a coordinate system rotating with the disc angular
velocity
${\bf \Omega}= {\bfm e}_z \Omega_K(R)$ at radius $r=R$, the simplified MHD
equations for the flow velocity ${\bf u}$ and the normalised magnetic field
vector ${\bf b} = {\bf B} / (4 \pi \rho)^{1/2}$ are
\begin{eqnarray}
\dot{\bf u} = -{\bf \nabla}\Pi - 2 {\bf \Omega} &\times& {\bf u} -
\left( r \frac{d\Omega_K}{dt} \right)_R u_x {\bf e}_y \nonumber \\
&+& \nu \nabla^2 	{\bf u} + {\bf b} \cdot \nabla {\bf b}
- {\bf u} \cdot \nabla {\bf u}, \nonumber
\end{eqnarray}
\[
\dot{\bf b} = \left( r \frac{d\Omega_K}{dt} \right)_R b_x {\bf e}_y
	- \eta \nabla^2 {\bf b} + {\bf b} \cdot \nabla {\bf u}
	- {\bf u} \cdot \nabla {\bf b},
\]
where $x=r-R$ and $y=R\phi$ are the local radial and azimuthal
coordinates respectively, $\eta$ is the resistivity, and $\nu$ is the
viscosity coefficient. $\Pi = P/\rho + b^2/2$ is the total scalar
pressure, while the ${\bf b} \cdot \nabla {\bf b}$ term contains the
off--diagonal terms of the MHD stress tensor.

As these equations are still too complex for a numerical treatment over a
sufficiently large dynamic range in scales, they were simplified further
following a suggestion by Desnyanski \& Novikov (\cite{Desnyanski:Novikov}) in
the context of hydrodynamics, and later applied to MHD by Gloaguen et al.
(\cite{Gloaguen:Leorat}). In this approximation, the
space of wave vectors $\bf k$ is discretized into a finite set of $k_n$,
and the non-linear
interaction between different components of the Fourier transforms of
${\bf u}$ and ${\bf b}$ is described by the set of equations
\begin{eqnarray}
\frac{du_n}{dt} + \nu k_n^2 u_n &=&
	\alpha k_n \left( u_{n-1}^2 - 2 u_n u_{n+1} \right) \nonumber \\
&+&	\beta k_n  \left( u_{n-1} u_n - 2 u_{n+1}^2 \right) - ( u
\leftrightarrow b ),
\nonumber
\end{eqnarray}
\begin{eqnarray}
\frac{db_n}{dt} + \eta k_n^2 b_n &=&
\alpha k_n \left( u_{n-1} b_{n-1} - 2 u_n b_{n+1} \right) \nonumber \\
&+&
\beta k_n  \left( u_{n-1} b_n - 2 u_{n+1} b_{n+1} \right) -
( b \leftrightarrow u ). \nonumber
\end{eqnarray}
In the three-dimensional generalisation of this approximation scheme,
the discretisation is made
by dividing the ${\bf k}$-space into into spherical shells, thus
discarding the information  regarding the direction of
${\bf k}$. This allows a greatly simplified system of nonlinear equations
to be written down, analogous to that of Gloaguen et al.
(\cite{Gloaguen:Leorat}), which are supplemented with additional terms to
account for the effects of differential rotation.

These equations, which are taken to model the turbulent cascade of MHD,
were solved numerically by Geertsema \& Achterberg
(\cite{Geertsema:Achterberg}).  They found that the turbulent shear
stress can be very large, and has large, chaotic fluctuations on time
scales of a few rotation periods.  Perhaps the most striking feature of
the simulations is, however, that the dissipation of energy at the
smallest scales of the turbulent cascade is very intermittent.  Energy
is released in avalanches with a wide range of sizes, and the power
spectrum of the dissipation rate obeys a $1/f$ power law over nearly two
orders of magnitude in intensity, see their Fig 13.  This behaviour is,
as already stated, similar to that of simple mathematical sandpile
models.  While a power law is to be expected from any scale-free model,
in particular a Kolmogorov-type one in the inertial range, we note that
the similarities between this system and mathematical sandpiles
apparently extend further.  Both are fundamentally governed by
nearest-neighbour interactions between a discrete number of nodes.  In
the MHD accretion model, these reside in ${\bf k}$-space, so the
interaction is between wave modes with similar wavelengths rather than
between adjacent regions in real space.  The dissipation is provided by
viscosity and resistivity at large $k$, resembling the removal of
material from the edge of a sandpile.  The MHD accretion model is, of
course, much more complex than the simple sandpile algorithms considered
so far in the literature.  Two fields, $\bf u$ and $\bf b$, each with
three components, are involved rather than the single height parameter
of conventional sandpile algorithms, and the time evolution is governed
by differential equations rather than difference equations.
Nevertheless, both the MHD accretion disc and simple sandpile models
appear to exhibit similar self-organised, critical behaviour, supporting
the claim often made in the sandpile literature that SOC is universal
phenomenon shared by large classes of cellular automata.

\section{Discussion}

We have shown in this paper that the numerical simulations by Geertsema \&
Achterberg (\cite{Geertsema:Achterberg}) of viscous resistive MHD turbulence in
an accretion disc give rise to behaviour characteristic of self-organised
criticality in a sandpile.  This similarity may help to explain certain
observed properties of a range of accreting astrophysical systems, which we
have reviewed.  Furthermore, the result is of intrinsic scientific interest as
an ab initio demonstration of the emergence of SOC from a system of MHD-based
equations.  Further analysis of the numerical results, to establish how the
SOC-sandpile phenomenology arises, would be of great interest both for
accretion disc astrophysics and for fundamental plasma physics.  At the present
stage of development of realistic MHD simulations (for example of the link
between smallscale dissipation and the time behaviour of the flux of matter
onto the central object), direct numerical proof of the occurence of
SOC-sandpile phenomenology is still a remote objective.  However, the capacity
of the SOC-sandpile paradigm to circumvent complex ``full mathematical
descriptions" is the basis for its present attraction in many fields of physics.
If future work confirms the importance of SOC-sandpile phenomenology in
accreting systems that is suggested in the present paper, it would permit a
dissociation (at least at zeroth order) of the detailed physics of turbulence
in the disc from the global modelling of this class of astrophysical object.

\begin{acknowledgements}
This work was supported in part by the Commission of the European Communities
under Contracts ERBCHRXCT940604 and FMRX-CT98-0168. The authors are
grateful to Bram Achterberg for helpful comments as a referee.
 \end{acknowledgements}