#### Cellular Automaton experiments on local galactic structure. II. Numerical simulations

Astron. Astrophys. Suppl. Ser.
Cellular Automaton experiments on local galactic structure. II. Numerical simulations
A. Lejeune 1
J. Perdang 0
0 Institute of Astronomy , Madingley Road, Cambridge CB3 OHA , UK , and Institut d'Astrophysique , 5, Avenue de Cointe, B
1 Institut de Physique , Sart
2 4000 Liege , Belgium
| This paper is a step towards demonstrating that the multi{parameter Cellular Automaton framework designed for the simulation of local galactic structure, developed in a companion paper (Perdang & Lejeune 1996, Paper I), is capable of duplicating the local irregularities observed in the structure of flocculent spiral galaxies. The numerical simulations exhibit the development of fractured galactic arms and the formation of fractal geometries associated with the matter distribution (fractal structure of the arms, of bulk dimension 1.7 and border dimension 1.3; distribution of di erent stellar components on fractal supports, of dimension 1.6, for reasonable estimates of the free model parameters). The prediction of fractional values for the di erent dimensions specifying the simulated structures can be exploited as a qualitative test of adequacy of the proposed model. The precise quantitative values of the observed dimensions, in conjunction with the observable global mass fractions of the di erent galactic components, play the parts of constraints for the free model parameters. We show that the currently theoretically inaccessible values of the free parameters of the formulation can be recovered from observation.
galaxies; structure | star formation | methods; numerical | galaxies; evolution
1. Introduction
In this Paper we report on an extended series of
numerical experiments carried out in the framework of a
2{dimensional Cellular Automaton (CA) model for local
galactic structure described in the companion paper
(Paper I). The model takes account of stellar evolutionary
e ects not dealt with in current CA models, as well as of
quasi-stationary motions more complex than usual di
erential rotation. Our results are indicative that the
floccular character of galaxies can be simulated in the CA
framework.
For technical details of the implementation of the
galactic physics we refer to Paper I. Our results are given
in a relative reference frame F in uniform rotation with
respect to the galaxy to be simulated.
2. Initial conditions. di usion
Remarks on infall and
In all of our experiments we start out with an initial di use
non{uniform gas covering the part of the galaxy we are
interested in (region Z over which the CA computations
Send o print requests to: A. Lejeune
?Permanent address
are carried out). Any cell is either in the empty state E,
or in the gas state D; with the notations of Paper I
F (0; i; j) = 0 or 1; for any cell (i; j):
(2:1)
The choice between an E or a D state is made on
the basis of a method reminiscent of a `mole's labyrinth'
(Herrmann 1983)
, already adopted in previous CA
simulations
(Lejeune & Perdang 1991)
: The cells (in the
auxiliary lattice L, cf. Paper I) are all provisionally populated
with gas; a number of b = 1, 2, ... empty `bubbles' are then
generated by T {step random walks starting at b randomly
chosen initial cells; for a selected number of di erent walks
b, the number of steps T per walk is chosen such as to
generate a preassigned fraction f (E) of empty space in the
lattice area; the remaining fraction of the area which is
lled with gas is denoted f (D) = 1 − f (E). Through the
coordinate transformation Eq. (5.1a) of Paper I we obtain
the initial state in the `physical' zone Z (the lattice L0).
Figure 1 provides a realisation of an initial con guration
(corresponding to the velocity eld of Fig. 3 of Paper I).1
1The purpose of this operation is to create an irregular cloudy
shape, a locally non{uniform initial matter distribution
de
A special simulation should be included in our
program to handle the internal dynamics of molecular clouds
reponsible for the formation of the complexity of the
actual cloud shapes. The random walk procedure as adopted
here should be viewed as a cheap substitute for simulating
the observed fragmented gas cloud border. Incidentally,
the numerical estimate of the fractal dimension ( 1.3)
of the border of our simulated cloud in Fig. 1 compares
favourably with the reported estimates of the fractal
dimensions of cloud contours.
We note also that with our algorithm the gas area is
either connected or it is disconnected, depending on whether
or not the `empty bubble' extends from the left to the
right border of our zone Z; provided that infall and
diffusion are weak enough, connectivity is preserved at any
later time. Under the rst alternative a signal in one cell
of the gas zone can propagate over the whole zone, i.e.
the state of any cell (i; j) at a given timestep can
influence the state of any other cell (i0; j0) at some later time
0; under the second alternative disconnected areas evolve
independently.
Infall, as handled in the present simulation (Reaction
(6.12) of Paper I) may generate a large number of D cells
within the empty bubbles during the initial stages (Tinf
100 timesteps) over which we let it operate. If Nc (
ny2, ny number of cells in the y direction) denotes the total
number of cells of our lattice ( 106), a rough estimate
of the number of these new gas cells N (D) created by the
implementation of infall is
N (D)
Nc f (E) P (E; D) Tinf :
(2:2)
In a T {step random walk the distance r(T ) a particle
traverses on average, as estimated by the expectation E of
r(T )2 is
E(r(T )2)
pdi T
(2:3)
in units of the cellsize (pdi , di usion probability
factor as implemented in our program). Hence, even for the
strongest possible di usion, (pdi = 1), over a timescale
T ( 103, for a typical evolutionary run), the average
distance (2.3), in units of the total size of the CA network
is
[E(r(T )2)]1=2=ny
(T =Nc)1=2;
(2:3a)
3 10−2 1 for our experiments. Di usion,
therefore, although not completely negligible, cannot play a
dominating part in our model. In the hydrodynamic
representations (Fig. 4a) the e ect of di usion is obliterated
for an averaging radius obeying p . This is indeed
consistent with our experiments.3
Table 1 summarises the evolutionary reactions
simulating the stellar evolutionary processes analysed in
Paper I. The table also lists the corresponding independent
transition probabilities per CA timestep (denoted here
P1; P2; :::; P15) as adopted in our numerical experiments.
Column (a) gives these parameter values for the galactic
evolution from which we have constructed synthetic
photographs. Column (b) corresponds to the reference
`average' model around which our calibration models for the
coe cients in representations (3.3) and (4.4) are
generated. Column (c) refers to one representative of the latter
collection of models; it provides an indication of the range
of model parameters explored in our experiments.
With the numerical values f (E) 0:25, P (E; D) = 1
10−5, we have N (D)=Nc 2:5 10−4 1, indicating that
this mechanism should not play a signi cant evolutionary
role.2
The action of di usion in the framework of our
simulation is more complex. On the computational level of the 3The di usion algorithm adopted in our program models
normal di usion, conforming to relation (2.3). For the purposes
individual cellstates F ( ; i; j) (exhibited in Fig. 1a for an of simulating motions of real individual galactic components
initial state, and in Fig. 3 for a later evolutionary phase), (proper motions of stars, for instance) normal di usion
apthe di usive mechanism detaches border gas cells, thereby pears to constrain too strongly the allowed velocities to
consticreating new disconnected patches of matter. However, on tute an acceptable model for the actual individual motions. We
a macroscopic, observationally meaningful `hydrodynamic believe that a model of anomalous di usion (or overdi usion)
scale', simulated in our experiments by spatial smoothing could capture more adequately the physical situation.
Anoma(Fig. 1b; compare also Fig. 4), the averaged patches re- lous di usion could be implemented in principle (probabilistic
main connected to the main clouds. displacement of a cell by a distance r obeying a Levy
distribution; for the latter the asymptotic probability density for large
void of any particular symmetry
(cf. also the procedure of displacements is of the form 1=r +1, 0 < < 2; cf. Feller 1971)
.
Heisenberg & von Weizs¨acker 1948). The random walk method We have not attempted to incorporate this mechanism in our
was chosen in our experiments since it is particularly easy to galactic model. We point out that for anomalous di usion the
implement. A more evolved CA approach capable of generating variance of the displacement does not exist, as a consequence
a variety of growth patterns, including realistic cloud shapes, of the contribution of the large displacements per timestep.
has been developed recently by one of us
(Perdang 1996)
. For this very reason anomalous di usion secures a more e
2In order to lead to the formation of a single connected gas cient transport: the exponent 1=2 in the RHS of (2.3a) is to
zone, and hence, in order to play a dominant part in the overall be replaced by 1= . It follows that provided that is taken
evolution, the characteristic time over which the infall should small enough, anomalous di usion can play a non{negligible
be operative, Tinf , should obey Tinf f(E)P1(E;D) : part over timescales as dealt with in our experiments.
3. The global evolution of the galactic species
In Figs. 2 we exhibit several evolutionary sequences of
the global mass fraction of our 6 galactic species. For the
parameter ranges adopted, 3 phases can be distinguished:
(i) An initial short transient phase, characterised by a
rapid building up of the protostar population and a
corresponding decline of the amount of gas. This phase
subdivides into a trigger{phase of star formation, of duration
determined by the probability of the cloud collapse (P1),
followed in turn by a more violent star formation process of
timescale xed by the catalytic reaction probability (P3).
In Frame (a), corresponding to a large probability P3, we
observe a prominent pulse of star formation of maximum
at step 20; in Frame (b) the probability P3 is reduced by
a factor 2 and the pulse of star formation is virtually
suppressed; in Frame (c), where P3 is further reduced by a
factor 2, the transient phase consists in a gentle rise of
the protostar population over about 150 steps.
(ii) A subsequent quasi{stationary phase in the
protostar, main sequence, and red giant population, during
which the formation rates of these species are essentially
balanced by their decay rates. The leakage of active
galactic matter towards inert matter (and low mass stars, which
play the part of inert matter as well within our
numerical approximations) ultimately implies that the `source' of
protostars is progressively depleted; on average the
protostar population is therefore bound to decrease. Frames
(a) and (b) indicate that quasi{stationarity is practically
realised beyond timestep 100. In all cases the slow
decrease due to loss of active matter can be approximated
by a linear time behaviour.
(iii) Finally, our model conceals an asymptotic
equilibrium state, for ! 1, in which all galactic matter is
in the form of inert or low mass stars (if the transition
probability P (L; I) is set equal to zero). Practically this
asymptotic state would be reached for evolutionary times
as such that as P (D; L) > 1, or as P (A; B) > 1, where
P (A; B) is the smallest transition probability in the global
transformation sequence D ! I via massive stars; in our
numerical experiments this stage is never attained.4
We mention at this stage that depending on the
precise parameter values, the quasi{stationary phase (ii) may
exhibit di erent time patterns. Besides the uniform
decrease, an oscillatory time-behaviour has been observed
if the intermediate steps of the global reaction D ! I
are arti cially set approximately equal. The low
amplitude fluctuations observable in the D and P components
of Fig. 2a may be a scar of such a periodic behaviour.
Moreover, when instead of examining the space averaged
global behaviour of the di erent species we examine the
detailed spatial behaviour of the species, spatial
oscilla4A semi{analytic discussion of the general time behaviour
based on an examination of the structure of the transition
matrix is given in Perdang (1993 Sect. 5) for a more schematic
version of the present model.
tions and intermittency phenomena are observed during
the globally oscillatory phase. We have not followed up
these types of complex behaviour in this paper.
In a mean eld approximation, in which the space
dependence is discarded, the global evolution equations
reduce to simple iterative schemes (di erence equations in
time). Approximating the di erences by di erentials, we
obtain ODEs for the evolution of the global population of
the species, thereby reducing our original CA formulation
to an elementary model of class (1) (Paper I). The
resulting ODEs can be interpreted as evolution equations of
chemical kinetics in a homogeneous medium; the kinetic
equations are read o from the stoichiometric equations
of Table 1. As already pointed out, the existence of
oscillating solutions in these kinetic equations, for some range
of the rate parameters, is well known: Standard
chemical kinetics can give rise to stable oscillatory behaviour,
provided that the reaction schemes involve a nonlinear
feedback, in the form of catalytic processes.5
In the context of the CA model of GSS oscillatory be
haviour is also encountered
(cf. in particular Seiden &
Schulman 1990 Fig. 16)
; the GSS oscillations appear in
some range of a refractory time, or period of immunity of
the gas against star formation; the latter plays the part of
a formal substitute for the details of the `chemical kinetics'
we are adopting in our model.
In principle, the relative populations of the species
(mass fractions), N (Y ; ; G)=PY N (Y ; ; G) = a(Y ;
; G), Y = D; P; M; R; I; L, over some large enough xed
area G of a galaxy, can be determined observationally at
the present age of a galaxy gal. Assuming that the galaxy
is in a stable quasi{stationary phase (ii), we nd from our
numerical experiments that the relative populations are
accurately represented in the form
5A triangular linear reaction scheme A ! B; B ! C; C ! A
shows an oscillatory behaviour if the 3 rates are identical (or
su ciently close); however, such an oscillation rapidly dies out
and the system settles in a stable stationary state. The
occurrence of a sustained oscillation around a stationary state
requires this state to be unstable. Instability demands that the
stationary state be on a `non{thermodynamic branch', which
is realised, in the chemical context, by catalytic and auto{
catalytic steps
(cf. Glansdor & Prigogine 1971)
. The various
catalytic steps in Table 1 were incorporated in our
simulation precisely to investigate the possibility of `chemical
instabilities' within the reaction network of the galactic species.
In our scheme, the quasi{stationary phase (ii) plays the part
of a non{thermodynamic branch parametrised by the `secular
time'. Asymptotically this branch merges with the
thermodynamic equilibrium state, namely state (iii). Small deviations
from the thermodynamic state inherit the stability of the
latter, so that the oscillations can occur su ciently far away from
the asymptotic state (iii) only; this thermodynamic property
implies that the oscillatory behaviour can survive over a nite
`secular' time only.
R
( a)
where the subscript o refers to the onset of phase (ii),
or more generally to any time su ciently exceeding the
lifetime of phase (i) of the galaxy. In our numerical
experiments we have set o = 100; except for a few experiments
involving small induced probabilitites this value of o
corresponds to a reference time within the stationary phase.
The parameter dependence of the two expansion
coefcients ao(Y ; G) and a0(Y ; G) is determined from a series
of test experiments. We represent this dependence by a
linear ansatz written in the form
ao(Y ; G) = aor(Y ; G) +
a0(Y ; G) =
a0r(Y ; G) +
X(Pk −Pkr) k(Y ; G); (3:2a)
k
X(Pk − Pkr) 0k(Y ; G); (3:2b)
k
where Pk is the transition probability of the kth
independent reaction step (notations of Table 1). A subscript r to
a model quantity aor, a0r indicates the value taken by this
quantity in the reference model de ned by the transition
probabilities (Pk = Pkr) given in Col. (b) of Table 1 (initial
amount of dust 72%). The coe cients of the expansions
(3.2), aor(Y ; G), a0r(Y ; G), and k(Y ; G), 0k(Y ; G), listed
in Table 2, are calibrated by a least squares t to a grid of
24 experiments in which each parameter ranges from 1/3
to 3 times its reference value.
To test for the goodness of the t we compared the
actual values ao and a0 of the 24 calibration models with
their linear estimates (Eqs. 3.2a,b). As a rule, provided
that ao is not too small (> 0:02), the relative agreement
is better than 20%. The linear ansatz (3.2a) is thus
satisfactory. Inaccuracies occur in runs in which the quasi{
stationary phase is not yet properly established at our
reference timestep o = 100 (which occurs if P1 or P2 are
too small). We observe also that for species L the
parameter ao is almost always well approximated (precision
better than 10%). On the other hand, the rate
parameters a0 are accurately given by the linear estimate only for
the inactive species I and L. For the remaining species
the estimate may di er by as much as a factor of 2 or 3;
within the accuracy of the linear ansatz (Eqs. 3.2a,b) it
is therefore legitimate to disregard the Pk{dependence in
the rates of the relative populations of the active species
A = P; M; R: a0(A; G) a0r(A; G).
The coe cient k(A; G) indicates how a variation
around the reference model of the kth reaction
probability, Pk, a ects the quasi{stationary population of the
active species A (= P; M; R); the induced variation of
the relative population obeys a(A; G) = k(A; G) Pk;
hence, if k(A; G) is positive [negative], then species A
is created [destroyed] | either directly or indirectly |
in process k. We note that the magnitude of the
sensitivity itself, j k(A; G) j, is not only large in a direct
creation or destruction process; it may remain so also in
indirect processes. As an example, the direct sensitivity ther remains a red giant, or it becomes an inert star, or it
of the protostar population to the probability of the di- transforms into gas; 2 transition probabilities are
indepenrect star induced protostar formation process, 3(P ; G), is dent. Since our model involves 21 transitions, or reactions
0.27; the R population, although only indirectly a ected (cf. Table 1) among 6 species, 15 transition probabilities
by the same reaction, has a larger sensitivity, 3(R; G) are independent (P1; P2; :::; P15 in Table 1).
= 0.71. Note that the sensitivities to spontaneous (gas Observationally we can determine 5 relative
populaand dust induced) star formation, 1(P ; G) and 1(R; G), tions of species, a(Y ; gal; G)obs, Y = P; M; R; I; L.
Subare much larger, 2.7 and 3.4 respectively; however, for stitution of Eqs. (3.2a,b) into Eq. (3.1) then yields the
all our runs the spontaneous transition probability P1 it- following linear system in the model parameters Pk
self is small ( 0.0015) in comparison with the star
induced transition probability P3 ( 0.033); accordingly the
overall contribution of the spontaneous formation process X[ k(Y ; G) + gal 0k(Y ; G)] Pk = a(Y ; gal; G)obs −
(P1 −P1r) 1(P ; G) remains always negligible with respect k
to the star induced formation (P3 − P3r) 3(P ; G).
As already mentioned, an examination of the orders of [aor(Y ; G)+ gala0r(Y ; G)]+X[ k(Y ; G)+ gal 0k(Y ; G)]Pkr
magnitude of the fragments (Pk − Pkr) k(Y ; G) entering k
the de nition of the quasi{stationary relative population (3:3)
ao(Y ; G) shows that some of them are negligible when the (5 equations in 15 unknowns). Further observational
inlatter coe cient is recomputed for the calibration models. formation is available from the spatial distribution of the
It is then reasonable to discard those fragments also in the di erent species. Our purpose is to show that in principle
general representation (3.2a). The P population, ao(P ; G), the latter can be quanti ed in a way as to supply
addiis dominated by the contributions of star induced star for- tional constraining equations in the unknowns Pk.
mation, of direct and star induced decay processes of
protostars into gas; star induced transformation of protostars 4. The space structure
into the M species also contribute; spontaneous P
formation is negligible, and so are most of the indirect processes.
For the M and R populations, ao(M ; G) and ao(R; G), all
steps yield essentially comparable orders of contribution.
Negligible contributions (Pk − Pkr) k(Y ; G) are listed in
Table 2 in the form k(Y ; G) = v 0, where v is the
actual numerical value of the sensitivity as obtained in our
calibration.
At any timestep , the evolution of each one of the 6
formal species E; D; P; M; R; I; L is represented by a
probabilistic transition (Table 1). For instance, for species R
our model accounts for 3 transitions (regarded as mutually
exclusive) translating the physical fact that a red giant
eiThe simulation of a mere global evolution of the
galactic species, summarised in a small number of time series,
would not require a detailed CA model: The global
evolution is satisfactorily reproduced by models of class (1),
namely by the low order ODEs generated from a mean
eld approximation; the only extra information supplied
by the CA model is a superposition of statistical
fluctuations which mimic to some extent the genuine observable
fluctuations.
The real interest of a CA approach is that it produces
the detailed time behaviour of a complex space dependence
of the physics of the galaxy (encoded in CA eld variable
F ( ; i; j)). Due to the intrinsically highly involved spatial
patterns, there seems to be no fully satisfactory way of
exhibiting the information carried by this eld variable in
any data{compressed form. A proper appreciation of the
full model results requires an analysis of the sequence of
frames de ned by the eld matrices F (0; i; j), F (1; i; j),
F (2; i; j), ... by an animation technique; only this form of
visualisation can properly reveal the local details of the
dynamics.
In Fig. 3 we supply a rst frame (a) showing the
microscopic CA state F ( ; i; j) of the galactic region at
timestep = 200 in the auxiliary lattice L together with
a space{averaged macroscopically meaningful picture of
the same zone (Gaussian averaging over 2 pixels), and a
second frame (b), at step = 500 in the physical lattice
L0 together with a space{averaged picture (radius 2 and
4 pixels) (model parameters: column a of Table 1; initial
conditions shown in Fig. 1). The intricacy of the averaged
computed structure at any timestep leaves little hope
Reference state
o = 100; G: whole lattice
reaction
for a duplication of these results by a low order di
erential model, or even by a tractable PDE model.
Figure 4a is an attempt at generating a synthetic black Figure 4b is a real HST photograph of a part of M100 of
and white photograph from the computed CA state at roughly same size as the computed galactic zone of Fig. 3b.
timestep 500; the grey level at each CA cell is de ned Although both pictures are not comparable in detail | we
through a weighted average of the contribution of the lu- should keep in mind that the synthetic gure represents a
minous matter (cf. caption); the nal picture is a spatial single realisation of a stochastic model | they do show one
Gaussian average over the `microscopic' grey levels (radius small{scale characteristic structural similarity: individual
= 2 pixels). The pattern of white patches is seen to match lsihgahpteps.atches in both pictures have statistically equivalent
almost exactly the lightblue areas in the averaged picture
of Fig. 3b.
reaction 7: P + mM + rR !
In our model these patches correspond to regions of This conclusion cannot be extended to arbitrary
galachigh concentrations of protostars; the existence of con- tic dynamics. In the presence of a chaotic flow eld, we
nected patches is related to the SPSF mechanism (and should expect the latter to a ect, and increase the
fracthe
Goldreich & Lynden{Bell (1965)
chaotic mechanism). tal dimension of the arm structure. On the other hand,
Figure 4c gives the structure of the levels of constant even in the case of a smooth velocity eld, and hence a
brightness of the synthetic photograph. smooth transformation between the physical space and
From Figs. 3 and 4 we can infer several general conclu- the auxiliary space, there is no guarantee that the
acsions on the space structure as generated in our models. tual technique of estimating the dimension does yield the
If we accept the view that the most luminous objects (the same value in both cases. In fact, the ln l against ln r
population of protostars) are the tracers of the galactic curve shows in practice a linear part over some r range,
arms (collection of bright patches in the synthetic photo- and it is the slope of the latter which estimates F ; at
graph Fig. 4a), then our model generates arms exhibiting small scales, r < ro (= actual resolution of the
gan irregular fragmented local geometry (on scales of the ure: 5 to 10 cells), the ln l { ln r curve flattens out; at
order of less than 200 pc); the large scale geometry (> large scales, r > rL (order of macroscale of structure
400 pc), on the other hand, appears to be channelled by in the gure < size of gure), the local slope of ln l {
the smooth dynamics of the galaxy (gravitation induced ln r can take on arbitrary values. Under a nonlinear but
shearing e ect of the stationary velocity eld). The sta- smooth transformation of the gure, if the transformation
bilised arm pattern tends to establish itself in the quasi{ itself exhibits structure over the range ( ro; rL), then
stationary regime (ii) of the global evolution of the species. the slope of the curve ln l0 ln r0 of the transformed gure
will be a ected over the transformed range of ( ro; rL).
Although theoretically the fractal dimension as de ned
for r ! 0, remains the same in the original and
trans4.1. Small scale geometry formed variables, practical estimates are bound to
operate at nite resolutions > ro. For instance, if we have a
As a measure of the local fragmentation we have deter- set of hierarchically distributed points on a straight line
mined fractal dimensions associated with the arm struc- (so that 0 < F < 1), experimentally approximated by
ture. This can be done in several ways: The fractal di- a nite number of dots distributed according to a nite
mension F of the boundary of the area occupied by the number of hierarchical levels, then a numerical estimate
luminous objects (Fig. 3) in (a) the physical lattice L0; of the fractal dimension F on the experimental realisation
and in (a') the auxiliary lattice L; one could likewise di- in ( ro; rL) yields an acceptable theoretical estimate of
rectly measure the fractal dimension of synthetic isophotes F . However, suppose we perform a stretching of the
ex(Fig. 4c). The natural method for estimating the fractal perimental points, x0 = Sx (x, initial coordinate of a point
dimension F of a topologically linear set (border of the on the line; x0 new coordinate;S, stretching factor); if S
area occupied by a population A, or xed isophote re- is chosen such that the intrinsic resolution in the new
cospectively) embedded in a plane, consists in measuring ordinates, ro0 = S ro rL, and if we estimate again
the length l( r) of this set at di erent resolutions r; the dimension F' from the slope of the ln l0 { ln r0 curve,
the dimension F 2 [1; 2] is then related to the slope of the at small r0 (< S ro), then we manifestly nd F 0 < F ;
ln l( r) { ln r relation, which is equal to 1 − F at a ne note that F 0 ! 0 for S large enough. This example should
enough resolution (provided that this relation is linear for make it clear that an a priori knowledge of the actual
inr ! 0). It is practically more convenient to estimate F trinsic resolution of our experimental data is required to
via the correlation dimension
(which is equal to the frac- make a reliable numerical estimate of a fractal dimension.
tal dimension if the latter exists; cf. for instance Perdang
1990, 1991)
; this latter method has been adopted in our The discrepancies in fractal dimension estimates of
atnumerical work. Theoretically, since L0 is a smooth de- tractors in phase spaces referred to di erent variables
(cf.
formation of L, the estimate of the fractal dimension is Cannizzo et al. 1990 for stellar pulsation attractors)
are
independent of the grid L or L0. related to this e ect.
The fractal dimension F shown in Table 3 is a
dimension of the area occupied by the arms rather than of the
border of the arm region (or length of an isophote). The
fact that this dimension turns out to be non{integer
reflects here the fragmented and hierarchical topology of the
arms; the number of arm fragments increases as we
increase the resolution; this behaviour is in agreement with
the observational situation (cf. the role of resolution in
Fig. 1 of Paper I). The evolution of the dimension F
exhibited corresponds to the run chosen for constructing the
synthetic photographs (parameter values of Col. (a) in Ta- gas distribution (empty bubbles in our simulation); the
ble 1). The e ect of the comparatively low spontaneous resulting smooth geometry is then irregularly modulated
star formation probability P1 is to delay the onset of the by the randomly acting physics of the star formation and
catalytic star formation process. As a consequence, the evolution processes.
establishment of the quasi{stationary state is delayed as In our experiments we start with large size empty
well ( o 300 instead of 100 for the larger P1 values). It fragmented bubbles floating in an originally uniform gas
transpires from Table 3 that the bulk fractal dimension is con guration. Since, as indicated, the additional dynamic
reasonably stationary (F 1.8) beyond timestep 300. mechanism included in our treatment, namely di usion, is
The dimension given in Table 3 is the correlation di- ine cient over distances of the order of a fraction of the
mension of the area occupied by the arm
(cf. Perdang size of our lattice, any large scale deformation of a bubble
1991)
; Fig. 5 also illustrates the behaviour of the frac- is a consequence of a stretching by the inhomogeneous
vetal dimension of the border (estimated again through a locity eld only. Through the stretching the bubble
transcorrelation dimension) which, as expected, is lower than forms into a single empty lane tending to become parallel
the previous dimension. Other fractal dimensions attached to the shear direction; it transforms into a double lane
with the arms (mass dimension of arm area, mass dimen- once the stretched bubble has acquired a length
exceedsion of border) have also been estimated; the latter may ing the length of the ring it belongs to, etc. The action of
di er by 10% from the corresponding correlation dimen- the shear on several bubbles leads to a system of parallel
sions; beyond step 300 all dimensions are essentially sta- empty lanes. Since no activity can occur in the empty
retionary (cf. Fig. 5). gions, the bubbles de ne the interarm zone. Support of the
It seems natural to expect that the ne structure of observable arm, the zone lled with gas likewise acquires
an arm (understood as the locus of the luminous objects) a shape of parallel arcs covering our lattice area.
is directly related to the kinetics of stellar evolution, and This large scale scenario, in substance already
demore speci cally to the star{induced star formation (tran- scribed by Weizsa¨cker (1951) as a mechanism responsible
sition probability P3). Accordingly the fractal dimension for the formation of spiral structure in galaxies, is
indeof the arm structure should help constrain the range of pendent of any reaction kinetics.6
the model parameters. If we summarise again the results We convinced ourselves of the e ciency of this
mechaof a collection of numerical experiments in an interpolation nism in the CA framework by running our simulation
proformula of type (3.2a) (disregarding here the time depen- gram with the evolutionary kinetics and di usion turned
dence), we obtain a constraining relation in the form (cf. o : Voids and matter are then observed to arrange in arms.
Eq. (3.3) with the dashed parameters left out) Test runs further make it clear that the global shape of
the arms is determined by the analytic form of the
stationXk kPk = F obs − Fr + Xk kPkr (4:1) tahrye vfrealoctcaitlydi melden.sWioenmoefntthieonmiantptears{sivnogidthinattetrhfaecoereitsichaelrlye
time independent ( 1.30 in our experiments); it is just
(F obs, observational determination of the arm dimen- xed by the construction algorithm of the bubble;
physision for a real galaxy). However, a few test experiments cally this dimension should then reflect the formation
prohave suggested that the dimensions of the arms are only cesses of the molecular clouds. In the actual simulation,
marginally sensitive to the individual stellar processes. It however, the stretching transforms the initially irregular
would seem that once the stationary arm pattern is estab- lines of constant density into almost straight lines along
lished, the arm structure acquires fairly universal fractal the shear direction. This is a consequence of a loss of
spadimensions which are independent of the precise physics tial resolution accompanying the stretching
transformaprevailing in the formation process of the arm. Such a pic- tion which, in principle, should reveal ner and ner
deture is not inconsistent with the percolation phase transi- tails (cf. our previous remark). Hence, since the CA model
tion describing the arm formation as proposed in GSS (cf. is working with a nite resolution, a numerical estimate
further comments in the next subsection). yields an apparent fractal dimension decreasing with time
In any event, our experiments are indicative that towards 1.
the knowledge of these dimensions does not provide (ii) Carriers of the collection of luminous objects, the
any signi cant constraints on the transition probabilities lanes of matter then naturally trace out an underlying
parametrising stellar evolution. large scale structure of progressively smoother geometry
on which the observable brightness patterns of the galaxy
4.2. Large scale structure
We can further infer from our simulations that the large
scale pattern traced out by the arms has a twofold origin:
(i) It is due to the action of the galactic dynamics (the
shear of the velocity eld) on an initial inhomogeneous
6In a paper with Heisenberg, the validity of this idea was
substantiated by a simulation in which an elongated irregular 2D
cloud of dots is rotated about an inner point of the cloud,
according to Kepler's third law. The cloud then develops into a
regular two{armed spiral
(Heisenberg & von Weizs¨acker 1948)
.
are superimposed (Fig. 4a). As already recalled, GSS have
interpreted the onset of the formation of a bright arm as
a percolation phase transition leading to the occurrence of
a `connected in nite cluster' of active stars. Their model
makes use of a single control parameter, p (probability of
induced star formation, i.e. counterpart of our parameter
P3); the phase transition occurs for p > pc, c a critical
value.
The analogue of the GSS percolation transition is
observed in our more elaborate model, the cluster being
identi ed by the protostar population. We mention that
at xed reaction probabilities the initial fraction of gas,
f (D), can likewise be chosen as a relevant control
parameter for the phase transition: An `in nite' cluster develops
as f (D) exceeds a critical value fcrit(D)
( 0.66 for the
choice of reaction parameters in Lejeune 1993)
.
The value of the critical control parameter (p, P3 or
f (D)) eventually depends on an arbitrary convention of
what we understand by `connection' of 2 neighbouring
cells. To illustrate this point, de ne two cells as being c{
connected if they are less than c cells apart; a `cluster' is
then a collection of c{connected cells. It is manifest that
the onset of the transition and the critical value of the
control parameter depend on the speci c arbitrary choice
of c.
This remark is indicative that a practically relevant
concept of `percolation' for the galactic context must be
predicated on the observational resolution (of which the
parameter c is a caricature). The formal percolation
concept as introduced in GSS has no directly accessible
quantitative observational counterpart. The resolution
dependence of the `connected in nite cluster' can be
demonstrated by image analysis of synthetic photographs (Figs.
4a,b) and real photographs (cf. Fig.1 in Paper I): As we
increase the smoothing e ect and generate lower
resolution images, we observe more extended connected bright
areas.
The GSS interpretation has the virtue, in our opinion,
of supplying a qualitative explanation of the formation
process of the bright arm structure in the galaxy by relating
it to the statistical mechanics of percolation. The
luminous arm structure (i.e. the extended connected patches
covered by protostars and luminous objects in general)
appear or disappear, as some xed combination of the
parameters regulating the evolution of galactic components
goes through a critical value. On the account of the
arbitrariness involved in the de nition of the `percolation
cluster', it is hard to see, however, how the GSS
percolation concept could help extract quantitative information
on the stellar processes from galactic observations. In fact,
if the transition is a genuine percolation phase transition,
then the relevant information it does supply is that it
carries no quantitative information on the speci c physical
mechanisms of the underlying model: The numerical
coefcients characterising a real percolation phase transition
(critical exponents and γ; cf. Eqs. (1.3) of Paper I), are
universal, i.e. determined by global topological and
geometrical factors (dimensionality and lattice symmetry);
they are independent of the precise physics responsible for
the formation of the cluster (exact evolutionary steps and
corresponding transition probabilities). In view of this
remark global parameters quantifying the arm structure are
not used as constraints for the physics of our model. The
fractal dimensions F already mentioned seemingly belong
into this category as well.
Since our model involves several species of luminous
objects (P , M and R), a phase transition is possible in
principle for each of these species. We have traced such a
situation under the conditions of an oscillatory regime.
Extended areas of the galaxy then undergo cyclical changes:
P ! M ! R ! P ! :::. If realisable in genuine galaxies
| the regime requires ranges of the transition
probabilities probably not consistent with realistic galactic
conditions | the physical nature of the luminous arm would
cyclically change, from a predominantly protostar
composition, to a main sequence and red giant composition.
Whether the transition isolated by GSS is to be viewed
as a genuine percolation phase transition remains to be
substantiated, even though the power{law behaviour of
the cluster size characteristic for percolation has been
veri ed under various conditions
(cf. Lejeune 1993)
. In fact,
a percolation phase transition, in its standard acceptance,
describes a transition between static states under a quasi{
static change of a control parameter (transition from
absence of a percolation cluster to presence of a cluster); the
notion of percolation involves no dynamics
(cf. Stau er
1985)
. In the present physical problem, we rather observe a
propagation, i.e. a dynamic phenomenon, with the in nite
cluster developing in time. Accordingly, the transition
appears to be more closely related to what has been referred
to as self{organised criticality
(Bak et al. 1987-1989)
, a
transition towards a dynamical state in which very large
numbers of degrees of freedom are excited. Near the onset
of self{organised criticality power{law dependences hold
for various magnitudes attached with the clusters of
excited degrees of freedom (distribution of cluster sizes,
distribution of lifetimes of clusters, etc). To the best of our
knowledge, and in contrast with the power{laws related
with percolation, there has been no proof so far that the
exponents of these laws are universal. Accordingly, if the
identi cation of the in nite cluster formation with a
transition to self{organised criticality applies in our case, then
the exponents could carry additional information on the
physics.
4.3. Space distribution of the luminous objects
An inspection of the picture of the cellstates (Fig. 3)
indicates that the overall space distribution of the
luminous objects, P , M and R is not homogeneous over
the area investigated. The luminous objects appear to be
concentrated over lower{dimensional subsets, and this
type of inhomogeneity is inherited by the inert objects I.
As a consequence of their formation mechanism (cf. Table
1), the low mass stars L, on the other hand, are regularly
scattered over the area originally covered by gas; but since
the latter itself is not uniform, the L distribution inherits
the lack of uniformity of the gas.
Physically the spatial distribution of a given species is
determined by a combination of three e ects: (a) the
initial con guration; (b) the mechanisms of formation and
destruction of this species; (c) the galactic kinematics
(overall velocity eld, and di usion). On the assumption
that the initial conditions are essentially the same for the
class of galaxies under consideration, and under the
proviso that we concentrate on small enough zones of these
galaxies (over which the kinematic stretching e ects are
ine cient), we are entitled to expect that information on
the transition probabilities Pk can be extracted from
conveniently chosen parameters describing the spatial
distribution of the species.
We have attempted to quantify the following single
aspect of the space pattern traced out by the species P ,
M , R and I. It is a well known fact that in CA
models individual states are often uniformly distributed over
self{similar, or statistically self{similar sets
(cf. Perdang
1993 Sect. 4 for an analytic discussion of several
examples)
. Therefore, we have addressed the somewhat broader
question of whether a mass dimension of the set covered
by species Z (= P; M; R; I) does exist; if so, this fractal
dimension can serve the purpose of a constraint for the
transition probabilities Pk. To this end we test whether a
scaling relation of the form
N (Z; i; j; G(r))
=
C rf(Z);
for r small; (4:2)
is applicable to our experiments; here N (Z; i; j; G(r)) is
the number of cells in state Z in a disc G(r) centred at an
arbitrary cell (i; j) of the auxiliary lattice L, or of the
physical lattice L0; the mass dimension f (Z) ( 2 [0; 2]) exists
provided that the interpolation coe cients C and f (Z)
of the experimental data are independent of the reference
cell (i; j). In the standard numerical method of estimating
f (Z) in our CA experiments di erences in lattice L and
L0 may arise as a consequence of the stretching and
accompanying loss of space resolution in L0 (cf. our previous
remarks). As already argued, under the loss of accuracy
due to a large enough stretching in one direction a
numerical estimate of f (Z) in L0 will be arti cially reduced as
compared to the estimate in L.7
7Test experiments have con rmed this point. For the model
corresponding to our synthetic photographs we have in L0:
f (P ) =1.60; f (M ) = 0.87; f (R) = 0.36 and f (I) = 0.10 at step
200 (compare with Table 4); the large di erences in the case of
the M and R species are due to an additional inaccuray related
to the small population of these objects. We have computed
Table 4 exhibits the conclusion of our test. It shows
that numerical values of f do indeed exist (plots of ln N
against ln r indicate a broad linear interval, and the slope
is independent of the central cell). The f values supplied
in the Table are all computed in the auxiliary lattice L in
order to avoid the loss of local accuracy accompanying the
stretching in L0. We believe that this procedure provides
the closest theoretical approximation to an observational
mass dimension as given in the frame of the observer,
provided only that the observational estimate is made over
a small part of a high resolution photograph (with G(r)
larger than the actual resolution).
The mass dimension for the di erent species does show
a signi cant dependence on the model parameters. In the
quasi{stationary phase (ii) this dimension is practically
time independent. Again, the apparently exceptional
behaviour of the model on which our synthetic photographs
are based, which exhibits a noticeable increase in f
between timesteps 200 and 500, is due to the fact that the
onset of quasi{stationarity is delayed. We approximate the
mass dimensions over phase (ii) by the linear ansatz
f (Z; ) '
fo(Z)
=
for(Z) +
(4:3)
(Z = P; M; R; I). The reference values for(Z) correspond
to the model whose parameters are listed in Col. b of Table
1; the collection of expansion coe cients k(Z) supplied
in Table 5 have been obtained from 15 calibration models
whose parameters Pk scatter around the reference model.
The counterpart of Eq. (3.3) becomes here
X(Pk − Pkr) k(Z);
k
X
k
k(Z) Pk
=
f (Z)obs
− for(Z) + X
k(Z)Pkr;
k
(4:4)
also the correlation dimension
(cf. Perdang 1991)
of the
distribution of the di erent species; we nd that as a rule, the latter
is smaller than the mass dimension, obeying approximately the
relation f (Z)−0.1. Since in the correlation dimension we sum
over all pairs of identical states, this measure includes here the
e ect of the input inhomogeneity in the matter distribution
(presence of empty bubbles); the mass dimension is estimated
at positions where matter is actually present. In our context
both measures refer to di erent aspects of the distributions.
(Z = P; M; R; I). Hence, given the observational estimates
of the current fractal dimensions of the distribution of the
di erent species P; M; R; I in a given galaxy (assumed to
be in the stable phase ii), f (Z)obs, relations (4.4) provide
4 additional linear constraining equations for the model
parameters Pk.8
5. Estimate of the free parameters
Our analysis has isolated 9 observational constraints in
the form of algebraic equations (Eqs. 3.3 and 4.4) to be
satis ed by our 15 a priori free model parameters Pk.
If 6 additional constraints were available, the proposed
CA formulation would have the status of a strictly
phenomenological or empirical model, analogous to a Fourier
representation for the simulation, the interpolation and
the extrapolation of a time series. On the theoretical side,
since a subset of the model parameters are recoverable
from stellar evolution theory, our CA model could serve
the purpose of testing current stellar evolution theory.
With 6 constraints being missing, we regard the
latter subset, namely the transition probabilities of isolated
objects, (P [P ; M ] = P4; P [M ; R] = P8; P [R; I] = P13;
P [L; I] = P15), as being given. The timescales of these
processes flow from evolution theory and are converted
into transition probabilities provided only that the stellar
mass function be known.
If we consider timescales less than a few Gy, then the
decay probabilities of the low mass stars become
negligible, and we may set further P [L; D] = P14 0 (together
with P [L; I] = P15 0). Moreover, as already pointed
out, the spontaneous star formation rate, P [D; P ] = P1,
has only a minor influence on the abundance of the
protostar species (except for the delay e ect mentioned). In
all cases the star{induced star formation rate is the
dominant mechanism. We believe therefore that we can
reasonably x the speci c parameter P1. The numerical value
8We have pointed out above that the critical exponents
referring to the percolation phase transitions are universal and
therefore they do not carry any information on the physics.
Since critical exponents are fractal dimensions in disguise, the
question arises whether conversely the fractal dimensions
attached with the geometric sets of the di erent species are not
intrinsically universal themselves. In fact, as observed
elsewhere
(Perdang 1993)
, in analytically solvable CA models over
square lattices only speci c dimensionalities turn up for the
carriers of individual states. It is not impossible then, that the
parameter and time dependences we observe in the
numerically determined fractal dimensions of the galactic simulation
merely reflect the fact that a true fractal structure is not yet
achieved. However, even if it were not a genuine fractal
dimension, f (Z) then would remain a quantitative index
characterising the space distribution of species Z. A valid model for a
given galaxy must duplicate the observational value of this
index. We believe, therefore, that even under those circumstances
the observational estimates of the exponents f (Z) continue to
carry information about the free parameters of our model.
D
+ dD
!
P
+ dD
D + dD ! L + dD
D + pP + mM + rR ! P + pP + mM + rR
P
!
R
L
not
not
!
parametrised
!
parametrised
I
I
we adopt is the value of the reference model (Col. (b) of
Table 1). Within this framework of approximations we are
then left with the following list of 9 parameters
(P2; P3; P5; P6; P7; P9; P10; P11; P12):
(5:1)
The constraints we have discussed in the previous sections
enable us then to derive these parameters directly from
observation, namely (i) from relative mass fractions of the
galactic components, and (ii) from fractal dimensions of
the space distributions of the components.
for(P )
for(R)
1(P )
1(M)
1(R)
1(I)
3(P )
3(M)
3(R)
3(I)
5(P )
5(M)
5(R)
5(I)
7(P )
7(M)
7(R)
7(I)
9(P )
9(M)
9(R)
9(I)
11(P )
11(M)
11(R)
11(I)
6. Conclusion and outlook
We believe that our experiments provide arguments that
the galactic arm structure results from a combination
of two physical ingredients which have already been
isolated separately in models of class (2)
(cf. the Heisenberg{
Weizs¨acker simulation, 1948)
, and in the simplest models
of class (3) (3{state CA experiments of MA and GSS).
The two ingredients, rather than being mechanisms
associated with what seemed to be regarded as rival theories,
appear here not only as complementary, but as
simultaneously necessary for the development of a realistic structure
of a luminous arm pattern:
On the one hand, the overall large scale pattern of
the matter distribution | the Grand Design | is
determined by the gravitation induced general flow eld (the
di erential rotation in a lowest order approximation). This
structure is derivable from classical N {body experiments.
On the other hand, the brightness distribution, as
directly accessible by photographs, reflecting the
distribution of the bright components, is the immediate outcome
of star formation and evolution processes. Since the latter
processes occur in regions of high matter density, the high
density regions are the carriers of the brightness
distribution. The visible arm pattern, while then globally shaped
by the matter distribution, possesses a structure of its own
conditioned by the kinetics of star formation and
evolution, which is superposed on the overall density
distribution. The small scale structure of the brightness
distribution is chiefly determined by this `chemical kinetics'.
We believe that this is the main conclusion of our
analysis. It is borne out by a comparison of our synthetic
photographs (Fig. 4a) with the pictures of real galaxies (HST
photograph of NGC 4321, Fig. 4c, or enhanced AG
photograph of NGC 5457, Fig. 1 of Paper I).
We are well aware that our treatment of the galactic
dynamics in terms of a stationary velocity eld in a
rotating reference frame is not fully satisfying. It ignores not
only a feedback of stellar evolution on the dynamics, but
above all it discards the alternative of a chaotic flow in
the region under investigation. It is indeed obvious that
chaos in the dynamics must a ect both qualitatively and
quantitatively the large scale and the small scale pattern
of the arms. Dynamical chaos is expected to increase the
fragmented character of the visible arm structure, and
presumably could become measurable via fractal dimensions
of the arms. The ideal model of galactic structure should
therefore not divorce the dynamics from star formation
and evolution. In principle, to the extent that gravitation
can be treated as an external force
(as it is essentially
done in standard self{consistent dynamic models; cf. for
instance the model by Patsis et al. 1991)
, the
gravitational e ects could be incorporated in the CA model in a
formally straighforward way. Unfortunately, such a
procedure is highly time consuming.
The second quantitatively relevant aspect of our model
is that it demonstrates that global observational galactic
data contain easily accessible information on stellar
evolution rates (including mass loss). We nd that the 9 model
parameters (5.1) are determined by a system of 9
equations which can be approximated by linear equations (Eqs.
3.3, 4.4) whose coe cients are observable. Four among the
latter are fractal dimensions associated with the spatial
distribution of the galactic species; while not yet
available, there should be no major di culty in deriving them
from surveys of nearby galaxies. We hope to address this
problem in the future. The remaining parameters are mass
fractions of the di erent galactic components which can be
found in the literature. Tables 2 and 3 provide a list of
tentative input coe cients to carry out the inversion problem.
It would be useful to recalibrate these coe cients, using
a more extended array of test models covering also wider
ranges of transition probabilities.
Finally, we should keep in mind that the arti cial
construction of the initial state as incorporated in our
program should be replaced by a simulation of a physically
more realistic process of molecular cloud formation
(aggregation model, for instance). Detailed physics of star
formation may be implemented as well, such as a
simulation of a hierarchical fragmentation process of the cloud
complex. Standard hydrodynamic simulations indicate in
fact that the isothermal collapse of rotating dense clouds
may lead to a succession of levels of condensations in
condensations
(Boss 1991)
. Such a mechanism could enhance
local irregularities in the arm structure, thereby leading
to a higher dimension F of the arm pattern. We plan to
address these questions in a later work.
Acknowledgements. The authors are indebted to G. Barbaro
for vital information on the parameters of stellar evolution,
to H. Robe for his careful reading of a rst version of this
paper, and to J.P. Swings for providing them with an early
HST photograph of M100. They wish to thank the referee Dr.
Comins for his comments which helped improve the
presentation of this paper. J.P. gratefully acknowledges a Royal Society
{ FNRS European Exchange Fellowship and nancial support
from the FNRS (Belgium).
Bak P. , Tang C. , Wiesenfeld K. , 1987 , Phys. Rev. Lett . 59 , 381
Bak P. , Tang C. , Wiesenfeld K. , 1988 , Phys. Rev. A38 , 364
Bak P. , Chen K. , Creutz M. , 1989 , Nat 342 , 780
Boss A.P. , 1991 , Nat 351 , 298
Cannizzo J. , Goodings D. , Mattei J.A. , 1990 , ApJ 357 , 235
Feller W. , 1971 , An Introduction to Probability Theory and its Applications , Vol. II. Wiley, New York
Gerola H. , Seiden P.E. , 1978 , ApJ 223 , 129 (quoted as GSS)
Gerola H. , Seiden P.E. , Schulman L.S. , 1980 , ApJ 242 , 517 (quoted as GSS)
Glansdor P. , Prigogine I. , 1971 , Thermodynamic Theory of Structure, Stability and Fluctuations. Wiley, London, Chapts. XIV and XV
Goldreich P. , Lynden{Bell D., 1965 , MNRAS 130 , 125
Heisenberg W. , von Weizs¨acker C.F. , 1948 , Z. Phys . 125 , 290
Herrmann H.J. , 1983 , J. Phys. A16 , L611
Lejeune A. , 1993 , in: Cellular Automata: Prospects in Astrophysical Applications . In: Perdang J.M. , Lejeune A . (eds.), World Scienti c, Singapore, p. 323
Lejeune A. , Perdang J. , 1991 , in: Computational Physics. In: Tenner A. (ed.), World Scienti c, Singapore, p. 407
Mueller M.W. , Arnett W.D., 1976 , ApJ 210 , 670 (quoted as MA)
Patsis P.A. , Contopoulos G. , Grosbol P. , 1991 , A &A 243 , 373
Perdang J. , 1990 , Vistas Astron . 33 , 249
Perdang J. , 1991 , in: Applying Fractals in Astrophysics, Lecture Notes in Physics m3 . Springer, Berlin. p. 1
Perdang J. , 1993 , in: Cellular Automata: Prospects in Astrophysical Applications . In: Perdang J.M. , Lejeune A . (eds.), World Scienti c, Singapore, p. 3
Perdang J. , 1996 , A Family of Endogenous Growth Models (in preparation)
Perdang J. , Lejeune A. , 1996 , (Paper I)
Sandage A. , Tammann G.A. , 1981 , A Revised Shapley{ Ames Catalog of Bright Galaxies , Carnegie Institute of Washington Publication 635, Washington D.C. (quoted as RSAC)
Sandage A. , Bedke J. , 1988 , Atlas of Galaxies, NASA SP{496 , Washington D.C. (quoted as AG)
Schulman L.S. , Seiden P.E. , 1983 , Ann. Israel Phys. Soc . 5 , 251 (quoted as GSS)
Schulman L.S. , Seiden P.E. , 1986 , Sci 233 , 425 (quoted as GSS)
Schulman L.S. , 1993 , in: Cellular Automata: Prospects in Astrophysical Applications . In: Perdang J.M. , Lejeune A . (eds.), World Scienti c, Singapore, p. 294 ( quoted as GSS )
Seiden P.E. , 1985 , in: The Milky Way Galaxy, IAU Symposium 106 . In: van Woerden H., Allen R.J. , Burton W.B. (eds.). Reidel, Dordrecht (quoted as GSS) , p. 551
Seiden P.E. , Gerola H. , 1979 , ApJ 233 , 56 (quoted as GSS)
Seiden P.E. , Gerola H. , 1982 , Fundam. Cosmic Phys. 7 , 241 (quoted as GSS)
Seiden P.E. , Schulman L.S. , 1990 , Adv. Phys. 39 , 1 (quoted as GSS)
Seiden P.E. , Schulman L.S. , Gerola H. , 1979 , ApJ 232 , 702
Seiden P.E. , Schulman L.S. , Feitzinger J.V. , 1982 , ApJ 253 , 91 (quoted as GSS)
Stau er D. , 1985 , Introduction to Percolation Theory, Taylor and Francis, London
von Weizs¨acker C.F. , 1951 in: Festschrift zur Feier des zweihundertj¨ahrigen Bestehens der Akademie der Wissenschaften in G¨ottingen , Mathematisch{Physikalische Klasse. Springer, Berlin, p. 86