text:
# Observations of Disk Galaxy Evolution

## 1. Scope of This Review

This review is an attempt to summarize current observational constraints on evolution in disk galaxies. Because of space limitations, we will mostly touch upon recent highlights, and do so in a simplified way that treats morphological, luminosity/size, and dynamical evolution independently (even though they are probably deeply interlinked). Throughout this review we will espouse the philosophy that the best way to understand how galaxies evolve is to watch them evolve in situ at high redshifts. Since our remit is to focus on evolution in disks, it is necessary to direct most of our attention to observations capable of distinguishing between disk galaxies and other objects in the field at high redshifts. We therefore unapologetically emphasize a morphological point of view, and focus mostly on recent work undertaken with the Hubble Space Telescope (HST).

## 2. A User’s Guide to Interpreting High-z Disk Galaxy Observations

Present observations of galaxy evolution extend to redshifts of at least $`z4`$ (Steidel et al. 2001), and possibly much higher (Dickinson et al. 2000). This can be something of a double-edged sword when it comes to the interpretation of observations. The tremendous benefits which come from being able to observe evolution directly at high redshifts come at the price of an extra layer of complexity in the interpretation of these observations, since it necessary to disentangle true evolution from a myriad of possible observational biases. These biases ricochet through all studies of high-$`z`$ galaxy evolution, and how they are handled conditions the interpretation of a large body of work. For example, consider Figure 1, which illustrates many of the evolving generic characteristics of high-redshift disks. How much of the spectacular evolution in the appearance of the typical disk systems as a function of redshift shown this figure is due to observational selection effects?

We now know enough about high-$`z`$ galaxy morphology to let us really come to grips with this question. In order of increasing severity, the three effects which most impact high-redshift galaxy observations are:

1. Resolution/pixellation effects. In most cases these are only a minor source of bias. HST’s resolving power and relatively small pixels, coupled with the inflection point in angular diameter distance as a function of redshift, mean than in many cases HST images actually have higher resolution than do calibration images of nearby galaxies. For example, under typical seeing conditions, only those galaxies at $`z<0.05`$ in the Sloan Digital Sky Survey have resolution equal to or better than that of HST at $`z=1`$.
2. Bandshifting of the rest wavelength of observation. Since most morphological studies of high-redshift galaxies have been done using HST’s F814W filter (approximately $`I`$-band), these studies correspond to rest-frame $`B`$-band (where local morphology is best understood) at $`z0.7`$. Therefore bandshifting effects are probably quite innocuous at $`z<1`$, but grow more worrisome at higher redshifts. As described below, there is already ample evidence for strong evolutionary effects by $`z=0.7`$, ruling out the extreme possibility (considered seriously early on) that most of the morphological changes witnessed on HST images are due to bandshifting effects. Furthermore, even the limited amount of imaging made with the HST Near-Infrared Camera Multi-Object Spectrometer (NICMOS) shows quite clearly that most $`z>1`$ objects classified as peculiar in $`I`$-band remain peculiar when imaged at rest optical wavelengths (Dickinson 1999). Finally, it is worth mentioning that the effects of bandshifting can be calibrated out of imaging data using pixel-by-pixel K-corrections provided deep data are available in several filters so that an estimate of an individual pixel’s color can be made and a template spectrum assigned (Abraham 1997; Bouwens, Broadhurst & Silk 1998).
3. The interplay between cosmological surface brightness dimming and evolution in stellar populations. This is certainly the most poorly constrained effect. Cosmological surface brightness diminution goes as $`(1+z)^3`$ (Tolman dimming), but this dimming in apparent surface brightness is partially offset by an (uncertain) evolutionary brightening in intrinsic surface brightness which operates in the opposite direction. As we will now show, evolutionary synthesis predictions can go some way toward making predictions for higher redshifts.

Figures 2 and 3 are an attempt to come to grips simultaneously with the effects of bandshifting, cosmological surface brightness dimming, and evolution. These figures plot the signal-to-noise/pixel as a function of age and redshift for a projected stellar mass density corresponding to that of the Milky Way at the solar radius. The region shown in gray is invisible even after a modest amount of binning. We assume a stellar population with an exponentially declining star-formation history characterized by an e-folding timescale of 1 Gyr111Note that when modeling individual pixels in resolved galaxies, a single short-timescale exponential can be used to provide a quite general description of star-formation characteristics. Longer exponential timescale (or even constant) star-formation models used to model the integrated properties of disks are well-approximated by the spatial superposition of multiple resolved short timescale exponentials with a range of ages. . In Figure 2 we model 15ks (i.e. approximately five orbit) HST integrations typical of moderately deep HST imaging in the Medium Deep Survey (Griffiths et al. 1994), the Hawaii Deep Survey (Cowie, Hu & Songaila 1995), and the LDSS/CFRS (Brinchmann et al. 1998) survey. We use the known characteristics of the Wide-Field Planetary Camera 2 (WFPC2) in the F814W filter, as well as the projected characteristics of the Advanced Camera for Surveys (ACS) in the SDSS $`z^{}`$ band, and the refurbished NICMOS in F160W (approximate $`H`$-band)222Both ACS and NICMOS are scheduled for installation aboard HST during Servicing Mission 3, currently scheduled for Feb/March 2002.. The bottom line is that Figure 2 shows that at $`z<1`$, there is relatively little bias in optical observations undertaken with HST. At this redshift, projected mass densities corresponding the intermediate regions of large disk galaxies are observable for around 70% of the maximum possible age for the galaxy (defined as the age of the Universe at the epoch of observation). However, the second panel in the figure shows that by $`z=2`$ biases with WFPC2 are severe: structural information at optical wavelengths comes only from the densest regions in galaxies (such as bulges), or from the youngest parts of galaxies (such as bright star-formation complexes with ages of less 1 Gyr). By $`z=2`$ it is simply no longer possible to use WFPC2 optical data to compare fairly fundamental characteristics, such as the size and morphology of galaxies, with those of local objects, at least without much deeper imaging than is possible with 5 orbits on HST. It is interesting to note that NICMOS allows fair mass densities to be probed at $`z=2`$, albeit at the cost of limited field of view.

What is the corresponding situation for the much deeper Hubble Deep Field (HDF) observations? Figure 2 shows that imaging at HDF-like depths with WFPC2 probes solar-neighborhood-like masses out to $`z=2`$ about as well as 5-orbit HST imaging does at $`z=1`$. At these depths ACS allows somewhat unbiased observations to be made out to $`z=3`$, while NICMOS is completely unbiased at this redshift.

Armed with this machinery, we can answer with some confidence the rhetorical question posed at the beginning of this section. Because the systems shown in Figure 1 are seen at $`z<1`$ and these images are from the HDF, it is clear that their morphological peculiarities are intrinsic, and thus a fair reflection of the underlying state of their stellar populations. If these systems had been at somewhat higher redshifts, but still at $`z<2`$, the interpretation of the structural changes would depend on whether or not they were consistent with a simple relative brightening (or fading) of known morphological subcomponents. This might be a difficult determination to make in some cases, but at least we could be reasonably confident that most of the galaxy was not hidden below the surface brightness detection threshold, which might easily be the case at higher redshifts.

## 3. Luminosity/Size Evolution

Because most studies of luminosity and size evolution are based on diagnostic diagrams showing central surface brightness plotted as a function of scale size, distinguishing between pure luminosity evolution and pure size evolution is difficult. For simplicity I will treat these two effects together.

Based on ground-based imaging obtained with the CFHT, Schade et al. (1996a) found that galaxy disks at $`0.5<z<1.1`$ have a mean rest-frame, inclination-corrected central surface brightness of $`\mu (B)_{AB}=19.8\pm 0.1`$ mag arcsec-2, which is $`1.6`$ mag brighter than the Freeman (1970) value. At low redshifts ($`0.2<z<0.5`$) these authors found a mean surface brightness of $`\mu (B)_{AB}=21.3\pm 0.25`$ mag arcsec-2, which is consistent with the Freeman value. Assuming most of this effect is due to luminosity (rather than size) evolution, the broad-brush picture painted is one of $`B`$-band evolution proceeding roughly linearly (in magnitude units) with redshift. This remains the most popular prescription for the luminosity evolution of disks, and is quite consistent with most recent studies, the most glaring exception being Simard et al. (1999), as described below.

The environmental dependance of disk surface brightness evolution has been explored in another ground-based sample studied by Schade et al. (1996b). These authors studied the differential properties of cluster and field galaxies from the CNOC survey. At redshifts of (0.23, 0.43, 0.55) the disk surface brightness in cluster late-type galaxies is higher in the $`B`$ band by $`\mathrm{\Delta }\mu (B)_0=(0.58\pm 0.12,1.22\pm 0.17,0.97\pm 0.2)`$mag, respectively, relative to Freeman’s law, which is consistent with the disk brightening seen in the field. However, Schade et al. emphasize that their sample of cluster galaxies is dominated by objects at large distances (up to 3 Mpc) from the dense cluster cores, so the implications of these findings for studies of the Butcher-Oemler effect and the morphology-density relation really must await investigations that thoroughly probe galaxy properties as a function of clustercentric radius.

The best HST-based follow-up investigations of luminosity/size evolution are the studies by Lilly et al. (1998) and Simard et al. (1999), which are based on samples with known redshifts obtained from the LDSS/CFRS and DEEP surveys, respectively. Lilly et al. (1999) examined the two-dimensional surface brightness profiles for 341 objects at $`0<z<1.3`$. The size function of disk scale lengths in disk-dominated galaxies was found to stay roughly constant to $`z1`$, at least for large disks where the sample was most complete and where the disk and bulge decompositions were most reliable. Lilly et al. define these objects to be those with exponential scale lengths $`>3.2h_{50}^1`$ kpc, and find good agreement between the scale lengths measured in their sample and a local calibration sample. This important result suggests that typical disks cannot have grown substantially with cosmic epoch since $`z=1`$, unless this growth is balanced by a corresponding destruction in the number of large disks through merging. It is probably worth emphasizing that while the disks in this sample appear to be the same size at high-$`z`$, they show significant differences from local disks: they have higher $`B`$-band disk surface brightnesses, bluer overall $`(UV)`$ colors, higher \[O II\] equivalent widths, and less regular morphologies than local counterparts. Lilly et al. conclude this is consistent with an increase in the star formation rate by a factor of about three to $`z0.7`$, in agreement with the expectations from recent models for the evolution of the Milky Way disk.

The measurements presented by Simard et al. (1999) are seemingly in close agreement with those given by Lilly et al. (1998), but these authors differ greatly in their interpretation. The Simard et al. study is based on a total sample of 190 field galaxies ($`I_{814}23.5`$ mag and $`0.1<z<1.1`$), and these authors find that if selection effects are ignored the mean disk surface brightness (averaged over all galaxies) increases by $`1.3`$ mag from z=0.1 to 0.9, in fairly good agreement with the Lilly et al. (1998) and Schade et al. (1995a,b) results. However, by creating mock catalogs and simulating the appearance of galaxies, these authors are able to quantify the surface brightness selection function for the sample as a whole, and conclude that most of this change could be due to comparing low-luminosity galaxies in nearby redshift bins to high-luminosity galaxies in distant bins. After allowing for this effect, these authors find that no discernible evolution remains in the disk surface brightness of bright ($`M_B<19`$) disk-dominated galaxies, and that at all redshifts to $`z1`$ the magnitude-size envelope for the bulk of galaxies is consistent with that of local galaxies. While the derivation of a formal selection function for the dataset based on simulations is laudable, we remain suspicious of the final result because the accessible area of size-surface brightness parameter space in a deep survey would seem to be larger than the region of parameter space occupied by real galaxies. For example, Driver (1999) uses the Hubble Deep Field to derive the “local” ($`0.3<z<0.5`$) bivariate brightness distribution of field galaxies within a modest (326 Mpc3, 47 galaxies) volume-limited sample. The sample probes the underlying galaxy population within the limits $`0.3<z<0.5`$, $`21.3<M_B<13.7`$ mag, $`18.0<\mu _B<24.55`$ mag arcsec-2. Galaxies occupy a relatively narrow locus in this parameter space, and Driver concludes that luminous low-surface brightness galaxies account for less than 10% of the $`L^{}`$ population. We conclude that it is unlikely that surface brightness selection plays a significant role in conditioning the numbers for photometric disk evolution, and that “what you see is what you get” when interpreting the size-magnitude diagram in the present generation of deep surveys which reach $`L^{}`$ at $`z<1`$.

Photometric redshifts can also be employed to constrain the sizes of high-redshift disk galaxies at higher redshifts than is possible spectroscopically, though with correspondingly larger uncertainties. In an interesting paper, Giallongo et al. (2000) compare the intrinsic sizes of the field galaxies with $`I<26`$ mag in the Hubble Deep Field and the NTT Deep Field to the predictions from a Cold Dark Matter (CDM) model, finding that the distribution of the galaxy disk sizes at different cosmic epochs is within the range predicted by typical CDM models, modulo an excess of small-size disks at $`z>0.5`$. This result is similar to earlier results obtained by Bouwens et al. (1998a,b) who adopted the novel (and, in our view, excellent) technique of “cloning” foreground galaxies to model their high-redshift counterparts in order to determine that faint HDF galaxies ($`I>24`$ mag) are much smaller, more numerous, and less regular than predicted by no-evolution models. These authors find that pure luminosity evolution in a flat $`\mathrm{\Lambda }`$-dominated cosmology can account for the properties of the brightest and biggest galaxies in the HDF, but that the excess of small faint galaxies remains difficult to explain. In any case, when speculating on the nature of the small, faint galaxy population on HST images we caution that it may not be too meaningful to try to understand such systems in the context of local well-ordered disks. As described in the next two sections, disk warping and morphological peculiarity are seen to become more important with increasing redshift, which makes one suspicious that the notion of an isolated well-ordered disk may disappear at $`z>1.5`$.

## 4. Dynamical Evolution

Dynamical investigations of disk galaxy evolution are in their infancy. Ambitious programs (such as the DEEP project) are in progress, but the published results from these comprehensive investigations have yet to appear. Vogt et al. (1997) reports initial results from a Keck-based study of the internal kinematics of eight objects (spanning the range $`21.8<M_B<19.0`$ mag for $`H_0=75`$ km s-1 Mpc-1 and $`q_0=0.05`$, i.e. predeominantly sub-$`L^{}`$). A small offset of $`<0.4`$ mag in $`B`$-band is seen, but no obvious change in the shape or slope of the relation with respect to the local Tully-Fisher relation is detected. The offset seems entirely consistent with the passive evolution described in the previous section. The DEEP project has reported obtaining kinematical information for over 100 high-redshift galaxies, and the preliminary results from this study have already been reported several times at meetings (Vogt 2000, 2001; Faber et al 2001). These data form a high redshift Tully-Fisher relation that spans five magnitudes and extends to well below $`L^{}`$, and the basic result reported in Vogt et al. (1997) seems unchanged: no obvious change in shape or slope relative to the local Tully-Fisher relation is seen, while a smooth but modest increase in disk surface brightness as a function of redshift is detected.

Since there appears to be strong kinematical evidence for luminous, large disk systems at high redshifts, it is interesting to enquire as to the highest redshift at which a bona-fide luminous disk galaxy has been seen. The current record holder is the object reported by van Dokkum and Sanford (2001), who discovered a rapidly rotating disk galaxy at $`z=1.34`$ with maximum rotation velocity $`V_{\mathrm{max}}=290\pm 20`$ km $`s^1`$, on the basis of resolved \[OII\] emission. The galaxy has $`M_V=22.4\pm 0.2`$ ($`3L^{}`$) in the rest frame, and its spectral energy distribution is very well fitted by that of a redshifted Sb/c galaxy. The observed kinematics, morphology (both HST optical and ground-based infrared), and spectral energy distribution are consistent with a massive spiral galaxy. Interestingly, the lower limit on $`V_{\mathrm{max}}`$ gives an upper limit on the offset from the present-day Tully-Fisher relation, and the galaxy is overluminous by less than $`0.7\pm 0.4`$ mag in rest-frame $`V`$-band. Of course, when it becomes possible (in the not-too-distant future) to construct samples of $`z>1`$ galaxies to probe the high-$`z`$ Tully-Fisher relation one will need to be very careful, since if disks are in their infancy at these redshifts they may be so tilted and chaotic that it is not clear what the physical nature of the “rotation velocity” really is. Resolved spectroscopy with Integral Field Units (see the paper by Bershady in this volume) may be essential for the success of dynamical programs at high-$`z`$.

## 5. Morphological Evolution

The number counts of disk galaxies are in good agreement with the predictions of passive-evolution models to at least $`I_{AB}=25.5`$ mag, the effective limit for crude morphological classification (into very broad early/disk/peculiar bins) in the Hubble Deep Fields. At this magnitude limit around 30% of all galaxies in the field are classified as “peculiar/merging”, i.e. they fall outside the conventional framework of the Hubble System altogether. As of this writing this basic result has been derived independently using at least three different methods for automatic classification (Abraham et al. 1996; Driver et al. 1998; Marleau & Simard 1998) and at least five sets of human eyeballs (Abraham et al. 1996; van den Bergh et al. 1996; Driver et al. 1998; Ellis, Abraham & Dickinson 2001), and the basic result therefore seems quite solid. However only about 5% of the galaxies in the HDF have known redshifts. Since the deepest existing redshift survey (the Caltech Faint Galaxy Redshift Survey, CFGRS; Cohen et al. 1996) begins to become seriously incomplete around $`R=24`$ in the HDF ($`R=23`$ in the HDF flanking fields), existing HST imaging data go significantly deeper than do spectroscopic redshifts. Several photometric redshift-based studies lend support to the passive evolution picture for disks, but there really is nothing like a proper spectroscopic redshift to lend credibility to this kind of work, so it came as something of a relief when Brinchmann et al. (1998) was able to show (on the basis of HST images for 341 galaxies with $`I_{AB}<22.5`$ mag of known redshift from the CFRS/LDSS redshift survey) that the proportion of spectroscopically-confirmed disks is indeed in good agreement with passive evolution predictions.

More fine-grained “precision morphology” (which can subdivide spiral galaxies into various classes, for example) still requires a trained human eye, and also requires complete redshift information and multi-filter data in order to synchronize data to a uniform rest wavelength. Using the CFGRS data, van den Bergh et al. 2000 were able to show that the fraction of intermediate/late- type (Sb - Ir) galaxies in the HDF the drops by a factor of $`2`$ out to $`z1`$, presumably because many late-type galaxies viewed at a look- back times of $`8`$ Gyr have been classified as “peculiar” or “merger”, while the fraction of such peculiar or merging galaxies exhibits a monotonic increase from 5% at $`z0`$, to 10% at $`z0.4`$, 19% at $`z0.7`$ and 30% at $`z0.95`$. It should be once again emphasized that the numbers of galaxies in this investigation is small, and the conclusions which can drawn from these data are correspondingly uncertain. It is also worth noting that the merger rate estimated from the HDF on the basis of morphological considerations is in reasonable agreement with the merger rate inferred from HST imaging of the CFRS fields, and from pair counts at lower redshifts.

Beyond $`z1`$, one must observe galaxies in the infrared in order to study them at familiar optical wavelengths in their rest frames. The NICMOS camera on HST had this capability, although its small field of view and limited lifetime meant that relatively few $`z>1`$ field galaxies were observed in the near infrared before the instrument ceased functioning. However, over the short lifetime of this camera it was possible to demonstrate that most galaxies with peculiar optical morphologies remain peculiar when viewed in the near infrared, at least out to $`z=3`$ (Dickinson 1999), and as shown in Figure 2 this really ought to be a fair test of the nature of these objects. It follows that the strange appearance of the peculiar objects at these redshifts is not merely the consequence of “morphological K-corrections” (i.e. of observing only irregularly distributed sites of star formation at rest ultraviolet wavelengths, while missing the bulk of the galaxy). The peculiar appearance of these objects reflects a genuinely irregular structural state in these galaxies. What is unclear, however, is what proportion of these objects are bona fide disks. Furthermore, in the most distant objects (i.e. at $`z>3`$), it may be the case that we really are only seeing the very tip of the baryonic mass iceberg (i.e. only the youngest stars in the galaxies, grouped together in irregular star formation complexes).

As described by Merrifield in this volume, a number of studies have also claimed that the percentage of barred spirals seems to decrease towards larger look-back times (van den Bergh et al. 1996; Abraham et al. 1999; van den Bergh et al. 2000). However, an alternative interpretation of the absence of barred and grand design spirals at high redshifts is that the features that define these categories of objects locally are undetectable outside of the nearby Universe. A useful tool in testing this is what Abraham & Merrifield (2000) have dubbed “Hubble space”, a quantitative two-parameter description of galactic structure that maps closely on to Hubble’s original tuning fork. By “artificially redshifting” samples of local galaxies and mapping them onto Hubble space it is possible to assess the extent to which characteristics such as bar strength should remain visible at high redshifts. Van den Bergh et al. (2002) present such an analysis for a statistically fair set of local spirals taken from the Ohio State University Bright Spiral Galaxy Survey (Frogel, Quillen & Pogge 1996). Results from this study are shown in Figure 4, where it is seen about $`2/3`$ of all luminous local SB galaxies would still classified as SB or S(B) at the resolution and noise level of the Hubble Deep Field, dropping to around 50% in the flanking fields. Therefore image degradation at high redshift can only account for a fairly small part of the observed decrease in the fraction of SB galaxies with increasing redshift, although we caution that the total the number of suitably bright (and suitably inclined) galaxies in the HDF data is small and more observations are needed to better establish the putative demise of barred spirals at high redshifts.

Explaining why the Hubble classification scheme appears to “evaporate” at redshifts $`z>0.3`$ is a challenge that can only be met by a greater understanding of galaxy dynamics at high redshifts. The gas fraction in typical spiral galaxies may have decreased by a large fraction between $`z=1`$ and $`z=0`$, and encounters with giant molecular clouds, spiral arms and perhaps clumps in the dark matter halo may have increased the velocity dispersion within the disk at high redshifts. Fuchs & von Linden (1998) suggest that either of these effects would impact the stability of spiral disks in the distant past, but the data simply do not yet exist to allow astronomers to determine whether either of these effects plays a dominant role. At this stage all that is clear is that the morphology of spiral galaxies is evolving rapidly, and systematically, even at quite low redshifts. Familiar types of galaxies, such as barred and grand-design spirals, appear to be relatively recent additions to the extragalactic zoo. The nature of the many morphologically peculiar galaxies at high redshift and their relationship to low-redshift disk systems remains a complete mystery. These objects might be mergers, proto-galaxies, new classes of evolved systems, or a combination of all three. Clearly we are only beginning to glimpse the morphology of galaxies in the distant Universe, and are just at the stage where we are groping to fit them into a physically meaningful organizational structure.

## 6. Conclusions and Speculations

It seems that the observable universe falls into two domains. At $`z<1`$, 3/4 of luminous galaxies are disks. Only a few of the remaining nearby ellipticals are presently forming stars at a significant rate. So almost all star formation at $`z<1`$ is taking place in disks. In these objects star formation is occuring in a manner consistent with passive evolution and with the slow build-up of stellar populations. On the other hand most star formation at higher redshifts is taking place in mergers, or in rather chaotic objects. Is there a connection between the transition from a merger-driven mode of star formation and the change in the slope of the star-formation history of the Universe plot (the “Madau plot”) at $`z11.5`$? It is not impossible that the changing mode of star formation is the underlying cause of the inflection in the Madau plot. If early disks are being destroyed at $`z>1`$ in this process then one would not expect the Galaxy to have an old metal-poor thin disk. In any case, we are only now entering the phase where optical and near-IR data will be deep enough to allow fair samples of disks at $`z>1`$ to be studied. The clear expectation is that studies of such objects will show that disk sizes decrease with increasing $`z`$ because of tidal truncation and destruction, and because gaseous disks are more prone to destruction by the winds induced by violent central bursts of star formation. In any case there is so much evidence that our own Galaxy’s disk grew from inside out that it would be surprising if this was not seen in other systems at high redshifts.

By $`z=1`$ it appears that most of the big disks have already settled into place, although their morphologies remain outside the formal Hubble sequence. At $`z0.7`$ about 2/3 of all Sbc-Sc galaxies are peculiar. Clearly late-type disks have mostly not yet had time to revert to their Hubble prototype form. On the other hand, the fraction of peculiar Sa’s at $`z0.7`$ appears to be minute. Some of this is no doubt due to the fact that it might be difficult to see some forms of peculiarity in compact early-type galaxies, but it seems important to determine whether Sa and Sab galaxies at $`z0.7`$ are basically similar in overall structure to nearby ones, or whether the nature of spiral structure changes with look-back time in a different way in Sa and Sc galaxies. The best preliminary evidence would seem to suggest that in Sc galaxies the spiral arms become patchier and more chaotic with increasing look-back time. On the other hand spiral arms appear to become less prominent in early-type galaxies at large look-back times. In any case, it is tempting to speculate that some of these effects might be related to the apparent paucity of barred spirals at $`z>0.7`$. Perhaps high-redshift disks are too dynamically hot to become unstable in a global bar forming mode.

## References

Abraham, R. G. & Merrifield, M. R. 2000, AJ, 120, 2835

Abraham, R. G., Merrifield, M. R., Ellis, R. S., Tanvir, N. R. & Brinchmann, J. 1999, MNRAS, 308, 569

Abraham, R. G. 1997, In “The Ultraviolet Universe at Low and High Redshift: Probing the Progress of Galaxy Evolution”, AIP Conference Proceedings, v.408., p.195

Abraham, R. G., Tanvir, N. R., Santiago, B. X., Ellis, R. S., Glazebrook, K., & van den Bergh, S. 1996, MNRAS, 279, L47

Bershady, M. 2002, this volume

Bouwens, R., Broadhurst, T., & Silk, J. 1998, ApJ, 506, 557

Brinchmann, J. et al. 1998, ApJ, 499, 112

Cohen, J. G., et al. 1996, ApJ, 471, L5

Cowie, L. L., Hu, E. M., & Songaila, A. 1995, AJ, 110, 1576

Dickinson, M. 1999. In After the Dark Ages: When Galaxies were Young (the Universe at $`2<z<5`$). Edited by S. Holt and E. Smith. American Institute of Physics Press, p. 122.

Dickinson, M. et al. 2000, ApJ, 531, 624

Driver, S. P., Windhorst, R. A., & Griffiths, R. E. 1995, ApJ, 453, 48

Driver, S. P., Fernandez-Soto, A., Couch, W. J., Odewahn, S. C., Windhorst, R. A., Phillips, S., Lanzetta, K., & Yahil, A. 1998, ApJ, 496, L93

Driver, S. P. 1999, ApJ, 526, L69

Ellis, R. S., Abraham, R. G., & Dickinson, M. 2001, ApJ, 551, 111

Faber, S. M., Phillips, A. C., Simard, L., Vogt, N. P., & Somerville, R. S. 2001, ASP Conf. Ser. 230: Galaxy Disks and Disk Galaxies, 517

Frogel, J. A., Quillen, A. C., & Pogge, R. W. 1996, ASSL Vol. 209: New Extragalactic Perspectives in the New South Africa, 65

Giallongo, E et al. 2000, ApJ, 530L, 73.

Griffiths, R. E. et al. 1994, ApJ, 437, 67

Le Fèvre, O. et al. 2000, MNRAS, 311, 565

Lilly, S. J., Le Fevre, O., Crampton, D., Hammer, F., & Tresse, L. 1995, ApJ, 455, 50

Lilly, S. et al. 1998, ApJ, 500, 75

Marleau, F. R. & Simard, L. 1998, ApJ, 507, 585

Merrifield M. R. 2002, this volume

Schade, D., Carlberg, R. G., Yee, H. K. C., Lopez-Cruz, O., & Ellingson, E. 1996, ApJ, 465, L103

Schade, D., Lilly, S. J., Le Fevre, O., Hammer, F., & Crampton, D. 1996, ApJ, 464, 79

Simard, L. et al. 1999, ApJ, 519, 563

Steidel, C. C., Pettini, M., & Adelberger, K. L. 2001, ApJ, 546, 665

van den Bergh, S. et al. 2002. (astro-ph/0202444).

van den Bergh S., Abraham R.G., Ellis R.S., Tanvir N.R. & Santiago B.X. 1996, AJ, 112, 359

van den Bergh, S., Cohen, J. G., & Crabbe, C. 2001, AJ, 122, 611

van den Bergh, S., Cohen, J. G., Hogg, D. W., & Blandford, R. 2000, AJ, 120, 2190

van Dokkum, P. G. & Stanford, S. A. 2001, ApJ, 562, L35

Vogt, N. P. 2001, ASP Conf. Ser. 240: Gas and Galaxy Evolution, 89

Vogt, N. P. 2000, ASP Conf. Ser. 197: Dynamics of Galaxies: from the Early Universe to the Present, 435

Vogt, N. P. et al. 1997, ApJ, 479, L121


format: text

text:
# A comparison of different methods in the study of dynamical fluctuations in high energy e+e-collisions

## I Introduction

Since the observation of exotic multiparticle events in Cosmic Ray experiments CR  and especially the discovery of unexpectedly large local fluctuations recorded by the JACEE collaboration JACEE , the investigation of nonlinear phenomena (NLP) in high energy collisions has attracted much attention Kittel . One of the signals of these NLP is the fractal property Mandelbrote  of the multiparticle final states in high energy collisions. This property can be characterized by the anomalous scaling of the probability moments defined as

$`C_q(M)={\displaystyle \frac{1}{M}}{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{p_m^q}{p_m^q}},`$ (1)
where a certain phase space region $`\mathrm{\Delta }`$ is divided in a proper way (isotropically for a self-similar fractal while anisotropically for a self-affine fractal FLPRL ) into $`M`$ sub-cells, $`p_m`$ is the probability for a particle to fall into the $`m`$th sub-cell, $`\mathrm{}`$ denotes the average over an event sample. If the $`C_q(M)`$’s have a power law behaviour with the diminishing of phase space scale:

$`C_q(M)M^{\varphi _q}(M\mathrm{}),`$ (2)
then it is said to be anomalous scaling and the system is a fractal.

In real experiments the probability $`p_m`$ is unobservable and the corresponding moments $`C_q`$ is unaccesible. This problem has been solved by Białas and Peschanski BP , who are able to show that the factorial moments

$`F_q(M)={\displaystyle \frac{1}{M}}{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{n_m(n_m1)\mathrm{}(n_mq+1)}{n_m^q}}`$ (3)
are equal to the probability moments $`C_q`$ provided the statistical fluctuations are Poissonian. Thus the scaling property of factorial moments, sometimes called intermittency, becomes a central problem in the study of nonlinear phenomena in high energy collisions.

Various methods have been developed in this study.

Firstly, people noticed that in the definition Eq.(3) of factorial moments the average over event sample (vertical average) is carried out first and then is the average over the $`M`$ sub-cells (horizontal average). It was proposed to exchange the order of these two averages and define the horizontal factorial moments (HFM) HFM  as

$`F_q^{(\mathrm{H})}(M)={\displaystyle \frac{M^1\underset{m=1}{\overset{M}{}}n_m(n_m1)\mathrm{}(n_mq+1)}{M^1\underset{m=1}{\overset{M}{}}n_m^q}}.`$ (4)
Accordingly, the $`F_q`$ defined in Eq.(3) is called vertical factorial moments (VFM).

Note that the equality of factorial moments $`F_q`$ and probability moments $`C_q`$ has been proved only for the VFM. Therefore, in the study of the nonlinear phenomena — fractal property of multiparticle system, the HFM is appropriate only when it is equal to VFM. We will see in the following that this equality holds in some cases but does not hold in some other cases.

Secondly, various methods have been proposed to correct for the unflatness of the phase-space variable distributions. One is to divide the factorial moments by a factor $`R_q`$ Fialkowski

$`F_q^\mathrm{C}={\displaystyle \frac{F_q}{R_q}},R_q={\displaystyle \frac{M^1\underset{m=1}{\overset{M}{}}n_m^q}{M^1\underset{m=1}{\overset{M}{}}n_m^q}}.`$ (5)
Another one is to change the phase space variable $`x`$ into the corresponding cumulant variable $`x_\mathrm{c}`$ before calculating the factorial moments. The cumulant variable is defined as Cummulant

$`x_\mathrm{c}={\displaystyle \frac{_{x_{\mathrm{min}}}^x\rho (x)dx}{_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}\rho (x)dx}},`$ (6)
where $`x_{\mathrm{min}}`$ and $`x_{\mathrm{max}}`$ are the two boundaries of the $`x`$ region, respectively.

Another problem arises while carrying on this kind of study in e+e-collisions. In these collisions the thrust (or sphericity) axis is chosen as the $`z`$ axis (longitudinal axis) to define the phase space variables: rapidity $`y`$, transverse momentum $`p_\mathrm{t}`$ and azimuthal angle $`\phi `$. Different frames could be used to define the azimuthal angle $`\phi `$. The first one is to choose the minor axis of thrust (or sphericity) analysis as the $`x`$ axis, and use it as the starting point for counting the azimuthal angle $`\phi `$. The second one is to put the $`z`$ axis still on the major thrust axis, but turn the coordinate system around it and let the new $`x`$ axis lie on the $`z_0`$-$`z`$ planLLL , where $`x_0,y_0,z_0`$ denote the axes of the lab system and $`x,y,z`$ those of the turned system, as shown in Fig.1. In the following this frame will be referred to as the “rotated frame”. The third method is to rotate the frame in each event for a random angle around the $`z`$ axis CGDatong . This is called the random frame.

All these methods have been used in the literature for studying the anomalous scaling of factorial moments, making the results hard to be compared. In the present paper we will examine these methods in some detail and discuss their applicability in physical problems. We will take e+e-collisions at the Z0 energy $`\sqrt{s}=91.2`$ GeV as example and use JETSET7.4 Jetset  Monte Carlo code to generate 500 000 multihadron events as the event sample.

## II Average distributions of phase space variables

In Fig.2 ($`a`$), ($`b`$) and ($`d`$) are shown the average distributions of $`y`$, $`p_\mathrm{t}`$ and $`\phi `$, where the rapidity is defined as $`y=0.5\mathrm{ln}[(E+p_z)/(Ep_z)]`$ with $`z`$ along the thrust axis; the azimuthal angle $`\phi `$ is defined in the plane perpendicular to the thrust axis, calculated with respect to the minor axis.

It can be seen from the figures that all the distributions are unflat. Especially, the distribution of $`p_\mathrm{t}`$ is exponential and is highly concentrated in low $`p_\mathrm{t}`$. A simple variable transformation $`p_\mathrm{t}\mathrm{ln}p_\mathrm{t}`$ OPAL  can make it a litter flatter as shown in Fig.2 ($`c`$).

## III The correction factor and cumulant variables

The unflat average distribution will cause additional variation of factorial moments with the diminishing of phase space scale and make the scaling property of factorial moments unequal to that of the probability moments even for the VFM. This effect has to be corrected.

Fialkowski proposed a factor Fialkowski , cf. Eq.(5), to correct for this effect. This method works good when the distribution of the phase space variable is not far from flat, e.g. the distribution of rapidity $`y`$ in a restricted central region $`|y|<Y_\mathrm{c}`$ with $`Y_\mathrm{c}=2`$ as shown in Fig.2 ($`a`$), and is not good when the distribution is far from flat. This is especially the case for the distribution of $`p_\mathrm{t}`$, cf. Fig.2 ($`b`$). Therefore, people sometimes transform $`p_\mathrm{t}\mathrm{ln}p_\mathrm{t}`$ first OPAL  and calculate $`F_q(\mathrm{ln}p_\mathrm{t})`$ instead of $`F_q(p_\mathrm{t})`$, and then correct the result by the factor $`R_q`$ given in Eq.(5). Note that the highly asymetric region $`0.1p_\mathrm{t}2`$ is transformed to the region $`2.3\mathrm{ln}p_\mathrm{t}0.69`$, which is distributed symmetrically around the pick of distribution, cf. the two arrows in Fig’s.2 ($`b`$) and ($`c`$).

Since a transformation to a flatter distribution is necessary before calculating factorial moments, it is evident that the best way is to transform all the phase space variables to a flat distribution first. This could be established through the transformation to cumulant variables Cummulant , cf. Eq.(6). The corresponding distributions are shown in Fig.3.

## IV Vertical and horizontal factorial moments

Now let us turn to the comparison of vertical and horizontal factorial moments (VFM and HFM).

As noticed in the Introduction, our aim is to explore the anomalous scaling of probability moments as shown in Eq.(2) but the equivalence of factorial and probability moments has been proved only for the VFM. So, the HFM is appropriate only when it is equal to VFM.

In Figures 4 and 5 $`(a),(b),(d),(e),(f)`$ are shown the 1-D and 3-D VFM (solid circles) and HFM (solid curves) for the moment orders $`q=2,3,4,5`$ calculated using the correction factor, Eq.(5), method and the cumulant variables Eq.(6), respectively.

It can be seen from the figures that the VFM and HFM are equal only when using together with the cumulant variables and are unequal, especially for 1-D $`F_q(\phi )`$ and $`F_q(p_\mathrm{t})`$, when using the correction factor method.

The results from VFM + Correction Factor method are about the same as that from VFM + Cummulant Variables for 1-D $`F_q(y)`$ and 3-D $`F_q(y,p_\mathrm{t},\phi )`$ but are not the case for 1-D $`F_q(\phi )`$ and $`F_q(p_\mathrm{t})`$, cf. Fig’s. 4 and 5.

## V The influence of phase space reduction

In hadron-hadron and nucleus-nucleus collisions the central rapidity regions with the rapidity $`y`$ restricted to $`|y|Y_\mathrm{c}`$ for some value of $`Y_\mathrm{c}`$ is commonly used. This is physically meaningful, because in these collisions the final state particles are mainly produced in the central region, while the particles in the regions $`|y|>Y_\mathrm{c}`$ are mainly come from the fragmentation of incident particles (leading particle effect).

On the contrary, in e+e-collisions, the multihadron final state is produced from a point source — virtual photon or Z0, and no leading particle effect is present. The rapidity $`y`$ is usually defined with respect to the thrust or sphericity axis. In this case, to carry on the study in a “central rapidity region” $`|y|Y_\mathrm{c}`$ is physically doubtful. This is especially evident when the collision energy is so high that 3, 4 or even more jets can be produced. In a 2-jet event the restriction $`|y|Y_\mathrm{c}`$ will cut out the most energetic particles from the two opposite jets symmetrically, but in a 3-jet event the same cut will cut out the most energetic particles only from one jet while retain almost all the particles in the other two jets. This asymmetric cut will results in unexpected phenomena, and the physical meaning may be difficult to interpret. Therefore, the study of multiparticle dynamics in e+e-collisions is better to be carried out in the full phase space. However, the central rapidity region is sometimes also used in the literature for the study of e+e-collisions OPAL . Therefore, to investigate the influence of rapidity cut is worthwhile.

In Fig.6 are shown the results of VFM and HFM in nearly the full phase space — $`5<y<5,0<\phi <2\pi ,0.1<p_\mathrm{t}<3`$ GeV, to be compared with the results shown in Fig.5 for a restricted phase space — $`2<y<2,0<\phi <2\pi ,0.1<p_\mathrm{t}<2`$ GeV.

It can be seen from Fig.6 that the first point in 3-D $`F_q(y,p_\mathrm{t},\phi )`$ and 1-D $`F_q(y)`$ and the first 3 points in 1-D $`F_q(\phi )`$ do not lie on a scaling curve together with the other points. This is due to the momentum conservation effect MMCN . In the anomalous-scaling study, these points should be omitted.

The momentum conservation effect will also be reduced in a restricted phase space region, which was first pointed out in Ref.MMCN  and has been proposed as a second method for eliminating the influence of this effect. This explains the reason why the first points in Fig’s.5 lie on the scaling curves.

## VI The fluctuations in azimuthal angle

The fluctuations in azimuthal angle are worthwhile special investigation. It is commonly expected that there should have cylindrical symmetry around the $`z`$ axis. If that is the case, then the fluctuations should have no correlation with the $`x`$ axis chosen for counting the azimuthal angle $`\phi `$. In Fig’s. 4, 5, 6 ($`b`$) are shown the $`F_q(\phi )`$ with $`x`$ axis along the minor of thrust analysis, while in the corresponding figures ($`c`$) are shown the results after rotating the $`x`$ axis to let it lie on the $`z`$-$`z_0`$ plane LLL , cf. Fig.1. It can be seen that in the rotated frame, Fig’s. 4, 5, 6 ($`c`$), the $`F_q(\phi )`$ increases much faster as the diminishing of phase space scale than that in the thrust-minor frame, Fig’s. 4, 5, 6 ($`b`$). As discussed in Ref. LLL  this is because the thrust-minor axis is basically determined by the first hard gluon emission and taking this axis as $`x`$ axis to count the azimuthal angle $`\phi `$ will highly reduce the fluctuation of the direction of first hard gluon emission. After rotation, the correlation between $`x`$ axis and the direction of first hard gluon emission is relaxed and the full dynamical fluctuations are exhibited.

We could also rotate the $`x`$ axis around $`z`$ for a random angle in each event CGDatong . The resulting $`F_q(\phi )`$ turn out to be the same as those in the rotated frame with $`x`$ on the $`z`$-$`z_0`$ plane shown in Fig’s.4, 5, 6 ($`c`$). This confirms the cylindrical symmetry of the fluctuation in $`\phi `$ after the correlation with the thrust-minor is relaxed.

## VII Conclusions

The following conclusions could be drawn from the above investigation:

1) The horizontal factorial moments (HFM) are equivalent to the vertical ones (VFM) only after the cumulant-variable transformation. Therefore, in the study of non-linear phenomena (intermittency or fractal) in high energy collisions the HFM could be used only in combination with the cumulant variables. On the other hand, the HFM is in its own right useful in single-event anaylsis. It can be seen from Eq.(4) that HFM is the average of the so called single-event factorial moments $`f_q^{(\mathrm{e})}`$ HWA

$`F_q^{(\mathrm{H})}(M)`$ $`=`$ $`f_q^{(\mathrm{e})}/f_1^{(\mathrm{e})}^q.`$ (7)
$`f_q^{(\mathrm{e})}(M)`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{m=1}{\overset{M}{}}}n_m^{(\mathrm{e})}(n_m^{(\mathrm{e})}1)\mathrm{}(n_m^{(\mathrm{e})}q+1)`$ (8)
where $`n_m^{(\mathrm{e})}`$ is the number of particles of a single event in the $`m`$th sub-cell. The fluctuation of the single-event factorial moments around its average (HFM) is a characteristic of single-event fluctuations HWA .

2) The scaling properties of factorial moments in transverse directions ($`\phi ,p_\mathrm{t}`$) are very sensitive to the correction method used. They are unstable when using the correction factor method, Eq.(5), even after the transformation $`p_\mathrm{t}\mathrm{ln}p_\mathrm{t}`$ has been made. Using this method, the VFM in $`\phi `$ and the HFM in $`p_\mathrm{t}`$ fall down instead of increase with the diminishing of phase space scale, while at the same time the HFM in $`\phi `$ and the VFM in $`p_\mathrm{t}`$ do increase with the diminishing of scale, cf. Fig’s. 4 ($`b`$) and ($`d`$).

3) In the full phase space, the first few points of factorial moments do not lie on the scaling curve with the other points, due to the momentum conservation effect. This effect can be eliminated either through neglecting these points or through a cut in phase space.

4) The thrust (or sphericity) major-minor frame is inappropriate for the study of the scaling property of the azimuthal angle $`\phi `$, because this frame is strongly correlated with the direction of first hard gluon emission. Rotate the $`x`$ axis to let it lie on the $`z`$-$`z_0`$ plane or rotate it randomly for each event can relax this correlation and exhibit the full dynamical fluctuations in $`\phi `$.

Therefore, the cumulant variables together with a frame rotated around the thrust (or sphericity) axis is the best for the investigation of the nonlinear phenomena (anomalous scaling of probability moments) in e+e-collisions. The VFM and HFM are equivalent in this case.

###### Acknowledgements.

The authors thank Wolfram Kittel and Hu Yuan for helpful discussions.


text:
# Unpredictability of wave function’s evolution in nonintegrable quantum systems

Abstract. It is shown that evolution of wave functions in nonintegrable quantum systems is unpredictable for a long time $`T`$ because of rapid growth of number of elementary computational operations $`𝒪(T)T^\alpha `$. On the other hand, the evolution of wave functions in integrable systems can be predicted by the fast algorithms $`𝒪(T)(log_2T)^\beta `$ for logarithmically short time and thus there is an algorithmic ”compressibility” of their dynamics. The difference between integrable and nonintegrable systems in our approach looks identically for classical and quantum systems. Therefore the minimal number of bit operations $`𝒪(T)`$ needed to predict a state of system for time interval $`T`$ can be used as universal sign of chaos.

Chaos as universal phenomenon exists in various systems (for example, in human society) and from this point of view the general approach to chaos should not be based on particular properties of a system. It is well known  that motion of nonintegrable classical systems has all attributes of chaos: complexity, unpredictability and randomness. However in nonintegrable quantum systems the apparent signs of chaos seems to be absent  and the main direction of studies in the field of ”quantum chaos” is semiclassical analysis of various quantum ”signatures” of classical chaos . Nevertheless, taking into account the fundamental correspondence principle it is reasonble to suppose that evolution of nonintegrable quantum systems should be also unpredictable. As a basis tool to analyze predictability of dynamics we use a number of elementary computational operations $`𝒪(T)`$ needed to determine a state of the system for time interval $`T`$.

Let us consider firstly how a number of elementary computational operations needed for prediction of system’s evolution depends on time $`T`$ and accuracy $`\mathrm{\Delta }`$ in classical mechanics. It is well known  that a distance $`\delta x(t)`$ between initially close phase space points grows with time as

$$\delta x(t)=\delta x(0)f(t),$$
(1)
where $`f(t)e^{\lambda t}`$ for nonintegrable systems and $`f(t)t`$ for integrable ones. The exponential growth of inevitable computational errors leads to unpredictability of long time evolution in nonintegrable classical systems. To predict a state (trajectory) of the system by the moment $`T`$ with accuracy $`\mathrm{\Delta }`$ we must make computations with accuracy at least

$$\delta =\frac{\mathrm{\Delta }}{f(T)}$$

and length of mantissa

$$nlog_2\delta =log_2f(T)log_2\mathrm{\Delta }.$$

It is reasonable to suppose that an algorithm which generates a trajectory has a power dependence of number of elementary bit operations $`𝒪(n)`$ on the length of mantissa (for example, $`𝒪(n)n`$ for addition and $`𝒪(n)n^2`$ for multiplication):

$$𝒪(T)n^\alpha (log_2f(T)log_2\mathrm{\Delta })^\alpha ,$$
(2)
where $`\alpha >1`$ is some number. Inserting now $`f(t)`$ in (2) we obtain for nonintegrable systems ($`\alpha >1`$)

$$𝒪(T)T^\alpha ,$$
(3)
while for integrable ones

$$𝒪(T)(log_2T)^\alpha .$$
(4)
In the first case a number of operations and hence a computational time needed for prediction grows as a power of time $`T`$. Therefore the prediction for a long time is impossible and evolution is not algorithmically ”compressible”. In regular systems there is a completely different situation because computational time grows logaritmically with time $`T`$ and prediction is always possible (”compressibility” of evolution). It is interesting to note that algorithmic ”uncompressibility” of some numerical sequence is the main sign of its randomness: such sequence we can only observe.

Let us consider now a time evolution of some state $`\psi (q,t)`$ in a stationary quantum system $`H(q,p)`$

$$i\mathrm{}_t\psi (q,t)=H\psi (q,t),$$
(5)
$$\psi (q,0)=\psi _0(q),$$

Any solution of (5) at $`t=T`$ can be represented as

$$\psi (q,T)=\underset{\mu }{}c_\mu e^{iE_\mu T/\mathrm{}}\varphi _\mu (q),$$
(6)
where $`\left\{E_\mu \right\}`$ and $`\left\{\varphi _\mu \right\}`$ – are exact eigenenergies and eigenfunctions of Hamiltonian $`H`$ and $`c_\mu =\varphi _\mu |\psi _0`$. So if quantum chaos exists it hides in eq. (6) which is essence of quantum unitary evolution. In analogy with classical case it would be reasonable to investigate a time dependence of two initially close wave functions $`\psi _1`$ and $`\psi _2`$. It is natural to characterize a difference between wave functions by the quantity

$$\delta \psi =\psi _2\psi _1=\sqrt{2(1Re\psi _1|\psi _2)}.$$
(7)
If at $`t=0`$ we have $`\delta \psi (0)<\epsilon `$ then at $`t=T`$ the difference (7) is also small $`\delta \psi (T)<\epsilon `$ because from unitarity it follows that $`\psi _1|\psi _2`$ does not vary with time. Therefore in quantum systems there is no divergence of initially close states which is the main sign of chaos in classical systems.

How then unpredictability of evolution can arise in quantum systems? To predict a state of the system $`\psi (q,T)`$ at moment $`t=T`$ it is necessary to know all eigenvalues and eigenfunctions in (6) exactly. Really we use in (6) approximate eigenvalues and eigenfunctions and therefore by some moment exact wave function and our prediction will be completely different. It is obvious that to increase a maximal time of accurate prediction we must improve accuracy of eigenvalues. The spectrum of integrable systems is easy computable and computational resources grow slowly with increasing of accuracy while for nonintegrable systems considerable efforts are needed to improve the spectrum and thus to increase the time of correct prediction. Therefore we can assume that quantum integrable systems differ significantly from nonintegrable ones in a sense of number of bit operations $`𝒪(T)`$ needed to predict a future. In the remaining part of this Letter we consider some details of the suggested mechanism of unpredictability and hence of chaos in quantum unitary evolution.

Let us assume that we have found approximate eigenenergies $`\{\stackrel{~}{E}_\mu \}`$ and eigenfunctions $`\{\stackrel{~}{\varphi }_\mu \}`$. Then we can predict by means (6) that at moment $`T`$ wave function will be

$$\stackrel{~}{\psi }(q,T)=\underset{\mu }{}\stackrel{~}{c}_\mu e^{i\stackrel{~}{E}_\mu T/\mathrm{}}\stackrel{~}{\varphi }_\mu (q),$$
(8)
where $`\stackrel{~}{c}_\mu =\stackrel{~}{\varphi }_\mu |\psi _0`$. Our main goal now is to find a time dependence of the difference between exact $`\psi (q,T)`$ and approximate $`\stackrel{~}{\psi }(q,T)`$ wave functions

$$\delta \psi (T)=\stackrel{~}{\psi }(T)\psi (T)=\sqrt{2(1Re(\psi (T)|\stackrel{~}{\psi }(T)))}.$$
(9)
Making use of (6) and (8) we have:

$$\psi (T)|\stackrel{~}{\psi }(T)=\underset{\mu ,\nu }{}c_\mu ^{}\stackrel{~}{c}_\nu e^{i(\stackrel{~}{E}_\nu E_\mu )T/\mathrm{}}\varphi _\mu |\stackrel{~}{\varphi }_\nu .$$
(10)
Now we need some assumptions about accuracy of spectrum obtained. If a maximal error of eigenfunctions $`\stackrel{~}{\varphi }_\mu `$ does not exceed $`\epsilon `$

$$\delta \varphi _\mu =\stackrel{~}{\varphi }_\mu \varphi _\mu ,$$

$$\delta \varphi _\mu <\epsilon .$$

then the error of $`E_\mu `$ is

$$\delta E_\mu =\stackrel{~}{E}_\mu E_\mu =2E_\mu Re(\varphi _\mu |\delta \varphi _\mu )+\delta \varphi _\mu |H|\delta \varphi _\mu ,$$
(11)
$$|\delta E_\mu |<2|E_\mu |\epsilon +\epsilon ^2H,$$

and error of $`c_\mu `$ is

$$\delta c_\mu =\stackrel{~}{c}_\mu c_\mu =\delta \varphi _\mu |\psi _0,$$

$$|\delta c_\mu |<\epsilon .$$

Scalar product $`\varphi _\mu |\stackrel{~}{\varphi }_\nu `$ we can write as

$$\varphi _\mu |\stackrel{~}{\varphi }_\nu =\delta _{\mu \nu }+\varphi _\mu |\delta \varphi _\nu =\delta _{\mu \nu }+R_{\mu \nu },$$

$$|R_{\mu \nu }|<\epsilon .$$

Now (10) becomes

$$\psi (T)|\stackrel{~}{\psi }(T)=\underset{\mu }{}c_\mu ^{}\stackrel{~}{c}_\mu (1+R_{\mu \mu })e^{i\delta E_\mu T/\mathrm{}}+\underset{\mu \nu }{}c_\mu ^{}\stackrel{~}{c}_\nu R_{\mu \nu }e^{i(E_\nu E_\mu +\delta E_\nu )T/\mathrm{}}.$$
(12)
It is clear that first sum carries the main contribution in $`\psi (T)|\stackrel{~}{\psi }(T)`$ because second sum has a small multiplier $`R_{\mu \nu }`$. To understand the main features of (10) it is sufficient to investigate time dependence of

$$\psi (T)|\stackrel{~}{\psi }(T)\underset{\mu }{}|c_\mu |^2e^{i\delta E_\mu T/\mathrm{}}.$$
(13)
At $`T=0`$ a value of $`\psi (0)|\stackrel{~}{\psi }(0)`$ equals $`1`$, then $`\psi (T)|\stackrel{~}{\psi }(T)`$ decreases and at some moment $`T_p`$ reaches values nearly zero. After this $`\psi (T)|\stackrel{~}{\psi }(T)`$ starts to oscillate about zero with some small amplitude. For times $`T>T_p`$ the value of $`\delta \psi (T)`$ oscillates about mean value $`\sqrt{2}`$ which corresponds to completely independent states. This typical behaviour of $`\psi (T)|\stackrel{~}{\psi }(T)`$ and $`\delta \psi (T)`$ is shown in Fig. 1. This means that for $`T>T_p`$ wave function $`\stackrel{~}{\psi }(T)`$ completely differs from exact function $`\psi (T)`$ and our further prediction is impossible.

We want now to find how $`T_p`$ and average amplitude of small oscillations $`\overline{A}`$ depend on errors $`\delta E_\mu `$ and coefficients $`c_\mu `$. Our intuition and numerical experience say that $`T_p`$ depends on dispersion of errors $`dE`$ while $`\overline{A}`$ is a function of number $`dim`$ of the most important addendums in (13). It is obvious that $`\psi (T)|\stackrel{~}{\psi }(T)`$ becomes nearly zero when arguments of exponents in (13) fill interval $`[0,2\pi ]`$. From this condition we have $`T_p`$:

$$T_p\frac{2\pi }{2dE/\mathrm{}}=\frac{\pi \mathrm{}}{dE}.$$
(14)
To estimate $`\overline{A}`$ let us make further simplification that all $`c_\mu `$ in (13) are equal $`1/\sqrt{dim}`$. In such a case (13) becomes

$$P(T)=Re(\psi (T)|\stackrel{~}{\psi }(T))=\frac{1}{dim}\underset{\mu =1}{\overset{dim}{}}\mathrm{cos}(\delta E_\mu T/\mathrm{})$$
(15)
and $`P^2(T)`$ is:

$$P^2(T)=\frac{1}{dim^2}\underset{\mu ,\nu }{}\mathrm{cos}(\delta E_\mu T/\mathrm{})\mathrm{cos}(\delta E_\nu T/\mathrm{})=$$

$$\frac{1}{dim^2}\underset{\mu }{}\mathrm{cos}^2(\delta E_\mu T/\mathrm{})+\frac{1}{dim^2}\underset{\mu \nu }{}\mathrm{cos}(\delta E_\mu T/\mathrm{})\mathrm{cos}(\delta E_\nu T/\mathrm{}).$$

It is easy to see that if we make time averaging of $`P^2(T)`$ then the first sum equals $`dim/2`$ while the second sum is zero because arguments of cosines are almost independent for different $`E_\mu `$ and $`E_\nu `$. In the result average amplitude of oscillations of $`Re(\psi (T)|\stackrel{~}{\psi }(T))`$ can be estimated as

$$\overline{A}\sqrt{\overline{P^2(T)}}=\frac{1}{\sqrt{2dim}}.$$
(16)
We tested our theoretical estimates (14) and (16) for wide range of parameters and found good agreement with numerical computations. So we have obtained a quite reasonable result that uncertainty in our knowledge of spectrum $`E\pm dE`$ leads to limitation in prediction time of $`\psi (T)`$ evolution $`T_p\mathrm{}/dE`$.

The next step consists of estimating of number of operations needed to obtain the spectrum with some accuracy $`dE`$. If the spectrum of a system can be computed by some formulae or by effective algorithms (as in the case of integrable systems), the number of operations grows slowly with decreasing of error $`dE`$. To explain this let us assume for simplicity that to compute the spectrum we must perform some series of successive actions which do not depend on accuracy required (for example $`E_N=\mathrm{}\omega (N+1/2)`$). Obviously the necessary number of bit operations $`𝒪(dE)`$ depends on mantissa length $`nlog_2dE`$ only

$$𝒪(dE)n^\gamma (log_2dE)^\gamma ,$$
(17)
where $`\gamma >1`$ – some coefficient characterizing the algorithm of computations. So we have obtained that in integrable systems a number of bit operations depends logarithmically on accuracy $`dE`$.

Fortunately the variety of systems in nature is not exhausted by easy computable and predictable systems. Description of any nonseparable system with strong coupling and especially of its excited states is nontrivial problem even for two degrees of freedom. The only reliable way to calculate such systems is to use various variation methods, i.e. to minimize the energy functional over some space of trial functions. Let us consider how the number of operations depends on spectrum accuracy for variational method of Ritz. It is well known  that Ritz’s method has a power convergence, i.e. approximate eigenfunctions $`\stackrel{~}{\varphi }_\mu `$ tend to exact ones $`\varphi _\mu `$ as

$$\stackrel{~}{\varphi }_\mu \varphi _\mu D^\alpha ,$$
(18)
where $`D`$ is dimension of space of trial functions and $`\alpha 1`$ – some coefficient. Obviously (see (11)) the rate of convergence of eigenvalues does not exceed the power law:

$$|\stackrel{~}{E}_\mu E_\mu |D^\alpha .$$
(19)
To make use of Ritz’s method it is necessary to calculate firstly a Hamiltonian matrix with number of elements $`D^2`$, and then to find its eigenvalues. The total number of operations to do this can be estimated as $`𝒪(D)D^\beta ,\beta >2`$. To achieve an accuracy $`\epsilon `$ the basis dimension $`D\epsilon ^{1/\alpha }`$ is required (see (19)) and hence the total number of operations is $`𝒪(\epsilon )\epsilon ^{\beta /\alpha }`$. Keeping in mind the relation between time of prediction and accuracy of eigenvalues (14) we obtain the number of bit operations which are needed to predict the wave function’s evolution over time interval $`T`$:

$$𝒪(T)T^{\beta /\alpha }.$$
(20)
For regular systems (17) the number of operations is equal to

$$𝒪(T)(log_2T)^\gamma .$$
(21)
It should be noted that our estimate (20) is valid for any algorithm with power convergence.

So we see that difference between integrable and nonintegrable systems in our approach looks identically in classical and quantum systems. The evolution of integrable systems can be predicted by fast algorithms for logarithmically short time (”compressibility” of evolution) while in nonintegrable systems such ”compressibility” is absent. The reason of unpredictability is universal in classical and quantum systems — rapid growth of number of elementary computational operations needed for prediction, but the mechanism of chaos is completely different.


text:
# 1 Introduction

## 1 Introduction

An important part of the humoral immune response is the germinal center (GC) reaction. GCs are responsible for an optimization process of antibodies with respect to a specific antigen. This process is called affinity maturation: During the GC reaction new plasma cells are generated which secrete antibodies of considerably higher affinity to the antigen.

The GC reaction is initiated by antigen-activated B-cells that migrate into the follicle system. Here, they start to proliferate in the environment of follicular dendritic cells (FDCs). The initiation is believed to be of oligoclonal character, i.e. the number of seeder B-cells is small and of the order of three . After three days of fast monoclonal expansion – the total number of proliferating B-cells (centroblasts) reaches about $`12000`$ – a phase of somatic hypermutation is started . The diversity of encoded antibodies is enhanced in this way. The centroblasts differentiate into antibody-presenting centrocytes  and an apoptotic process is initiated. However, they have the possibility to get into interaction with the antigen-presenting FDCs and with T-helper cells. It is believed that this interaction depends on the affinity of antibody and antigen, and that these centrocytes which successfully bind the antigen are rescued from apoptosis . This provides a more-step selection process  of these B-cells with high affinity to the antigen. Positively selected B-cells further differentiate into plasma- and memory-cells (shortly denoted as output cells). In this way the answer of the immune system is optimized with respect to the antigen.

GCs develop a very specific morphology . The proliferating and mutating centroblasts are collected in the dark zone while centrocytes and FDCs build up the light zone. It is unclear how long the dark zone remains present during a GC reaction. The total duration of a GC reaction is about $`21`$ days . Dark zones have been observed to appear at day $`4`$ and to vanish at day $`8`$ . However, there also exists evidence for dark zones of longer duration .

In the present article a possible correlation of the GC morphology and affinity maturation is investigated. in the framework of a mathematical model  which considers the spatial distribution of cells in the GC for the first time. The assumptions made in this model will be compared to other model architecture. The postulates of the model are summarized and shortly described in Sec. 2. Note, that all parameters of the model are quantitatively determined in narrow connection to experiment. An analysis of the robustness of the results has been performed previously . Here, we aim to report one major outcome of the model (Sec. 3). An interpretation and possible implications of the results as well as methodological aspects are discussed in Sec. 4.

## 2 The scheme of the model

A short description of a previously introduced mathematical model for the morphological organisation and cell dynamics of the GC  is provided in this section. The GC is simulated on an equidistant two-dimensional lattice with lattice constant $`10\mu m`$. This corresponds to the average cell diameter of B-cells in GCs. The radius of the lattice is $`220\mu m`$, corresponding to a typical radius of a GC. Each lattice point can be occupied by exactly one centroblast, centrocyte, or output cell. All cells actively and isotropically move on the lattice. The diffusion constants are adapted corresponding to the different diameters of centroblasts and centrocytes . FDCs are represented by a soma at one lattice point and four (in 2 dimensions) dendritic arms of $`30\mu m`$ length. FDCs are assumed to be immobile reflecting a rather stable FDC network observed in experiment.

It has been previously shown that the development of dark zones can be explained on the basis of non-local cell-cell interactions . Such an interaction may be provided by a chemotaxis gradient which acts on the motility of centrocytes and stems from FDC and/or T-cells and/or naive B-cells in the mantle zone . Another possibility (which will be used in this work) is a diffusing signal molecule which is produced by FDCs and bound by centroblasts . Note, that this implies a separation of signals acting on proliferation and differentiation of centroblasts, as has been proposed in corresponding experiments . The signal molecules are clustered in quanta that diffuse on the lattice according to a classical diffusion equation. The diffusion is not influenced by the presence of cells at the same lattice point. One quantum corresponds to the signal concentration that is necessary to initiate the centroblast differentiation process into centrocytes. Using this non-local concept an intermediate dark zone is produced . The duration of the dark zone basically is dependent on the amount of secreted signal molecules and its diffusion constant. The ratio of centroblast differentiation and proliferation rates changes the duration of the dark zone as well. However, this ratio also has influence on the total life time of the GC as a whole and, therefore, is determined independently.

The affinity of the encoded antibodies to the antigen is formulated with the well known shape space concept . Each type of antibody is represented on a four-dimensional lattice which is ordered in such a way, that neighboring points have similar affinity to the antigen. A hypermutation is represented by a jump to a neighbor point. The affinity between the antibodies on a centrocyte and the antigen on an FDC is modeled by a gaussian affinity weight function centered at the optimal antibody type .

The dynamical properties of the different cell types on the lattice are summarized in the following. The parameter values have been quantitatively determined using experimental constraints. In many cases, indeed, the parameters were directly accessible in experiments. Others had to be determined indirectly using experimental observations of the general GC properties (for more details we refer to ):
Centroblasts proliferate with a (constant) rate of $`1/6hr`$ . At each division a somatic hypermutation occurs with the probability $`1/2`$ . They differentiate in dependence on a differentiation signal that is secreted by the FDCs, and diffuses over the lattice. The differentiation process is activated when a centroblast meets a threshold quantum of differentiation signal at the same lattice point. Activated centroblasts differentiate with a rate of $`1/6hr`$ into centrocytes . A finite life time of centroblasts is not imposed. However, the effective life time is shorter than $`1`$ day due to centroblast differentiation.
Centrocytes die with a rate of $`1/7hr`$ . They bind to the FDCs according to the affinity to the antigen according to an affinity weight function. Bound centrocytes remain bound for $`2hr`$ . They are thought to be rescued from apoptosis during this time. Positively selected centrocytes further differentiate with a rate of $`1/7hr`$ into either re-proliferating centroblasts (with probability $`80\%`$) , or into output cells (with probability $`20\%`$). The differentiation into output cells is delayed by $`48hr`$ with respect to the starting time of hypermutation , i.e. it starts at day 6 of the GC reaction.
Output cells leave the GC by diffusion and do not further interact with other cells in the GC.
Dead cells are eliminated from the GC.

The simulations are started with $`3`$ randomly distributed seeder B-cells and $`20`$ FDCs. The fact that centroblasts proliferate at least in parts outside the FDC network during their proliferation phase turned out to be a necessary requirement for the development of dark zones . This is ensured by a random distribution of the FDCs on $`70\%`$ of the (maximum) GC volume. The seeder cells are of low but non-vanishing affinity to the antigen. They can reach the optimal antibody-type with $`5`$ to $`10`$ mutations . The simulations are insensitive to a change of the time-step-width which is $`0.004hr`$ for the presented results. In a stochastic model the outcome of the simulation depends on the used generator of random numbers and on its initialization. Therefore, the results are given with a standard deviation corresponding to this uncertainty.

## 3 Optimisation of affinity maturation

At first some basic properties of the GC reaction are reported as generated by the model simulation. Assuming the already introduced centroblast differentiation signal (secreted by FDCs and bound by centroblasts) a dark zone develops. It appears at day $`4`$ of the reaction, and remains stable for a duration that depends on the production rate of the signal molecules. In order to simulate different durations of dark zones, the signal production rate is varied. For each production rate the differentiation rate of centroblast is adjusted correspondingly (within physiological constraints) in order to ensure a comparable final state of the GC after 21 days of the reaction. The development of the dark zone as well as its depletion are not principally affected by the variation of other parameters within their physiological constraints. For example a smaller proliferation rate basically scales the whole GC reaction without changing the general behavior.

The time course of the total GC volume is in accordance with experimental observations  provided that the dark zone vanishes between day $`6.5`$ and $`10`$ of the GC reaction (i.e. that the signal production rate is chosen correspondingly): After an exponential increase of the total cell population, a maximum is reached after $`4`$ days of the reaction. The total cell population then is diminished steadily until the end of the reaction after $`21`$ days. At this time only about $`50`$ proliferating B-cells remain in the GC . Taking these results together, the general GC morphology is well described by the model results.

The average affinity of B-cells is enhanced in four phases during the whole GC reaction. This is best illustrated by a typical example with a dark zone present until day $`8.3`$. In Fig. 1

the time course of the fraction of high affinity centroblasts and centrocytes is shown . High affinity denotes those B-cells which bind the antigen with a probability of more than $`30\%`$. As can be seen, such a cell does not exist at the beginning of the GC reaction. They develop after the start of somatic hypermutations after $`3`$ days. Still one observes a short delay because some mutations have to occur before the first cells appear that have an above threshold affinity to the antigen. Then the relative number of high affinity cells steadily grows. In accordance with experiment, good cells already dominate around day $`10`$ of the GC reaction . One observes an intermediate steep increase, that approximately starts when the dark zone is depleted. This correlation does not seem to be obvious in view of this single example. However, it has been observed in all simulations, especially considering different dark zone durations. In the late phase of the GC reaction the curve reaches a plateau on a high level, and the large majority of B-cells are high affinity cells.

In the same figure the fraction of high affinity output cells is shown. This curve does not show the value at each moment of the GC reaction but the sum of all output cells that has been produced until time $`t`$. This value is a measure for the total quality of the produced output cells. After day $`5`$ the quality of the output cells is steadily increased during the GC reaction. The part of high affinity output cells reaches $`77\%`$ in the present example. The average over all simulation with dark zones that vanish between day $`8.5`$ and $`10.5`$ of the reaction is $`75.2\%\pm 4.2\%`$ (the error is one standard deviation). In view of the fact that at the beginning of the GC reaction no high affinity B-cell existed at all, this affinity enhancement is remarkable.

The main task of the present article is to analyze a possible correlation of the duration of the dark zone and the achieved affinity maturation, i.e. the total output quality at the end of the reaction. A statistical analysis has to be based on comparable GC simulation. By changing the production rate of the centroblast differentiation signal molecule not only the duration of the dark zone is varied. Also the total duration of the reaction is changed. As stated at the beginning of this section, the centroblast differentiation rate is adjusted correspondingly so that the final numbers of B-cells $`N`$ after $`21`$ days of the reaction are of the same range. This is especially important as the affinity maturation process depends on the value of $`N`$ . The average number of B-cells for all $`226`$ simulations is $`\overline{N}=49\pm \mathrm{\hspace{0.17em}23}`$. The error denotes one standard deviation. Only those simulations are taken into account that generate a final number of B-cells $`N`$ within one standard deviation of this average value.

The output quality is plotted against the duration of the dark zone (see Fig. 2).

The quality of the produced output cells becomes optimized for dark zones that vanish between day $`7`$ and $`11`$ of the GC reaction. For shorter and longer dark zones, the resulting quality of the output cells is reduced. The same holds true for the total number of output cells (data not shown). It is worth pointing out, that intermediate dark zones vanishing around day 9 optimize affinity maturation on the level of quantity and quality of output cells.

## 4 Discussion

We used a previously developed model to elucidate a possible correlation between the duration of the dark zone and the efficiency of the affinity maturation process. In a first step the simulated GCs were compared to GCs observed in experiment. The essential properties of real GC reactions were correctly reproduced. This includes the appearance of the dark zone, the time course of the total volume, as well as the reached affinity maturation.

It has been previously found  that dark zones do not appear in simulations that are based on local cell interactions only. Therefore, it was necessary to introduce an non-local cell interaction into the model in order to understand the intermediate appearance of dark zones as observed in experiment .

A centroblast differentiation signal molecule has been assumed that is secreted by FDCs and consumed by the centroblasts. It has been previously shown that such a signal is a promising candidate to explain appearance and depletion of dark zones in GCs . However, the signal molecule has to be understood as a hypothesis, because it is not known how the centroblast differentiation process is initiated in real GC reaction. There are experimental hints that centroblast, indeed, differentiate due to an interaction with FDCs or T-helper cells .

Most parameters in the model are strongly constrained by experimental data and consequently the variation of parameters is restricted. However, the production of the hypothesized differentiation signal does not underly such a restriction. It turned out that the duration of dark zones strongly depends on the production rate of the differentiation signal and that other model parameters are less important. This situation opens the possibility to vary the signal production rate and in this way to test the affinity maturation process of GCs with dark zones of different duration. Note, that a corresponding test is not possible in experiments.

A statistical analysis of simulations with different durations of dark zones leads to the conclusion that the quality of the output cells averaged over all produced output cells during the whole GC reaction is optimized for intermediately appearing dark zones that vanish between day $`7`$ and $`11`$ of the GC reaction. Note that also the quantity of produced output cells is optimized for these dark zone durations. In addition, the time course of the total GC volume is in agreement with observed time courses for the same dark zone durations .

This result suggests a relation of two, at first sight, very different categories: The GC morphology and affinity maturation of B-cells. The morphology of the GC is basically determined by a non-local interaction with other cell-types. However, the function of the specific GC morphology is not restricted to a spatial arrangement of interacting cells. One may suspect that specific cell arrangements are advantageous for the success of the GC reaction. The duration of the dark zone does not only determine the time course of the cell distribution in GCs but also to be a critical parameter for the quality and amount of output of the GC reaction. This leads to the question of how the spatial distribution of centroblasts and centrocytes influences the affinity maturation process.

On one hand the existence of the dark zone is necessary to produce a large pool of different B-cell types, i.e. a high diversity of encoded antibody-types. This diversification is optimally realized with a fast proliferation of centroblasts without major interaction to antigens. The mutations are likely to occur randomly  and it is supposed that starting from low (but non-vanishing) affinity seeder cells better B-cells are found during this process, eventhough most clones will be of even less quality from a probabilistic point of view. The diversification process continues in the dark zone when the selection process has already started. Note that recycled B-cells that return into the dark zone don’t have been observed in the simulations . Therefore, the B-cells that proliferate in the dark zone are not recycled B-cells but directly stem from the original seeder cells.

After about $`8`$ days some high affinity cells have been found and identified in the selection process that takes place in the light zone. These roughly optimized B-cells are mostly recycled cells  and restart to proliferate in the light zone. A further enhancement of affinity to the antigen is based on these preselected B-cells. How does the ongoing process of diversification in the dark zone develop? The cells produced in this random process cannot overcome a certain B-cell quality. Therefore, a re-proliferation of already selected B-cells is more promising for the further development of the GC reaction. This process resembles a directioned selection process that replaces the random process in the dark zone. Ongoing proliferation and mutation in the dark zone does not only become senseless (as the quality of the cells in the dark zone are below the average B-cell quality in the light zone). They would also hinder the optimization process in the light zone because B-cells of relatively low affinity would take part in the selection process. Consequently, the selection of high affinity cells would be inhibited by a repetition of a first selection process with B-cells stemming from the dark zone. An early depletion of the dark zone eliminates the low affinty B-cells from the GC and in this way allows a fine-tuning of preselected B-cells. We conclude that a dark zone is needed in an early stage of the GC reaction to give centroblasts the possibility to proliferate and mutate independently of the FDC network. At later stages a fine-tuning of already preselected B-cells in the direct neighborhood of the antigen presenting FDCs becomes more important and a persisting dark zone inhibits affinity maturation.

One should be aware that this interpretation only provides a possible explanation of the correlation between morphology and affinity maturation that has been found in the framework of a spatial model for GCs. The basic concepts used here are in accordance with a widely accepted picture of GCs. Nevertheless, the used method has to be critically reviewed. The simulation are based on a regular lattice. We can widely exclude an effect of lattice symmetries on the results by comparing results from different lattice architectures . We have presented a statistical evaluation of two-dimensional simulations. Indeed, the general behaviour of the simulations does not change in three dimensions .

A more relevant problem is the representation of cells of different size on equal lattice nodes. We have shown that such a description leads to wrong results if considering chemotaxis as the driving mechanism for the development of the dark zone . In such a scenario typical cell velocities differ in two orders of magnitude and a more reliable description of cell volumes becomes unavoidable. We have therfore tested if the results are altered using this more sophisticated cell volume concept. For cells as slow as in the present simulations the more simple volume concept turns out to be completely sufficient (data not shown).

In regular lattice simulations the number of neighboring lattice nodes is fixed. In a one-cell-one-node description this determines the number of possible interaction partners of a cell. The effect of this restriction is not easily estimated. The simulation should be compared to a lattice free description. To this end we have implemented the first three dimensional dynamical Delaunay triangulation code with variable number of vertices . The cells are identified with the vertices (in a continuous space) and the Delaunay triangulation is used in order to calculate the neighbors of each cell. The triangulation has to be maintained after the movement of cells (vertices). In order to allow for proliferation and apoptosis the creation and deletion of vertices has to be included. The application of this method to GC reaction will provide another test for the stability of the results presented here.

### Acknowledgments

I thank Tilo Beyer, Andreas Deutsch and Gernot Schaller for intense discussions and valuable comments.