text:
# On the size difference between red and blue globular clusters 1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract No. NAS5-26555.

## 1 Introduction

Since the work of Kinman (1959) and Zinn (1985), it has been known that globular clusters in the Milky Way can be divided into (at least) two sub-populations with distinct kinematical and chemical properties. Undoubtedly one of the most important developments in research on globular clusters (GCs) within the last decade is the discovery that similar substructure is seen in the color (and hence, presumably, metallicity) distributions of GC systems in many *early-type* galaxies (Zepf & Ashman, 1993; Gebhardt & Kissler-Patig, 1999; Kundu & Whitmore, 2001; Larsen et al., 2001). There is increasing evidence that GC sub-populations in early-type galaxies share many properties with those in spirals (Forbes, Brodie, & Larsen, 2001), and that the GCs in different galaxy types may be closely related. Characterizing and understanding the properties (spatial and metallicity distributions, kinematics etc.) of GC sub-populations is currently a subject of much investigation and holds the promise of revealing important information about events in the evolution of their host galaxies.

One piece of the puzzle is to establish just how similar are GCs in different galaxies. The Hubble Space Telescope can resolve the spatial profiles of globular clusters well beyond the Local Group, although careful modeling of the undersampled point spread function (PSF) of the WFPC2 camera is necessary in order to derive reliable size information for typical extragalactic GCs. Kundu & Whitmore (1998) were among the first to measure the sizes of extragalactic GCs and found that GCs in the lenticular galaxy NGC 3115 had half-light radii of about 2 pc, similar to or perhaps slightly smaller than those in the Milky Way. When measuring sizes for the blue (metal-poor) and red (metal-rich) GC subpopulations separately, Kundu & Whitmore (1998) noted a size difference of about 20% with the red GCs being systematically smaller. A similar size difference was found between red and blue GCs in NGC 4486 (Kundu et al., 1999). Subsequently, size differences between GC subpopulations have been found in many other galaxies including NGC 4472 (Puzia et al., 1999), NGC 4594 (the “Sombrero”; Larsen, Forbes & Brodie, 2001), M31 (Barmby, Holland, & Huchra, 2002) and other early-type galaxies (Larsen et al., 2001). However, it should also be noted that Harris et al. (2002) found *no* size-color relation for a sample of 27 GCs in NGC 5128.

Although the results on NGC 5128 suggest that the size difference may not be *universal*, there seems to be compelling evidence that it is at least a wide-spread phenomenon, observed in spirals as well as in early-type galaxies. Understanding the origin of this size difference is a high priority, since any intrinsic correlation between the sizes and metallicities of star clusters might hold clues to their formation mechanisms. However, before one attempts to explain the size difference in terms of e.g. the properties of the proto-cluster clouds, variations in external pressure etc., other more straight-forward explanations need to be ruled out.

One possibility is that the observed size difference is a projection effect, resulting from a correlation between GC size ($`r`$) and Galactocentric distance ($`R`$) combined with different radial distributions of the GC subpopulations (throughout this paper, we will adopt the following convention: a small $`r`$ refers to the size (radius) of an individual globular cluster, while a capital $`R`$ refers to the distance of that cluster from the galaxy center. We will use subscripts ($`R_{\mathrm{proj}}`$, $`R_{3\mathrm{D}}`$) to specify projected and 3-D quantities). van den Bergh, Morbey, & Pazder (1991) found that the size-$`R`$ relation in the Milky Way can be approximated by a square-root relation, $`r\sqrt{R}`$, and they also noted that similar relations exist in M31, NGC 5128 and the LMC. However, because of the limited spatial coverage of HST, little is known about size-$`R`$ trends in other more distant galaxies.

Concerning the wide-field spatial distributions of extragalactic GC systems, most of the available information comes from ground-based studies which generally provide the only way to obtain complete coverage at large radii. In addition to the Milky Way and M31, galaxies where the radial distributions of GC sub-populations have been studied include NGC 1380 (Kissler-Patig et al., 1997), NGC 1399 (Dirsch et al., 2003), NGC 4472 (Lee, Kim, & Geisler, 1998) and NGC 4486 (Côté et al., 2001). In all of these cases, the red GC sub-population is more centrally concentrated than the blue one. Therefore, red GCs observed at a given *projected* radius will on average be physically closer to the galaxy center than the blue ones. Clearly, there is reason to suspect that projection effects might account for at least some of the differences in the mean sizes of GC samples. Our goal in this paper is to investigate and quantify this idea in more detail.

## 2 Observations: sizes of GC subpopulations

While size information is now available for many extragalactic GC systems, the spatial coverage is in most cases limited to a single WFPC2 pointing, usually centered on the galaxy. Therefore, information about variations with galactocentric distance is limited. Even with several pointings, the surface density of GCs drops off rapidly at large radii, and in all but the richest systems it is difficult to obtain reliable statistics based on the small number of clusters contained within an off-center WFPC2 field. In practice, then, the number of suitable datasets that can serve as a basis for our analysis is limited. We have selected three well-studied galaxies with multiple HST/WFPC2 pointings and rich GC systems (NGC 4472, NGC 4486 and NGC 4594). The WFPC2 data are the same as those used in Larsen et al. (2001) and we refer to that paper for details about the data reduction. Briefly, the globular cluster sizes were measured by convolving a series of analytic models with the WFPC2 PSF, adjusting the FWHM of the model until the best match to the observed profile was obtained. For this purpose, we used the ISHAPE code described in Larsen (1999). The main systematic errors involved in this process are related to che choice of model PSF and are difficult to quantify, but comparison of size measurements on different frames have indicated that this method yields FWHM values consistent to about 0.1 WF pixel for individual clusters (Larsen, Forbes & Brodie, 2001; Larsen et al., 2001). In addition, we include data for the Milky Way from Harris (1996) and M31 (size information kindly provided by P. Barmby).

As an illustration of the typical size differences between GC sub-populations, Fig. 1 shows the size distributions for GCs in each of the WFPC2 pointings for the three early-type galaxies. Here we are not concerned with the absolute values of the GC sizes so for the HST data we simply use the FWHM of the cluster profiles, measured in pixels and corrected for the WFPC2 PSF. To avoid any systematic differences due to the different pixel scales of the PC and WF chips, we only use sizes measured on the WF chips. We divide between “red” and “blue” clusters at $`VI=1.05`$, corresponding roughly to a metallicity of $`\text{[Fe/H]}1`$. Table 1 lists the mean sizes of red and blue globular clusters in each of the WFPC2 fields, along with the probability $`P`$ (from a Kolmogorov-Smirnoff test) that the size distributions are drawn from the same parent distribution. For all of the *central* WFPC2 pointings, Table 1 and Fig. 1 confirm a significant difference between the size distributions of blue and red GCs ($`P<0.03`$), amounting to a difference of 15–30% in the mean sizes. However, at larger radii the situation is less clear: While the red GCs are smaller than the blue GCs in all of the NGC 4486 pointings except NGC 4486-E, the difference is only significant in the three innermost ones. In NGC 4472, the two outer pointings (B and C) cover roughly the same radial range, but a size difference is seen in only one of them and is only significant at the $``$89% confidence level. In NGC 4594, only the central pointing has enough clusters to obtain useful information about size distributions. In this central pointing there is a $`15`$% difference between the mean sizes of red and blue GCs and the distributions are different at the 97% confidence level.

The size distributions for Milky Way and M31 GCs are shown in Fig. 2. In the Milky Way we have adopted a division at $`\text{[Fe/H]}=1.0`$, while the classification of clusters in M31 as metal-poor or metal-rich was based on a variety of photometric and spectroscopic criteria (Barmby et al., 2000). As previously noted (Larsen et al., 2001; Barmby, Holland, & Huchra, 2002), metal-rich and metal-poor GCs in the Milky Way and M31 display a size difference reminiscent of that observed in the early-type galaxies, but at least in the Milky Way it is likely that the difference is at least partly due to the different radial distributions of the two GC subpopulations. In Larsen et al. (2001) we tried to reduce this problem by comparing clusters in smaller radial bins, but the modest number of GCs in the Milky Way limits the extent to which such sub-division is feasible. In Sec. 3.2 below we return to this issue.

## 3 Modeling projection effects

### 3.1 GC size versus luminosity

Before we proceed, we briefly discuss one other effect which could potentially affect the comparison of mean sizes for GC subpopulations. Because the mass-to-light ratio of a single-aged stellar population (such as a star cluster) is a function of metallicity, the mass limit of a magnitude-limited sample of GCs is, in principle, also metallicity-dependent. Thus, a size-mass relation could introduce spurious size differences between magnitude-limited samples of different metallicities. Furthermore, age differences between GC sub-populations could also lead to different mass-to-light ratios. However, in the Milky Way there is a striking *absence* of any size-luminosity relation for GCs (van den Bergh, Morbey, & Pazder, 1991), and the same seems to be the case in early-type galaxies (Kundu & Whitmore, 2001) and even for young clusters in mergers (Zepf et al., 1999). This suggests that the magnitude (or mass) limit of the sample is largely irrelevant when comparing cluster sizes.

In Fig. 3 we plot the size measurements for globular clusters in NGC 4472, NGC 4486 and NGC 4594 as a function of their magnitude. The lines are least squares linear fits to the data, and they all have very small slopes that differ from 0 by less than $`1\sigma `$. Again, this suggests that the size difference cannot be explained as a result of different M/L ratios of the GC sub-populations. We note in passing that this lack of a size-luminosity relation is a very interesting and puzzling result in its own right, with the obvious implication that the mean *densities* of stellar clusters are directly proportional to their masses.

### 3.2 Revisiting the Milky Way

In order to calculate the projected size distributions of GC sub-populations, two basic ingredients are needed: 1) a relation between GC size and galactocentric distance, assumed to be *the same* for both (or all) GC subpopulations, and 2) density profiles, assumed to be *different* for each population. The size-$`R`$ relation is not easily constrained for extragalactic GC systems, partly because the GC sizes are close to the resolution limit of WFPC2 and small differences in telescope tracking, focusing etc. could introduce spurious effects when comparing measurements from different pointings. A more fundamental problem is that only the projected size-$`R`$ relations can be directly observed, while we need the run of GC size versus physical (3-D) galactocentric distance. For the following discussion we will therefore start out by *assuming* a relation of the form

$$r=a\sqrt{R_{3\mathrm{D}}}+b,$$
(1)
i.e. similar to that observed in the Milky Way, but note that we have allowed for addition of a constant term to the simple square-root law suggested by van den Bergh, Morbey, & Pazder (1991). We emphasize that we do not intend to imply any physical significance of this particular relation but simply use it as a convenient fitting function. Note, also, that any intrinsic scatter in the relation (which might depend on $`R_{3\mathrm{D}}`$ in non-trivial ways) is ignored.

Fig. 4 shows $`r_{\mathrm{eff}}`$ versus $`R_{3\mathrm{D}}`$ for Milky Way GCs. The upper panel shows a least-squares fit of Eq. (1) to all clusters with size information. The best fit is

$$r_{\mathrm{eff}}=(0.87\pm 0.16\mathrm{pc})\sqrt{R_{3\mathrm{D}}/3\mathrm{k}\mathrm{p}\mathrm{c}}+(1.79\pm 0.26\mathrm{pc}),$$
(2)
where $`r_{\mathrm{eff}}`$ is the effective (half-light) radius and we have normalized the relation to a scale radius of 3 kpc (roughly the median galactocentric distance of the metal-rich GCs). Of course, choosing a different scale radius is equivalent to changing the coefficient $`a`$. The two lower panels show the metal-poor (“blue”) and metal-rich GCs separately. Especially for the metal-rich GCs, the slope of the size-$`R`$ relation is very poorly constrained, but to check whether any significant size differences exist between the two GC populations we adopt the same slope as for the combined sample and just fit the constant term. For the metal-poor and metal-rich GCs this term is $`b=1.76\pm 0.13`$ pc and $`b=1.86\pm 0.23`$ pc, respectively. Formally, the constant term is now slightly larger for the *metal-rich* clusters, contrary to the trend seen in the projected size distributions, but the difference is not statistically significant. In practice, then, there is no evidence for a significant size difference between the two Milky Way GC populations when the radial dependence is averaged out. We therefore conclude that the size difference in Fig. 2, at least for the Milky Way, can be fully explained as a result of the different radial distributions of the two GC subpopulations.

### 3.3 Extragalactic GC systems

While ground-based data offer superior coverage at large radii compared to most HST datasets, completeness issues are more severe near the center. However, in order to model the GC size distributions it is important to know the behavior of the spatial distributions of GCs near the center, as well as at larger distances. Of the HST datasets used here, NGC 4486 has the best radial coverage, extending out to 9$`\mathrm{}`$. In Fig. 5 we show the GC surface densities as a function of projected distance from the galaxy center (in arcmin) for the two GC subpopulations. It is evident that the red GCs are much more centrally concentrated. We have adopted a magnitude limit of $`V=24`$, which is 1.5 mag brighter than the 50% completeness limit even in the central pointing (Larsen et al., 2001), so we do not expect completeness effects to be significant except perhaps in the central 20$`\mathrm{}`$–30$`\mathrm{}`$. Also superimposed on the plots are de Vaucouleurs $`R^{1/4}`$ laws and King (1966) profiles. For the King profiles we have arbitrarily fixed the concentration parameter at $`R_{\mathrm{tidal}}/R_{\mathrm{core}}=20`$. The two best-fitting de Vaucouleurs profiles have effective radii of $`2\mathrm{}`$ (red clusters) and $`24\mathrm{}`$ (blue clusters), although these numbers are rather sensitive to the background correction. The two King profiles have effective radii of $`2\stackrel{}{\mathrm{.}}05`$ and $`4\stackrel{}{\mathrm{.}}56`$, a much smaller difference than for the de Vaucouleurs profiles, but note that our radial coverage only extends out to about 1/3 of the half-light radius of the de Vaucouleurs fit to the blue GC distribution.

In the outer, poorly constrained regions, the King profiles decline more rapidly than the $`R^{1/4}`$ law, resulting in the very different half-light radii for the two types of profiles. In fact, a fundamental difference between the King and de Vaucouleurs profiles is that the former has a well-defined tidal radius, while the latter does not. However, our primary objective here is not to obtain a correct fit at very large galactocentric distances (which contribute with only a few clusters in projection), but merely to find appropriate “fitting functions” that allow us to model the spatial distributions of GCs and deproject the observed surface densities to obtain estimates of 3-D space densities. Fig. 5 suggests that the King profiles provide a fit that is at least as good as an $`R^{1/4}`$ law in the inner regions, and probably better.

One remaining issue is to convert the observed surface density profiles to space densities, $`\rho (R_{3\mathrm{D}})`$. For the de Vaucouleurs profiles we use the series expansion tabulated by Bendinelli, Ciotto & Parmeggiani (1993). For the King profiles, the space density is automatically obtained during the process of computing the King model.

### 3.4 Projected mean properties of GC sub-populations

Near the center, the King profiles flatten out whereas the de Vaucouleurs profiles diverge. This difference in behavior has important implications for the projected properties of GC systems distributed according to either of the two profiles. To illustrate this point, consider the mean 3-D distance $`R_{3\mathrm{D}}`$ of GCs observed at a given projected distance $`R_{\mathrm{proj}}`$. These are related as

$$R_{3D}_{R_{\mathrm{proj}}}=\frac{_{R_{\mathrm{proj}}}^{R_{\mathrm{max}}}R_{3\mathrm{D}}n(R_{3\mathrm{D}})𝑑R_{3\mathrm{D}}}{_{R_{\mathrm{proj}}}^{R_{\mathrm{max}}}n(R_{3\mathrm{D}})𝑑R_{3\mathrm{D}}}$$
(3)
where

$$n(R_{3D})\frac{R_{3D}\rho (R_{3D})}{\sqrt{R_{3D}^2R_{\mathrm{proj}}^2}}$$
(4)
The upper integration limit, $`R_{\mathrm{max}}`$, should in principle be set to the outer radius of the GC system. This is a poorly constrained, but relatively uncritical number as the density in the outer parts is very low. We will generally adopt $`R_{\mathrm{max}}=20R_{\mathrm{eff}}^{}{}_{R}{}^{}`$, where $`R_{\mathrm{eff}}^{}{}_{R}{}^{}`$ is the effective radius of the red GC subpopulation. This is well beyond the limit where reliable data exists for most GC systems.

Fig. 6 shows $`R_{3D}`$ as a function of $`R_{\mathrm{proj}}/R_{\mathrm{eff}}^{}{}_{R}{}^{}`$ for de Vaucouleurs and King profiles corresponding to the fits in Fig. 5. The straight dotted line in each figure is simply a 1:1 relation, plotted for reference. Fig. 6 shows that, at any $`R_{\mathrm{proj}}`$, the “blue” GCs are indeed located at significantly larger mean distances from the galaxy center than the red ones. For the de Vaucouleurs profile, most of the clusters observed near the center are also located at small physical distances from the center. For the King profiles, which flatten out near the center, the mean 3D distance from the center at $`R_{\mathrm{proj}}=0`$ is substantially larger than for the de Vaucouleurs profiles.

In a similar fashion to Eq. (3), the mean cluster size $`r`$ at $`R_{\mathrm{proj}}`$ is

$$r_{R_{\mathrm{proj}}}=\frac{_{R_{\mathrm{proj}}}^{R_{\mathrm{max}}}r(R_{3\mathrm{D}})n(R_{3\mathrm{D}})𝑑R_{3\mathrm{D}}}{_{R_{\mathrm{proj}}}^{R_{\mathrm{max}}}n(R_{3\mathrm{D}})𝑑R_{3\mathrm{D}}}$$
(5)
### 3.5 Application to the M31 GC system

In Fig. 7 we apply Eq. (5) to the M31 GC system, adopting the same size-$`R`$ relation as in the Milky Way. The top and bottom panels show GC size vs. projected galactocentric distance for metal-poor and metal-rich GCs, respectively. As noted by other authors, there is a clear trend of GC size increasing with $`R`$ for the M31 GCs (Crampton et al., 1985; Barmby, Holland, & Huchra, 2002). Note, however, that the samples of M31 GCs are likely to suffer from severe selection effects, making the spatial distributions uncertain (Barmby et al., 2000).

In each panel of Fig. 7, the solid and dashed lines represent the size-$`R`$ relations obtained by projecting the Milky Way relation (Eq. 2) using King and de Vaucouleurs GC density profiles with effective radii of 17$`\mathrm{}`$ and 30$`\mathrm{}`$. These effective radii are based on the data in Barmby et al. (2000) and are rather similar to those for the Milky Way GC system. The first impression is that the projected Milky Way relation appears to overestimate the GC sizes in M31, but this may be partly due to the fact that the Milky Way contains a number of GCs with large sizes which shift the mean relation upwards. Such clusters are not detected in M31, but could well have evaded detection especially near the center of the galaxy where the bright background from the bulge and disk makes it difficult to see low-surface brightness, extended objects.

As indicated in the Figure, the mean differences between observed and model sizes are $`0.63\pm 0.18`$ pc and $`1.06\pm 0.13`$ pc for the metal-poor and metal-rich GCs, respectively. Thus, in addition to the systematic offset between modeled and observed relations, there is a 0.4 pc difference between metal-rich and metal-poor clusters which is formally unaccounted for by simple projection of the Milky Way size-$`R`$ relation. However, there are several outlying datapoints, especially in the upper panel. The median is less sensitive to such outliers, and if we use the median instead of mean as an indicator of the difference between model and observations (labeled $`\mathrm{Obs}\mathrm{Model}_{\mathrm{MED}}`$ in Fig. 7) then the difference is reduced to 0.14 pc. In summary, it seems quite likely that the size difference between the GC sub-populations in M31 may again be explained largely as a result of differences in the radial distributions, as in the Milky Way, although it would desirable to confirm this with size measurements for a larger sample of clusters.

### 3.6 Projected GC size vs. radius trends in the HST datasets

NGC 4472 and NGC 4486 are both located in the Virgo cluster, for which we will assume a distance of 15 Mpc. NGC 4594 is somewhat closer, at a distance of about 9 Mpc (Ford et al., 1996). At 15 Mpc, 1 WF pixel corresponds to 7.3 pc. With a conversion factor of 1.48 between the FWHM and half-light radii of the King models used to fit the cluster profiles, the Milky Way relation (Eq. 2) then translates to $`\mathrm{FWHM}(0.08\mathrm{pixels})\sqrt{R_{3\mathrm{D}}/0\stackrel{}{\mathrm{.}}7}+\mathrm{\hspace{0.17em}0.17}`$ pixels (where the scale radius of 3 kpc corresponds to $`0\stackrel{}{\mathrm{.}}7`$ at the assumed distance). It is, however, not obvious just how to apply the Milky Way size-$`R_{3\mathrm{D}}`$ relation to other galaxies. In Fig. 8 we display the size measurements for red and blue GCs in NGC 4486 as a function of $`R_{\mathrm{proj}}`$. The radial range, $`10\mathrm{}`$ or 44 kpc, is roughly similar to those for the Milky Way and M31 plots in Fig. 4 and Fig. 7. In each panel, the mean size vs. $`R_{\mathrm{proj}}`$ relation obtained by projecting the Milky Way relation is shown with dashed lines, using Eq. (5) and the King and de Vaucouleurs profiles from Fig. 5. At the scale and radial range of this figure, the two model GC surface density profiles give very similar mean size vs. $`R_{\mathrm{proj}}`$ relations.

While the sizes are measured on several WFPC2 pointings and systematic differences might be present from one pointing to another, it is undeniable that the Milky Way relation applied directly provides a rather poor fit to the NGC 4486 data. Part of the mismatch might be due to systematic errors in the PSF modeling, which could shift all datapoints up or down, but application of the Milky Way relation also gives a too steep slope. There is, of course, no particular reason why the size-$`R`$ relations should be the same in the Milky Way and in NGC 4486. The solid lines in Fig. 8 show the size-$`R_{\mathrm{proj}}`$ relation obtained by projecting a size-$`R_{3\mathrm{D}}`$ relation with half the slope of that in the Milky Way and shifted downwards by 0.05 pixel. This provides a much better match to the observations, although the match remains less than satisfactory in the outer regions. Note that, instead of explicitly applying a shallower slope to the size-$`R_{3\mathrm{D}}`$ relation, we could have normalized $`R_{3\mathrm{D}}`$ to a 4 times larger scale radius. Although the physical basis for the size-$`R`$ relation is poorly understood, it is interesting to note that half-light radius of the red GC population in NGC 4486, $`2\mathrm{}`$ or 9 kpc, is also roughly 3–4 times larger than that of the metal-rich Milky Way GCs.

The other HST datasets do not provide additional constraints on any size-$`R`$ relation for GCs, as illustrated in Fig. 9 where we compare GC size vs. $`R_{\mathrm{proj}}`$ for NGC 4472, NGC 4486 and NGC 4594. For the other two galaxies the radial range is more limited and in NGC 4594 the number of clusters in the outer parts is very small (note that the NGC 4594 data cover a smaller radial range in kpc, due to the smaller distance of this galaxy). However, the behavior of GC size vs. $`R_{\mathrm{proj}}`$ in these galaxies does not appear to be significantly different from that in NGC 4486. Though the fit to the NGC 4486 data may be further improved by additional iterations on the input size-$`R_{3\mathrm{D}}`$ relation, we will adopt the half-slope, zero-point shifted relation in Fig. 8 for the following discussion, given as

$$\text{FWHM}=(0.12\mathrm{pixels})+(0.04\mathrm{pixels})\sqrt{R_{3\mathrm{D}}/0\stackrel{}{\mathrm{.}}7}$$
(6)
Fig. 10 shows the ratio of the mean sizes of “blue” and “red” GCs as a function of projected galactocentric distance, calculated according to Eq. (6) and using the GC surface density profile fits in Fig. 5. The solid and dashed lines are for the de Vaucouleurs and King profiles, respectively. In addition, the dotted line shows the ratio for King profiles for which the effective radii of the blue and red GC surface density distributions differ by a factor of 4 (instead of $`2.3`$). The projected radius $`R_{\mathrm{proj}}`$ is plotted in units of the effective radius of the red GC population. Near the center, a difference of about 12% in the mean sizes is indeed produced for the King profiles, while the size difference remains below 10% for the de Vaucouleurs profiles. At larger radii the size difference decreases regardless of the choice of model. While the details of Fig. 10 clearly depend on many poorly constrained parameters, it provides another hint that projection effects might account, at least partially, for the differences between the projected size distributions of metal-rich and metal-poor GCs.

### 3.7 Modeling the size distributions

A lot of useful information is lost by simply averaging over the size distributions and comparing mean sizes. Even if the average trends are reproduced, the detailed behavior of the observed distributions might not be properly accounted for. Instead, it is more illustrative to directly compare the predicted size *distributions* with the observations.

At a given projected distance $`R_{\mathrm{proj}}`$, the number of clusters $`N(R_{3\mathrm{D}})dR_{3\mathrm{D}}`$ in a small radial range $`dR_{3\mathrm{D}}`$ around $`R_{3\mathrm{D}}`$ is proportional to the right-hand side of Eq. (4). It follows that the size distribution, expressed as the number of clusters $`N(r)dr`$ in a small size bin $`dr`$ centered on $`r`$ can be expressed as

$$N(r)dr=\frac{\rho R_{3\mathrm{D}}}{\sqrt{R_{3\mathrm{D}}^2R_{\mathrm{proj}}^2}}\frac{dR_{3\mathrm{D}}}{dr}dr$$
(7)
where it is assumed that $`r`$ is a monotonic function of $`R_{3\mathrm{D}}`$ so that the inverse function $`R_{3\mathrm{D}}(r)`$ exists. Note also that the space density of clusters, $`\rho `$, must be evaluated at the appropriate radius, $`R_{3\mathrm{D}}(r)`$.

The distribution functions for $`R_{3\mathrm{D}}`$ (Eq. 4) at various projected radii are shown in Fig. 11 for King profiles (a) and de Vaucouleurs laws (b). Each panel shows the $`R_{3\mathrm{D}}`$ distribution functions for projected radii of $`R_{\mathrm{proj}}=0.0`$, 0.5 and 1.0. The solid lines are for models with effective radii of 1.0 (both panels). Following the fits to the NGC 4486 GC system (Fig 5), the dashed lines show King models with $`R_{\mathrm{eff}}=2.3`$ (a) and de Vaucouleurs models with $`R_{\mathrm{eff}}=10`$ (b). For both models the $`R_{3\mathrm{D}}`$ distributions are strongly peaked near $`R_{\mathrm{proj}}`$ and quite similar for models with different effective radii, except for the King models near the center. Hence, the most significant differences in the size distributions of blue and red GCs are expected near the center, and should be much more noticeable for King density profiles.

In Fig. 12 we show the projected size distributions at the same $`R_{\mathrm{proj}}`$ as in Fig. 11, using the size-$`R_{3\mathrm{D}}`$ relation for NGC 4486 from Fig. 8 (Eq. 6). As suspected, the King and de Vaucouleurs laws produce very different size distributions near the center. Even though the de Vaucouleurs law does produce a difference in the mean sizes (though smaller than for the King profiles), this is mainly due to a more extended “tail” in the size distribution of the blue clusters. For the King profiles, on the other hand, the simulated size distributions show two clearly separated peaks near the center, more akin to the observed distributions. Changing the constant term ($`b`$) in Eq. (1) simply shifts both distributions horizontally, while changing the scale factor ($`a`$) changes the width of the two distributions and their separation.

In Fig. 13 we compare the observed size distributions in NGC 4486 with our simple model. Panel (a) shows the size distributions obtained by projecting the square-root size-$`R_{3\mathrm{D}}`$ law (Eq. 6) at 15–30 discrete $`R_{\mathrm{proj}}`$ (hence the jagged appearance of some of the curves), spanning the radial range covered by each WFPC2 field. The individual contributions were weighted by the azimuthal coverage of the WFPC2 field at the corresponding $`R_{\mathrm{proj}}`$. We have not attempted to normalize the model distributions to the observed numbers of clusters in each bin, although the relative numbers of red and blue clusters do change according to the assumed (King) surface density profiles. The observed and modeled number ratios N(red)/N(blue) are listed for each for each WFPC2 pointing in Table 2 and agree within the uncertainties, except perhaps for the outermost pointing where deviations from the King profile may begin to become apparent.

While the simulated size distributions in Fig. 13 (a) do show a significant difference between red and blue GCs, the offset between the two distributions is clearly not as pronounced as for the actual observed distributions. Thus, at least for the simple square-root law, it is difficult to fully explain the size difference as a projection effect. In order to obtain a wider separation between the peaks of the two GC size distributions, we need a steeper relation between GC size and galactocentric distance in the central parts of the GC system. However, adopting a steeper slope for the square-root law would be incompatible with the relatively flat size-$`R_{3\mathrm{D}}`$ trend observed at larger radii (Fig. 9). We therefore need a curve which is steeper at small radii compared to the square-root law, but flattens out at large radii. One possibility is to use a relation consisting of two linear segments:

$$\text{FWHM}=\{\begin{array}{cc}a_1R_{3\mathrm{D}}+b_1\hfill & \hfill \text{for }R_{3\mathrm{D}}<R_1\\ a_2R_{3\mathrm{D}}+b_2\hfill & \hfill \text{for }R_{3\mathrm{D}}>R_2\end{array}$$
(8)
with $`R_2>R_1`$ and using a cubic spline to interpolate between the two linear segments in the interval $`R_1<R_{3\mathrm{D}}<R_2`$. Just piecing a relation together directly from two linear segments with an abrupt change in slope at a certain radius will not work, because the term $`dR_{3\mathrm{D}}/dr`$ in Eq. (7) will then be discontinuous, leading to a “jump” in the size distribution. However, by choosing an adequate separation between $`R_1`$ and $`R_2`$ the change in slope will be gentle enough to avoid spurious effects in the resulting size distribution. After some experimentation, we found that a reasonable match to the observed size distributions was obtained for $`a_1=0.07`$ pixels arcmin<sup>-1</sup>, $`b_1=0.05`$ pixels, $`R_1=1\stackrel{}{\mathrm{.}}5`$, $`a_2=0.005`$ pixels arcmin<sup>-1</sup>, $`b_2=0.21`$ pixels and $`R_2=4\stackrel{}{\mathrm{.}}5`$. The curve corresponding to these coefficients is compared with the square-root law in Fig. 14 and in panel (b) of Fig. 13 we show the corresponding size distributions. By construction, the match to the observed size distributions is now much better in the inner pointing (uppermost panel), though the model now appears to actually overpredict the size difference in the range $`0\stackrel{}{\mathrm{.}}10<R<3\stackrel{}{\mathrm{.}}10`$. The “tail” extending up to sizes $`>0.5`$ pixels in the observed distributions is not reproduced in any of the model distributions, but may well be due to intrinsic and/or measurement scatter in the GC size distributions, which is not taken into account in our simple model. Probably for the same reason, the simulated distributions in the outer pointings are much more peaked than the observed ones. Regardless, Fig. 13 (b) shows that projection effects can indeed produce size differences that are quite similar to the observed ones.

Finally, we caution against over-interpreting Fig. 13. It should be emphasized that the measured GC sizes are only a fraction of a WFPC pixel and any differences in mean size from one panel to another are probably only marginally significant. For example, Table 1 indicates a decrease in the mean cluster size from field NGC 4486-C to NGC 4486-D of about 0.06 pixels, which is almost certainly an artifact. The random rms error on the size measurements is about 0.1 pixel (Larsen et al., 2001) and may account for a significant fraction of the widths of the observed size distributions, although the Milky Way data suggests that some intrinsic scatter is also present.

## 4 Summary and concluding remarks

We have re-examined the size distributions of globular cluster sub-populations in three early-type galaxies with HST/WFPC2 imaging, as well as those in the Milky Way and M31. As in previous studies, we find that the blue clusters are, on average, larger than the red ones. The difference is, however, most pronounced in the central parts of the galaxies, and its statistical significance decreases strongly at larger radii.

In order to test whether the size difference might be a projection effect, resulting from a correlation between cluster size and galactocentric distance combined with different radial distributions of the GC subpopulations, we first fitted a square-root law to size vs. galactocentric distance ($`R_{3\mathrm{D}}`$) for metal-rich and metal-poor GCs in the Milky Way. Within the uncertainties, the two GC sub-populations in the Milky Way appear to follow the same size-$`R_{3\mathrm{D}}`$ relation. For extragalactic GC systems, only projected size distributions can be observed. In general, the size distributions observed at a given projected galactocentric distance will depend on the relation between GC size and 3-D galactocentric distance, as well as the 3-D density profiles of the two GC populations. These relations are poorly constrained by current data, but we have explored a limited number of analytical approximations that seem to be adequate for the NGC 4486 GC system (for which the best data are available). Specifically, we modeled the projected size distributions of GC subpopulations in NGC 4486 by assuming de Vaucouleurs $`R^{1/4}`$ and King profiles for the GC density profiles and projecting plausible three-dimensional size-$`R_{3\mathrm{D}}`$ relations for the GC subpopulations. For the de Vaucouleurs profiles we fail to reproduce the observed size differences, but for the King profiles we find that, with some fine-tuning of the input size–$`R_{3\mathrm{D}}`$ relation, we are able to produce size distributions for the red and blue GC subpopulations that are quite similar to the observations. In our model, the size difference between red and blue GCs is expected to decrease strongly beyond 1 effective radius of the GC system. This is consistent with the lack of any size-color relation for GCs in NGC 5128 noted by Harris et al. (2002), as all of the clusters in their sample are located beyond 6 kpc (compared to an effective radius of the underlying stellar halo of 5.2 kpc (Peng et al., 2002)).

Our approach clearly involves many simplifying assumptions, such as ignoring any intrinsic scatter in the GC sizes at any given galactocentric distance and using weakly constrained input relations for GC size vs. $`R_{3\mathrm{D}}`$. It nevertheless seems likely that projection effects might indeed account for a significant fraction, if not all, of the observed differences between the size distributions of GC sub-populations. Thus, at this point it may be unnecessary to resort to explanations involving different formation and/or destruction mechanisms. Since size differences have now been observed in many galaxies and the radial distributions of blue and red GCs generally appear to be different, the globular systems in those galaxies probably follow similar size-$`R`$ relations. It cannot be ruled out that intrinsic size differences exist in some cases. Indeed, we know of some cases of clusters with unusually large sizes, such as the faint Palomar-type clusters in the Milky Way halo and the “faint fuzzy” clusters in the disks of the nearby lenticular galaxies NGC 1023 and NGC 3384 (Larsen & Brodie, 2000; Brodie & Larsen, 2002). These objects have effective radii that are a factor of 3–4 larger than those of “normal” star clusters, and might have formed under different conditions.

This work has benefited from discussions and correspondence with many of our colleagues, including D. Forbes, J. Strader and J. Cohen. The presentation has also benefitted from the useful comments of an anonymous referee. We acknowledge support by National Science Foundation grants AST9900732 and AST0206139.


format: text

text:
# References

LOW FRICTION FLOWS OF LIQUIDS AT NANOPATTERNED INTERFACES

Cécile Cottin-Bizonne, Jean-Louis Barrat111Author for correspondance: barrat@dpm.univ-lyon1.fr; Tel: 33 4 72 44 85 65 Fax: 33 4 72 43 26 48,Lydéric Bocquet, Elisabeth Charlaix

Département de Physique des Matériaux, CNRS and Université de Lyon, Bâtiment Léon Brillouin, 43 Boulevard du 11 Novembre, 69622 Villeurbanne Cedex

With the recent important development of microfluidic systems, miniaturization of flow devices has become a real challenge. Microchannels, however, are characterized by a large surface to volume ratio, so that surface properties strongly affect flow resistance in submicrometric devices. We present here results showing that the concerted effect of wetting properties and surface roughness may considerably reduce friction of the fluid past the boundaries. The slippage of the fluid at the channel boundaries is shown to be drastically increased by using surfaces that are patterned at the nanometer scale. This effect occurs in the regime where the surface pattern is partially dewetted, in the spirit of the ’superhydrophobic’ effects that have been recently discovered at the macroscopic scales. Our results show for the first time that, in contrast to the common belief, surface friction may be reduced by surface roughness. They also open the possibility of a controlled realization of the ’nanobubbles’  that have long been suspected to play a role in interfacial slippage .

The nature of the boundary condition for fluid flows past solid surfaces is a subject of ancient interest , which has been revived recently by a large number of experiments and new theoretical approaches. The possibility of investigating flows at small scales in a quantitative manner, as opened by the development of nanoscale measurements (Surface force apparatus or Atomic Force Microscope) has allowed a number of experimental determinations of the ’slip length’ $`\delta `$ (see figure1) that is used to characterize this boundary condition . Optical techniques, such as fluorescence correlation or recovery methods  have also shown evidence for the existence of a nonzero slip length.

From a theoretical point of view, the parameters controlling the magnitude of the slip length are still largely unknown. At a macroscopic scale, the strength of the interaction between a solid and a liquid is most obviously characterized through the wetting behaviour. Weak interactions result in nonwetting behaviour, with large contact angles for a drop of liquid resting on the solid substrate. This characterization is purely thermodynamic, and has in principle no direct influence on nature of fluid flow past the interface. Molecular dynamics studies and mode coupling calculations , however, have shown that wettability of a perfect surface can be correlated to the magnitude of the hydrodynamic slippage. Qualitatively, a non wettable substrate is only weakly coupled to the liquid, with a depleted surface layer that can be seen as an atomic scale ’air cushion’. Therefore momentum transfer parallel to the interface is inefficient, and a large ’slip length’ results.

Experimentally, it appears that although wettability is an important parameter, different results can be obtained for substrate/liquid combinations with similar wetting properties . Another parameter of obvious importance, which may explain such variability, and has not been taken into account previously in molecular simulations, is surface roughness. In fact, it was shown by Richardson  that roughness suppresses slippage on a macroscopic scale, for any type of microscopic conditions.

Zhu’s experiments  have indeed shown that, in a situation where slippage is observed on smooth surfaces, it can be suppressed by increasing surface roughness. Other experiments carried out by Watanabe et al.  on highly water-repellent walls (i.e. walls bearing a pattern of narrow parallel grooves) resulted on the other hand on important slippage at the wall. In this work, we present a numerical study of slippage at interfaces bearing a model roughness, which takes the form of a nanoscale periodic pattern. The effect of such patterns on wetting properties has been studied extensively, both experimentally and theoretically, but their influence on dynamics has received much less attention.

The configuration considered in our study is a fluid slab confined between two parallel solid walls. The bottom wall is decorated with a periodic array of square shaped dots of height $`h`$ and width $`a`$ (see figure 2). The typical lateral size of the cell is $`L_x=L_y=t=20\sigma `$. Periodic boundary conditions in the directions $`x`$ and $`y`$ parallel to the wall are used. In our simulations, all the interactions are of the Lennard-Jones type

$$v_{ij}=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}c_{ij}\left(\frac{\sigma }{r}\right)^6\right]$$
(1)
The fluid and the solid atoms have the same molecular diameter $`\sigma `$ and interaction energies $`ϵ`$. The variation in the $`c_{ij}`$ (the index $`i,j=F,S`$ refers to the fluid or solid phase) is a convenient control parameters that can be varied to adjust the surface tensions. The solid substrate is described by atoms fixed on a the (100) plane of an FCC lattice. We have worked with two values of $`c_{FS}`$, $`c_{FS}=0.5`$ or $`c_{FS}=0.8`$, which correspond to contact angles (deduced from Young’s law) of $`\theta =137^o`$ or $`\theta =110^o`$, respectively. The simulations are carried out at constant temperature $`k_BT/ϵ`$=1. In flow experiments, the velocity component in the directions orthogonal to the flow was thermostatted, in order to avoid viscous heating within the fluid slab. All results reported are obtained within a linear response regime. Most numerical results will be given in Lennard-Jones units (L.J.u.), taking the diameter $`\sigma `$ as the unit of length and the energy $`ϵ`$ as the unit of energy.

We first briefly describe the static properties of the system, which were obtained using the following procedure. For a fixed number of liquid atoms $`N_L`$ and distance between the planes, the normal pressure is obtained from the average force on the substrates along the $`z`$ direction. This normal pressure can be varied either by changing the distance between the two walls at a fixed surface density, or by modifying the number of particles at fixed distance.

Figure 2 shows the pressure-thickness curve obtained when using the first procedure. Two branches separated by a typical ’van der Waals loop’ are clearly visible, indicating the existence of a phase transition between two possible situations. At higher normal pressures, the liquid occupies all the cell, including the grooves separating the square dots. At lower pressures, partial dewetting is observed and a composite interface is formed, the space between the dots being essentially free from liquid atoms (see figure 3). Note that the range of pressure in figure 2 is rather small, so that the properties of the bulk fluid over this range are essentially constant (the relative density change being typically less than 0.2%)

In spite of the very small sizes involved, a qualitative interpretation of the observed behaviour can easily be given in terms of macroscopic capillarity. If we consider a system at fixed normal pressure $`P_N`$, the difference in Gibbs free energies between wetted configuration (case $`A`$) and the formation of a composite interface (case $`B`$) can be written as:

$$G_BG_A=(t^2+4aha^2).(\gamma _{LV}\mathrm{cos}\theta )+(t^2a^2).\gamma _{LV}+P_N(t^2a^2)h$$
(2)
The composite interface is therefore favoured when

$$P_N<P_{composite}=\frac{\gamma _{LV}.(\mathrm{cos}\theta +1)}{h}\frac{4a\gamma _{LV}\mathrm{cos}\theta }{t^2a^2}$$
(3)
Although a quantitative agreement can hardly be expected in view of the small sizes in our system, we have checked that equation 3 correctly predicts the general trends observed in our simulation, when the height or width or the square dots are varied. It can be used, e.g., to understand the influence of an increase in the corrugation wavelength ($`t`$ and $`a`$). At fixed amplitude $`h`$, such an increase results in a decrease of $`P_{composite}`$, down to unphysically negative pressures. However dewetting may persist up to large, say micrometric, scales for adequately chosen asperity sizes .

We now consider the essential objective of our study, namely the influence of roughness on dynamical properties of the confined fluid layers. Our study involved parallel Couette flow, with the upper wall is moved with velocity $`U`$ and the lower wall is moved with a velocity $`U`$ (typically $`U=0.3`$ in reduced Lennard-Jones units). Thermostatting through velocity rescaling in the direction perpendicular to the flow is used to keep the temperature constant. For flat walls , the slip length $`\delta `$ is defined as the distance between the wall position and the depth at which the extrapolated velocity profile reaches the nominal wall velocity, $`v=U`$. In the presence of square dots, the same definition is used. Obviously the presence of dots makes the choice of the wall position somewhat arbitrary. We choose to define the wall as the position of the bottom layer of substrate.

In the following, we will discuss the particular case of a fluid-solid interaction $`c_{FS}=0.5`$, which corresponds to a contact angle $`\theta =137^o`$ on a flat substrate. For this particular value of the interaction, and in the range of pressures we have investigated, the fluid displays a moderate amount of slip at a flat interface, that can be characterized by a slip length $`\delta `$ in the range $`2025\sigma `$ (the actual value being slightly pressure dependent). Taking a molecular size of $`0.5`$ to $`1`$nm, this corresponds to a value of $`1025`$nm.

In the presence of the square dot pattern, two very different situations have to be distinguished. The first case is that of a completely wetted substrate (figure 3). A typical velocity profile for this situation is shown in figure 4. First, it is seen from this profile that the introduction of the patterned substrate does not modify the flow in the vicinity of the upper, structureless wall. This boundary can still be characterized by a slip length $`\delta =22\sigma `$, identical to what would be obtained with two flat substrates at the same pressure. In contrast, the slip is strongly suppressed at the lower wall, where it can be characterized by a value $`\delta 2\sigma `$. This value would be slightly higher ($`7\sigma `$) if the corrugation was taken at the crest of the pattern, but would still be much smaller than the value obtained for a perfectly flat surface.

A completely different result is obtained in the case where the pressure is low enough that a composite interface is formed. A typical velocity profile corresponding to this situation is shown in the right panel of figure 4. As in the wetted situation, the pattern on the lower wall does not modify the slip effect on the upper one. The slip effect on the lower wall is, on the other hand, strongly enhanced by the presence of the composite interface. Numerically, the slip length increases by a factor of about $`2.5`$, reaching a value $`\delta 57\sigma `$.

Qualitatively, the increase in the slip length should be associated with the absence of friction at the vapour-liquid interface, which can be described by a zero stress boundary condition . Hence, increasing the liquid vapour interfacial area by using a more ’spiky’ pattern should result in larger slip length. Using a smaller value of $`a`$ ($`a=4.9\sigma `$) we find indeed that the slip length can reach a value $`\delta =130\sigma `$. It is also interesting to compare these simulations results with what could be inferred from approximate hydrodynamic calculations such as that of Hocking . Hocking considered in particular flow past a composite interface, and showed that this results in a partial slip condition when the fluid filling the corrugation (in our case vapour) is of lower viscosity. Unfortunately, a direct comparison with our results is not possible, due to the use of rather different corrugation models and to the ’no-slip’ condition used by Hocking at the molecular level. We nevertheless expect that our results should stimulate further work using continuum approaches.

Our results confirm that mesoscopic roughness at the solid liquid interface can drastically modify the interfacial flow properties, in the same manner as it affects static wetting properties. Experimentally, very few results involving surfaces of controlled roughness are available, and, as mentioned above, these experiments yield contrasting results. We believe this variability may find its roots in different wetting situations realized at the mesoscale. It has also been advocated that the existence of ’nanobubbles’ at the liquid solid interface is an important factor in slippage phenomena. Our simulations, also they do not, strictly speaking, confirm the existence of such bubbles, show nevertheless that a composite interface can indeed enhance slippage considerably. In our simulations, this enhancement corresponds to an equilibrium situation, which is rapidly achieved at such small scales. Regarding this point, the situation in experiments, often done with surfaces bearing a micrometric pattern, is less clear cut since metastable trapping of bubbles (possibly made of dissolved gases) may occur into submicrometric channels. However, the essential ingredient of the effect studied here, i.e. dewetting, is known to be present on patterned surfaces, from the nanometer up to the micrometric scale. Therefore, although our simulations are performed at nanometer scales, a friction reduction is expected over a large range of length scale, up to microscale patterns. The effect is however expected to be stronger for the more ”spiky” nanometer pattern. In fact, very spectacular drops in the flow resistance of droplets have been reported  on surfaces that where decorated with a spiky nanoscale pattern, while the effect of a similar microscale pattern was less pronounced.

In summary, use of patterned surfaces, and treated to produce a ’water repellent’ like effect, appears to be a promising way towards devices that would allow flow of liquids in small size channels with very small flow resistance.

Acknowledgments It is a pleasure to thank Prof. H.A. Stone for interesting discussions. We thank the DGA for its financial support, and the PSMN (ENS-Lyon) and CDCSP (University of Lyon) for the use of their computational facilities.


text:
# Randomness Relative to Cantor Expansions

## 1 Varying Alphabets and the Cantor Expansion

Algorithmic Information Theory (see ) deals with random sequences over a finite (not necessarily binary) alphabet. A real number is random if its binary expansion is a binary random sequence; the choice of base is irrelevant (see  for various proofs).

Instead of working with a fixed alphabet we can imagine that the letters of a sequence are taken from a fixed sequence of alphabets. This construction was introduced by Cantor as a generalization of the $`b`$–ary expansion of reals. More precisely, let

$$b_1,b_2,\mathrm{}b_n,\mathrm{}$$

be a fixed infinite sequence of positive integers greater than 1. Using a point we form the finite or infinite sequence

$$0.x_1x_2\mathrm{}$$
(1)
such that $`0x_nb_n1`$, for all $`n1`$. Consider the set of rationals

$$s_1=\frac{x_1}{b_1}\text{,}s_2=\frac{x_1}{b_1}+\frac{x_2}{b_1b_2}\text{,}\mathrm{}\text{,}s_n=s_{n1}+\frac{x_n}{b_1b_2\mathrm{}b_n}\text{,}\mathrm{}$$
(2)
The above sum is bounded from above by 1,

$$0s_n\underset{i=1}{\overset{n}{}}\frac{b_i1}{b_1b_2\mathrm{}b_i}=1\frac{1}{b_1b_2\mathrm{}b_n}<1,$$

so there is a unique real number $`\alpha `$ that is the least upper bound of all partial sums (2). The sequence (1) is called the Cantor expansion of the real $`\alpha [0,1]`$.

If $`x_n=b_n1`$, for all $`n1`$, then $`s_n=11/(b_1b_2\mathrm{}b_n),`$ so $`\alpha =1`$. If $`b_n=b`$, for all $`n1`$, then the Cantor expansion becomes the classical $`b`$–ary expansion. If $`x_n=1`$ and $`b_n=n+1`$, for all $`n1`$, then $`\alpha =e`$.

The genuine strength of the Cantor expansion unfolds when various choices and interactions on different scales are considered.

The main result regarding Cantor expansions is:

###### Theorem 1

Fix an infinite sequence of scales $`b_1,b_2,\mathrm{}`$. Assume that we exclude Cantor expansions in which starting from some place after the point all the consecutive digits are $`x_n=b_n1`$. Then, every real number $`\alpha [0,1]`$ has a unique Cantor expansion (relative to $`b_1,b_2,\mathrm{}`$) and its digits are determined by the following relations:

$$\rho _1=\alpha ,x_1=b_1\rho _1,\rho _{n+1}=b_n\rho _nx_n,x_{n+1}=b_{n+1}\rho _{n+1}.$$

Consequently, if we exclude Cantor expansions in which starting from some place after the point all the consecutive digits are $`x_n=b_n1`$. then given $`\alpha [0,1]`$ there is a unique sequence $`𝐱^\alpha X^{(f)}`$ whose Cantor expansion is exactly $`\alpha `$. If $`𝐱X^{(f)}`$, then we denote by $`\alpha ^𝐱`$ the real whose Cantor digits are given by the sequence $`𝐱`$, hence $`𝐱^{\alpha ^𝐱}=𝐱`$ and $`\alpha ^{𝐱^\alpha }=\alpha `$.

For more details regarding the Cantor expansion see .

## 2 Examples

First, following  we consider the British system in which length can be measured in miles, furlongs, chains, yards, feet, hands, inches, lines. These scales relate in the following way: 1 mile = 8 furlongs = $`810`$ chains = $`81022`$ yards = $`810223`$ feet = $`8102234`$ hands = $`81022343`$ inches = $`8102234312`$ lines. Hence, the sequence of scales starts with $`b_1=10,b_2=8,b_3=10,b_4=22,b_5=3,b_6=4,b_7=3,b_8=12`$ and can be continued ad infinitum. For example, the number $`0.963(11)232(10)00\mathrm{}0\mathrm{}`$ represents a length of 9 miles, 6 furlongs, 3 chains, 11 yards, 2 feet, 3 hands, 2 inches and 10 lines.

For our second example we consider a ball in gravitational fall impinging onto a board of nails with different numbers $`b_n+1`$ of nails at different horizontal levels (here, $`n`$ stands for the $`n`$th horizontal level and $`b_n`$ is the basis corresponding to the position $`n`$). Let us assume that the layers are “sufficiently far apart” (and that there are periodic boundary conditions realizable by elastic mirrors). Then, depending on which one of the $`b_n`$ openings the ball takes, one identifies the associated number (counted from $`0`$ to $`b_n1`$) with the $`n`$th position $`x_n\{0,\mathrm{},b_n1\}`$ after the point. The resulting sequence leads to the real number whose Cantor expansion is $`0.x_1x_2\mathrm{}x_n\mathrm{}`$.

As a third example we consider a quantum correspondent to the board of nails harnessing irreducible complementarity and the randomness in the outcome of measurements on single particles. Take a quantized system with at least two complementary observables $`\widehat{A},\widehat{B}`$, each one associated with $`N`$ different outcomes $`a_i,b_j`$, $`i,j\{0,\mathrm{},N1\}`$, respectively. Notice that, in principle, $`N`$ could be a large (but finite) number. Suppose further that $`\widehat{A},\widehat{B}`$ are “maximally” complementary in the sense that measurement of $`\widehat{A}`$ totally randomizes the outcome of $`\widehat{B}`$ and vice versa (this should not be confused with optimal mutually unbiased measurements ).

A real number $`0.x_1x_2\mathrm{}x_n\mathrm{}`$ in the Cantor expansion can be constructed from successive measurements of $`\widehat{A}`$ and $`\widehat{B}`$ as follows. Since all bases $`b_n`$ used for the Cantor expansion are assumed to be bounded, choose $`N`$ to be the least common multiple of all bases $`b_n`$. Then partition the $`N`$ outcomes into even partitions, one per different base, containing as many elements as are required for associating different elements of the $`n`$th partition with numbers from the set $`\{0,\mathrm{},b_n1\}`$. Then, by measuring

$$\widehat{A},\widehat{B},\widehat{A},\widehat{B},\widehat{A},\widehat{B},\mathrm{}$$

successively, the $`n`$th position $`x_n\{0,\mathrm{},b_n1\}`$ can be identified with the number associated with the element of the partition which contains the measurement outcome.

As an example, consider the Cantor expansion of a number in the bases 2, 6, and 9. As the least common multiple is 18, we choose two observables with 18 different outcomes; e.g., angular momentum components in two perpendicular directions of a particle of total angular momentum $`\frac{9}{2}\mathrm{}`$ with outcomes in (units are in $`\mathrm{}`$)

$$\{\frac{9}{2}\text{,}4,\frac{7}{2}\text{,}\text{.}\text{.}\text{.,}+\frac{7}{2}\text{,}+4,+\frac{9}{2}\}\text{.}$$

Associate with the outcomes the set $`\{0,1,2,\mathrm{},17\}`$ and form the even partitions

$`\{\{0,1,2,3,4,5,6,7,8\},\{9,10,11,12,13,14,15,16,17\}\},`$
$`\{\{0,1,2\},\{3,4,5\},\{6,7,8\},\{9,10,11\},\{12,13,14\},\{15,16,17\}\},`$
$`\{\{0,1\},\{2,3\},\{4,5\},\{6,7\},\{8,9\},\{10,11\},\{12,13\},\{14,15\},\{16,17\}\},`$
(or any partition obtained by permutating the elements of $`\{0,1,2,\mathrm{},17\}`$) associated with the bases 2, 6, and 9, respectively.

Then, upon successive measurements of angular momentum components in the two perpendicular directions, the outcomes are translated into random digits in the bases 2, 6, and 9, accordingly.

As the above quantum example may appear “cooked up”, since the coding is based on a uniform radix $`N`$ expansion, one might consider successive measurements of the location and the velocity of a single particle. In such a case, the value $`x_n`$ is obtained by associating with it the click in a particular detector (or a range thereof) associated with spatial or momentum measurements. Any such arrangements are not very different in principle, since every measurement of a quantized system corresponds to registering a discrete event associated with a detector click .

## 3 Notation and Basic Results

We consider $`\mathrm{I}\mathrm{N}`$ to be the set of non-negative integers. The cardinality of the set $`A`$ is denoted by $`\mathrm{card}(A)`$. The base 2 logarithm is denoted by $`\mathrm{log}`$.

If $`X`$ is a set, then $`X^{}`$ denotes the free monoid (under concatenation) generated by $`X`$ with $`e`$ standing for the empty string. The length of a string $`wX^{}`$ is denoted by $`|w|`$. We consider the space $`X^\omega `$ of infinite sequences ($`\omega `$-words) over $`X`$. If $`𝐱=x_1x_2\mathrm{}x_n\mathrm{}X^\omega `$, then $`𝐱(n)=x_1x_2\mathrm{}x_n`$ is the prefix of length $`n`$ of $`𝐱`$. Strings and sequences will be denoted respectively by $`x,u,v,v,w,\mathrm{}`$ and $`𝐱,𝐲,\mathrm{}`$. For $`w,vX^{}`$ and $`𝐱X^\omega `$ let $`wv,w𝐱`$ be the concatenation between $`w`$ and $`v,𝐱`$, respectively.

By “$``$” we denote the prefix relation between strings: $`wv`$ if there is a $`v^{}`$ such that $`wv^{}=v`$. The relation “$``$” is similarly defined for $`wX^{}`$ and $`𝐱X^\omega `$: $`w𝐱`$ if there is a sequence $`𝐱^{}`$ such that $`w𝐱^{}=𝐱`$. The sets $`\mathrm{𝐩𝐫𝐞𝐟}(𝐱)=\{w:wX^{},w𝐱\}`$ and $`\mathrm{𝐩𝐫𝐞𝐟}(B)=_{𝐱B}\mathrm{𝐩𝐫𝐞𝐟}(𝐱)`$ are the languages of prefixes of $`𝐱X^\omega `$ and $`BX^\omega `$, respectively. Finally, $`wX^\omega =\{𝐱X^\omega :w\mathrm{𝐩𝐫𝐞𝐟}(𝐱)\}.`$ The sets $`(wX^\omega )_{wX^{}}`$ define the natural topology on $`X^\omega `$.

Assume now that $`X`$ is finite and has $`r`$ elements. The unbiased discrete measure on $`X`$ is the probabilistic measure $`h(A)=\mathrm{card}(A)/r`$, for every subset of $`X`$. It induces the product measure $`\mu `$ defined on all Borel subsets of $`X^\omega `$. This measure coincides with the Lebesgue measure on the unit interval, it is computable and $`\mu (wX^\omega )=r^{|w|}`$, for every $`wX^{}`$. For more details see .

In dealing with Cantor expansions we assume that the sequence of bases $`b_1,b_2,\mathrm{}b_n,\mathrm{}`$ is computable, i.e. given by a computable function $`f:\mathrm{I}\mathrm{N}\mathrm{I}\mathrm{N}\{0,1\}`$. Let $`X_i=\{0,\mathrm{},f(i)1\}`$, for $`i2`$, and define the space

$$X^{(f)}=\underset{i=1}{\overset{\mathrm{}}{}}X_i\mathrm{I}\mathrm{N}^\omega .$$

The set

$$\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})=\{w:w=w_1w_2\mathrm{}w_n,w_iX_i,1in\}$$

plays for $`X^{(f)}`$ the role played by $`X^{}`$ for $`X^\omega `$.

Prefixes of a sequence $`𝐱X^{(f)}`$ are defined in a natural way and the set of all (admissible) prefixes will be denoted by $`\mathrm{𝐩𝐫𝐞𝐟}(𝐱)`$. As we will report any coding to binary, the length of $`w=w_1w_2\mathrm{}w_n\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})`$ is $`w=\mathrm{log}(_{i=1}^nf(i));|w|=n`$. In $`X^{(f)}`$ the topology is induced by the sets $`[w]_f=\{𝐱X^f:w\mathrm{𝐩𝐫𝐞𝐟}(𝐱)\}`$ and the corresponding measure is defined by

$$\mu ([w]_f)=\underset{i=1}{\overset{|w|}{}}(f(i)^1),$$

for every $`w\mathrm{𝐩𝐫𝐞𝐟}(X^f)`$. An open set is of the form $`[A]_f=\{𝐱:n(𝐱(n)A)\}`$, for some set $`A\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})`$. The open set $`[A]_f`$ is computably enumerable if $`A`$ is computably enumerable. Only the equivalence between the notions of Cantor–randomness and weakly Chaitin-Cantor–randomness will be proven

The following two lemmas will be useful:

###### Lemma 2

Let $`0a<2^m`$ and let $`\alpha ,\beta `$ be two reals in the interval $`[a2^m,(a+1)2^m]`$. Then, the first $`m`$ bits of $`\alpha `$ and $`\beta `$ coincide, i.e., if $`\alpha =_{i=1}^{\mathrm{}}x_i2^i`$ and $`\beta =_{i=1}^{\mathrm{}}y_i2^i`$, then $`x_i=y_i`$, for all $`i=1,2,\mathrm{},m`$.

###### Lemma 3

Let $`b_1,b_2,\mathrm{}`$ be an infinite sequence of scales and $`a=j/(b_1b_2\mathrm{}b_m)[0,1].`$ Let $`\alpha ,\beta `$ be two reals in the interval $`[a,a+1/(b_1b_2\mathrm{}b_m)]`$. Then, the first $`m`$ digits of the Cantor expansions (relative to $`b_1,b_2,\mathrm{}`$) of $`\alpha `$ and $`\beta `$ coincide, i.e., if $`\alpha =_{i=1}^{\mathrm{}}x_i/(b_1b_2\mathrm{}b_i)`$ and $`\beta =_{i=1}^{\mathrm{}}y_i/(b_1b_2\mathrm{}b_i)`$, then $`x_i=y_i`$, for all $`i=1,2,\mathrm{},m`$.

## 4 Definitions of a Random Sequence Relative to the Cantor Expansion

In this section we propose five definitions for random sequences relative to their Cantor expansions and we prove that all definitions are mutually equivalent. We will fix a computable sequence of scales $`f`$.

We say that the sequence $`𝐱X^f`$ is Cantor–random if the real number $`\alpha ^𝐱`$ is random (in the sense of Algorithmic Information Theory). e.g., the sequence corresponding to the binary expansion of $`\alpha `$ is random.

Next we define the notion of weakly Chaitin–Cantor random sequence. To this aim we introduce the Cantor self-delimiting Turing machine (shortly, a machine), which is a Turing machine $`C`$ processing binary strings and producing elements of $`\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})`$ such that its program set (domain) $`\mathrm{𝑃𝑅𝑂𝐺}_C=\{x\{0,1\}^{}:C(x)\text{ halts}\}`$ is a prefix-free set of strings. Sometimes we will write $`C(x)<\mathrm{}`$ when $`C`$ halts on $`x`$ and $`C(x)=\mathrm{}`$ in the opposite case.

The program-size complexity of the string $`w\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})`$ (relative to $`C`$) is defined by $`H_C(w)=\mathrm{min}\{|v|:v\mathrm{\Sigma }^{},C(y)=w\}`$, where $`\mathrm{min}\mathrm{}=\mathrm{}`$. As in the classical situation the set of Cantor self-delimiting Turing machines is computably enumerable, so we can effectively construct a machine $`U`$ (called universal ) such that for every machine $`C`$, $`H_U(x)H_C(x)+O(1)`$. In what follows we will fix a universal machine $`U`$ and denote $`H_U`$ simply by $`H`$.

The sequence $`𝐱X^f`$ is weakly Chaitin-Cantor–random if there exists a positive constant $`c`$ such that for all $`n\mathrm{I}\mathrm{N}`$, $`H(𝐱(n))xc`$.

The sequence $`𝐱X^f`$ is strongly Chaitin-Cantor–random if the following relation holds true: $`lim_n\mathrm{}(H(𝐱(n))x)=\mathrm{}`$.

The sequence $`𝐱X^f`$ is Martin-Löf-Cantor–random if for every computably enumerable collection of computably enumerable open sets $`(O_n)`$ in $`X^{(f)}`$ such that for every $`n\mathrm{I}\mathrm{N}`$, $`\mu (O_n)2^n`$ we have $`𝐱_{n=1}^{\mathrm{}}O_n`$.

The sequence $`𝐱X^f`$ is Solovay-Cantor–random if for every computably enumerable collection of computably enumerable open sets $`(O_n)`$ in $`X^{(f)}`$ such that $`_{n=1}^{\mathrm{}}\mu (O_n)<\mathrm{}`$ the relation $`𝐱O_n`$ is true only for finitely many $`n\mathrm{I}\mathrm{N}`$.

###### Theorem 4

Let $`𝐱X^{(f)}`$. Then, the following statements are equivalent:

1. The sequence $`𝐱`$ is weakly Chaitin-Cantor–random.
2. The sequence $`𝐱`$ is strongly Chaitin-Cantor–random.
3. The sequence $`𝐱`$ is Martin-Löf-Cantor–random.
4. The sequence $`𝐱`$ is Solovay-Cantor–random.

These equivalences are direct translations of the classical proofs (see, for example, ).

Moreover, we have the following relations.

###### Theorem 5

Let $`𝐱X^{(f)}`$. Then, the sequence $`𝐱`$ is weakly Chaitin-Cantor–random if $`𝐱`$ is Cantor–random. If the function $`f`$ is bounded, then every weakly Chaitin-Cantor–random $`𝐱`$ is also Cantor–random sequence.

Proof.  The argument is modification of the proof idea of Theorem 3 in .

Assume first that $`𝐱X^{(f)}`$ is not Cantor–random and let $`\alpha =\alpha ^𝐱`$. Let $`𝐲=y_1y_2\mathrm{}`$ be the bits of the binary expansion of $`\alpha `$. We shall show that $`𝐲`$ is not a binary random sequence.

Fix an integer $`m1`$ and consider the rational

$$\alpha (m)=\underset{i=1}{\overset{m}{}}\frac{x_i}{b_1b_2\mathrm{}b_m}\text{.}$$

We note that $`w=x_1x_2\mathrm{}x_m`$ is in $`\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})`$ and $`w=\mathrm{log}(b_1b_2\mathrm{}b_m)`$. Further on, $`0<\alpha (m)<\alpha `$ and

$$\alpha \alpha (m)\underset{t=m+1}{\overset{\mathrm{}}{}}\frac{x_t}{b_1b_2\mathrm{}b_t}\underset{t=m+1}{\overset{\mathrm{}}{}}\frac{b_t1}{b_1b_2\mathrm{}b_t}=\frac{1}{b_1b_2\mathrm{}b_m}\text{.}$$

Next we define the following parameters:

$$M_m=\mathrm{log}(b_1b_2\mathrm{}b_m),$$
(3)
$$a_m=\alpha (m)2^{M_m}.$$
(4)
and we note that

$$\alpha \alpha (m)\frac{1}{b_1b_2\mathrm{}b_m}2^{M_m}.$$
(5)
We are now in a position to prove the relation: for every integer $`m1`$,

$$[\alpha (m),\alpha ][a_m2^{M_m},(a_m+2)2^{M_m}).$$
(6)
Indeed, in view of (5) and (4) we have $`\alpha <(a_m+2)2^{M_m}`$ as:

$$\alpha 2^{M_m}\alpha (m)2^{M_m}+1<a_m+2.$$

Again from (4), $`a_m\alpha (m)2^{M_m}`$.

Using (6), from $`w=x_1x_2\mathrm{}x_m`$ plus two more bits we can determine $`y_1y_2\mathrm{}y_{M_m}`$, that is, from the first $`m`$ digits of the Cantor expansion of $`\alpha `$ and two additional bits we can compute the first $`M_m`$ binary digits of $`\alpha `$. In view of Lemma 2 we obtain a computable function $`h`$ which on an input consisting of a binary string $`v`$ of length 2 and $`w`$ produces as output $`𝐲(M_m)`$.

We are ready to use the assumption that $`𝐲`$ is random but $`𝐱`$ is not Cantor–random, that is, there is a universal self-delimiting Turing machine $`U^2`$ working on binary strings and there is a positive constant $`c`$ such that for all $`n1`$,

$$H_{U^2}(𝐲(n))nc,$$
(7)
and for every positive $`d`$ there exists a positive integer $`l_d`$ (depending upon $`d`$) such that

$$H(𝐱(l_d))𝐱(l_d)d.$$
(8)
We construct a binary self-delimiting Turing machine $`C^2`$ such that for every $`d>0`$, there exist two strings $`l_d`$ and $`v,s_{l_d}\{0,1\}^{}`$, such that $`|v|=2,`$ $`|s_{l_d}|𝐱(l_d)d=\mathrm{log}(b_1b_2\mathrm{}b_{l_d})d`$ and $`C^2(v,s_{l_d})=𝐲(M_{l_d})`$.

Consequently, in view of (7) and (8), for every $`d`$ we have:

$`M_{l_d}c`$ $``$ $`H_{U^2}(𝐲(M_{l_d}))`$
$``$ $`H_{C^2}(𝐲(M_{l_d}))+O(1)`$
$``$ $`|s_{l_d}|+2+O(1)`$
$``$ $`\mathrm{log}(b_1b_2\mathrm{}b_{l_d})+O(1)`$
$`=`$ $`M_{l_d}+O(1)d,`$
a contradiction.

Recall that $`\alpha =_{i=1}^{\mathrm{}}x_i/(b_1b_2\mathrm{}b_i)=_{i=1}^{\mathrm{}}y_i2^i`$. Now we prove that $`𝐱`$ is Cantor–random whenever $`𝐲`$ is random. Let $`m1`$ be an integer and let $`\alpha _2(m)=_{i=1}^my_i2^i`$. Given a large enough $`m`$ we effectively compute the integer $`t_m`$ to be the maximum integer $`L1`$ such that

$$2^m\frac{1}{b_1b_2\mathrm{}b_L}\text{.}$$
(9)
We continue by proving that for all large enough $`m1`$:

$$[\alpha _2(m),\alpha ][\alpha (t_m)\frac{1}{b_1b_2\mathrm{}b_{t_m}},\alpha (t_m)+\frac{1}{b_1b_2\mathrm{}b_{t_m}}]\text{.}$$
(10)
We note that $`\alpha _2(m)<\alpha `$ and

$$\alpha =\underset{i=1}{\overset{\mathrm{}}{}}\frac{x_i}{b_1b_2\mathrm{}b_i}\alpha (t_m)+\underset{i=t_l}{\overset{\mathrm{}}{}}\frac{x_i}{b_1b_2\mathrm{}b_i}\alpha (t_m)+\underset{i=t_l}{\overset{\mathrm{}}{}}\frac{b_i1}{b_1b_2\mathrm{}b_i}\alpha (t_m)+\frac{1}{b_1b_2\mathrm{}b_{t_m}}\text{.}$$

As $`\alpha \alpha (t_m)+1/(b_1b_2\mathrm{}b_{t_m})`$ we only need to show that $`\alpha (t_m)\alpha _2(m)+1/(b_1b_2\mathrm{}b_{t_m})`$. This is the case as otherwise, by (9), we would have:

$$\alpha (t_m)>\alpha _2(m)+\frac{1}{b_1b_2\mathrm{}b_{t_m}}\alpha _2(m)+2^m\alpha ,$$

a contradiction.

In case when $`f`$ is bounded, assume by contradiction that $`𝐱`$ is Cantor–random but $`𝐲`$ is not random, that is there exists a positive constant $`c`$ such that for all $`n1`$ we have:

$$H(𝐱(n))\mathrm{log}(b_1b_2\mathrm{}b_n)c,$$
(11)
and for every $`d>0`$ there exists an integer $`n_d>0`$ such that

$$H_{U^2}(𝐲(n_d))<n_dd.$$
(12)
In view of Lemma 3 and (10) there is a computable function $`F`$ depending upon two binary strings such that $`|v|=2`$, $`F(𝐲(n_d),v)=𝐱(t_{n_d})`$, so the partially computable function $`FU^2`$ which maps binary strings in elements of $`\mathrm{𝐩𝐫𝐞𝐟}(X^{(f)})`$ is a Cantor self-delimiting Turing machine such that for every $`d>0`$ there exists a binary string $`s_{n_d}`$ of length less than $`n_dd`$ and a binary string $`v`$ of length 2 such that $`F(U^2(s_{n_d}),v)=𝐱(t_{n_d})`$.

As $`f`$ is bounded, the difference $`t_{m+1}t_m`$ is bounded. In view of (9), for large $`m1`$, $`b_1b_2\mathrm{}b_{t_m}>m1`$, so we can write:

$`n_dc1`$ $``$ $`\mathrm{log}(b_1b_2\mathrm{}b_{t_{n_d}})c`$
$``$ $`H_U(𝐱(t_{n_d}))`$
$``$ $`H_{FU^2}(𝐱(t_{n_d}))+O(1)`$
$``$ $`|s_{n_d}|+2+O(1)`$
$``$ $`n_dd+O(1),`$
a contradiction. q.e.d.

###### Open Question 6

It is an open question whether the above result holds true for unbounded functions $`f`$.

Consider the following statement:

Let $`𝐱`$ be a binary sequence. If there exists a computable infinite set $`M`$ of positive integers and $`c>0`$ such that for every $`mM`$, $`H_{U^2}(𝐱(m))mc`$, then $`𝐱`$ is random.

Note that if the above statement would be true, then the answer to the Open Question would be affirmative.

It is interesting to note that in case of unbounded functions $`f`$ we may have Cantor–random sequences $`𝐱X^{(f)}`$ which do not contain a certain letter, e.g. $`0X_i`$.

###### Example 7

Let $`f(i)=2^{i+2}`$. Then the measure of the set $`F=_{i=1}^{\mathrm{}}X_i^{}`$, where $`X_i^{}=X_i\{0\}`$ satisfies $`\mu (F)=_{i=1}^{\mathrm{}}(12^{i1})>0`$. Thus $`F`$ contains a Cantor–random sequence $`𝐱`$.

However, by construction, $`𝐱`$ does not contain the letter $`0`$ which is in every $`X_i`$.

## 5 On the Meaning of Randomness in Cantor’s Setting

So far, a great number of investigations have concentrated on the meaning and definition of randomness in the standard context, in which bases remain the same at all scales. That is, if one for instance “zooms into” a number by considering the next place in its expansion, it is always taken for granted that the same base is associated with different places.

From a physical viewpoint, if one looks into a physical property encoded into a real in, say, fixed decimal notation, then by taking the next digit amounts to specifying that physical property more precisely by a factor of ten. A fixed “zoom” factor may be the right choice if all physical properties such as forces and symmetries and boundary conditions remain the same at all scales. But this is hardly to be expected. Take, for instance, a “fractal” coastline. How is it generated? The origins of its geometry are the forces of the tidal and other forces on the land and coastal soil. That is, water moving back and forth, forming eddies, washing out little bays, and little bays within little bays, and little bays within little bays within little bays, … and so on. There may be some structural components of this flow which results in scale dependence. Maybe the soil-water system forming the landscape will be “softer” at smaller scales, making bays relatively larger that their macroscopic counterparts. Indeed, eventually, at least at subatomic scales, the formation of currents and eddies responsible for the creation of ever smaller bays will break down.

In such cases, the base of the expansion might have to be modified in order to be able to maintain a proper relation between the coding of the geometric object formed by the physical system and the meaning of its number representation in terms of “zooming”. All such processes are naturally stochastic, and therefore deserve a proper and precise formalization in terms of random sequences in Cantor representations.


text:
# Equilibrium Distribution of Mutators in the Single Fitness Peak Model

## Abstract

This paper develops an analytically tractable model for determining the equilibrium distribution of mismatch repair deficient strains in unicellular populations. The approach is based on the single fitness peak (SFP) model, which has been used in Eigen’s quasispecies equations in order to understand various aspects of evolutionary dynamics. As with the quasispecies model, our model for mutator-nonmutator equilibrium undergoes a phase transition in the limit of infinite sequence length. This “repair catastrophe” occurs at a critical repair error probability of $`ϵ_r=L_{via}/L`$, where $`L_{via}`$ denotes the length of the genome controlling viability, while $`L`$ denotes the overall length of the genome. The repair catastrophe therefore occurs when the repair error probability exceeds the fraction of deleterious mutations. Our model also gives a quantitative estimate for the equilibrium fraction of mutators in Escherichia coli.

Mutator, quasispecies, repair catastrophe

In order to preserve the integrity of their genomes, living systems have evolved sophisticated mechanisms for correcting errors in their DNA sequences Voet and Voet (1995). Otherwise, genetic damage due to environmental factors such as radiation, metabolic free radicals, and mutagens, combined with replication errors, would lead to unviable organisms due to the unrecoverable loss of genetic information. This phenomenon, which was first characterized by Eigen in Eigen (1971), is known as the “error catastrophe” Eigen (1971); Eigen et al. (1989). It has since been studied in a number of theoretical papers Gallucio (1997); Nilsson and Snoad (2002, 2000) (and references therein), and has also been observed experimentally Crotty et al. (2001).

Some of the error-correcting ability in living systems is already built into the DNA polymerases themselves. In Escherichia coli, the proofreading ability of the DNA polymerases Pol I and Pol II results in an error probability of $`10^610^7`$ per base pair Voet and Voet (1995). Additional enzymes continuously scan the DNA molecule, repairing lesions and mismatches that occur due to environmental damage.

A key error-repair mechanism is known as mismatch repair, and occurs immediately following DNA replication. The mismatch repair system scans the DNA molecule, identifies, and then corrects mismatched base pairs. Mismatch repair in E. coli reduces the error probability in DNA replication to $`10^810^{10}`$ per base pair Voet and Voet (1995). Cells with inactivated mismatch repair consequently have mutation rates which are $`10`$ to $`10,000`$ times higher than cells whose mismatch repair system is functioning. Because of their higher than wild-type mutation rates, these “mutator” strains are believed to play an important role in the emergence of antibiotic resistance, and cancer in multicellular organisms Sniegowski et al. (1997); LeClerc et al. (1996); Taddei et al. (1997); Tenaillon et al. (1999); Giraud et al. (2001); Matic et al. (1997).

To develop a model for the equilibrium distribution of mutators versus nonmutators in a unicellular population, we consider a genome of alphabet size $`S`$ (“bases” $`0,1,\mathrm{},S1`$), consisting of two genes. The first gene consists of $`L_{via}`$ bases, and controls the viability of the genome. The second gene consists of $`L_{rep}`$ bases, and codes for the enzymatic machinery involved in mismatch repair. If we let $`\sigma `$ denote an arbitrary gene sequence, then we may write, $`\sigma =\sigma _{via}\sigma _{rep}`$.

We assume a single fitness peak (SFP) model for both genes. Thus, there is a unique, “fit” sequence $`\sigma _0=\sigma _{via,0}\sigma _{rep,0}`$. A cell with genome $`\sigma `$ has a first-order growth rate constant $`k>>1`$ if $`\sigma _{via}=\sigma _{via,0}`$, and $`1`$ otherwise. Mismatch repair has an error probability of $`ϵ_r`$ per mismatched base pair, and is functioning only when $`\sigma _{rep}=\sigma _{rep,0}`$.

While somewhat artificial, the SFP model has been successfully applied in Kamp and Bornholdt (2002) toward understanding the correlations between antibody and viral mutation rates. Furthermore, because proteins generally have a key set of conserved residues which more or less dictate their final structure and function, the corresponding gene has a subsequence of conserved bases required for its proper function Mirny and Shakhnovich (2001); Voet and Voet (1995). Thus, by summing over the unconserved bases, it is possible to reduce the fitness landscape to an SFP in the conserved subsequence. Therefore, there is reason to believe that many of the phenomenological aspects of our system can be captured by an SFP-based approach, and that such an approach can also be semiquantitative in a number of cases.

The basic equation governing the dynamics on the genome space has the form of Eigen’s quasispecies equations Eigen (1971); Eigen et al. (1989),

$$\frac{dx_\sigma }{dt}=(\kappa _\sigma \overline{\kappa }(t))x_\sigma +\underset{\sigma ^{}\sigma }{}[\kappa _m(\sigma ^{},\sigma )x_\sigma ^{}\kappa _m(\sigma ,\sigma ^{})x_\sigma ]$$
(1)
where $`x_\sigma `$ denotes the fraction of the population with genome $`\sigma `$, $`\kappa _\sigma `$ is the growth rate constant of $`\sigma `$, $`\kappa _m(\sigma ,\sigma ^{})`$ denotes the mutation rate constant from $`\sigma `$ to $`\sigma ^{}`$, and $`\overline{\kappa }(t)_\sigma \kappa _\sigma x_\sigma (t)`$.

We assume that replication errors are sufficiently small so that we need only worry about point mutations, and that since $`k>>1`$, any flow off of the viability peak is unidirectional. Furthermore, since we are only interested in the relative distribution of viable mutators and non-mutators, we only focus on a subspace of sequences given by $`\sigma =\sigma _{via,0}\sigma _{rep}`$. Thus, from now on, we shall simplify matters and redenote $`\sigma _{rep}`$ as $`\sigma `$. The full genome is $`\sigma _{via,0}\sigma `$ by implication. On this subspace, the effective growth rate constant becomes $`k(1L_{via}ϵ_\sigma )`$, due to leakage off of the fitness peak. Here $`ϵ_\sigma `$ denotes the per base pair replication error probability, and is equal to $`ϵϵ_r`$ if $`\sigma =\sigma _0`$, and $`ϵ`$ otherwise. Finally, by the symmetry of our system, we may make the further assumption that $`x_\sigma `$ depends only on the Hamming distance $`HD(\sigma ,\sigma _0)`$ from $`\sigma _0`$. Thus, defining $`\mathrm{\Omega }_l(\sigma )=\{\sigma ^{}|HD(\sigma ^{},\sigma )=l\}`$, we may then also define $`x_l=x_\sigma `$, where $`\sigma \mathrm{\Omega }_l(\sigma _0)`$. Note that point mutations between $`x_l`$ and some $`x_\sigma `$ may only occur if $`x_\sigma \mathrm{\Omega }_{l,l\pm 1}(\sigma _0)`$. However, we may negect intra-$`\mathrm{\Omega }_l(\sigma _0)`$ couplings, due to cancellation of mutational inflows and outflows.

A $`\sigma \mathrm{\Omega }_l(\sigma _0)`$ may be connected via a point mutation to a $`\sigma ^{}\mathrm{\Omega }_{l1}(\sigma _0)`$ by changing any one of the $`l`$ bases distinct from the corresponding bases in $`\sigma _0`$ back to the corresponding base in $`\sigma _0`$. Thus, there are $`l`$ possible connections. A $`\sigma \mathrm{\Omega }_l(\sigma _0)`$ may be connected via a point mutation to a $`\sigma ^{}\mathrm{\Omega }_{l+1}(\sigma _0)`$ by changing any one of the $`L_{rep}l`$ bases equal to the corresponding bases in $`\sigma _0`$. Since there are $`S1`$ possibilities per base, the result is $`(L_{rep}l)(S1)`$ connections. The net mutational flow is then,

$$\underset{\sigma ^{}\sigma }{}[\kappa _m(\sigma ^{},\sigma )x_\sigma ^{}\kappa _m(\sigma ,\sigma ^{})x_\sigma ]=\frac{kl}{S1}(ϵ_{l1}x_{l1}ϵ_lx_l)+k(L_{rep}l)(ϵ_{l+1}x_{l+1}ϵ_lx_l)$$
(2)
where $`ϵ_0=ϵϵ_r`$, and $`ϵ_l=ϵ`$ for $`l1`$. We divide the $`ϵ`$’s by $`S1`$ because a point mutation can occur to any one of the $`S1`$ bases distinct from the changed base.

We also have, $`\kappa (t)=k(1L_{via}ϵ)+kL_{via}ϵ(1ϵ_r)x_0`$. Now, define $`C_l=\left(\genfrac{}{}{0pt}{}{L_{rep}}{l}\right)(S1)^l`$, the number of elements in $`\mathrm{\Omega }_l(\sigma _0)`$, and set $`z_l=C_lx_l`$. If we reexpress the dynamical equations in terms of $`z_l`$, then at equilibrium we obtain the system of equations,

$`0`$ $`=`$ $`{\displaystyle \frac{L_{via}}{L_{rep}}}(1ϵ_r)z_0(1z_0)+{\displaystyle \frac{z_1}{L_{rep}(S1)}}ϵ_rz_0`$
$`0`$ $`=`$ $`{\displaystyle \frac{L_{via}}{L_{rep}}}(1ϵ_r)z_0z_1+ϵ_rz_0(1{\displaystyle \frac{1}{L_{rep}}}+{\displaystyle \frac{1}{L_{rep}(S1)}})z_1+{\displaystyle \frac{2}{L_{rep}(S1)}}z_2`$
$`\mathrm{}`$
$`0`$ $`=`$ $`{\displaystyle \frac{L_{via}}{L_{rep}}}(1ϵ_r)z_0z_l+(1+{\displaystyle \frac{1}{L_{rep}}}{\displaystyle \frac{l}{L_{rep}}})z_{l1}(1{\displaystyle \frac{l}{L_{rep}}}+{\displaystyle \frac{l}{L_{rep}(S1)}})z_l+{\displaystyle \frac{l+1}{L_{rep}(S1)}}z_{l+1}`$
$`\mathrm{}`$
$`0`$ $`=`$ $`{\displaystyle \frac{L_{via}}{L_{rep}}}(1ϵ_r)z_0z_{L_{rep}}+{\displaystyle \frac{z_{L_{rep}1}}{L_{rep}}}{\displaystyle \frac{z_{L_{rep}}}{S1}}`$ (3)
Except for the last equation, the $`(l+1)^{\text{st}}`$ equation has a mutational contribution from $`z_{l+1}`$ which scales as $`1/L_{rep}`$. This means that the contribution to $`z_l`$ due to backmutation from $`z_{l+1}`$ becomes negligible for large sequence lengths. This makes sense, since for finite $`l`$, the ratio $`C_l^{}/C_l\mathrm{}`$ as $`L_{rep}\mathrm{}`$ for $`l^{}>l`$, so the probability of mutating to lower values of $`l`$ vanishes in the limit of infinite sequence length.

The term, $`(L_{via}/L_{rep})(1ϵ_r)z_0(1z_0)`$, and the corresponding terms in the other equations arise from the $`\overline{\kappa }(t)`$ term in the original quasispecies equations (Eq. (1)) of our model. Because the nonmutator sequence has a lower level of leakage off of the viability peak than the mutator sequences, it has a higher effective growth rate constant. Thus, when dealing with population fractions as opposed to absolute populations, the result is an effective positive flow into the nonmutator sequence given by the term $`(L_{via}/L_{rep})(1ϵ_r)z_0(1z_0)`$. By a similar argument, it can be seen why the corresponding terms in the remaining equations are negative.

We wish to solve these equations for a fixed value of $`\alpha L_{via}/L_{rep}`$ in the limit of infinite sequence length $`L`$. Let us focus first on the behavior of $`z_0`$ in this limit. In the first equilibrium equation, the $`z_1`$ term drops out as $`L_{rep}\mathrm{}`$, giving,

$$0=z_0(\alpha (1ϵ_r)(1z_0)ϵ_r)$$
(4)
which has the solutions $`z_0=0,1ϵ_r/(\alpha (1ϵ_r))`$. The first solution is inconsistent with the requirement that $`z_0=1`$ when $`ϵ_r=0`$. However, the second solution only holds as long as $`z_0[0,1]`$. Clearly, $`z_01`$ $``$ $`ϵ_r`$. The other requirement that $`z_00`$ gives,

$$ϵ_r\frac{\alpha }{1+\alpha }=\frac{L_{via}}{L}$$
(5)
Defining $`ϵ_{r,crit}=\frac{\alpha }{1+\alpha }`$, we see that in the limit of infinite genome length, our system has two “phases.” For $`ϵ_r<ϵ_{r,crit}`$ the system is in a “non-mutator,” or, equivalently, “repairer” phase, in which the fraction of non-mutators is a quantity which depends only on $`\alpha `$ and $`ϵ_r`$. At $`ϵ_r=ϵ_{r,crit}`$ the system undergoes a “phase” transition, which we term the “repair catastrophe,” after which the system is in a “mutator” (“non-repairer”) phase. In this phase, there is essentially no preference for being a non-mutator, and the fraction of non-mutators becomes inversely proportional to the total number of gene sequences.

A key parameter to study in the phase behavior of our model is the localization length, given by,

$$l=\underset{l=1}{\overset{L_{rep}}{}}lz_l$$
(6)
This quantity measures the mean Hamming distance of the population from the nonmutator sequence. To compute $`l`$ below the phase transition in the limit of $`L_{rep}\mathrm{}`$, we may note that for finite $`l`$ our equilibrium equations become,

$`0`$ $`=`$ $`\alpha (1ϵ_r)z_0(1z_0)ϵ_rz_0`$
$`0`$ $`=`$ $`\alpha (1ϵ_r)z_0z_1+ϵ_rz_0z_1`$
$`\mathrm{}`$
$`0`$ $`=`$ $`\alpha (1ϵ_r)z_0z_l+z_{l1}z_l`$
$`\mathrm{}`$ (7)
We have already solved the first equation. The next two equations can be solved together to give, for $`l1`$,

$$z_l=ϵ_r(1+\alpha (1ϵ_r)z_0)^lz_0$$
(8)
It should be noted that, while each $`z_l`$ converges to the corresponding formula given above as $`L_{rep}\mathrm{}`$, the convergence is not uniform, since the larger the $`l`$, the larger $`L_{rep}`$ must be made to get $`z_l`$ within some specified cutoff of its $`L_{rep}=\mathrm{}`$ value.

Define $`z_{l,\mathrm{}}=lim_{L_{rep}\mathrm{}}z_l`$. It may be readily checked that $`_{l=0}^{\mathrm{}}z_{l,\mathrm{}}=1`$, so total population is conserved in this limiting process. The localization length is given by,

$$l=\underset{l=1}{\overset{\mathrm{}}{}}lz_{l,\mathrm{}}=\frac{1ϵ_{r,crit}}{ϵ_{r,crit}}\frac{ϵ_r}{ϵ_{r,crit}ϵ_r}$$
(9)
Note that, as expected, the localization length is finite for $`ϵ_r<ϵ_{r,crit}`$, but diverges at the phase transition.

It is also useful to solve the equilibrium distribution exactly for the case $`\alpha =0`$. This corresponds to $`L_{rep}=L`$, that is, the entire genome consists of the repair gene. Note that $`ϵ_{r,crit}=0`$, so that the system is always in the mutator phase. In this case, it may be shown that the equilibrium solution is given by, $`z_l=\left(\genfrac{}{}{0pt}{}{L_{rep}}{l}\right)(S1)^lϵ_rz_0`$ for $`l>0`$, and so the requirement that $`_{l=0}^{L_{rep}}z_l=1`$ gives, $`z_0=1/(1+ϵ_r(S^{L_{rep}}1))`$. It is readily shown that for large $`L_{rep}`$, the localization length $`l(11/S)L_{rep}`$. This result is equivalent to the case of a uniform distribution, which makes sense since for $`\alpha =0`$ there is no preference for being a nonmutator in the limit of large $`L_{rep}`$.

The phase behavior which emerges from this model may be understood as follows: For highly efficient repair, the selective advantage for being a nonmutator is sufficiently large to cause the population to equilibrate in a localized cluster about the nonmutator sequence. When repair is inefficient, the accumulation of deleterious mutations in both mutators and nonmutators is comparable, and hence the mutators, which are entropically strongly favored, dominate the population. The selective advantage for being a nonmutator is dictated by $`\alpha `$, since for low $`\alpha `$ there is relatively little leakage off of the fitness peak, while for high $`\alpha `$ there is a large amount of leakage off of the fitness peak. Thus, for low $`\alpha `$, repair has to be highly efficient to give the nonmutators a sufficient selective advantage to be in the nonmutator phase, while for high $`\alpha `$, nonmutators have a significant advantage even for relatively inefficient repair.

One of the main features to note regarding the mutator-nonmutator equilibrium is that it is independent of the background error probability $`ϵ`$. This feature is interesting because as $`ϵ0`$, the difference in viability between the mutators and the nonmutators disappears. Thus, one might naively expect $`ϵ_{r,crit}`$ to be a function of $`ϵ`$, but in the limit of small $`ϵ`$ (so that only point mutations are important), this is not the case.

Figure 1 shows a plot of $`z_0`$ versus $`ϵ_r`$ for $`\alpha =1/9,1,`$ and $`9`$. We used $`S=2`$, and took a value of $`L_{rep}=1,000`$ in order to sufficiently converge the calculations. The equilibrium equations were solved by using fixed-point iteration at every $`ϵ_r`$. Note that the phase transition does indeed occur at the predicted values of $`ϵ_{r,crit}`$. The analytical, $`L_{rep}=\mathrm{}`$ curves lie essentially on top of our numerical results, and were therefore not plotted here.

Figure 2 shows the corresponding plots of $`l`$ versus $`ϵ_r`$. Note that the localization lengths settle at the value of $`L_{rep}/2`$ at the transition, which is the expected finite $`L_{rep}`$ behavior for the mutator phase.

Finally, we may use our model to estimate the equilibrium fraction of mismatch repair deficient strains in E. coli. The E. coli genome has $`4,639,221`$ base pairs, comprising $`4,403`$ genes et al. (1997). Based on calculations for Saccharomyces cerivisiae, or Baker’s yeast, we esimate that between $`1830\%`$ of these genes are “viability” genes, i.e, required for E. coli survival et al. (2002, 1999) (unfortunately, similarly detailed data is not currently available for E. coli, so we had to make an estimate based on available information). Thus, we assume approximately $`1,000`$ “viability” genes, which we gather into the viability peak of our model. The mismatch repair system is controlled by the MutH, MutL, MutS, and UvrD (or MutU) proteins, giving $`4`$ repair genes. If we simply use the average gene length, and assume that the same fraction of base pairs must be conserved in both the viability genes and in the mismatch repair genes, then in our model we obtain $`\alpha 1000/4=250`$. Since mismatch repair has a failure probability of $`10^410^1`$ per mismatched base pair, we estimate an equilibrium fraction of mutators in the range of $`4\times 10^74\times 10^4`$. For E. coli, the observed equilibrium fraction of mutators is on the order of $`10^510^3`$ LeClerc et al. (1996). While encouraging, our result is nevertheless based on a number of simplifying assumptions. The strongest evidence in support of our model would be the experimental observation of the repair catastrophe itself. While it is not clear how to selectively control the efficiency of the mismatch repair system, if possible this would allow a direct experimental test of our model.

As a concluding remark, we should note that our prediction of a repair catastrophe in mutator-nonmutator equilibrium suggests that phase transitions may underlie the behavior of a variety of biological systems. A classification of the phase behaviors inherent in various biological networks will greatly increase our understanding of the underlying dynamics governing such systems.

###### Acknowledgements.

This research was supported by an NIH postdoctoral research fellowship, and by the Howard-Hughes Medical Institute graduate research fellowship. The authors would like to thank Dr. Andrew Murray and Ethan Perlstein for helpful conversations regarding this work.