text:
# The local Hubble flow: A manifestation of dark energy

## 1 Introduction

The Hubble diagram for high redshift Type Ia Supernovae (Riess et al. 1998, Perlmutter et al. 1999) and the Boomerang and MAXIMA-1 measurements of the first acoustic peak location in the angular power spectrum of the cosmic microwave background (CMB) (de Bernardis 2000; Jaffe et al. 2000) strongly favor a standard cosmological model having the critical density ($`\mathrm{\Omega }=\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) and a dominant $`\mathrm{\Lambda }`$-like, “dark energy’ component at the present epoch ($`\mathrm{\Omega }_\mathrm{\Lambda }^00.7`$). In this Letter we argue that the cosmological $`\mathrm{\Lambda }`$ component is also important for understanding the very local Hubble diagram. Namely, there is a recognized riddle in observational cosmology – the linear and quiet Hubble flow in the very clumpy local universe (Sandage et al. 1972; Weinberg 1977; Peebles 1992; Baryshev 1994; Karachentsev & Makarov 1996; Baryshev et al. 1998; Sandage 1999).

We show that the Hubble flow emerges there, at a distance of about 1.5 Mpc, where the antigravity of the dark energy starts to dominate over the gravity of clumpy matter. We suppose that this is not a coincidence but the key to the riddle of the local Hubble law.

## 2 The enigma of the local Hubble flow

### 2.1 Formulation of the problem

In the standard cosmology the dynamics described by the Hubble law, $`V=H\times r`$, is a strict consequence of a uniform distribution of self-gravitating matter (Robertson 1955; Peebles 1993). Hence the predicted linear velocity field is only valid for scales where the universe is uniform. However, a very puzzling fact, long ago noted by Sandage, Tammann & Hardy (1972; STH), is that Hubble discovered his law in the distance interval $`120`$ Mpc, where the galaxies are very clumpy distributed. Indeed, this was deep inside a cell of uniformity, the size of which is still debated, but there are modern estimates from 30 Mpc up to 200 Mpc (Sylos Labini, Montuori & Pietronero 1998; Wu, Lahav & Rees 1999) Inside such cells the galaxies are fractally distributed, being clumpy at all smaller scales. This means that ”in fact, we would not expect any neat relation of proportionality between velocity and distance for these \[closeby\] …galaxies” ( Weinberg 1977).

On the other hand, recent studies of the local volume ($`<10`$ Mpc), based on accurate distances, confirm a ”cold” linear Hubble flow, with about the global Hubble constant $`H_0`$ (Sandage 1986; Karachentsev & Makarov 1996; Teerikorpi 1997; Giovanelli et al. 1999: Ekholm et al. 1999). Sandage (1999) recently expressed this surprising situation, as the ”extremely local rate \[of expansion\] is the same as the global rate to better than $`10\%`$”, and an ”explanation of why the local expansion field is so noiseless remains a mystery”. Indeed, the velocity dispersion around the local Hubble law is very low ($`\sigma _V5070`$ km/s) (Sandage 1986; Karachentsev & Makarov 1996; Ekholm et al. 2000; in the last work Cepheid distances gave $`40`$ km/s). That here is a real problem, has been also shown with N-body simulations (Governato et al. 1997): the expected velocity scatter close to galactic clumps like our LG is $`150700`$ km/s, depending on the CDM model.

One may thus divide the enigma of the local Hubble flow into the following Problems:

(1) The linear velocity law appears at distances (1.5 Mpc) which are about 1 percent of the scale where the Universe may finally appear uniform. Why does our highly non-uniform closeby environment expand as if it were uniform? Why is the rate of expansion the same locally and globally, i.e. deep inside and well outside of the cell of uniformity?

(2) What makes the local Hubble flow so quiet? Not only does the Hubble law exist, but it has a remarkably small scatter.

### 2.2 Previous attempts to explain the Hubble enigma

In their seminal paper, STH discussed the hierarchical galaxy distribution suggested by de Vaucouleurs (1971). Wertz (1971) and Haggerty & Wertz (1971) had made calculations on how the Hubble law deviates from linearity inside such a structure. STH did not find the predicted non-linearity and rejected the hierarchical model. However, they emphasized that the co-existence of the linear Hubble law and the large local clouds of inhomogeneous matter still is a dilemma. STH proposed two possible solutions. The linearly expanding universe may be inhomogeneous if its mean density is very low, $`<<\rho _{crit}`$. Or it may be filled by a totally uniform dark matter.

Modern redshift surveys, together with accurate measurements of the local velocity field, have sharpened this ”Hubble - de Vaucouleurs” paradox (Baryshev et al. 1998). Indeed, the 3-dimensional galaxy maps at scales from 0.1 Mpc up to 20-200 Mpc are well described by a scale invariant fractal distribution (Davis & Peebles 1983; Pietronero 1987; Peebles 1993; Sylos Labini, Montuori & Pietronero 1998; Teerikorpi et al. 1998; Wu, Lahav & Rees 1999). 111Note that Mandelbrot (1982) replaced the old concept of hierarchical galaxy distribution by the more adequate concept of stochastic fractals. Using new data on large scale structure and on the Hubble law, the solutions mentioned by STH were studied in Baryshev et al. (1998). ¿From the linear theory of the gravitational growth of density fluctuations (without dark energy), they concluded that either $`\mathrm{\Omega }_m0.01`$ or for the uniform dark matter $`\mathrm{\Omega }_{dark}0.99`$.

Nevertheless, these explanations are not sufficient, if considered together with recent cosmological observations which give $`\mathrm{\Omega }_m0.3`$. Also there is no direct evidence for a high-density uniform dark matter. Hence the riddle still exists.

## 3 Dark energy and dynamical model

### 3.1 Vacuum, quintessence, dark energy

Recent detection of the dark energy component provide a natural candidate for a high-density truly uniform background, giving a novel possibility for resolving the Hubble enigma.

In the cosmological theory there are several kinds of dark energy with a positive energy density ($`ϵ_Q=\rho _Qc^2>0`$) and negative pressure ($`p_Q<0`$). The dark energy or quintessence is a common name for: the classical Einstein’s cosmological constant $`\mathrm{\Lambda }`$, time-variable cosmological constant $`\mathrm{\Lambda }(t)`$, cosmological vacuum, effective scalar field and other forms of exotic matter with negative pressure (see reviews by Bahcall et al. 1999; Sahni & Starobinsky 2000). Dark energy is now much discussed and used for interpretation of observational cosmological tests (Wang et al. 2000; Podariu & Ratra 2000). The main reason to consider evolving DE is the “cosmic coincidence”: Why is the rapidly decreasing energy density of matter at the present epoch so close to that of the unchanging vacuum density?

The most striking property of the dark energy having the equation of state $`p_Q=wϵ_Q`$, with $`w[1,0)`$, is that its gravitating mass

$$M_Q=\frac{4\pi }{3}(1+3w)\rho _Qr^3$$
(1)
is negative for $`w<1/3`$. This produces cosmological ”antigravity” or a repulsion which accelerates the expansion of the universe. E.g. for Einstein’s $`\mathrm{\Lambda }=8\pi G\rho _\mathrm{\Lambda }/c^2`$ (cosmological vacuum with $`w=1`$), the gravitating mass is $`M_\mathrm{\Lambda }=8\pi \rho _\mathrm{\Lambda }r^3/3`$.

### 3.2 Dynamical model and the critical distance $`r_Q`$

We consider the dark energy component uniformly spreading everywhere, and existing alongside with luminous and dark matter in the immediate environment of our Galaxy. Hence, the local dynamics is determined by the competition of the gravity of dark (and luminous) matter and the antigravity of the DE. The mass-to-luminosity ratios in systems of different scales suggest that $`M/L`$ remains constant for $`r>0.5`$ Mpc (Bahcall, Lubin, & Dorman 1995), hence dark matter should be distributed like luminous matter on such scales. The correlation analysis of the space distribution of galaxies shows that the mean density of luminous matter (and hence, of dark matter) decreases with increasing scale so that the matter mass increases as $`M_m(r)r^D`$, where $`D`$, the fractal dimension, is between 1 and 2 at least on scales up to 20 Mpc (Wu, Lahav & Rees 1999 ). As its density is constant ($`M_Q(r)r^3`$), the DE component starts to dominate after some distance $`r_Q`$.

For a simple estimate of the scale beyond which vacuum or quintessence dominates, we consider two models that may imitate the environment of the Local Group, where the local Hubble law is accurately observed. The first model is a generalization of Sandage’s (1986) point-mass model, where a mass $`M`$ is placed on the dark energy background with density $`\rho _Q`$. In the second model, a spherical distribution of mass $`M_m(r)`$ is placed on the DE background.

The dynamics of expansion is described by the equation of motion which is the (1,1)-component of Einstein’s field equations under the assumptions of spherical symmetry, dust-like matter, and cosmological dark energy given by the energy-momentum (EM) tensor $`T_k^i=\rho _Qc^2diag(1,w,w,w)`$:

$$\ddot{r}=GM_{eff}/r^2;M_{eff}=M_m(r)+M_Q(r).$$
(2)
Here $`M_m(r)=M`$, or $`M_m(r)=M_{}(r/r_{})^D`$, respectively for the first and second model, $`M_Q`$ is given by Eq.(1).

For the point-mass the critical distance $`r_Q`$ corresponding to $`\ddot{r}=0`$, i.e. the distance where DE antigravity compensates matter gravity, is defined by:

$$r>r_Q=(3M/(4\pi \stackrel{~}{w}\rho _Q))^{1/3}$$
(3)
where $`\stackrel{~}{w}3w1`$. For $`w1/3`$ there is no repulsive force and $`r_Q`$ does not exist.

For the second model the dark energy term dominates dynamically at distances

$$r>r_Q=(r_{})^{\frac{D}{3D}}(M_{})^{\frac{1}{3D}}(\frac{3}{4\pi \stackrel{~}{w}\rho _Q})^{\frac{1}{3D}}$$
(4)
### 3.3 The temporal behavior of $`r_Q`$ and critical time $`t_Q`$

In the case of cosmological constant a galaxy which now is on the border of the vacuum dominated region around a mass clump, previously lied inside the gravity dominated sphere. In terms of comoving coordinates, $`r_Q`$ shifts outwards for increasing redshift. This prompts us to ask how $`r_Q`$ behaves if the DE density changes, as in quintessence models, where $`r_Q`$ may still shift outwards in comoving coordinates, but slower than for the constant vacuum.

We illustrate the shift of $`r_Q`$ in time for three types of DE models: 1) Einstein’s cosmological constant $`\mathrm{\Lambda }`$, 2) the quintessence with $`w=2/3`$, and 3) the case of “coherent” evolution of a quintessence ($`\rho _Q=k\rho _{hm}`$, where $`\rho _{hm}`$ is the uniform matter component) during the late history of the universe, with $`w=2/3`$, $`k=1`$. The condition that the covariant divergence of the total EM tensor (quintessence + matter) is zero implies that $`\dot{\rho }=3(\rho +\frac{p}{c^2})\dot{a}/a`$, where $`\rho =\rho _Q+\rho _m`$ and $`p=p_Q`$ for dust-like matter ($`p_m=0`$). Hence the DE density behaves as

$$\rho _Qa^{3(\frac{1+k+wk}{1+k})}$$
(5)
For $`1/k=0`$ (no homogeneous matter) Eq.(5) gives the usual (non-coherent) quintessence behavior $`\rho _Qa^{3(1+w)}`$.

Let us consider the ratio between two metric distances $`r_{gal}/r_Q`$, where $`r_{gal}`$ is the distance to a galaxy which now is in the DE dominated region and takes part in the Hubble flow. The distance $`r_{gal}`$ is simply proportional to $`a`$. The behavior of the critical distance $`r_Q`$ depends on the DE and matter models. E.g. for the point-mass model the ratio is $`r_{gal}/r_Q=a^{(w/(1+1/k))}`$, which for $`1/k=0`$ gives $`a^w`$, corresponding either to Einstein’s $`\mathrm{\Lambda }`$ ($`w=1`$) or to the quintessence (for $`w1/3`$).

We define a characteristic time for the antigravity dominance. It is counted from the critical epoch $`t_Q`$, corresponding to the critical scale factor $`a_Q=a(t=t_Q)`$ when DE antigravity starts to dominate over the gravity of matter concentration component for a galaxy which presently ($`t=t_0`$) is at the distance $`r=r_{gal}=2r_Q^0`$, i.e. well within the DE dominated region. E.g. for Sandage’s point-mass model thus defined critical scale factor is

$$a_Q=2^{(1+1/k)/w}$$
(6)
With $`1/k=0`$ this gives $`a_Q=0.5`$ for $`w=1`$ and $`a_Q=0.354`$ for $`w=2/3`$, while with $`k=1`$,$`w=2/3`$ the critical epoch is at $`a_Q=0.125`$. This illustrates the difference in critical times for different quintessence models.

## 4  Specific properties of dark energy dominated dynamics

### 4.1 Structure evolution in DE dominated regions

The theory of structure formation in the general case of an evolving DE component is still under construction, though there are some results for particular cases (Sahni & Starobinsky 2000; Fabris & Goncalves 2000). If during the structure formation period most of the energy of the universe produces antigravity and hence resists gravitational collapse, it is impossible for the structure to grow at all. On the other hand, if antigravity dominates only a short period, as in the case of the cosmological constant, then even for $`\mathrm{\Omega }_\mathrm{\Lambda }1`$ it has only a slight effect on the infall velocities around growing structures (Peebles 1984; Lahav et al. 1991). This is clear from Eq.(2) which gives for the ratio of DE and point-mass accelerations $`g_Q/g_ma^{3w}`$, hence for $`w=1`$ the dynamical influence of the constant $`\mathrm{\Lambda }`$-component decreases very fast, while for $`w=2/3`$ the antigravity acts longer.

Another important consequence of the DE component is that within $`\mathrm{\Lambda }`$ dominated regions initial peculiar velocities decrease due to adiabatic cooling. The linear analysis of gravitational instability, following Zeldovich (1965), shows that only decreasing modes exist in the vacuum-dominated region (Chernin, Teerikorpi & Baryshev 2000). This property of the accelerated universe is important for understanding the local Hubble flow. Our simple dynamical models led to the notion of the distance $`r_Q`$ to the border between the matter and DE dominated regions. If a galaxy now beyond $`r=r_Q`$ also in the past was long enough in the DE dominated region, this would readily explain the smooth Hubble law at such distances. Then initial peculiar velocities induced by the mass of the Local Group (and similar more distant groups) have decreased because of the adiabatic cooling ($`\delta V(now)=\delta V(z)a(z)`$).

Therefore a compromise situation is possible when the time-variable DE is dominant during a period, sufficient for the structure growth in early times and for the cooling of velocity dispersion in late times. This means that such DE dominated regions between mass concentrations are ”pacific oceans” where the Hubble law appears and the global Hubble constant may be measured even locally.

### 4.2 Restoration of the Hubble flow in volumes with bulk motion

We referred to two properties of the cosmological vacuum, namely uniformity and antigravity. Let us now consider an important third property of the vacuum. It is clear that the medium with the equation of state $`p=\rho c^2`$ described by the cosmological constant has EM tensor which is proportional to the metric tensor $`T_i^k=\alpha g_i^k`$. Hence it has the remarkable mechanical property of vacuum: rest and motion can not be discriminated relative to it. In other words, any matter motion is co-moving to this medium.

When the vacuum dominates over the self-gravity of matter, matter masses like galaxies move as test particles in the antigravity of the vacuum. The vacuum accelerates the motions tending to form a regular kinematic pattern with a linear velocity field. The expansion rate in such a flow depends only on the value of the vacuum density.

These considerations can be applied to a local vacuum-dominated area, no matter how fast is the bulk motion of the area against the CMB. They are valid as well for the large-scale matter distribution; in both cases, the flow will be linear and the rate of expansion will be the same on any scale. And the flow will also be stable against velocity or density perturbations.

## 5 The critical distance $`r_Q`$ for the Local Group

Let us calculate the critical distance for the Local Group in the case of the point-mass model and the cosmological vacuum ($`w=1`$ and $`\rho _Q=\rho _\mathrm{\Lambda }`$ is constant in time).

The most recent estimate for the mass of the LG gives $`M_{LG}=2\times 10^{12}M_{}`$ (van den Bergh 1999). The vacuum density from SNIa observations is $`\rho _\mathrm{\Lambda }0.7\rho _{crit}=4.7\times 10^{30}h_{60}^2`$ g/cm3. Then the distance $`r_\mathrm{\Lambda }`$, where the vacuum starts to dominate, is $`1.5`$ Mpc in the point-mass model. If one changes $`M_{LG}`$ or $`\rho _\mathrm{\Lambda }`$ by a factor of two, the critical distance is changed by 26 percent. If one chooses the mass distribution model with fractal dimension $`D=1`$ and $`M_{}=2\times 10^{12}M_{}`$ at $`r_{}=1`$ Mpc, then Eq. (4) gives $`r_\mathrm{\Lambda }=1.8`$ Mpc.

Hence, $`r_\mathrm{\Lambda }`$ is robustly put into the range from 1 to 2 Mpc. This value is surprisingly close to the distance where the Hubble law emerges, which is 1.5 Mpc (Sandage 1986). Is this just a coincidence? Not at all, we suggest that this is rather a major feature of the local matter flow. Indeed, the dynamical dominance of the dark energy at 1-2 Mpc and beyond means that the homogeneous vacuum provides the dynamical background for the Hubble flow of matter, on these distances.

## 6 The local and global values of the Hubble constant

Now we can understand Problem I – the near-equality of the local and global Hubble constants – if the undistorted Hubble law is the signature of the dynamical dominance of the uniform DE component. In regions, which have for a sufficient time been DE-dominated, the rate of expansion – in the first (and main!) approximation – is expressed via the DE density: $`H_0(8\pi G\rho _Q/3)^{1/2}`$. The lumpy, fractally distributed matter gives only a small correction to this formula for the global Hubble constant, because its average density is much less than the DE density. The very interesting aspect of the dominating uniform DE component is that it allows us to predict the value of the Hubble constant from a local estimate of the dark energy density.

If we identify the distance $`r_{start}`$ where the Hubble law starts, with the critical radius $`r_Q`$ for our Local Group, then the predicted value of $`H_0`$ is

$$H(predicted)(\frac{2GM_{LG}}{\stackrel{~}{w}})^{1/2}(r_{start})^{3/2}$$
(7)
where $`M_{LG}`$ is the mass of the Local Group. With the ”standard” values $`M=2\times 10^{12}M_{}`$ and $`r_{start}=1.5`$ Mpc, the Eq.(7) gives for $`H_0`$ values from 53 km/s/Mpc ($`w=1`$) to 74 km/s/Mpc ($`w=2/3`$), covering the values obtained for $`H_0`$ by different research groups. More conservatively, these are upper limits, if in fact $`r_Q<r_{start}`$. We regard this agreement with direct measurements of $`H_0`$ as evidence for the cosmological vacuum and quintessence models. Note that this reasoning for the DE component is based on locally measured quantities ($`M_{LG}`$ and $`r_{start}`$) and is independent of the supernovae determination of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ (which did not depend on the value of $`H_0`$).

## 7 The quietness of the local Hubble flow

The possibility to solve Problem II may be illustrated by the following simplified reasoning. At any epoch $`a(t)`$ and on the scale $`r`$ one may regard $`\delta VHr`$ as an upper limit to the scatter around the Hubble law. For models considered in Sec.3.3 the curvature is zero and hence from Eq.(5) it follows that the Hubble constant $`H(a)`$ depends on the scale factor as $`a^{(3/2)\beta }`$, where $`\beta =(1+k+wk)/(1+k)`$. For example in the case of vacuum, i.e. $`w=1`$ and $`1/k=0`$ we have $`H=const`$, and in the case of dust matter $`k=o`$ hence $`Ha^{3/2}`$.

So the scatter $`\delta V(a)<H_0r_0a^{1(3/2)\beta }`$, where $`r_0`$ is the considered scale at the present epoch. E.g. for $`r_0=5`$ Mpc and $`H_0=60`$ km/s/Mpc, $`\delta V(a)<300a^{1(3/2)\beta }`$ km/s. This is also valid at the critical epoch $`a_Q`$ after which the peculiar velocities within a DE dominated region start to cool down adiabatically, so that $`\delta V(now)=\delta V(a_Q)a_Q`$. Hence one expects now $`\delta V(now)<300a_Q^{2(3/2)\beta }`$ km/s.

For simple numerical examples we use the critical epoch $`a_Q`$ calculated in Sec.3.3. In the case of the cosmological constant we have $`a_Q=0.500`$, hence $`\delta V(now)<75`$ km/s. For the quintessence with $`w=2/3`$, $`a_Q=0.354`$ and $`\delta V(now)<63`$ km/s, and for the case of coherent evolution $`w=2/3`$, $`k=1`$, the critical scale factor $`a_Q=0.125`$ and now one expects to observe $`\delta V(now)<38`$ km/s.

Of course, these numbers just illuminate the trend for different DE models. A complete treatment of the dispersion problem needs high-resolution large-volume simulations of structure formation on the evolving DE background.

## 8 Conclusions

We conclude that a cosmological dark energy component allows one to understand the emergence of the linear and quiet Hubble flow at the border of the Local Group. The coincidence of the $`r_Q`$ with the starting distance of the Hubble law is seen to have the same roots as the “cosmic coincidence” of $`\mathrm{\Lambda }`$ and matter densities. A summary of our results is:

* We introduce the concept of the critical distance $`r_Q`$ where the repulsive force starts to dominate over the gravity of matter. For the Local Group $`r_Q1.5`$ Mpc in the case of $`\rho _Q=0.7\rho _{crit}`$. This is close to the distance where the local Hubble flow emerges. We make this new cosmic coincidence as a starting point for understanding the cold local Hubble flow as reflecting the dark energy dominated region around us.

* In dark energy dominated regions of the universe (”pacific oceans”) the linear Hubble law exists due to the homogeneity of the DE component, with the same Hubble constant at first approximation determined by the dark energy density. This resolves Problem I that the Hubble law not only starts at the border of the DE dominated region, but it appears with the $`H_0`$ as expected from the dark energy density.

* The solution of Problem II follows from the fact that a galaxy spending a sufficient time in the DE dominated region looses its initial peculiar velocity which cools down adiabatically. The mysterious quietness of the local Hubble flow does not look at all dramatic in the context of the dark energy driven expansion.

* Our results support a dark energy component which presently dominates the dynamics of the Universe on scales from the Local Group to the Hubble radius. Even the closeby galaxy universe is now seen as a cosmic laboratory where all the physical ingredients of the universe: luminous matter, dark matter, and the dark energy, may be detected. Hence, the very local volume is extremely important for the study of global properties of the universe. Cosmology starts immediately beyond the border of the Local Group!

In other environments the DE dominance may start from smaller or larger distances. So, around the Coma cluster, the vacuum starts dominating from about 20 Mpc. For a cosmologist in Coma, cosmology begins around such a distance. In such DE dominated regions the quiet Hubble flow is produced by the same vacuum density and hence the Hubble constant is the same.

One now realizes that when Hubble found his linear redshift law in the closeby lumpy environment, he actually saw the influence of the mysterious uniform dark energy. Though the value of the vacuum or DE density is still poorly understood in cosmological theory, the observations make one agree with the old words by Lemaître who was the first to associate $`\mathrm{\Lambda }`$ with the vacuum having positive energy density and negative pressure leading to a cosmological repulsion: ”Everything happens as though the energy in vacuo were different from zero.” (Seitter & Duemmler 1989).

Acknowledgements

We are grateful to Chris Flynn for helpful comments. This study has been supported by The Academy of Finland (project “Cosmology in the local galaxy universe”).

References

Bahcall, N.A., Lubin, L.M., Dorman, V. 1995, Astrophys.J., 447, L81

Bahcall N., Ostriker J., Perlmutter S., Steinhardt P. 1999, Science 284, 1481

Baryshev Yu. 1994, Astron.Astrophys.Trans.,5, 15

Baryshev Yu., Sylos Labini F., Montuori M., Pietronero L., Teerikorpi P. 1998, Fractals, 6, 231

Chernin A., Teerikorpi P., Baryshev Yu. 2000, Adv. Space Res.(in press)

Davis M., Peebles P.J.E. 1983, ApJ, 267, 465

de Bernardis et al. 2000, Nature, 404, 955

de Vaucouleurs G. 1971, Science, 167, 1203

Ekholm T., Lanoix P., Teerikorpi P., Paturel G., Fouqué P. 1999, A&A, 351, 827

Ekholm T., Baryshev Yu., Teerikorpi P., Hanski M., Paturel G. 2000, A&A (submitted)

Fabris, J.C., Goncalves, S.V.B., 2000, gr-gc/0010046

Giovanelli R., Dale D., Haynes M., Hardy E., Campusano L. 1999, ApJ, 525, 25

Governato F. et al. 1997, New Astr., 2, 91

Haggerty M., Wertz R.J. 1971, MNRAS, 155, 495

Jaffe A.H. et al. 2000, astro-ph/0007333

Karachentsev, I., Makarov, D. 1996, AJ, 111, 794

Lahav O., Lilje P., Primack J., Rees M. 1991, MNRAS, 251, 126

Mandelbrot B. 1982 The Fractal Geometry of Nature (W.H. Freeman )

Peebles P.J.E. 1984, ApJ, 284, 439

Peebles P.J.E., 1992, Lecture in the symposium ”The Cosmic Microwave Background Radiation and the Large-Scale Structure of the Universe after COBE” (Stockholm 1992)

Peebles, P.J.E. 1993 Principles of Cosmology (Princeton: Princeton Univ. Press)

Perlmuter, S. et al. 1999, ApJ, 517, 565

Pietronero L. 1987, Physica, A144, 257

Podariu S., Ratra B. 2000, ApJ, 532, L109

Riess, A.G. et al. 1998, AJ, 116, 1009

Robertson, H.P. 1955, PASP, 67, 82

Sahni V., Starobinsky A., 2000, astro-ph/9904398

Sandage, A. 1986, ApJ, 307, 1

Sandage, A. 1999, ApJ, 527, 479

Sandage A., Tammann G., Hardy E. 1972, ApJ, 172, 253

Seitter, W.C., Duemmler, R. 1989, in Morphological Cosmology (eds. P. Flin, H.W. Dverbeck), Springer, Berlin, pp. 377-387

Sylos Labini, F., Montuori, M., Pietronero, L. 1998, Phys.Rep., 293, 61

Teerikorpi, P. 1997, ARA&A, 35, 101

Teerikorpi, P. et al. 1998, A&A, 334, 395

van den Bergh, S. 1999, A&A Rev., 9, 273

Wang L., Caldwell R.R., Ostriker J.P., Steinhardt P.J. 2000, ApJ, 530, 17

Weinberg, S. 1977 The First Three Minutes (Basic Books, New York), p.26

Wertz J.R. 1971, ApJ, 164, 227

Wu, K.K.S., Lahav, O., Rees, M.J. 1999, Nature 397, 225

Zeldovich, Ya.B. 1965, Advan. Astron. Ap., 3, 241


format: text

text:
# Adaptive binning of X-ray galaxy cluster images

## 1 Introduction

A problem often encountered with the processing of images, especially those from X-ray telescopes, is that there is a large dynamic range present in intensity. In order to examine the structure of emission in regions of an image where there are few photons, it is necessary to bin the data using bins of large angular size. An alternative is to smooth the data with a Gaussian kernel of large angular size. Conversely, to examine structure in areas of high emission, a small bin-size or kernel is sufficient. If one wishes to examine regions with a wide variation in count-rate, then some procedure which has a variable bin or kernel is required to create statistically optimal maps.

Another related problem, which motivated the present work, is the creation of optimal colour maps. The ratio of photon counts in a low intensity region will have a large associated statistical error, so the counts need to be averaged over a large region to form the ratio. The algorithm needs to account for the statistical accuracy in the separate bands, not just the total image.

One technique which allows structure on a wide range of scales to be revealed is adaptive kernel smoothing (Ebeling, White & Rangarajan 2000; Huang & Sarazin 1996). The asmooth algorithm of Ebeling et al. (2000) convolves an image with a Gaussian kernel. If the signal within the kernel, applied to part of an image, is of a chosen minimum significance above the local background, then the convolved signal is added to the output image. The kernel is increased in size until all the counts in the input image are added to the output image.

asmooth is very good at finding low surface brightness features, such as filaments. It also especially useful for producing images for display, and for identifying bright features against a background. It is not a universal smoothing tool, however, and care needs to be taken before using it in a particular situation. For instance, it is designed to identify positive features against a background, but not to identify ‘holes’ in emission. Due to the way it identifies significance against the local background, it may produce spurious features in areas where there is bright flat emission, for instance the cores of galaxy clusters, where Poisson noise may become significant. Used with care, asmooth is a useful routine in many situations.

The analysis method outlined in this paper, *adaptive binning*, was originally developed to produce temperature, column density and colour maps of the Perseus cluster in Fabian et al. (2000). We required a simple algorithm to compare the relative fluxes in different X-ray bands to a theoretical plasma emission model. In order to be able to spatially resolve the ratios, it is necessary to bin or smooth the data. Binning, rather than smoothing, is a better method for producing spatial estimates of counts for analysis, because it does not spread counts around an image beyond the boundaries of the bin.

In this paper we use data analysed from the 24 ks ACIS-S *Chandra* X-ray observation of the Perseus Cluster, Abell 426, from Fabian et al. (2000), to demonstrate the adaptive binning method. The Perseus Cluster is the brightest X-ray cluster in the sky, and contains a large cooling flow (Fabian et al. 1994) with a mass deposition rate of approximately $`300\mathrm{M}_{}\mathrm{yr}^1`$. Interesting substructure in the cluster due to its central radio source makes it an ideal object for our purpose. We will not attempt here to analyse the physics behind the observations, which is left to Fabian et al. (2000) and further work. We also study the method by applying it to simulated data.

The algorithm is presented in a logically consistent way, from intensity adaptive binning to colour adaptive binning. Colour binning, as mentioned above, was our original motivation, developed to analyse X-ray colour cluster images.

## 2 Intensity binning

Adaptive intensity binning is the simplest case of the algorithm. It attempts to adaptively bin a single image based on the number of photons in each region. The basic method is to bin pixels in two-dimensions by a factor of two, until the fractional Poisson error of the count in each bin becomes less than or equal to a threshold value. When the error is below this value, those pixels are not binned any further. The algorithm, in detail, is as follows.

1. Each pixel in the image is put into its own ‘bin’, the term we use for a collection of pixels. A pixel here means one of the individual picture elements which form the input image. Essentially the image is initially divided into imaginary $`1\times 1`$ pixel bins.
2. If there are $`n_i`$ pixels in bin $`i`$, the total count in the bin is $`c_i`$, and the background per pixel is $`b`$, the net count in the bin is simply defined by

$$s_ic_in_ib.$$
(1)
3. The fractional error on the net count in the bin is

$$\frac{\sigma _{\left(s_i\right)}}{s_i}=\frac{\sqrt{c_i+n_ib}}{c_in_ib}.$$
(2)
This is also the error on the average count in each pixel.
4. If the fractional error is less than or equal to a threshold value, then the pixels in the output image which correspond to the pixels in the input bin are set to the average mean count, $`s_i/n`$. The fractional error of the net count in the bin is also stored in the pixels in an ‘error image’. Additionally the bin is marked as having been processed.
5. Each bin is merged into its neighbouring three bins, to make new bins containing $`2\times 2`$ of the previous bins. The four bins with the lowest $`x`$ and $`y`$ coordinates (lowest declination and highest right-ascension) are merged, as is each consecutive set of four bins. Any bins which have already been processed are ignored in the merging. It is useful to remember we are considering a bin as a list of pixels. Pixels which have already been set in the output image are ignored in future iterations.
6. The process is repeated from 1 until there is only a single bin remaining.
7. A ‘bin-map’ is also produced by the algorithm, giving an identification number for each processed bin in terms of the pixels which it contains. Using the bin-map, any image of that size can be binned using the same bins.

### 2.1 Demonstration

A simple demonstration of the algorithm operating on a $`4\times 4`$ pixel image is shown in Fig. 1. We demonstrate the process using a threshold fractional error of 0.1 and no background counts. (a) shows the image before binning. We first look for individual pixels (bins of $`1\times 1`$ pixels) with a Gaussian fractional error in the count less than or equal to the threshold. Only one pixel does; it is shown ‘painted’ in (b). We then bin the remaining pixels by a factor or 2. We examine the binned pixels to see whether any have errors less than the threshold. The three pixels in the top-left corner do, so they are averaged together and painted in (c). We then bin again with a $`4\times 4`$ pixel bin. The remaining pixels have an error less than the threshold, so they are averaged and painted in (d), the final output image. Had the error on the final bin been larger than the threshold, it would have been binned anyway.

### 2.2 Real example

We show in Fig. 2 a ‘raw’ image of the Perseus cluster in the 0.5–7 keV band from the *Chandra* observation, binned with 2 arcsecond pixels. The image has been exposure-corrected. It also has been corrected for the readout node lines of the ACIS-S detector (this smears out the point source at the centre). Contours were placed on the image showing 8 levels spaced equally in square-root space between 20 and 300 counts.

The Perseus Cluster is very bright, so we see much structure in the data without doing more than simple binning. The contours shown have been smoothed, but even they break up at large radii due to the count rate being swamped by Poisson noise.

Fig. 3 shows the raw image of Perseus adaptively binned with a threshold fractional error of 0.06. The X-ray background was ignored since it was low, even near the image edges. The contours on the image near the centre match the contours in the raw data well, showing the algorithm works in this regime. At the sides, the bins are quite large and blocky, but otherwise there do not appear to be any edge-effects. There is still noise present in the image, but at a level much reduced from the raw data.

Fig. 4 shows the cluster again adaptively binned, but now with a threshold error of 0.04. Note how the bins are larger, but the level of noise is significantly less. A real feature has been lost from the image, however. A point source present in the upper left radio lobe has disappeared. Its counts were merged into rest of the emission from that region.

In Fig. 5 is shown the error map for the adaptively binned map above. In it are easily visible the changes in bin size as the count rate decreases towards the outside. As the count rate decreases, the error of the bins increases until it reaches the threshold, and then the bin-size doubles.

For comparison, Fig. 6 shows an adaptively smoothed (AS) image of the cluster calculated from the data in Fig. 2. It was smoothed with a minimum significance of 4–$`\sigma `$ using the asmooth algorithm of Ebeling et al. The contours are at the same levels as the raw image.

The AS image contains sharp positive features, such as the edges of the radio lobes. It does not perform as well in the determination of the number of counts in the negative features, such as the radio lobes themselves, where the number of counts per pixel is half that of the raw data. However, the algorithm is designed to find positive features, so using it to look at holes is a misapplication. Also apparent are some edge effects, where asmooth does not find enough counts to place a high significance on the generated features. The edge effects are avoidable by smoothing a larger area of sky than required, and cropping the image thereafter. One disadvantage of this is that asmooth running time increases quickly with image size.

In Fig. 7 we show another adaptively binned image of the Perseus cluster. This, however, was created from data from an observation using the ACIS-I detector on *Chandra*, with an exposure of 18.6 ks. We present it because the number of counts per pixel is lower than the ACIS-S image, and so it is a good demonstration of the algorithm used on data with more noise.

Adaptive binning of the ACIS-I image brings out an interesting feature of the cluster; it clearly demonstrates that the southern rim of the northern radio lobe (hole) lies south of the nucleus. There is a linear diagonal structure present to the south of the nucleus, running from the north-west to the south-east. This structure is also present in the raw data, but is not clear on the ACIS-S image, despite the longer exposure, as the ACIS-S raw image contains a dark strip due to the node-line of the detector, where the effective exposure is short.

The dark bin to the south-east of the nucleus shows a ‘stranded bin’, an occasional problem with the algorithm, which we discuss later.

### 2.3 Notes on the algorithm

#### 2.3.1 General

1. The algorithm allows bins of width $`n`$ to contain pixels with coordinates $`(x,y)`$ in the range $`nix<n(i+1)`$, $`njy<n(j+1)`$, where $`i`$ and $`j`$ are integers, and $`i0`$, $`j0`$. The coordinate origin is $`(0,0)`$. $`x`$ and $`y`$ increase in the directions of decreasing right-ascension and increasing declination, respectively.
2. If the unbinned image does not have sides which are an equal length of pixels, or are not a power of two, then the bins on one or both sides will be truncated. This can lead to bins with a small numbers of counts and a large error on the final pass.
3. As is the case with conventional binning, sources with small spatial extent may be split between two bins if they cross a bin boundary. This can be solved by allowing sub-bin positioning, described as follows. For each pass, bins are allowed to be placed at positions displaced by a multiple of some fraction of their length from their conventional positions. The errors are calculated for each possible bin, and the bins are sorted into order of ascending error, discarding those with an error larger than the threshold. The bins are painted in that order, ignoring those bins which had already been partially painted.

The disadvantage of this technique is that the output image appears ‘dithered’. Bins are painted close to one another, but without enough space to put another bin between them. The result is a loose cluster of small bins inside a larger one. The original algorithm appears to produce better results and is less dependent on small variations of errors.
4. The bins are doubled in size between each pass, rather than slowly increased in size, as this ensures a whole number of smaller bins fit inside a larger bin. If this is not the case, then gaps are left between bins in the binning, leading to dithering.
5. Allowing rectangular bins is not useful. It is ambiguous in which direction a rectangle could lie, or be extended. Their use could introduce spurious linear structures into the output image.
6. The use of non-square bins is a possible technique. Any two dimensional object that tiles together could be used, but it is not clear what is the optimal shape. Square bins are a good choice for a number of reasons, including the fact that pixels are naturally square in most detectors; square bins act as scaled pixels. However, there is probably some room for additional work to find better shapes. An alternative idea is to use contour levels in an AS image to define ‘bins’, as demonstrated by Sanders, Fabian & Allen (2000).

#### 2.3.2 Overlayed colour images

The adaptive-binning algorithm is useful for producing images demonstrating in real colour where the soft, medium and hard areas of emission are. First an intensity image can be binned to produce a bin-map. Then images in three energy bands can be binned themselves using the bin-map. A software package such as gimp can then be used to combine the three images as, for example, the red, green and blue layers in a single image. The contrast of the image can then be increased to highlight the areas of hard and soft emission.

The advantage of doing this rather than using raw data to make the image, is that the low-intensity regions are not dominated by noise. The technique was used to create the colour image of the Perseus cluster in Fabian et al. (2000).

#### 2.3.3 Stranded pixels and contiguous regions

One ‘feature’ of the algorithm is that pixels may become ‘stranded’. If a bin has a reduced count relative to its neighbours, it may be left out on a pass and have no nearby pixels to be binned with on the next. It may be left until the final pass, and that area will be binned with the remaining pixels. Such pixels could be cosmetically removed by replacing them with their original value.

A modified version of the algorithm avoids stranded pixels. If a bin consists of more than one non-contiguous sets of pixels, then the bin is split into several contiguous regions, and each region treated by itself. A stranded set of pixels will remain until the end, where it will be binned by itself, as it is isolated. Here, two pixels are said to be connected if one pixels is one of the neighbouring eight pixels of the other pixel.

This modification has a couple of disadvantages. Firstly, finding contiguous sets of pixels is fairly slow. Also, parts of the image which would be connected were the image larger, are treated by themselves, and therefore have a larger statistical error on their count. With these problems in mind, it is probably a useful modification for many applications of the algorithm. For simplicity, we will use the unmodified version of the algorithm in the following sections.

### 2.4 Simulated cluster

In order to properly test the properties of the algorithm, we attempted to simulate a simple image of a cluster. We used a $`\beta `$-model to simulate the cluster with the following form (Sarazin 1988):

$$S_b(r)=S_0\left[1+\left(\frac{r}{r_c}\right)^2\right]^{0.53\beta },$$
(3)
where $`S_b`$ is the surface brightness at a radius $`r`$, $`r_c`$ is a critical radius and $`\beta `$ is a parameter. We created an image of $`512\times 512`$ pixels, where $`r_c`$ was 128 pixels, $`S_0`$ was 100 counts/pixel, and $`\beta `$ was 0.67. The meaning of the word count, as used here for surface brightness, is the expectation value of an observation, and is allowed fractional values. The cluster was positioned in the centre of the image. It is shown in Fig. 8(a). The contours show constant surface brightness from 10-90 counts in linear 10 count intervals.

To simulate how such a cluster would appear, if observed, we took the surface brightness at a particular pixel, and randomly generated a value from a Poisson distribution with that surface brightness as the expectation value, to which we added a background with an expectation value of 20 counts/pixel. The Poisson image is shown in Fig. 8(b), with contours between 30 and 100 with the 10 count intervals. The contours are further out because of the background. The data was smoothed in contouring, as it was too noisy otherwise.

We processed the Poisson simulated data with our adaptive binning algorithm to produce Fig. 8(c). Here a fractional error threshold of 0.02 was used. Fig. 8(d) shows the output of the method using a threshold of 0.04, i.e. twice as large a threshold as (c). The contours are at the same levels as the original surface brightness image.

For comparison, we show the output from the asmooth program using a minimum feature significance of 4-$`\sigma `$ in Fig. 8(e). asmooth, by default calculates the local background and does not take a global value. To create this image we smoothed the simulated data and subtracted the background manually. The count per pixel becomes negative in the outer regions, as asmooth does not know about the global background, which also means features shown may have a lower significance than 4-$`\sigma `$.

#### 2.4.1 Spatial distributions of differences

We constructed in Fig. 9 images displaying the absolute fractional differences between the reconstructed and model surface brightness images. Spatial correlations in errors should be visible on the images.

Fig. 9(a) shows the absolute differences between an adaptively binned image with 0.02 fractional errors (Fig. 8(c)), and the original surface brightness image binned with the same bins. (b) shows the same except using bins constructed for 0.06 fractional errors (Fig. 8(d)). There do not appear to be any correlations between the two images, in terms of where the differences are large and small. The differences appear noise-like and random. Note, however, that if the adaptively binned image were binned with inappropriate bins, then this would not show in (a) and (b) due to the surface brightness image being binned with the same bins.

Fig. 9(c) is an image showing the absolute differences between adaptively binned image with 0.02 fractional errors and the unbinned surface brightness image. Due to the count rate varying across the bin, only the centres have a low absolute fractional difference. This would be the same for any binning process.

Fig. 9(d) shows the differences between the AS image in 8(e) and the surface brightness image for comparison. Visible are the edge effects of the algorithm. Also present are some peaks at the cluster core. These are due to the local background estimation of the algorithm measuring significant detections against the flat core of the cluster centre.

#### 2.4.2 Histograms of differences

A more useful and quantitative analysis of the performance of the adaptive binning algorithm can be made by plotting the fractional differences (the absolute value of which is in the above images) as histograms, shown in Fig. 10. (a) displays the difference between the 0.02 error adaptively binned image and the surface brightness image binned using the same bins *for each pixel*. (b) is the same as (a), except using 0.06 error binned images. (c) shows a histogram of the differences between the 0.02 error adaptively binned image and the unbinned surface brightness image. The difference between an image produced by asmooth (Fig. 8(d)) and the surface brightness image is shown in (d) for comparison.

Fig. 10(a) and (b) approximate Gaussian distributions. The widths of the distributions are close to the threshold fractional error values used by the algorithm, 0.02 and 0.06. (c) shows a wider distribution, which is obtained due to the count varying across the bin in the surface brightness image, but not in the adaptively binned image. The results in (a) and (b) give confidence in the effectiveness of the algorithm, as the distributions are symmetric and Gaussian.

The histogram for the AS image in (d) has a very narrow peak, but there are long tails, mainly due to the edge-effects. Excluding the outer region almost removes the tail.

We also analysed the differences between the input surface brightness image and the adaptively binned image using the $`\chi ^2`$ statistic. If there are $`S_i`$ total counts in the surface brightness image in bin $`i`$, and $`s_i`$ in the adaptively binned image, and $`n_i`$ pixels in the bin, then the statistic is

$$\chi _i^2=\frac{(S_is_i)^2}{S_i+n_ib}.$$
(4)
We plotted the distribution of $`\chi ^2`$ for the bins in the adaptively binned image. It showed reasonable agreement with the distribution predicted for one degree of freedom.

## 3 Colour binning

A simple way to investigate the inner structure of a cluster is to create X-ray images in different energy bands. The images can then be divided to show the areas of hard and soft X-ray emission, known as a colour image. Additionally plasma codes, such as the mekal thermal model, can be fitted to the relative counts, enabling one to map temperature, metallicity and absorption.

One faces a similar problem to the simple intensity binning problem. How can ratios of counts be formed so that the error on the result is accurate enough? A binning algorithm must be able to use the same size bins on the images to be divided, plus it must adapt its bin size to the counts in each of the bands, rather than just one.

The colour adaptive binning algorithm is similar to the intensity algorithm, except it uses an input image for each band, to produce an output image for each band. Each input image is binned using the same bins. The bins are defined by the error on a ‘combined colour’, which folds the ratios of all the bands, not the error on the intensity. If $`s_k^i`$ is the background-subtracted count in band $`i`$ for bin $`k`$, the combined colour is defined as

$$C_k=\left[\left(s_k^1/s_k^2\right)/s_k^3\right]/\mathrm{}.$$
(5)
The fractional error on the combined colour, if there are $`N`$ bands, is

$$\frac{\sigma _{\left(C_k\right)}}{C_k}=\left[\underset{i=1}{\overset{N}{}}\frac{c_k^i+n_kb^i}{\left(c_k^in_kb^i\right)^2}\right]^{1/2}.$$
(6)
This error is a useful quantity because it takes into account the count in each band in a symmetric way. To make an actual colour map, two binned output images are divided.

An alternative approach is to adaptively bin a total intensity image. The resultant bin-map can then be applied individually to each of the X-ray bands. Ratios can be formed by dividing the output images. This approach has its advantages when several bands are being used. If a sharp feature is only present in one band, then the source will be binned with a large bin using the combined-colour method, due to the large statistical uncertainty on the colour. However, binning using the total intensity leads to colours with larger range of errors than the combined-colour approach.

### 3.1 Real example

The combined colour adaptive binning method has been applied to the Perseus data shown earlier. Three images of the cluster were created in different bands, labelled from 0–2, in the energy ranges 0.5–1, 1–2 and 2–7 keV.

The colour adaptive binning algorithm above was applied to the data using a fractional errors on the combined colour, $`\left(s^0/s^1\right)/s^2`$, of 0.12 and 0.2. Again the background counts in each image were small enough to be neglected.

Fig. 11 shows the $`s^0/s^1`$ colour image produced with a fractional error of 0.2, which is sensitive to X-ray absorption; darker shades indicate more absorption. The absorption shadow of an infalling dwarf galaxy, discussed by Fabian et al. (2000), is clearly seen close to the centre of the image. The rest of the image has approximately uniform colour, indicating there are no further strong absorption features.

Fig. 12 shows the $`s^1/s^2`$ colour image produced with a fractional error of 0.12, which is sensitive to temperature; lighter shades indicate lower temperature. From this simple analysis, there appears to be a temperature gradient in the cluster, with cooler temperatures towards the centre and the brighter regions. This agrees with cooling flow models. The analysis emphasises a spiral structure, which is present right down to the cluster centre (compare with Fig. 7).

Using the combined colour method to generate the bins, the fractional errors on the bins for the individual colours, both range from 0.05 to 0.1 using a fractional error of 0.12. We have shown the colours with different fractional errors due to the relative strengths of features in the colours.

### 3.2 Spectral fitting

As mentioned above, one can fit the relative counts in different X-ray bands against the relative counts predicted by a plasma code, to estimate physical properties, such as temperature and metallicity. To fit many parameters, several bands are required. It is useful to choose bands with approximately the same number of counts, the choice of which can be derived from the spectrum.

For a more detailed analysis, full spectral fitting is necessary. We have investigated taking an intensity-binned image of the object, automatically extracting spectra for each of the bins, automating xspec to fit each spectrum, and then plotting the fitted parameters spatially. As a side issue, the data must be binned so that there are enough counts per pixel to assume Gaussian errors. Also, spatial variations of the response matrix of the detector and telescope may need to be considered.

## 4 Availability

An implementation of the adaptive binning algorithm is available written in the C++ programming language. See the URL http://www-xray.ast.cam.ac.uk/~jss/adbin/ for instructions for download.

## 5 Summary

We have presented an algorithm to adaptively bin data. We have demonstrated its use on X-ray intensity and colour images of the Perseus cluster. We also applied the algorithm to simulated data. The method has also been used to create multi-layer colour images, where the colours represent the intensity in different bands. The fractional error threshold can be varied depending on how much statistical noise is acceptable in the analysed image.

The method is not limited to binning on intensity or colour. Any sort of data may be binned, providing there is a method to compute the statistical fractional error, or weighting, of a group of pixels. The algorithm may be useful in presenting the output from numerical simulations.

The method is also not limited to two-dimensional data. We have used it already on one-dimensional cuts through images. It may also be simply extended to work on $`N`$-dimensional data sets.

## Acknowledgements

ACF and JSS thank the Royal Society and PPARC for support, respectively.


text:
# Simulated emergence of cyclic sexual-asexual reproduction

## 1 Introduction

The question of why sexual reproduction prevails in the immense majority of metazoan species still defies a complete and satisfactory answer - for a recent review we direct the reader to Reference . A number of different theories have been put forth to try to explain this puzzle, but all have so far fallen short of becoming a paradigm in the field. One of the major difficulties one finds in this research is the lack of convincing observational data with which to compare one’s theoretical effort. This is especially critical when addressing the onset of sexual reproduction itself, and in this case we can only speculate about this origin, but cannot test these speculations. A more promising line of research lies in analyzing the maintenance of sexual reproduction against the establishment of an apparently more efficient asexual variety. This efficiency stems from the absence of the necessity to generate males, which by itself alone should give the asexual variety the upper hand against a competing sexual alternative . Explaining the evolutionary choice favoring sexual reproduction in spite of its shortcomings has been the long-term goal for a number of researchers in this area.

It is not surprising that computer simulations have become a major tool in this field over the last few years. With the introduction of a number of simple simulational models, theoretical ideas began to be put to a “live” test on virtual populations. These models encompass dynamic rules that mimic the action of natural evolution in a variety of ways, and create a manageable representation of the conflict between selection and mutation, believed to be the driving evolutionary force. Among the models that have been used, we make special mention of the bit-string class, using Monte Carlo techniques . They seem particularly well suited to represent the microscopic dynamics of genetic evolution and allow for a very efficient coding on personal workstations; recent and comprehensive reviews of a popular model of this class, the Penna model, are available in the literature .

Penna’s model has in fact been recently used to address the foregoing problem. The greater genetic diversity generated by sexual reproduction and its consequences, such as survival after catastrophes  and winning when in competition against an asexual variety in the presence of genetically coupled parasitism , were already considered in its context. A recent simulation including the action of pleiotropic genes, i.e. those that have multiple effects, showed that in a region of its parameter space sexual reproduction generates larger equilibrium populations than its asexual counterparts when evolving independently . We report here on the use of a modified Penna model in a different, though closely related, problem.

The problem in question here is that of the holocyclic behaviour of some species. By this is meant the capability shown by these species of seasonally alternating their reproductive regimes. Of particular interest in this class of species are the aphids, or green flies.

Aphids are a group of about $`\mathrm{4\; 400}`$ species of small insects that feed on the phloem fluid of plants . Many species of aphids attack important agricultural crops, and are therefore of major economic importance, and display complex life cycles. They are one of the few groups of animals that undergo cyclical parthenogenesis, that is, the alternation of a varying number of asexual generations (parthenogenesis) with a single generation of sexual reproduction. All asexual generations are entirely female. Species that produce both sexual and asexual morphs are called holocyclic, as opposed to those anholocyclic ones that undergo obligate parthenogenesis. In a holocyclic life cycle both asexual and sexual morphs are produced at different times of the year. Males are present only in one part of the life cycle. They are usually produced only in the autumn by parthenogenetic females. Females reproducing sexually produce just one egg, that in temperate regions overwinters and hatches only in the spring.

The purpose of the present work is to study the emergence of an holocyclic life cycle through a simple evolutionary model, with dynamics designed to mimic the main features of the relation between the aphids and their environment. To summarize, we show that this simple model has a first order holocyclic transition dependent on one parameter. On one side of the transition there is no alternation of reproductive regime and the species is asexual; on the other, the holocyclic order is chosen by evolution as a means of insuring the survival of the species.

This paper is organized as follows: in the next section we present the model that was used to simulate an aphid-like species; our results are discussed in Section III, and we conclude in Section IV.

## 2 The Model

One of the main features of the environmental impact on the aphids is their great susceptibility to winter conditions. The majority of the population cannot survive in very cold conditions. On the other hand, aphid eggs generated by sexual reproduction resist harsh weather conditions. They overwinter and hatch only when spring comes. These are the elements to be introduced into a model that intends to test the hypothesis that the onset of holocyclic order is a direct consequence of evolution dynamics acting on cold-susceptible species.

The model we worked with is an adaptation of the Penna model for age-structured populations, to which we added the features described above. Since the phenomenon in question does not depend on the age structure, this model is not the only possible choice. Nevertheless, it has a very simple representation of the conflict between mutations and selection, which is central to any evolutionary problem, together with a built-in mechanism that decreases the survival probability with age, that other models have to introduce in an ad hoc way. We will briefly describe the Penna model in its mixed sexual-meiotic parthenogenesis version. A similar version of this model - but without the specific aphid-type dynamics - was used in Ref. .

The genome of each (diploid) organism is represented by two computer words. In each word, a bit set to one at a position (“locus”) corresponds to a deleterious mutation - a “perfect” strand would be composed solely of zeros. The reader should note that we don’t allow in this context for a beneficial mutation in the wildtype genome, that would eventually lead us to the problem of speciation. The effect of a harmful mutation may be felt by the individual at all ages equal to or above the numerical order of that locus in the word. As an example, a bit set to one at the second position of one of the bit-strings means that a harmful effect may become present in the life history of the organism to which it corresponds after it has lived for two time periods. The diploid character of the genome is related to the effectiveness of the mutations. A mutation in a position of one of the strands is felt as harmful either because of homozygose or because of dominance. For the former, a mutation must be present in both strings at the same position to be effective. The concept of dominance on the other hand relates to loci in the genome in which a mutation in just one strand is enough to make it affect the life of the organism. The life span of an individual is controlled by the amount of effective mutations active at any instant in time. This number must be smaller than a specified threshold to keep the individual alive; it dies as soon as this limit is reached.

Reproduction is modeled by the introduction of new genomes in the population. Each female becomes reproductive after having reached a minimum age, after which she generates a fixed number of offspring at the completion of each period of life. The meiotic cycle is represented by the generation of a single-stranded cell out of the diploid genome. To do so, each string of the parent genome is cut at a randomly selected position, the same for both strings, and the left part of one is combined with the right part of the other, thus generating two new combinations of the original genes. The selection of one of these completes the formation of the haploid gamete coming from the mother.

The difference between sexual and parthenogenetic reproduction appears at this stage. For the first, a male is selected in the population and undergoes the same meiotic cycle, generating a second haploid gamete out of his genome. The two gametes, one from each parent, are now combined to form the genome of the offspring. Each of its strands was formed out of a different set of genes.

For meiotic parthenogenesis, all genetic information of the offspring comes from a single parent. Its gamete is generated through the meiotic cycle described above and then cloned to compose an homozygous genome for the offspring. For both regimes, the next stage of the reproduction process is the introduction of $`M`$ independent mutations in the newly generated genetic strands. In this kind of model it is normal to consider only the possibility of harmful mutations, because of their overwhelming majority in nature. If male offspring can be born, the gender of the newborn is then randomly selected, with equal probability for each sex.

A last ingredient of the model is a logistic factor, called the Verhulst factor, which accounts for the finite carrying capacity of the environment for this particular species. It introduces a mean-field probability of death for an individual, coming from causes with dynamics not included in the model, and for computer simulations has the benefit of limiting the size of populations to be dealt with. This factor is essentially the ratio between the actual population in any time step and a parameter of the model, traditionally and perhaps improperly called the carrying capacity of the environment. Because it is not connected to the quality of the individual’s genome, the usage of this factor in a model of evolutionary dynamics has its shortcomings , and one has to make sure that the outcome of the simulations is not biased by the particular choice of strategy for its implementation.

The passage of time is represented by the reading of a new locus in the genome of each individual in the population, and the increase of its age by one. After having accounted for the selection pressure of a limiting number of effective (harmful) mutations and the random action of the Verhulst dagger, females that have reached the minimum age for reproduction generate a number of offspring. The simulation runs for a pre-specified number of time steps, at the end of which averages are taken over the population(s).

In the mixed sexual-parthenogenesis version, each female in the population has a season-dependent reproductive mode, which is also genetically inherited. Four bits in the genome represent the reproductive mode for each season. A bit set to $`0`$ ($`1`$) in a position indicates that at the season associated to that particular position this female will generate offspring through sex (meiotic parthenogenesis). This reproductive pattern is passed on to the female offspring. Mutations can also occur on this pattern, with some probability. One out of the four seasons is randomly selected and the reproduction mode for that season is switched.

The duration of the seasons is set to some number of time periods. At each time step, the season to which it corresponds is calculated. In the fall, females that reproduce by parthenogenesis can generate both male and female offspring. In the winter, the carrying capacity of the environment is reduced by some factor, called here the compression, or burdening, factor. Individuals with age zero that were generated by sexual reproduction - and are protected by the eggshell - suffer only from the normal Verhulst factor, representing the action of predators, and do not age; they only hatch in the spring. All the others, irrespective of gender, age or reproductive pattern, suffer the action of the enhanced Verhulst factor. The fertility of females also depends on the reproductive regime, and we fixed it to $`1`$ offspring per time period for sexual reproduction.

Averages are taken of the total population, number of females with each reproductive pattern - both total and only the ones mature for reproduction - for each season. These numbers are normalized to the normal carrying capacity.

The runs from which our data was collected had as common parameters $`12`$ aphid time periods for season duration, thus fixing our time scale as roughly a week per time step, a probability for reproductive pattern switching of $`0.01`$, and an equal selection of each gender for asexually-generated offspring in the fall. In all cases, the initial population had a season-independent sexual reproductive pattern, and was composed of the same number of individuals in each gender.

For the values of the standard parameters of the Penna model we used:

1) Carrying capacity: $`C=\mathrm{800\; 000}`$;

2) Threshold of harmful mutations: $`T=3`$;

3) Minimum reproduction age: $`R=10`$;

4) Mutation rate: $`M=2`$;

5) Birth rate: $`b=10`$ per female per time step for asexual reproduction;

6) Initial population: $`\mathrm{10\; 000}`$ males and $`\mathrm{10\; 000}`$ females.

## 3 Simulation results

We begin by showing the total population for each season in Figure 1. The figure for the population is normalized to the normal carrying capacity of the environment, one of the parameters of the model. Data is collected as an average over the last $`960`$ steps of a run of $`\mathrm{50\; 000}`$ Monte Carlo steps, for winter compression factors of $`10`$ and $`100`$. In both cases, the population increases from spring to fall, and decreases, as expected, in the winter season. For a compression factor large enough, the harsh conditions in this last season cause the population to vanish almost completely. The species does not disappears thanks to the remaining offspring of females that reproduced sexually during the last time step in the fall.

Figure 2 shows the fraction of mature females - those that have an age greater that the minimum reproduction age - that reproduce sexually in each season. If the winter has a mild effect on the effective carrying capacity, the evolution pressure favors reproduction through asexual meiotic parthenogenesis. Sexually-reproducing females are still present in the population, as a result of a small probability of back mutations from the prevailing asexual mode. The situation changes dramatically when the effects of the winter get stronger. If the corresponding compression factor is capable of bringing down the population in winter time to a small enough fraction of the normal carrying capacity, the reproductive regime in the fall switches to a largely predominant sexual variety. The data shown in this histogram correspond to situations that are way out of the transition region, and must be compared with Figure 1. In particular, the data shown there for the winter season determine what we see in the present histogram for the fall season.

The data for the winter season in Figure 2 deserve a short comment. Since the mature female population in winter time has all but vanished for a compression factor of $`100`$, evolution exerts no pressure whatsoever over the reproductive regime in this season. As a consequence, the distribution of females among the two regimes is random, and the resulting fraction of sexually reproducing ones fluctuates around $`0.5`$.

We proceed to characterize the transition between the anholocyclic and holocyclic regimes. The order parameter of this transition is clearly the fraction of sexually reproducing females in the fall season, which is $`0`$ in one regime and $`1`$ in the other. Figure 3 shows the behaviour of this parameter as a function of the surviving population in winter. A very sharp switch in regimes can be seen when this parameter is $`0.0049`$. Some of the runs near the transition region also had an unusually long relaxation time. These two characteristics point to the identification of a first order transition, and the long relaxation times can be associated to long-lived meta-stable states equivalent to super-cooled states in fluid transitions. A sharp transition of this kind has also been found in the study of the evolution of reproductive regimes in rapidly mutating ecologies , and may well be an universal feature of living systems.

The natural candidate in the model to act as a measure of the distance to the transition point appears at first to be the compression factor. But the simulation results show very strong fluctuations in the order parameter when this compression parameter is in the range $`10.5`$ \- $`12.5`$. A closer examination of the outcome of these simulations clarifies this situation. In Figure 4 the final value of the population is shown as a function of this compression factor. Results are shown for $`10`$ different runs for each value of this parameter, lasting $`\mathrm{200\; 000}`$ Monte Carlo steps each. The dispersion of the final population for a compression factor in the above range causes an end result that can be in either side of the transition. Although not directly accessible as an input parameter to the simulations, the real control variable in this case is the surviving population in the winter, and Figure 3 is an eloquent support for this claim. Also shown in Figure 4 are the average values for the population in winter as a function of the compression factor, with error bars representing the standard deviation. The horizontal dotted line stands for the value of the population at the transition point, $`0.0049`$.

## 4 Conclusions

As was already shown in Ref. , in a stable environment and when no pleiotropic genes are present, the advantage of not having to carry males gives asexual reproduction the upper hand in a direct competition against sex. This is seen by the small fraction of females that carry a sexual reproductive pattern in spring and summer shown in Figure 1. And even when their offspring are given the advantage of the eggshell protection in winter time, the switching of the dominant mode of reproduction in the fall only occurs for a very small surviving population in the winter. For large values of this surviving population, the pattern of reproduction is asexual and season-independent, switching to a cyclic asexual-sexual season-dependent one when the surviving population is smaller than some threshold. The establishment of an holocyclic order in the population derives thus from the advantages of overwintering under eggshell protection, and has no direct relation to the genetic diversity generated by sex. We could identify it as a first-order process driven by the surviving population; this claim is supported by the sharp transition observed when the fraction of sexually reproducing females in the fall is considered as an order parameter for this transition.

## Acknowledgments

We are greatly indebted to F. Delmotte, who introduced us to specialized literature on the biology of aphids, and to our referee for a number of suggestions that improved the clarity of our presentation and updated our references. J.S.S.M.’s work is supported as a Visiting Fellow by CIRES, University of Colorado at Boulder, and A.R. acknowledges financial support by the Brazilian agency CNPq.


text:
# Sets of singular restrictions of symplectic forms

## 1 Introduction

Let $`V`$ be an $`n`$-dimensional vector space over some field $`F`$. We start with the following well-known definitions.

###### Definition 1

For an automorphism $`\sigma Aut(F)`$ and other vector space $`V^{}`$ over the field $`F`$ a mapping $`f:VV^{}`$ is called $`\sigma `$-linear if

$$f(x+y)=f(x)+f(y)x,yV$$

and

$$f(ax)=\sigma (a)f(x)xV,aF.$$

In the case when

$$Kerf=0\text{ and }f(V)=V^{}$$

$`f`$ is a bijection of $`V`$ onto $`V^{}`$ and we say that $`f`$ is a $`\sigma `$-linear isomorphism of $`V`$ onto $`V^{}`$. Moreover, if $`V=V^{}`$ then the isomorphism $`f`$ will be called a $`\sigma `$-linear transformation of $`V`$. We shall often exploited the term ”linear” in plase of $`\sigma `$-linear.

###### Definition 2

Let $`\sigma _1,\sigma _2Aut(F)`$. A mapping

$$\mathrm{\Omega }:V\times VF$$

is a $`(\sigma _1,\sigma _2)`$-bilinear form on $`V`$ if

$$\mathrm{\Omega }(,x)\text{ and }\mathrm{\Omega }(x,)$$

are $`\sigma _1`$-linear and $`\sigma _2`$-linear functionals on $`V`$. In what follows we shall often say that a form is ”bilinear”; the full name will be used only for some special cases.

###### Definition 3

A bilinear form $`\mathrm{\Omega }`$ defined on $`V`$ is called symplectic if the equality $`\mathrm{\Omega }(x,x)=0`$ holds for each $`xV`$.

###### Remark 1

Recall that a bilinear form $`\mathrm{\Omega }`$ on $`V`$ is called skew-symmetric if

$$\mathrm{\Omega }(x,y)=\mathrm{\Omega }(y,x)x,yV.$$

Each symplectic form is skew-symmetric. The inverse statement holds only if the characterystic of the field $`F`$ is not equal to $`2`$. In the case when $`charF=2`$ there exists a skew-symmetric forn which is not symplectic.

###### Remark 2

An immediate verification shows that for each skew-symmetric $`(\sigma _1,\sigma _2)`$-bilinear form $`\mathrm{\Omega }`$ we have $`\sigma _1=\sigma _2`$.

Now assume that our symplectic form $`\mathrm{\Omega }`$ is non-singular (then the number $`n`$ is even). Denote by $`𝔾_k^n`$ the Grassmannian manifold of $`k`$-dimensional planes contained in $`V`$ and consider the set $`S_k^n(\mathrm{\Omega })`$ of all planes $`s𝔾_k^n`$ such that the restriction of the form $`\mathrm{\Omega }`$ onto $`s`$ is singular.

###### Remark 3

The forms $`a\mathrm{\Omega }`$ and $`\sigma \mathrm{\Omega }`$ are symplectic for each $`aF`$ and $`\sigma Aut(F)`$. Moreover,

$$S_k^n(\mathrm{\Omega })=S_k^n(a\mathrm{\Omega })=S_k^n(\sigma \mathrm{\Omega }).$$

Any two non-singular symplectic form $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ on $`V`$ are equivalent; i.e. there exists a $`\sigma `$-linear transformation $`f`$ of $`V`$ such that

$$\mathrm{\Omega }_2(x,y)=\mathrm{\Omega }_1(f(x),f(y))x,yV.$$

This implies that the set $`S_k^n(\mathrm{\Omega }_2)`$ can be transfered to $`S_k^n(\mathrm{\Omega }_1)`$ by the transformation of $`𝔾_k^n`$ induced by $`f`$.

###### Remark 4

Each $`\sigma `$-linear transformation of $`V`$ ($`\sigma Aut(V)`$) induces some transformation of the Grassmannian manifold $`𝔾_k^n`$. This fact will be often exploited in what follows.

Any symplectic form defined on an odd-dimensional vector space is singular. Therefore, the sets $`S_k^n(\mathrm{\Omega })`$ and $`𝔾_k^n`$ are coincident if the number $`k`$ is odd. For the case when $`k`$ is even the structure of $`S_k^n(\mathrm{\Omega })`$ is more complicated. However, for the cases $`k=2,n2`$ there exists a simple geometrical characteristic of this set. We shall study it above.

Let $`s𝔾_m^n`$ and $`mk`$. Denote by $`𝔾_k^n(s)`$ the set of all planes $`l𝔾_k^n`$ satisfying the condition

$$ls\text{ if }m>k$$

or

$$sl\text{ if }m<k.$$

It must be pointed out that $`𝔾_k^n(s)`$ can be considered as some Grassmannian manifold.

###### Definition 4

Let $`X𝔾_k^n`$ and $`k=2`$ or $`n2`$. In the case $`k=n2`$ we say that $`X`$ satisfies condition $`𝒮`$ if for each plane $`s𝔾_{n1}^n`$ there exists a line $`F(s)`$ contained in $`s`$ and such that

$$𝔾_k^n(s)X=𝔾_k^n(s)𝔾_k^n(F(s));$$
(1)
in other words, the set $`𝔾_k^n(s)X`$ consists of all planes containing some line. In the case $`k=2`$ the set $`X`$ satisfies condition $`𝒮`$ if for each line $`s𝔾_1^n`$ there exists a plane $`F(s)𝔾_{n1}^n`$ containing $`s`$ and such that equality (1) holds true.

It is not difficult to prove (Section 3) that the set $`S_k^n(\mathrm{\Omega })`$ ($`k=2,n2`$) satisfies condition $`𝒮`$. This paper is devoted to prove the inverse statement.

###### Theorem 1

For a set $`X𝔾_k^n`$, $`k=2`$ or $`n2`$, the following two conditions are equivalent:

1. $`X`$ satisfies condition $`𝒮`$,
2. there exists a non-singular symplectic form $`\mathrm{\Omega }`$ on $`V`$ such that $`X=S_k^n(\mathrm{\Omega })`$.

For the case when $`k2,n2`$ (clearly, the number $`k`$ is even) the similar characteristic of the set $`S_k^n(\mathrm{\Omega })`$ is not found yet.

## 2 Bijections of $`𝔾_k^n`$ onto $`𝔾_{nk}^n`$ defined by bilinear forms

Let $`\mathrm{\Omega }`$ be a non-singular bilinear form on $`V`$ and $`s𝔾_k^n`$. Denote by $`s_\mathrm{\Omega }^{}`$ the $`\mathrm{\Omega }`$-orthogonal complement to $`s`$; i.e.

$$s_\mathrm{\Omega }^{}=\{yV|\mathrm{\Omega }(x,y)=0xs\}.$$

Then $`s_\mathrm{\Omega }^{}𝔾_{nk}^n`$ and the form $`\mathrm{\Omega }`$ defines the mapping of

$$F_{knk}^n(\mathrm{\Omega }):𝔾_k^n𝔾_{nk}^n$$

tranferring each plane to the $`\mathrm{\Omega }`$-orthogonal complement. It is a bijection onto.

We shall use the following properties of bijections defined by non-singular bilinear forms.

1. For each non-singular form $`\mathrm{\Omega }`$ on $`V`$ there exists a non-singular form $`\mathrm{\Omega }^{}`$ such that

$$(F_{knk}^n(\mathrm{\Omega }))^1=F_{nkk}^n(\mathrm{\Omega }^{}).$$

Define

$$\mathrm{\Omega }^{}(x,y)=\mathrm{\Omega }(y,x)x,yV.$$

It is trivial that the form $`\mathrm{\Omega }^{}`$ is non-singular and satisfies the required condition.
2. The equality

$$(F_{knk}^n(\mathrm{\Omega }))^1=F_{nkk}^n(\mathrm{\Omega })$$

holds for some non-singular bilinear forn $`\mathrm{\Omega }`$ if and only if this form is reflexive; i.e. the following two conditions

$$\mathrm{\Omega }(x,y)=0\text{ and }\mathrm{\Omega }(y,x)=0$$

are equivalent for any $`x,yV`$.
3. For any two non-singular bilinear forms $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ on $`V`$ there exist transformations $`f`$ and $`g`$ of $`𝔾_k^n`$ and $`𝔾_{nk}^n`$ (respectively) induced by linear transformations of $`V`$ and such that

$$F_{knk}^n(\mathrm{\Omega }_2)=F_{knk}^n(\mathrm{\Omega }_1)f=gF_{knk}^n(\mathrm{\Omega }_1)$$

(see ).
4. Properties A and C imply that for any two non-singular bilinear forms $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ on $`V`$ the composition

$$F_{nkk}^n(\mathrm{\Omega }_2)F_{knk}^n(\mathrm{\Omega }_1)$$

is a transformation of $`𝔾_k^n`$ induced by some $`\sigma `$-linear transformation of $`V`$.

## 3 Proof of implication $`(\mathrm{ii})(\mathrm{i})`$

Let $`\mathrm{\Omega }`$ be a non-singular symplectic form on $`V`$. It is trivial that a planes $`l𝔾_k^n`$ belongs to the set $`S_k^n(\mathrm{\Omega })`$ if and only if

$$dimll_\mathrm{\Omega }^{}1.$$
(2)
This implies the following inclusion

$$𝔾_{n2}^n(s)𝔾_{n2}^n(s_\mathrm{\Omega }^{})𝔾_{n2}^n(s)S_{n2}^n(\mathrm{\Omega })s𝔾_{n1}^n.$$

Assume that the inverse inclusion fails and there exists

$$l𝔾_{n2}^n(s)S_{n2}^n(\mathrm{\Omega })$$

which does not contain the line $`s_\mathrm{\Omega }^{}`$. The condition $`lS_{n2}^n(\mathrm{\Omega })`$ guarantees the existence of a non-zero vector $`xl`$ such that $`\mathrm{\Omega }(x,y)=0`$ for each vector $`yl`$. The plane $`s`$ is generated by the plane $`l`$ and the line $`s_\mathrm{\Omega }^{}`$. Therefore, the similar equality holds for each vector $`ys`$; i.e. $`xs_\mathrm{\Omega }^{}`$. These arguments dispruve our hypothesis and the requred inclusion is proved. In other words, we ontain the equality

$$𝔾_{n2}^n(s)S_{n2}^n(\mathrm{\Omega })=𝔾_{n2}^n(s)𝔾_{n2}^n(s_\mathrm{\Omega }^{})s𝔾_{n1}^n$$
(3)
showing that $`S_{n2}^n(\mathrm{\Omega })`$ satisfies condition $`𝒮`$.

It was noted above that a plane $`l𝔾_k^n`$ belongs to the set $`S_k^n(\mathrm{\Omega })`$ if and only if equation (2) holds. This implies that the bijection $`F_{knk}^n(\mathrm{\Omega })`$ transfers $`S_k^n(\mathrm{\Omega })`$ to $`S_{nk}^n(\mathrm{\Omega })`$. Moreover, $`F_{n\mathrm{2\hspace{0.17em}2}}^n(\mathrm{\Omega })`$ transfers equation (3) to the following equality

$$𝔾_2^n(s)S_2^n(\mathrm{\Omega })=𝔾_2^n(s)𝔾_2^n(s_\mathrm{\Omega }^{})s𝔾_1^n;$$

i.e. $`S_2^n(\mathrm{\Omega })`$ satisfies condition $`𝒮`$.

## 4 Reduction of the case $`k=2`$ to the case $`k=n2`$

Assume that Theorem 1 is proved for the case $`k=n2`$ and consiner a set $`X𝔾_2^n`$ satisfying condition $`𝒮`$. Let $`g`$ be a bijection of $`𝔾_2^n`$ onto $`𝔾_{n2}^n`$ defined by some bilinear form. An immediate verification shows that $`g(X)`$ is a subset of $`𝔾_{n2}^n`$ satisfying condition $`𝒮`$ and there exists a non-singular symplectic form $`\mathrm{\Omega }`$ such that

$$g(X)=S_{n2}^n(\mathrm{\Omega }).$$

Then

$$F_{n\mathrm{2\hspace{0.17em}2}}^n(\mathrm{\Omega })g(X)=S_2^n(\mathrm{\Omega }).$$

It was noted above that $`F_{n\mathrm{2\hspace{0.17em}2}}^n(\mathrm{\Omega })g`$ is a transfomation of $`𝔾_2^n`$ induced by a $`\sigma `$-linear trandformation $`f`$ of $`V`$ (property D). This implies that

$$X=S_2^n(\mathrm{\Omega }^{}),$$

where $`\mathrm{\Omega }^{}`$ is the non-singular symplectic form defined by the condition

$$\mathrm{\Omega }^{}(x,y)=\mathrm{\Omega }(f(x),f(y))x,yV.$$

## 5 Proof of implication $`(\mathrm{i})(\mathrm{ii})`$

Each set $`X𝔾_{n2}^n`$ satisfying condition $`𝒮`$ is defined by some mapping

$$F:𝔾_{n1}^n𝔾_1^n.$$

###### Remark 5

It is trivial that for the case when $`X=S_{n2}^n(\mathrm{\Omega })`$ the mapping $`F`$ coincides with $`F_{n\mathrm{1\hspace{0.17em}1}}^n(\mathrm{\Omega })`$.

We state that the mapping $`F`$ satisfies the following two conditions:

1. each plane $`s𝔾_{n1}^n`$ contains the line $`F(s)`$;
2. for any two planes $`s_1,s_2𝔾_{n1}^n`$ the line $`F(s_1)`$ is contained in $`s_2`$ if and only if $`F(s_2)`$ is contained in $`s_1`$.

First condition is trivial.

To prove second condition consider the plane

$$l=s_1s_2𝔾_{n2}^n.$$

If $`lX`$ then the lines $`F(s_1)`$ and $`F(s_2)`$ is not contained in $`l`$. Then condition (1) guarantees that $`F(s_1)`$ is not contained in $`s_2`$ and $`F(s_2)`$ is not contained in $`s_1`$. In the case when $`lX`$ the plane $`l`$ contains $`F(s_1)`$ and $`F(s_2)`$. The fulfilment of condition (2) is proved.

Let $`\mathrm{\Omega }`$ be a non-singular bilinear form on $`V`$. Recall that a vector $`yV`$ is called $`\mathrm{\Omega }`$-orthogonal to a vector $`xV`$ if $`\mathrm{\Omega }(x,y)=0`$. If our form is reflexive then $`y`$ is $`\mathrm{\Omega }`$-orthogonal to $`x`$ if and only if $`x`$ is $`\mathrm{\Omega }`$-orthogonal to $`y`$. In this case we can say that the vectors $`x`$ and $`y`$ are $`\mathrm{\Omega }`$-orthogonal.

Now fix some reflexive non-singular bilinear form $`\mathrm{\Omega }`$ on $`V`$ and consider the mapping

$$f=FF_{1n1}^n(\mathrm{\Omega }):𝔾_1^n𝔾_1^n.$$

It is not difficult to see that it satisfies the next conditions:

1. for each $`l𝔾_1^n`$ the lines $`l`$ and $`f(l)`$ are $`\mathrm{\Omega }`$-orthogonal;
2. for any two lines $`t_1,t_2𝔾_1^n`$ the line $`t_1`$ and $`f(t_2)`$ are $`\mathrm{\Omega }`$-orthogonal if and only if $`t_2`$ and $`f(t_1)`$ are $`\mathrm{\Omega }`$-orthogonal.

We want to show that the mapping $`f`$ is induced by some $`\sigma `$-linear transformation of $`V`$. Then

$$F=f(F_{1n1}^n(\mathrm{\Omega }))^1=fF_{n\mathrm{1\hspace{0.17em}1}}^n(\mathrm{\Omega })$$

(the form $`\mathrm{\Omega }`$ is reflexive) and property C implies the existence of a non-singular bilinear form $`\mathrm{\Omega }^{}`$ on $`V`$ such that

$$F=F_{n\mathrm{1\hspace{0.17em}1}}^n(\mathrm{\Omega }^{}).$$

Condition (1) guarantees that any plane $`s𝔾_{n1}^n`$ contains the $`\mathrm{\Omega }^{}`$-orthogonal complement $`s_\mathrm{\Omega }^{}^{}`$. It is not difficult to see that each form satisfying this condition is symplectic and the required statement will be proved.

First of all show that the mapping $`f`$ is surjective. Let $`l𝔾_1^n`$ and $`l_1,\mathrm{},l_{n1}`$ be linearly independent lines generating the plane $`l_\mathrm{\Omega }^{}`$. Then the lines $`f(l_1),\mathrm{},f(l_{n1})`$ generate a plane $`s`$ the dimension of which is not greater than $`n1`$. Consider a line $`l^{}`$ contained in $`s_\mathrm{\Omega }^{}`$. Condition (2’) shows that the line $`f(l^{})`$ is $`\mathrm{\Omega }`$-orthogonal to the lines $`l_1,\mathrm{},l_{n1}`$; therefore, $`f(l^{})=l`$.

We state that for each plane $`s𝔾_{n1}^n`$ there exist planes $`s_1,s_2𝔾_{n1}^n`$ such that

$$f(𝔾_1^n(s))=𝔾_1^n(s_1)$$

and

$$f^1(𝔾_1^n(s))=𝔾_1^n(s_2).$$

Let $`l𝔾_1^n`$ be the $`\mathrm{\Omega }`$-orthogonal complement to the plane $`s`$. The mapping $`f`$ is a surjection and there exists a line $`l_1𝔾_1^n`$ such that $`f(l_1)=l`$. Denote by $`s_1`$ the $`\mathrm{\Omega }`$-orthogonal complement to $`l_1`$. After that denote by $`s_2`$ the orthogonal complement to the line $`f(l)`$. Then condition (2’) implies the fulfilment of the requred equalities.

Now assume that for some line $`l𝔾_1^n`$ there exist two lines $`l_1,l_2𝔾_1^n`$ such that

$$f(l_i)=li=1,2.$$

The previous arguments show that

$$f(𝔾_1^n(l_\mathrm{\Omega }^{}))=𝔾_1^n((l_1)_\mathrm{\Omega }^{})$$

and

$$f(𝔾_1^n(l_\mathrm{\Omega }^{}))=𝔾_1^n((l_2)_\mathrm{\Omega }^{}).$$

Then the planes $`(l_1)_\mathrm{\Omega }^{}`$ and $`(l_2)_\mathrm{\Omega }^{}`$ are coincident; i.e. $`l_1=l_2`$. In other words, we have proved that the mapping $`f`$ is a injective.

Therefore, $`f`$ is a bijection onto. Moreover, $`f`$ and $`f^1`$ transfer each collection of linearly independent lines to a collection of linearly independent lines. The Fundamental Theorem of Projective Geometry \[1–3\] states that in this case $`f`$ is induced by some $`\sigma `$-linear transformation of $`V`$.

1. \] E. Artin, Geometric Algebra, Interscience, New York, 1957.
2. \] J. Dieudonne, La Geometrie des groupes classiques, Springer – Verlag, Berlin – New York, 1971.
3. \] O. T. O’Meara, Lectures on linear groups, Providence, Rhode Island, 1974.
4. \] M. A. Pankov, Irregular subsets of the Grassmannian manifolds, 1999 (book, math.AT/9910081).