text:
# Face-on dust disks in galaxies with optical jets1footnote 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, Incorporated, under NASA contract NAS5-26555.

## 1 Introduction

The reasons why radio galaxies project high energy jets of relativistic plasma over vast, galactic and even extragalactic, distances from a tiny volume at the nucleus of the host galaxy remain unknown. Essentially all radio galaxies seem to display jets when observed with sufficient sensitivity and resolution, suggesting this is a fundamental characteristic of radio galaxies and energy transport. A small fraction of these jets are seen in the optical. This represents a potentially profound insight into the nature of jets since the presence of optical synchrotron radiation demands relativistic particles of high energy and short lifetime, unlike the radio, and in-situ acceleration processes are probably required.

It is widely believed that a rotating black hole at the nucleus, coupled to an inner accretion disk is somehow responsible for the onset of a jet. The extraordinary collimation is likely to occur at the very nucleus, though galactic influences at larger distances may play a role in collimating the jet, Junor, Biretta & Livio (1999). It is also generally accepted that both orientation (the angle of the emerging jet and of the inner accretion disk to the line of sight) and relativistic beaming, play roles in the appearance of any particular source to the observer. By such means, FR-I (or low power) radio galaxies may, when viewed from the appropriate direction, appear as luminous, variable BL Lac objects, Urry & Padovani (1995).

Optical emission from jets, in the beaming scenario in which the jet moves relativistically towards the observer, is favoured because both the typical synchrotron spectral break is blueshifted closer to the optical window, and also because the intrinsic intensity “beam pattern” is strongly forward projected, with an opening angle of order $`1/\gamma `$ where $`\gamma `$ is the relativistic factor, $`\gamma =\sqrt{(1/(1v^2/c^2))}`$.

Other physical processes may play a role, however. The high-pressure environment close to a galaxy nucleus could quench the progression of the jet and the high pressure may translate to high synchrotron losses and a greater tendency for optical radiation. The case of 3C 264 offers evidence of such an influence, as there the optical jet terminates at the very edge of an inner dust lane, Baum et al. (1997). Age, or rather youth, may also play an important role. Best, Longair & Röttgering (1996) show the evolution of optical morphology of radio galaxies with increasing radio source size, for example, and discuss this in the context of jet-induced star formation and the alignment effect.

Here, we make a very simple observation that we believe strongly supports the beaming and orientation paradigm, which is that the majority of nearby radio galaxies showing optical jets have compact well-ordered dust disks that are almost circular and are hence likely to be face-on. This implies that jets emerge close to perpendicular to such disks and that they are pointing almost towards us. We also note that the nuclear point source luminosities are consistent with this general behaviour. In § 2 we describe the data, in § 3 we discuss and conclude.

## 2 HST imaging observations of 3CR galaxies

In a series of proposals, we acquired multiwavelength images of 3CR radio galaxies and quasars using the Wide Field and Planetary Camera 2 (WFPC2) on-board the Hubble Space Telescope (HST). The observations utilised the efficient snapshot mode of observing, whereby short exposures are used to fill scheduling gaps. The data and observations are described in de Koff et al. (1996), McCarthy et al. (1997), Martel et al. 1998 (b), Lehnert et al. (1999), Martel et al. (1999). An early goal of our observations was to use the high sensitivity and spatial resolution of HST to locate new optical synchrotron jets in radio galaxies and hence understand the systematics relating high energy synchrotron emission with the host properties. That effort has been successful, and is ongoing. New optical jets were found in 3C 78, Sparks et al. (1995), 3C 15, Martel et al. 1998 (a) and elsewhere, Sparks et al. (1996). To date, there are about 12 extragalactic radio sources with known optical jets, and a similar number of candidate optical jets.

The majority of nearby confirmed optical jets are in FR-I galaxies and for the purposes of this paper we define a sample of nearby FR-I galaxies by restricting attention to 3CR radio galaxies from Bennett (1962) whose redshift $`z<0.1`$ and $`|b^{II}|>10^{}`$, Spinrad et al (1985). Galaxies classified as either FR-I or FR-I/II in Zirbel and Baum (1995) were selected, resulting in a sample of 25 radio sources, counting 3C 75 as two, as listed in Table 1. Five have known optical jets: 3C 15, 3C 66B, 3C 78, 3C 264, and 3C 274 (M87).

Nuclear point sources are common in our high resolution images of radio galaxies, and may provide additional insight into the beaming hypothesis. We also therefore measured the nuclear core brightness where there was an obvious point source using aperture photometry and the magnitude zeropoint from the HST Data Handbook. The magnitudes are on the Vega system and are presented in Table 1.

## 3 Dust lanes in nearby 3CR FR-I radio galaxies

de Koff et al. (2000) investigated the systematics of dust lanes in 3CR galaxies from the snapshot observations and show that in FR-I galaxies, dust lanes tend to be regular and compact and that in such sources there is a strong tendency for the radio axis to be perpendicular to the dust disk. Typical scale lengths for the dust disks are hundreds of parsecs. Martel et al. (2000) present an in-depth analysis of the dust-lanes in several of these nearby radio galaxies.

Baum et al. (1997) presented a detailed study of 3C 264 and showed its remarkable circular dust disk, with the bright optical jet apparently terminating at the edge of the disk, Fig. 2. The morphology of this dust lane was quite unique, with no other dust lanes of the de Koff et al. (2000) sample having such a circular appearance. However, detection of round, low optical depth dust features can be difficult since they can be significantly harder to see than edge-on features that cut into the galaxy isophotes.

Motivated by the appearance of the 3C 264 disk, and a similar one in 3C 66B, Sparks et al. (2000), we inspected all the remaining low$`z`$ FR-I optical jet galaxies to search for the presence of subtle dust features, and in addition those of the well-defined sample of 25 (above) which had not already been investigated. This was done in case, by looking harder at the jet galaxies than the others, we were biassing our result. However, of the 25 galaxies, 24 have clearly detectable dust. Even the remaining one, the Southern component of 3C 75, shows a hint of a very faint feature in an archival HST image. The dust detection rate is 96% and we conclude therefore that we are not witnessing the effects of missing subtle dust features in galaxies: our dust detection rate is essentially complete. Fig. 1 shows dust images of those galaxies not shown in the references above.

Fig. 2 shows the dust in the optical jet galaxies shown in a variety of ways. All five of the optical jet objects have multicolor HST images available. For 3C 264 and 3C 66B no special processing is needed, and we simply show a single direct image for each. For the others, we fitted isophotes to the red F702W galaxy images with iterative sigma clipping to reject isophotal outliers (i.e. the jets). The isophotal profiles were then used to make a smooth, perfectly elliptical model of the red light in the galaxy (where the effects of dust are minimized), and the bluer data is divided by this red model as in Sparks et al. (1985) to generate an image of the dust. These are the $`(VR_{\text{model}})`$ images, or “dust images” referred to elsewhere in the text.

The center row of Fig. 2 shows the galaxy hosting 3C 15. To the left is a direct image, whose irregularities are attributable to dust. In the center is a pure $`(VR)`$ image, and on the right is the $`(VR_{\text{model}})`$ image. The large scale circular features in the $`(VR_{\text{model}})`$ image are most likely artifacts from the modelling, unlike the more compact feature towards the center which also shows in the $`(VR)`$ data image (i.e. the ratio of the $`V`$ image to the $`R`$ image, as opposed to the ratio of the $`V`$ image to a smooth, elliptical model of the $`R`$ image). The dust optical depth is low (as might be expected if the dust lane is close to face-on), but by comparing the different images, and perhaps most easily seen in the center one, there is a faint, almost complete, ring of absorption and reddening, rather deeper to the West. The jet is on a much larger scale than the dust, as it is in 3C 66B.

The lowest row of Fig. 2 shows the dust in 3C 78. To the left is a direct $`V`$ image, in the center is a direct $`R`$ image and on the right, the $`V`$ image divided by a smooth model of the $`R`$ image. Again, by comparing amongst the images, the absorption can be seen in the color image and in the direct $`V`$. It is almost absent in the $`R`$ image. Once more, absorption can be traced entirely around the nucleus in a compact ring.

Hence, remarkably, four of the five optical jet galaxies show circular dust lanes. The other one is M87 whose dust shows a more chaotic morphology, with distinct, roughly radial filaments that coincide with line-emitting gas, Sparks, Ford & Kinney (1993), Sparks (1999). None of the other 25 galaxies show similar circular dust lanes.

To obtain a better estimate of the axial ratios of each of the “round” dust lanes, four cuts were made at $`45^{}`$ angles to one another across the clearest images. The diameters were marked by eye, and these measurements were then used to determine the parameters of the least-squares best-fit ellipse. The resulting measurements are given in Table 1. Note that uncertainties tend to bias the measurements to smaller $`b/a`$ values if the disk is truly circular. The remaining dust lane axial ratios were taken either from the literature as cited, or measured from the direct images referenced using a major and minor axis cut.

Fig. 3 shows the axial ratio distributions of the dust lanes in radio galaxies for the 15 galaxies with reasonably well-defined dust disks, and it is clear that the optical-jet objects have a significantly rounder distribution. Formally, the difference is significant at the 0.1% level (i.e. a probability of chance occurrence of less than 0.001), based on the plotted axial ratio distributions and using the Fisher exact probability test. The diameters of the dust disks (or rings) are approximately 1300, 350, 430 and 620 pc for 3C 15, 66B, 78 and 264 respectively, assuming $`H_0=75`$ km/s/Mpc.

## 4 Discussion

### 4.1 Dust and jet properties

The simplest way to understand this result is that a round dust disk is a dust disk seen close to face-on. Since we have shown in de Koff et al. (2000), that radio jets tend to lie perpendicular to dust disks, especially in the settled FR-I objects, the natural inference therefore is that we are also seeing the optical jets almost pole on, that is pointing directly towards us. de Koff et al. (2000) found that for FR-I sources with dust disks of size $`<2.5`$ kpc, the radio jet lies within less than or about 15 degrees of the disk (perpendicular) axis (but with a small number of exceptions such as 3C31). The roundest disk in the sample without jets would imply an inclination $`\varphi `$ of order 40 degrees to the line-of-sight if it is intrinsically circular. Similarly, the flattest disk in the jet sample would require an angle $`\varphi `$ to the line of sight of $`\varphi 30`$ degrees, taken at face value. This suggests that the dispersion in angles $`\varphi `$ is dominated by the projection of the assembly of dust disk plus jet, rather than by an intrinsic scatter of jet angles to dust disk angles. It also implies that there is a critical angle, $`\varphi _c`$, such that if jets point within that angle we see optical emission. From these numbers we expect $`\varphi _c`$ to lie in the range 30—40 degrees.

There are some small asymmetries and irregularities apparent in the circular dust lanes, and these may be due to the effects of what inclination there is to the plane of the sky, or else to intrinsic irregularities within the disk. Detailed analysis of the optical depth around the disk becomes highly model dependent as there is nowhere to obtain the unobscured light profile.

Why jets and dust lanes are close to perpendicular remains a mystery. Since these are relaxed, orderly systems, perhaps the disk, galaxy and black hole spin axes have had time to co-align, Natarajan & Pringle (1998). Alternatively, the jet ejection process itself may react back on the disk and force it into a perpendicular orientation, Quillen and Bower (1999). Also, in these round dust disks, the outer edge of the disk is remarkably sharp. Perhaps this is an indication of a more settled and evolved disk than the irregular morphologies often seen, de Koff et al. (2000), which would in turn suggest that evolutionary processes within the jets are not responsible for the optical emission. One would expect galaxy dynamical timescales to be much longer than the evolutionary timescale for a radio jet (given its high speed and short synchrotron loss lifetime).

If the sample is a set of randomly oriented linear jets, we expect the distribution of $`\varphi `$ to be $`sin(\varphi )`$ for $`\varphi `$ in the range 0 to 90 degrees. Hence, statistically we would expect a fraction $`1cos(\varphi _c)`$ to have an angle to the line of sight $`\varphi `$ less than a critical angle $`\varphi _c`$. If $`\varphi _c40`$ degrees, then the fraction is 23% which for a sample size of 25 means six objects (three if $`\varphi _c=30`$ degrees). This is obviously consistent with the five that we find, if the optical jet galaxies represent those whose jets point within 30—40 degrees to the line of sight.

It is possible that more face-on disks tend to appear more irregular and hence are excluded from the statistical comparison. There are no compelling reasons to think that this is actually the case, however. One of the 10 galaxies with “irregular” dust has an optical jet; by chance we would expect 1.8 galaxies to be oriented with $`\varphi <35`$ degrees so there is nothing inconsistent in the statistics to suggest a large fraction of the irregular dust disks are actually face-on.

Hence, we are led to a picture where optical jets “appear” when they point to within (of order) $`\varphi <\varphi _c35\pm 5`$ degrees to our line-of-sight, based on the apparent axial ratios of their associated dust lanes.

This is reasonable qualitatively, given expected relativistic beaming and Doppler boosting of the underlying synchrotron break frequency. $`\gamma `$ values of order ‘a few’ are typically thought to occur jets in radio galaxies, and indeed superluminal proper motion is observed in the jet of M87 which is very direct evidence for $`\gamma `$ values of order $`\gamma 6`$, Biretta, Sparks & Macchetto (1999). If the beam angular width is $`1/\gamma `$, a critical angle of $`40`$ degrees implies a rather lower value, $`\gamma 1.4`$ or $`v0.7c`$. If $`\varphi _c=30`$ degrees, the corresponding $`\gamma =2`$ and $`v=0.85c`$. Such numbers should be taken as illustrative: clearly, in principle, beams are observable outside the primary beam width, at much lower surface brightness. Much deeper observations with fainter surface brightness limits may help to disentangle the effects of multiple velocity components and intrinsic beam pattern.

Laing et al. (1999) discuss the case for relativistic beaming in FR-I radio galaxies in detail. Based on comparison with radio data, primarily, they develop a consistent model in which jets are intrinsically symmetric, initially relativistic but decelerating, and faster on-axis than at the edges. This is quite consistent with our inference of jet $`\gamma >1.4`$ provided that either optical emission becomes unobservable at lower values of $`\gamma `$ (where the beam is broad enough for many more jets to be included within it) or else the jets never have velocity less than $`v=0.7c`$. Laing et al. (1999) finds a best fitting velocity for single-velocity models, of $`v_0=0.72c`$, but they also achieve good fits with lower minimum velocities when a range of velocities is allowed. It seems more plausible that the optical emission becomes undetectable (with current observations) at lower velocity.

### 4.2 Nuclear core sources

Nuclear point sources are common in radio galaxies, and are often apparently bright in galaxies with optical jets, Martel et al. 1998 (a). Chiaberge, Capetti & Celotti (1999) and Hardcastle & Worrall (2000) showed that the optical core source luminosity is very well correlated with the radio core flux, and hence is most likely also synchrotron emission arising from the base of the jet. In a pilot study of five objects Capetti & Celotti (1999) showed that there are indications of a correlation between gas disk inclination and core prominence, consistent with the expectations of the proposed unification of BL Lac and FR-I radio galaxies. To take this a step further, we plot the nuclear magnitude where there was obviously a point source against the axial ratio of the dust lanes, Fig. 4. There is a trend, with the brightest cores in the most face-on disks, and a correlation coefficient of $`0.83`$ yielding a formal significance level of better than 1% . However, many targets do not appear on the plot. Either they have cores but irregular dust, including two of the brightest, or else they do not have well-defined cores in our data. We consider the plot therefore primarily indicative and intriguing for future study.

To examine the FR-I/BL Lac unification idea, we include on Fig. 4 the average of the BL Lac objects from the comparison sample of Capetti & Celotti (1999), showing the $`\pm 1\sigma `$ dispersion in their magnitudes. We assume this representative BL Lac object is viewed directly along the line-of-sight. If the average spectral index is $`\alpha 1`$ in the optical, the relativistic boosting for continuous outflow is $`(1\beta cos\theta )^{2\alpha }`$, where $`\beta =v/c`$ and $`b/a=cos\theta `$. For illustration, the best-fit value for $`\gamma 5.6`$ is shown, with the zeropoint of the proportionality chosen so that the curve passes through the BL Lac point and assuming that all nuclei have the same intrinsic brightness. The derived $`\gamma `$ is sensitive to $`\alpha `$, with lower $`\alpha `$ yielding higher $`\gamma `$. While there are major uncertainties and selection effects (we only see the most luminous cores; there is no expectation that all cores have identical intrinsic flux) it is interesting to note that there is a good consistency with conventional beaming models and relativistic factors between 1 and 10. A more detailed investigation into these issues would be useful.

### 4.3 Caveats

Perhaps the most worrisome objection to the scenario described in § 4.1 is that in the case of 3C 264, the optical jet flares and essentially ends right at the edge of the dust disk. This leads one to suppose either (i) a causal connection, as discussed by Baum et al. (1997), or (ii) a coincidence. None of the other optical jets show any sign of interaction with the host disks: for the most part they are on a much larger spatial scale. In the scenario where the disk is face-on and the jet almost perpendicular to it, the coincidence of the jet flaring out at the edge of the disk must be accepted. In the contrary situation, where the jet ends at the edge of the disk due to a physical interaction with the disk, the alternate coincidence that of all galaxies, this one has the most face-on disk must be accepted. More exotic options, such as dust spheres in the centers of optical-jet galaxies rather than disks, can be excluded from detailed analysis of the color distribution in the case of 3C 264, Martel et al. (2000).

It may be possible to resolve this issue, either with increased sample sizes for the optical jet cases, although that seems unlikely given our 100% detection rate for nearby FR-I dust lanes, while more distant objects are much harder to resolve, or else by demonstrating perhaps through spectroscopic observations, that there is a genuine interaction between the 3C 264 dust lane and the outward propagating jet. In the latter case, we would favour a scenario in which more than one physical process can be relevant, however it would remain very curious as to why that one disk, of all disks, is so close to circular!

## 5 Conclusions

We have derived a simple but appealing result: those galaxies with optical jets that are sufficiently nearby to allow a detailed view of the nuclear environments and whose nuclei are not so luminous as to swamp that view, show in 4 out of 5 cases evidence for a circular dust disk. None of the other 21 galaxies in the sample show such a feature. We interpret this as support for a two-fold model in which jets emerge perpendicular to dust lanes and optical emission arises when those jets are pointing towards us, specifically within a critical angle of $`\varphi _c30`$$`40`$ degrees. Nuclear core luminosities show a trend with more luminous cores in the face-on disks, as expected from the beaming scenario.

Support for this work was provided by NASA through proposal number GO-5476.01-93A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.


format: text

text:
# Theory of periodic swarming of bacteria: application to Proteus mirabilis

## I Introduction

To gain an insight into the development and dynamics of various multicellular assemblies, we must understand how cellular interactions build up the structure and result in certain functions at the macroscopic, multicellular level. Microorganism colonies are one of the simplest systems consisting of many interacting cells and exhibiting a non-trivial macroscopic behavior. Therefore, a number of recent studies have focussed on experimental and theoretical aspects of colony formation and the related collective behavior of microorganisms .

The swarming cycles exhibited by many bacterial species, notably Proteus (P.) mirabilis, have been known for over a century . When Proteus cells are inoculated on the surface of a suitable hard agar medium, they grow as short “vegetative” rods. After a certain time, however, cells start to differentiate at the colony margin into long “swarmer” cells possessing up to 50 times more flagella per unit cell surface area. These swarmer cells migrate rapidly away from the colony until they stop and revert by a series of cell fissions into the vegetative cell form, in a process termed consolidation. The resulting vegetative cells grow normally for a time then swarmer cell differentiation is initiated in the outermost zone (terrace), and the process continues in periodic cycles resulting in a colony with concentric zonation depicted in Fig 1. Similar cyclic behavior has been observed in an increasing number of Gram-negative and Gram-positive genera including Proteus, Vibrio, Serratia, Bacillus and Clostridium (for reviews see ).

The reproducibility and regularity of swarming cycles together with the finding that its occurrence is not limited to a single species suggest that periodic swarming phenomena can be understood and quantitatively explained on the basis of mathematical models. In this manuscript we first give an overview of the relevant experimental findings related to the swarming of P. mirabilis, then construct a simple model with two limit densities. We then investigate the behavior of the model as a function of the control parameters and compare it to experimental results.

## II Overview of experimental findings

#### a Differentiation.

As reviewed in , the differentiation of vegetative cells is accompanied by specific biochemical changes. Swarmer cells enhance the synthesis of flagellar proteins, extracellular polysaccharides, proteases and virulence factors, while exhibit reduced overall protein and nucleic acid synthesis and oxygen uptake. These findings may be explained by arguing that the production and operation of flagella is expensive and may require the repression of non-essential biosynthetic pathways. The largely (10-30 fold) elongated swarmer cells develop by a specific inhibition of cell fission which seems not affecting the doubling time of the cell mass or DNA.

The differentiation process is initiated by a number of external stimuli including specific signalling molecules and physico-chemical parameters of the environment. As an example for the latter, the viscosity of the surrounding medium is presumably sensed by the hampered rotation of the flagella . Neither the signal molecules that initiate the differentiation nor the involved intracellular signalling pathways are identified yet, but a corresponding transmembrane receptor has been found recently . The structure of this receptor, together with other findings reviewed in , suggest a ‘quorum sensing’ regulatory pathway  characteristic for many, cell density dependent collective bacterial behaviors like sporulation, luminescence, production of antibiotics or virulence factors.

#### b Migration of swarmer cells.

It is well established  that swarmer cell migration does not require exogenous nutrient sources, since swarmer cells replated onto media devoid of nutrients continue to migrate normally. The ability of migration depends on the local swarmer cell density as isolated single swarmer cells cannot move, while a group of them can. It was also demonstrated  that the mechanism by which bacteria swarm outwards involved neither repulsive nor attractive chemotaxis. The typical swimming velocity (i.e., in liquid environment) of swarmer cells is $`100`$mm/h  which is also their maximal swarming speed and rate of colony expansion on soft agar plates . In the usual experimental conditions for investigating swarming colony formation the front advances with a speed of $`0.510`$ mm/h . Unfortunately, in these cases there is no information available on the velocity of individual swarmer cells, but it must be between the colony expansion speed and the swimming velocity.

#### c Consolidation.

The molecular mechanisms of consolidation, i.e., the downregulation of the gene activity responsible for swarming behavior  is even less known than that of differentiation. If the swarming motility utilizes intracellular energy reserves as has been suggested , then swarmer cells must have a finite lifetime. In addition to the septation of swarmer cells taking place at the outermost terrace, inside the colony the differentiation process, i.e., the supply of fresh swarmer cells must also be shut off. The cessation of swarmer cell production does not seem to be due to severe nutrient depletion since vegetative cells keep growing (although with decreasing growth rate) well inside the colony for many hours after the last swarmer cells were produced in that region .

#### d Colony formation.

The cycle time (total length of the migration and consolidation periods) has been found  to be rather stable ($`4h`$) for a wide range of nutrient and agar content of the medium. The size of the terraces and the duration of the migration phases were strongly influenced (up to an order of magnitude, and up to a factor of 3, respectively) by the agar hardness. The nutrient concentration did not have an observable effect on these quantities from $`0.01\%`$ up to $`1\%`$ glucose concentrations. There was, however a remarkable positive correlation between the cycle time and the doubling time (ranging from $`0.7`$h up to $`1.8`$h) of the cells.

Two interfacing colonies inoculated with a time difference of a few hours and therefore being in different phases of the migration-consolidation cycle, were found to maintain their characteristic phases . Thus, the control of the swarming cycle must be sufficiently local.

#### e Cycle rescheduling.

A few experiments investigated the cell density dependence of the duration of quiescent growth (lag phase) prior to the first migration phase. These studies clearly revealed that vegetative cells have to reach a threshold density to initiate swarming , in accord with the suggested quorum sensing molecular pathway of the initiation of swarmer cell differentiation.

Agar cutting experiments demonstrated that a cut inside the inner terraces does not influence the swarming activity . However, when the cut has been made just behind the swarming front, the duration of the swarming phase was shortened and consolidation was lengthened by up to $`40\%`$ .

Even more interestingly, mechanical mixing of the cell populations before the expected beginning of consolidation expands the duration of the swarming phase considerably, by up to $`50\%`$. This finding, together with replica-printing experiments demonstrates that at the beginning of the consolidation phase still a large pool of swarmer cells exists and seems to be “trapped” at the rear of the outermost terrace.

## III The model

Taking into account the above described experimental findings, here we construct a model which is capable of explaining most of the observed features of colony expansion through swarming cycles.

(i) The model is based on the vegetative and swarmer cell population densities only, denoted by $`\rho ^o`$ and $`\rho ^{}`$, respectively. These values are defined on the basis of cell mass instead of cell number, therefore one unit of swarmer cells is transformed into one unit of vegetative cells during consolidation.

(ii) Vegetative cells grow and divide with a constant rate $`r_01`$h-1 . This will later allow us to establish a direct correspondence between cell density increase and elapsed time.

(iii) Usually, swarmer cell differentiation is initiated when the local density of the vegetative cells exceeds a threshold value ($`\rho _{min}^o10^2`$ cells/$`\mu `$m2 ). (Prior the first swarming phase experiments indicate the presence of an extra time period $`t_{\mathrm{}}`$ which is probably associated with the biochemical changes required to develop the ability of the swarming transition. This effect is present only at the seeding of the colony, thus $`t_{\mathrm{}}=0`$ otherwise.) When $`\rho ^o=\rho _{min}^o`$ at time $`t_0`$, some of the vegetative cells enter the differentiation process, modeled by introducing a rate $`r`$. Since the biomass production rate is assumed to be unchanged during the differentiation process, the rate of producing new vegetative cells is $`r_0r`$ and the differentiating cells elongate with the normal growth rate $`r_0`$.

(iv) The full development of swarmer cells, i.e., a typical 20-fold increase in length needs a time ($`t_d\mathrm{ln}20/r_03`$h) comparable with, or even longer than the length of a consolidation period, hence can not be neglected. Therefore, the first “real” swarmer cells, which are able to move appear only at $`t_0+t_d`$.

(v) The production of swarmer cells is limited in time and the length $`\tau `$ of the time interval during which vegetative cells can enter the differentiation process is another phenomenological parameter of our model.

As we here focus on the periodicity of colony expansion, we do not consider what happens in the densely populated regions after $`t_0+\tau `$. Specifically, in our model any activity of the vegetative cells ceases in these parts of the system.

(vi) Swarmer cells can migrate only if their density exceeds a threshold density $`\rho _{min}^{}`$. Above that threshold, swarmer cells are assumed to move randomly with a diffusion constant $`D_0`$. The finite lifetime of swarmer cells is incorporated into the model through a constant rate ($`r^{}`$) decay.

Unfortunately there are no good estimates on these parameters in the literature. According to recent experimental observations , $`\rho _{min}^{}`$ is less than $`10\%`$ of the value of $`\rho ^o`$ prior the beginning of the swarming phase, i.e., $`0.1\rho _{min}^oe^{r_0t_d}10^2`$ cells/$`\mu `$m2. $`D_0`$ may be estimated as $`v_0^2t_p/2`$ with $`v_0`$ being the typical speed of swarming cells and $`t_p`$ being the persistence time of their motion. As $`v_0`$ is approximately $`30`$mm/h (see Sec II.b) and $`t_p`$ is in the order of minutes, $`D_0`$ is estimated to be in the order of 10 mm2/h.

The above considerations lead to the following set of equations

$`\dot{\rho }^o(t)`$ $`=`$ $`r_0\rho ^o(t)+\mathrm{\Gamma }(t)\mathrm{\Gamma }^{}(t)`$ (1)
$`\dot{\rho }^{}(t)`$ $`=`$ $`\mathrm{\Gamma }(t)+\mathrm{\Gamma }^{}(tt_d)e^{r_0t_d}+D(\rho ^{})\rho ^{}`$ (2)
where the consolidation ($`\mathrm{\Gamma }`$) and differentiation ($`\mathrm{\Gamma }^{}`$) terms are given by

$$\mathrm{\Gamma }=r^{}\rho ^{}\text{ and }\mathrm{\Gamma }^{}=r(\rho ^o,t)\rho ^o.$$
(3)
Eqs. (2) can be significantly further simplified by making use of the possibility to measure elapsed time with the increase in $`\rho ^o`$. In particular, neglecting correction terms related to consolidation in areas where $`\rho ^o>\rho _{min}^o`$, i.e., where primarily differentiation takes place, we can cast $`\mathrm{\Gamma }^{}(tt_d)e^{r_0t_d}`$ in the form

$`r\rho ^o(tt_d)e^{r_0t_d}=r\rho _{min}^oe^{(r_0r)(tt_dt_0)}e^{r_0t_d}=`$ (4)
$`r\rho _{min}^oe^{r_0(tt_0)}e^{r(tt_0t_d)}=r\rho ^o(t)e^{rt_d}`$ (5)
for $`t>t_0+t_dt_1`$. Let us introduce a transformed population density $`\rho `$ as

$$\rho (t)=\{\begin{array}{ccc}\rho ^o(t)& \text{for }\rho ^o<\rho _{min}^o& \text{i.e., for }t<t_0\hfill \\ \rho ^o(t)e^{r(tt_0)}& \text{for }\rho ^o>\rho _0^o& \text{and }t<t_1\hfill \\ \rho ^o(t)e^{rt_d}& \text{for }\rho ^o>\rho _0^o& \text{and }t>t_1\hfill \end{array}$$
(6)
which is in fact the total density of the vegetative and the differentiating, but not yet fully differentiated swarmer cells (see Fig. 2). With this notation, using (5) and similar considerations for $`t_0<t<t_0+t_d`$, Eqs. (2) can be written into a simple, not retarded form

$`\dot{\rho }`$ $`=`$ $`r_0\rho +r^{}\rho ^{}r(\rho )\rho `$ (7)
$`\dot{\rho }^{}`$ $`=`$ $`r^{}\rho ^{}+r(\rho )\rho +D(\rho ^{})\rho ^{},`$ (8)
where

$$r(\rho )=\{\begin{array}{cc}r& \text{for }\rho _{min}<\rho <\rho _{max}\hfill \\ 0& \text{otherwise}\hfill \end{array}$$
(9)
with $`\rho _{min}=e^{rt_d}\rho _{min}10^1`$ cells/$`\mu `$m2 and $`\rho _{max}(\tau )=e^{(r_0r)\tau }\rho _{min}`$. If the $`\mathrm{\Gamma }(r_0r)\rho `$ condition does not hold in the $`[t_0,t_1]`$ time interval, then $`\rho _{min}`$ must be also treated as a dynamical variable. This case will not be considered here. In the following we use $`\rho `$ for the characterization of vegetative cell density.

## IV Results

### A Numerical method in 1D

The model defined through Eqs. (8) has the following seven parameters: the rates $`r_0`$, $`r`$, $`r^{}`$, threshold densities $`\rho _{min}`$, $`\rho _{max}`$ (or $`\tau `$), $`\rho _{min}^{}`$ and the diffusivity $`D_0`$. However, this number can be reduced to four by casting the equations in a dimensionless form using $`1/r_01`$h as time unit, $`\rho _{min}0.1`$ cells/$`\mu `$m2 as density unit and $`x_0=\sqrt{D_0/r_0}3`$mm as the unit length. The resulting control parameters are $`r/r_0`$, $`r^{}/r_0`$, $`\rho _{max}/\rho _{min}=\mathrm{exp}[(r_0r)\tau ]`$ and $`\rho _{min}^{}/\rho _{min}`$. To obtain continuous density profiles, the step-function dependence of $`D`$ on $`\rho ^{}`$ was replaced by

$$D(\rho ^{})=\frac{D_0}{2}\left[1+\mathrm{tanh}2\alpha \frac{\rho ^{}\rho _{min}^{}}{\rho _{min}^{}}\right]$$
(10)
with $`\alpha =10`$ providing a rather steep, but continuous crossover.

Representative examples for the time development of the model are shown in Fig. 3. The production of swarmer cells is localized, and determined by the density profile of vegetative cells at the end of migration periods. In this particular model $`\rho (x)`$ is decreasing towards the colony edge, therefore in the migration phases the source of swarmer cells is moving outwards. The front of swarmer cells is expanding from the inside of the last terrace. Because of the decay term $`\mathrm{\Gamma }`$, cells become non-motile first at the colony edge.

### B Phase diagram

Each of the dimensionless control parameters can have an important effect on the dynamics of the system. As an example, if the duration $`\tau `$ of swarmer cell production is increased, then the consecutive swarming cycles are not separated and a continuous expansion takes place with damped oscillations (Fig. 3b). To map the behavior of the system as a function of the control parameters, the following procedure was applied. Migration periods were identifed by requiring $`\mathrm{max}_x\rho ^{}(x)>\rho _{min}^{}`$. For a given set of parameters we determined the lengths $`\{t_i\}`$ of the consecutive migration periods, and the system was classified as periodic if the three largest values of the set $`\{t_i\}`$ were the same within $`20\%`$. Otherwise, the expansion was classified as continuous as long as $`\mathrm{max}_i\{t_i\}`$ was large enough: comparable with the total duration of the simulated expansion.

The behavior of the model is summarized in Fig. 4, where the boundaries of the various regimes are plotted for three different values of $`r^{}/r_0`$. We found, that $`r`$ and $`\rho _{max}`$ can be combined into one relevant parameter, the swarmer cell production density, as

$$P=_{\mathrm{}}^{\mathrm{}}\mathrm{\Gamma }^{}(x,t)𝑑t=\frac{r}{r_0r}(\rho _{max}\rho _{min}),$$
(11)
which quantity does not depend on the choice of position $`x`$. As the insert demonstrates, for a given $`P`$, the actual values of $`r`$ or $`\rho _{max}`$ are irrelevant to this kind of classification in the parameter regime investigated. The general structure of the phase diagram was found to be similar for various values of $`r^{}`$. For large enough $`P`$ or low enough $`\rho _{min}^{}`$ a continuous expansion takes place, while for too small $`P`$ or large $`\rho _{min}^{}`$ the expansion of the system is finite. For intermediate values of these parameters an oscillating growth develops exhibiting well distinguishable consolidation and migration phases. As the lifetime of the swarmer cells is increased, the parameter regime, in which periodic behavior is exhibited, is shrinked and moved towards lower $`P`$ values.

One can easily estimate the position of the boundary of the non growing phase based on that (i) the width $`w`$ of the terraces is small (this assumption is justified later, in Fig. 7.), thus (ii) the time required for the diffusive expansion of the swarmer cells is much shorter than their lifetime, which, in turn, is (iii) shorter than the duration of a swarming cycle: $`r^{}/r_01`$. The amount of swarmer cells produced in one period is $`Pw`$. Neglecting the decay during expansion, the width $`w^{}`$ of the next, new terrace can be determined from the conservation of cell number as

$$2\rho _{min}^{}w^{}=w(P+\rho _r^{}\rho _{min}^{}),$$
(12)
where $`\rho _r^{}`$ denotes the swarmer cell density remaining from the previous swarming cycle and the symmetric expansion of the released swarmers was also taken into account. To achieve a sustainable growth $`w^{}w`$ is required, resulting in a condition $`3\rho _{min}^{}P+\rho _r^{}`$. If $`\rho _r^{}P`$, as one can expect for $`r^{}/r_01`$, we get for the boundary of the non growing phase

$$P=3\rho _{min}^{},$$
(13)
which, as Fig. 4 demonstrates, is indeed in good agreement with the numerical data.

### C Terrace formation.

The average length $`T`$ of a full swarming cycle was calculated by determining the position of the peak in the power spectrum of $`S(t)=_0^{\mathrm{}}\rho ^{}(x,t)𝑑x`$, the time dependence of the total number of swarmer cells in the system. As Fig. 5 demonstrates, the dimensionless cycle time values are widely spread between values of $`3`$ and $`12`$. However, $`T`$ is only sensitive to changes in $`r_0`$, $`r^{}`$ and $`\rho _{min}^{}/\rho _{min}`$, hence it does not depend on $`\rho _{max}`$ or $`\tau `$.

The average expansion speed $`v`$ and terrace size $`w`$ were also calculated in the parameter regime resulting oscillatory expansion of the colony. First we determined the time $`t_{1/3}`$ when the system reached $`1/3`$ of its maximal simulated expansion $`R_{max}=R(t_{max})`$, with $`R(t)`$ being the position of the expanding colony edge and $`t_{max}`$ is the total duration of the simulation. The average speed was then calculated for the time interval between $`t_{min}=\mathrm{max}(t_{1/3},t_{max}5T)`$ and $`t_{max}`$: for the last $`5T`$ long time interval, or for the last $`2/3`$rd of the total expansion, depending on which was smaller. After obtaining $`v`$ as $`[R_{max}R(t_{min})]/(t_{max}t_{min})`$, the average terrace width was calculated as $`w=vT`$. Fig. 6. shows the dependence of these parameters on the swarmer cell production density $`P`$ and migration density threshold $`\rho _{min}^{}`$. In general, decreasing $`\rho _{min}^{}`$ or increasing $`P`$ results in an increase in both $`w`$ and $`v`$. As Fig. 7. demonstrates, for a given $`r^{}`$, the relevant parameter controlling $`w`$ is $`P/\rho _{min}^{}`$.

The results on the cycle time $`T`$ (Fig. 5.) can be interpreted as follows. As $`\rho _{min}^{}`$ is a good estimate on the density of swarmer cells in the expanding front, at the end of migration phase the vegetative cell density within the new terrace is given by $`r^{}\rho _{min}^{}T^{}`$ with $`T^{}`$ being the duration of the migration phase. Now the length of the consolidation phase, $`TT^{}`$, is determined by the requirement that $`\rho `$ must reach $`\rho _{min}`$:

$$\rho _{min}=r^{}\rho _{min}^{}T^{}e^{r_0(TT^{})}=r^{}\rho _{min}^{}\frac{T^{}}{e^{r_0T^{}}}e^{r_0T}.$$
(14)
As $`T^{}1/r_0`$, the estimate (14) is simplified to

$$r_0T=\mathrm{ln}\frac{r_0\rho _{min}e}{r^{}\rho _{min}^{}},$$
(15)
which gives a rather accurate fit to the numerically determined data (Fig. 5).

### D Lag phase

Since the duration ($`T_L`$) of the lag phase (the time period before the first migration phase) has been in the focus of many recent experiments, now we turn our attention towards this quantity. At least four processes determine $`T_L`$. First, there is a time $`t_{\mathrm{}}`$ associated with the biochemical changes required to switch into swarming mode. As discussed in Sec. II.a, these processes take place only prior the first swarming phase, and are presumably related to sensing the altered environmental conditions. Second, the cell population must reach the threshold density $`\rho _{min}^o`$ (at time $`t_0`$). Third, $`t_d`$ time is required to produce fully differentiated swarmer cells (at time $`t_1`$), and finally, the density of the swarmer cells must reach the migration threshold $`\rho _{min}^{}`$. Let us investigate how these parameters depend on the initial inoculum density $`\rho _0`$.

As $`\rho `$ grows with a rate $`r_0`$ until the appearance of swarmer cells,

$$r_0t_1=\mathrm{max}[t_{\mathrm{}}+t_d,\mathrm{ln}(\rho _{min}/\rho _0)].$$
(16)
The time development of $`\rho ^{}`$ can be estimated by the integration of Eq.(8) with $`r^{}=0`$ (i.e., assuming $`\mathrm{\Gamma }\mathrm{\Gamma }^{}`$) yielding

$$\rho ^{}(t+t_1)=\frac{r}{r_0r}\rho _0e^{r_0t_1}[e^{(r_0r)t}1].$$
(17)
Therefore, $`t^{}=T_Lt_1`$ is given by

$$r_0t^{}=\frac{1}{1r/r_0}\mathrm{ln}\left(\frac{r_0r}{r}\frac{\rho _{min}^{}}{\rho _0e^{r_0t_1}}+1\right),$$
(18)
an expression usually giving a minor correction to $`t_1`$.

Fig. 8a. shows the above calculated $`r_0T_L`$ vs $`\rho _0`$ for $`t_{\mathrm{}}=0`$. The increase in length of the swarmer cells was assumed to be $`20`$-fold, thus $`r_0t_d=\mathrm{ln}203`$, which value can be seen for $`\rho _0\rho _{min}`$. In the opposite limit, when $`\rho _0\rho _{min}`$, we have $`r_0T_L\mathrm{ln}\rho _0+const`$. These relations allow the determination of both $`r_0`$ and $`r`$ (using the known value of $`\rho _{min}/\rho _{min}^{}`$) from the experimental data on $`T_L(\rho _0)`$.

### E Comparison with experiments

Most of the published experimental data are related to the average period length, $`T`$ and terrace size $`w`$. From these parameters the average expansion speed can be calculated as $`v=w/T`$, i.e., $`v`$ is not an independent quantity. As we could see in the previous paragraph, from the density dependence of the lag phase the parameters $`r_0`$, $`t_{\mathrm{}}+t_d`$ and $`r`$ can be estimated. Notice that this estimate on $`r_0`$ is in principle different from the value obtained by the usual methods based on densitometry in liquid cultures. Technically, $`\rho _{min}^{}`$ could be also determined , but such measurements are not published yet.

There are four Proteus strains studied systematically in experiments: the PRM1, PRM2, BB2000 and BB2235 strains (see Table I). To extract the values of the model’s parameters the following procedure was applied. (i) We estimated $`r_0`$ based on lag phase length measurements. (ii) From the calculated $`r_0T`$ values the $`\rho _{min}^{}/\rho _{min}`$ ratio was estimated (assuming $`r^{}/r_0=1`$) based on Eq.(15), see Fig. 5. (iii) Using Eqs. (17) and (18), by a nonlinear fitting procedure (Levenberg-Marquardt method, , see Fig. 8b.) $`\rho _{min}`$, $`r`$, and $`t_{\mathrm{}}+t_d`$ was determined. (The latter value is not relevant in respect the periodicity of the behavior.) (iv) Knowing $`\rho _{min}^{}`$ and $`\rho _{min}`$, from the experimental terrace width data $`x_0`$ and $`P`$ can be estimated using Fig. 6. (v) Finally, $`\tau `$ is given by

$$r_0\tau =\mathrm{ln}\left(\frac{Pr_0}{\rho _{min}r}+1\right)$$
(19)
The parameter values of the model are summarized in Table II, together with the predictions on $`T`$, $`v`$ and $`w`$. An excellent agreement can be achieved with biologically relevant parameter values.

Two classes of model parameters should be distinguished: (a) the ones which are related to the growth and differentiation of the cells ($`r_0`$, $`r^{}`$, $`\rho _{min}`$, and $`P(r,\tau )`$) and (b) those which depend on agar softness ($`D_0`$ and $`\rho _{min}^{}`$). For a given strain we expect that a change in the agar concentration influences only the latter group, while changes in temperature may affect both, but primarily $`r_0`$. In fact, as Table I. demonstrates, by changing $`D_0`$ and keeping all the other, growth-related parameters constant, we could quantitatively reproduce the colony behavior observed on various agar concentrations. Similar statement holds for the temperature effects as well, where the only parameter we changed was the growth rate $`r_0`$.

## V Discussion

Periodic bacterial growth patterns have been in the focus of research in the last few years. Since a colony can be viewed as a system where diffusing nutrients are converted into diffusing bacteria, one may not be surprised by the emergence of spatial structures . However, the periodic patterns of bacterial colonies are qualitatively different from the Liesegang rings (for a recent review see ) developing in reaction-diffusion systems: the spacing between the densely populated areas is uniform and independent of the concentration of the other diffusing species, i.e., the nutrients. The Turing instability is also well-known for producing spatial structures , but in that case the pattern emerges simultaneously in the whole system. It is also known, that bacteria can aggregate in steady concentric ring structures as a consequence of chemotactic interactions , but as we discussed in Sec. II., it is established that swarming of P. Mirabilis does not involve chemotaxis communication. Therefore, none of the well known generic pattern forming schemes can explain the colony structure of swarming bacteria.

As we mentioned in the introduction, oscillatory growth is also exhibited by other bacterial species. One of them, Bacillus subtilis, has been the subject of systematic studies on colony formation and a number of models have been constructed to explain the observed morphology diagram (for recent reviews see ). Only one model addressed the problem of migration and consolidation phases: Mimura et al  set up a reaction-diffusion system in which the decay rate of the bacteria was dependent both on their concentration and the locally available amount of nutrients. The periodic behavior is then a consequence of the following cycle: if nutrients are used up locally, then the bacterial density starts to decay preventing the further expansion of the colony. Nutrients diffuse to the colony and accumulate due to the reduced consumption of the already decreased population. The increased nutrient concentration gradually allows the increase in population density and the expansion of the colony, which starts the cycle from the beginning. While this can be a sound explanation for B. subtilis, as we discussed in Sec. II., the nutrient limitation clearly can not explain neither the differentiation nor the consolidation of P. Mirabilis swarmer cells.

Another recent study  focused on the swarming of Serratia liquefaciens. In that case the structure of the molecular feedback loops are better explored, and were resolved in the model. The production of a wetting agent was initiated by high concentrations of specific signalling molecules. The colony expansion was considered to be a direct consequence of the flow of the wetting fluid film, in which process the only effect of bacteria (besides the aforementioned production) was changing the effective viscosity of the fluid. The wetting agent production was downregulated through a negative feedback loop involving swarmer cell differentiation. This scenario is certainly not applicable to P. Mirabilis, where swarmer cells actively migrate outwards and their role is quite the opposite: enhancing the expansion of the colony.

The first theoretical analysis focusing on P. mirabilis was performed by Esipov and Shapiro (ES) in . Their model was constructed based on assumptions similar to ours, and could reproduce the alternating migration and consolidation phases during the colony expansion. However, the complexity of the ES model involves a rather large number of model parameters, which practically impedes both the full mapping of the parameter space and the quantitative comparison of the model results with experimental findings. The major differences between our and the ES model can be summarized as follows: (i) we do not resolve the age of the swarmer population. Instead, we have a density measure and a constant decay rate implying an exponential lifetime distribution on the (unresolved) level of individual cells. Since the available microbiological observations  suggest only that the lifetime is finite, there is no reason for preferring any specific distribution. (ii) We did not incorporate into our model an unspecified “memory field” with a built-in hysteresis. Instead, we implemented a density-dependent motility of the swarmer cells, which behavior has been indeed observed . (iii) In our model the fully differentiated swarmer cells do not grow, which assumption is probably not fundamental for the reported behavior, but it seems to be more realistic because of the repression of many biosynthetic pathways . Finally, (iv) we do not consider any specific interaction between the motility of swarmer cells and the non-motile vegetative cell population. Although such interactions probably exist, they are undocumented, and as we demonstrated, are not required for the formation of periodic swarming cycles. However, such effects can be important in the actual determination of the density profiles.

With these differences, which are not compromising the biological relevance of the model, we were able to map completely the phase diagram, establish approximate analytical formulas and estimate the value of all model parameters in the case of four different strains. In addition, experimental data measured under various conditions could be explained with one particular parameter setting in the case of the PRM1 strain indicating the predictive power of our approach. Our model is a minimal model in the sense that all of the explicitly considered effects (thresholds, diffusion, etc.) were required to produce the oscillatory behavior, thus, it can not be simplified further. Such minimal models can serve as a comparison baseline for later investigations of various specific interactions.

The values of the microscopic parameters of the model can be either measured directly (like $`r_0`$, $`\rho _{min}^{}`$, $`\rho _{min}`$, $`\rho _{max}`$ or $`r^{}`$) or can be determined indirectly from experimental data (as $`\rho _{min}^o`$ and $`r`$). Most of these measurements have not yet been performed, we hope that our work will motivate such experiments further examining the validity of our assumptions. In fact, one of the parameters, $`r^{}/r_0`$ was set to 1 during the fitting processes, as currently there is no available data to estimate its value. Our numerical results suggest that it is probably larger than $`0.3`$, and it is unlikely to be larger than 2 (meaning an average lifetime less than 30 minutes). Within this range our qualitative conclusions are valid, while the numerical values of the parameter estimates can change up to a factor of $`3`$.

The behavior of “precocious” swarming mutants reported in  deserves special attention. First, we would like to comment on the huge difference found in the value of the transition rate $`r`$ (see Table II). We emphasize that this is not an arbitrary output of a multiparameter fitting process. First, we have reasons to believe, that the motility thresholds of the two BB strains are rather similar. Knowing the growth rates and the cycle times, Eq. (15) shows us that the difference in the values of $`\rho _{min}`$ can not exceed one order of magnitude. Assuming then this maximal difference in $`\rho _{min}`$, $`r`$ remains the only free variable in Eqs. (16)-(18), and the fitting can be performed unambigously. Thus, we are quite confident that such a large difference exists in $`r`$ showing that the rsbA gene (in which these strains differ) influences not only the cell density threshold, but the rate of differentiation as well. It is also interesting to note that in Fig. 8b the behavior of the PRM2 and PRM1 strains reflect a relation very similar to that of the BB2000 and BB2235 strains. Finally, our calculations predicted a slightly longer cycle time for the precocious swarming mutant BB2235, which is also in accord with the actual experimental findings (see Fig 2. of ).

In our model the assumed functional form of the density-dependence of the diffusion coefficient is somewhat different from the most often considered one , namely

$$D(\rho )\rho ^k.$$
(20)
The advantage of (20) is that it allows analytic solutions for certain cases , however, it describes an unlimited, arbitrarily fast diffusion inside the colony where the density is high. In contrast, in real colonies the diffusion of cells is certainly bounded, and the expansion of the boundary can be often limited by the supply of cells from behind . Therefore we believe that our thresholded formulation (10) is a better approximation of what is taking place inside the real colonies.

Finally we would like to comment on the role of nutriens in the swarming behavior of P. mirabilis. In our model there is a phenomenological parameter ($`\tau `$) determining how long the swarmer cells are produced at a given position in the colony. When investigating the dependence of the cycle time on this parameter, as Fig. 5 demonstrates, we found an extremely weak effect. Thus, at least within the framework of this model there is no contradiction between the assumption that the swarmer cell production ceases due to nutrient (or accumulated waste) limitations, and the seemingly nutrient-independent cyclic behavior. In fact, this idea can be developed further. By increasing $`\tau `$ (or decreasing the motility threshold $`\rho _{min}^{}`$) we arrive into a regime where the migration/consolidation phases are not clearly separable as a motile swarmer cell population exists even when the expansion of the colony is slower. Experiments mapping the morphology diagram of P. mirabilis (Fig. 2 of ) showed that there are certain values of agar hardness and nutrient concentration, for which the expansion of the colony is still oscillating, but the periodic density changes are smeared out due to the presence of motile swarmer cells in the consolidation periods. If one associates the increasing agar hardness with increasing $`\rho _{min}^{}`$ and the nutrient concentration with $`\tau `$ then one can qualitatively reproduce those (i.e., the $`P_r`$ and $`P_h`$) regions of the morphology diagram.

###### Acknowledgements.

One of the authors (M.M.) is grateful to T. Matsuyama and H. Itoh for many stimulating discussions on experimental results. This work was supported by funds OTKA T019299, F026645; FKFP 0203/197 and by grants No 09640471 and 11214205 from the Ministry of Education, Science and Culture of Japan.


text:
# Exact Criterion for Determining Clustering vs. Reentrant Melting Behavior for Bounded Interaction Potentials

## I Introduction

The phase behavior of systems whose constituent particles interact by means of pair potentials diverging at the origin is a problem that has been extensively studied in the last few decades. The whole range of inverse power-law pair potentials have been examined, ranging from hard spheres (HS) to the one-component plasma (OCP) and it has been established that excluded volume effects are mainly responsible for bringing about the freezing transition. The crystal structure in which a liquid freezes is subsequently determined by the steepness of the repulsion, with hard repulsions favoring a face centered cubic (fcc) lattice and soft ones a body centered cubic (bcc) lattice mcconnell:etal:prl:93 . Power-law diverging potentials result into freezing at arbitrarily high temperatures. However, the divergence of the potential alone is not enough to cause such a phenomenon, as demonstrated recently by Watzlawek et al. watzlawek:etal:prl:99  who employed a logarithmically divergent pair potential, suitable to describe effective interactions between star polymers in good solvents. It was shown that the strength (prefactor) of the logarithmic potential, determined by the number of arms $`f`$ of the stars, is crucial in determining whether the system freezes or remains fluid at all densities. As a consequence, the phase diagram of star polymers was predicted to display reentrant melting and a critical freezing value $`f_c=34`$ of arms, such that for $`f<f_c`$ the system remains always fluid watzlawek:etal:prl:99 .

Another interesting class of interactions are those which do not diverge at the origin, i.e., they are bounded. Such potentials arise naturally as effective interactions between the centers of mass of soft, flexible macromolecules such as polymer chains ard:peter:00 , dendrimers likos:ballauff , polyelectrolytes etc. Indeed, the centers of mass of two macromolecules can coincide without violation of the excluded volume conditions, hence bringing about a bounded interaction. Moreover, the same mechanisms that exist for tuning the usual, diverging interactions between colloidal particles can be applied in order to tune the bounded interactions: the solvent quality, temperature, chain length, salt concentration etc. will all affect the effective potential. Thus, it appears to be useful to consider such potentials in some generality in order to be able to draw conclusions about the expected phase behavior of systems interacting by means of these.

Two model systems in this category have already been studied in some detail. One is the penetrable spheres model (PSM) marquest:witten:89 ; likos:pensph:98 , in which the pair potential is a positive constant $`\epsilon `$ for distances $`r<\sigma `$ and vanishes otherwise. The other is the Gaussian core model (GCM), introduced in the mid-seventies by Stillinger stillinger:76 . In the GCM, the pair potential $`v(r)`$ has the form $`v(r)=\epsilon \mathrm{exp}[(r/\sigma )^2]`$, with $`\epsilon `$ being an energy and $`\sigma `$ a length scale. It has been shown that the GCM models very accurately the effective interactions between the centers of mass of linear polymer chains olaj:lantschbauer:77 ; grosberg:82 ; schaefer:baumgaertner:86 ; krueger:etal:89 ; dautenhahn:hall:94 ; ard:finken:hansen:99 ; ard:peter:00 ; bolhuis:jcp:00 .

The PSM was studied by means of cell-model calculations and computer simulations likos:pensph:98 , liquid-state integral equation theories fernaud:jcp:00  and density-functional theory schmidt:cecam:99 ; the fluid structure of the PSM has been further studied recently by Rosenfeld et al. by using ideas based on the universality of the bridge functional rosenfeld:etal:pensph . It was found that no reentrant melting takes place because the solid always lowers its free energy by allowing for multiply occupied crystal sites, a mechanism that is called clustering likos:pensph:98 . The clustering mechanism stabilizes therefore the solid at all temperatures. Hence, the topology of the phase diagram of the PSM is similar to that of power-law diverging potentials, when details about the clustering structure of the solids are disregarded. On the other hand, the GCM has been studied even more extensively by means of molecular dynamics simulations stillinger:weber:78 ; stillinger:weber:80 , high-temperature expansions stillinger:jcp:79  and the discovery of exact duality relations in the crystalline state stillinger:prb:79 . Recently, a full statistical-mechanical study of the GCM was performed and it was established that the topology of the phase diagram of the GCM resembles that of star polymers. Freezing and reentrant melting accompanied by an upper freezing temperature were quantitatively calculated lang:etal:gcm . The question that arises, therefore, is the following. Given a nonattractive and bounded pair potential which satisfies the following requirements guaranteeing stability and the existence of the thermodynamic limit ruelle :

1. it is bounded;
2. it is positive definite;
3. it decays fast enough to zero at large separations, so that it is integrable and its Fourier transform exists;

to which topology belongs the phase diagram of the system? In this paper, we present an exact criterion which gives an answer to this question and show representative results for model systems which confirm its validity. The rest of the paper is organized as follows: in section II we present the physical arguments supporting the mean-field theory of the models and in section III we discuss the existence of a spinodal instability in this theory and its implications on the phase bahavior. We present a systematic comparison between theory and simulation in section IV and we draw the generic phase diagrams of such systems in section V. Finally, in section VI we summarize and conclude.

## II The model and the mean-field limit

We will work with a general interaction $`v(r)=\epsilon \varphi (r/\sigma )`$ satisfying the requirements put forward above. Here, $`\epsilon `$ and $`\sigma `$ are an energy and a length, respectively, and $`\varphi (x)`$ is some dimensionless function. The latter does not have to be analytic, i.e., discontinuities in the potential or its derivatives are allowed. Without loss of generality, we assume $`\varphi (0)=1`$. Let us call $`\widehat{\varphi }(Q)=\sigma ^3\stackrel{~}{\varphi }(Q)`$ the dimensionless Fourier transform (FT) of the interaction. For more concreteness (and for the purposes of demonstration) we introduce in addition the family of bounded potentials $`v_\xi (r)`$ depending on a tunable parameter $`\xi `$,

$$v_\xi (r)=\epsilon \frac{1+e^{\sigma /\xi }}{1+e^{(r\sigma )/\xi }},$$
(1)
where $`\xi `$ is a ‘smoothing parameter’ having dimensions of length. The case $`\xi =0`$ recovers the PSM whereas as $`\xi `$ grows the interaction becomes smoother. Due to its resemblance to the Fermi-Dirac distribution, we call this family the Fermi distribution model (FDM). The additional factor $`1+e^{\sigma /\xi }`$ in the numerator of the rhs of eq. (1) ensures that the potential varies from $`\epsilon `$ at $`r=0`$ to zero at $`r\mathrm{}`$, for all $`\xi `$.

We introduce dimensionless measures of temperature and density as

$`t`$ $`=`$ $`{\displaystyle \frac{k_BT}{\epsilon }}=\left(\beta \epsilon \right)^1;`$ (2)
$`\eta `$ $`=`$ $`{\displaystyle \frac{\pi }{6}}\rho \sigma ^3={\displaystyle \frac{\pi }{6}}\overline{\rho },`$ (3)
where $`k_B`$ is Boltzmann’s constant and $`\rho =N/V`$ is the density of a system of $`N`$ particles in the volume $`V`$. We will refer to $`\eta `$ as the ‘packing fraction’ of the system.

The key idea for examining the high temperature and/or high density limit of such model systems is the following. We consider in general a spatially modulated density profile $`\rho (𝐫)`$ which does not vary too rapidly on the scale $`\sigma `$ set by the interaction. At high densities, $`\rho \sigma ^31`$, the average interparticle distance $`a\rho ^{1/3}`$ becomes vanishingly small, and it holds $`a\sigma `$, i.e., the potential is extremely long-range. Every particle is simultaneously interacting with an enormous number of neighboring molecules and in the absence of short-range excluded volume interactions the excess free energy of the system evans:79  can be approximated by a simple mean-field term, equal to the internal energy of the system:

$$F_{\mathrm{ex}}[\rho (𝐫)]\frac{1}{2}𝑑𝐫𝑑𝐫^{}v(|𝐫𝐫^{}|)\rho (𝐫)\rho (𝐫^{}),$$
(4)
with the approximation becoming more accurate with increasing density. Then, eq. (4) immediately implies that in this limit the direct correlation function $`c(r)`$ of the system, defined as evans:79

$$c(|𝐫𝐫^{}|;\rho )=\underset{\rho (𝐫)\rho }{lim}\frac{\delta ^2\beta F_{\mathrm{ex}}[\rho (𝐫)]}{\delta \rho (𝐫)\delta \rho (𝐫^{})},$$
(5)
becomes independent of the density and is simply proportional to the interaction, namely

$$c(r)=\beta v(r).$$
(6)
Using the last equation, together with the Ornstein-Zernike relation hansen:mcdonald , we readily obtain an analytic expression for the structure factor $`S(Q)`$ of the system as

$$S(Q)=\frac{1}{1+\overline{\rho }t^1\widehat{\varphi }(Q)}.$$
(7)
This mean-field approximation (MFA) has been put forward and examined in detail in the context of the Gaussian core model independently by Lang et al. lang:etal:gcm  and by Louis et al. ard:mft:00 . The model is particularly relevant from the physics point of view, due to its connection to the theory of effective interactions between polymer chains ard:peter:00 ; bolhuis:jcp:00 . Here, we establish the validity of the MFA at high densities for bounded, positive-definite interactions in general and we examine its implications for the global phase behavior of such systems.

Bounded and positive-definite interactions have been studied in the late seventies by Grewe and Klein grewe:klein:jmpa:77 ; grewe:klein:jmpb:77 . The authors considered a slightly different model than the one considered here, namely a Kac potential of the form:

$$v(r)=\gamma ^d\psi (\gamma r),$$
(8)
where $`d`$ is the dimension of the space and $`\gamma 0`$ is a parameter controlling the range and strength of the potential. Moreover, $`\psi (x)`$ is a nonnegative, bounded and integrable function:

$$0\psi (x)A<\mathrm{},C=d^d𝐱\psi (x)<\mathrm{}.$$
(9)
Grewe and Klein were able to show rigorously that at the limit $`\gamma 0`$, the direct correlation function of a system interacting by means of the potential (8) is given by eq. (6) above. The connection with the case we are discussing here is straightforward: as there are no hard cores in the system or a lattice constant to impose a length scale, the only relevant length is set by the density and is equal to $`\rho ^{1/3}`$ in our model and by the parameter $`\gamma ^1`$ in model (8). In this respect, the limit $`\gamma 0`$ in the Kac model of Grewe and Klein is equivalent to the limit $`\rho \mathrm{}`$ considered here. However, in the Kac model the strength of the interaction goes to zero simultaneously with the increase in its range. Moreover, the validity of the mean-field expression at large but finite densities and at low temperatures has not been tested in detail.

## III Spinodal instability and freezing

We employ the MFA as a physically motivated working hypothesis for now and, by direct comparison with simulation results, we will show later that it is indeed valid. Within the framework of this theory, an exact criterion can be formulated, concerning the stability of the liquid phase at high temperatures and densities. The function $`\varphi (x)`$ was assumed to be decaying monotonically from unity at $`x=0`$ to zero at $`x\mathrm{}`$. For the function $`\widehat{\varphi }(Q)`$, there are two possibilities: (i) It has a monotonic decay from the value $`\widehat{\varphi }(Q=0)=\sigma ^3𝑑𝐱\varphi (x)>0`$ to the value $`\widehat{\varphi }(Q)=0`$ at $`Q\mathrm{}`$. We call such potentials $`Q^+`$-potentials. Obviously, the Gaussian interaction belongs to this class. (ii) It has oscillatory behavior at large $`Q`$, with the implication that it is a nonmonotonic function of $`Q`$, attaining necessarily negative values for certain ranges of the wavenumber. We call such potentials $`Q^\pm `$-potentials. Long-range oscillations in $`Q`$-space imply that $`\varphi (x)`$ changes more rapidly from unity at $`r=0`$ to zero at $`r\mathrm{}`$ in the $`Q^\pm `$-class than in the $`Q^+`$ one. Moreover, let us call $`Q_{}`$ the value of $`Q`$ at which $`\widehat{\varphi }(Q)`$ attains its minimum, negative value.

If we are dealing with a $`Q^\pm `$-potential eq. (7) implies that $`S(Q)`$ has a maximum at precisely the wavevector $`Q_{}`$ where $`\widehat{\varphi }(Q)`$ attains its negative minimum, $`|\widehat{\varphi }(Q_{})|`$ and this maximum becomes a singularity at the ‘spinodal line’ $`\overline{\rho }t^1|\widehat{\varphi }(Q_{})|=1`$, signaling the so-called Kirkwood instability of the system grewe:klein:jmpb:77 ; klein:grewe:jcp:80 ; klein:brown:jcp:81 ; klein:etal:physica:94 . The theory has a divergence, implying that the underlying assumption of a uniform liquid is not valid and the system must reach a crystalline state. Indeed, on the basis of the fluctuation-dissipation theorem, $`S(Q)`$ can be interpreted as a response function of the density to an infinitesimal external modulating field at wavenumber $`Q`$ hansen:mcdonald  and a diverging value of this response function clearly signals an instability. If the Fourier transform of $`\varphi (x)`$ has negative Fourier components, then an increase in temperature can be compensated by an increase in density in the denominator of eq. (7), so that $`S(Q_{})`$ will have a divergence at all $`t`$. We thus conclude that $`Q^\pm `$-systems freeze at all temperatures.

If we are dealing with a $`Q^+`$-potential ($`\widehat{\varphi }(Q)`$ monotonic), then eq. (7) implies that $`S(Q)`$ is also a monotonic function of $`Q`$ at high densities lang:etal:gcm . For such potentials, one can always find a temperature high enough, so that the assumptions of eq. (4) hold and then eq. (7) forces the conclusion that freezing of the system is impossible at such temperatures. This does not imply, of course, that such systems do not freeze at all; one simply has to go to a low enough temperature and density, so that the mean-field assumption does not hold and the interaction is much larger that the thermal energy. Then, the system will display a hard-sphere type of freezing, to be discussed more explicitly below. An upper freezing temperature $`t_\mathrm{u}`$ must exist for $`Q^+`$-potentials, implying that such systems must remelt at $`t<t_\mathrm{u}`$ upon increase of the density. Hence, we reach the conclusion that $`Q^+`$-systems display an upper freezing temperature and reentrant melting. The criterion says nothing about the crystal structure of the solid, however, which always depends on the details of the interaction as well as the density lang:etal:gcm ; stillinger:stillinger:97 .

For potentials is in the $`Q^+`$-class, the mean-field arguments presented above hold not only at high temperatures but also at low ones, provided that the requirement $`\rho \sigma ^31`$ is satisfied, because these are molten at high densities for all nonzero temperatures. The validity of the mean-field theory for $`Q^+`$-type systems, even at very low temperatures, was confirmed recently by direct comparison with simulation results for the particular case of the Gaussian potential lang:etal:gcm . If the potential is in the $`Q^\pm `$-class, the mean field approximation holds provided that the system is not already frozen, as we will confirm shortly. Moreover, both kinds of systems display an unusual kind of ‘high density ideal gas’ limit. Indeed, taking the expression (7) for $`S(Q)`$ and using the relation $`S(Q)=1+\overline{\rho }\widehat{h}(Q)`$ hansen:mcdonald , where $`\widehat{h}(Q)`$ is the dimensionless Fourier transform of the pair correlation function $`h(r)`$ of the uniform fluid, we obtain:

$$\widehat{h}(Q)=\frac{t^1\widehat{\varphi }(Q)}{1+\overline{\rho }t^1\widehat{\varphi }(Q)}.$$
(10)
At low $`Q`$’s, where $`\widehat{\varphi }(Q)`$ is of order unity, the term proportional to the density in the denominator dominates in the limit of high densities and $`\widehat{h}(Q)`$ scales as $`1/\overline{\rho }0`$. At high $`Q`$’s, the Fourier transform $`\widehat{\varphi }(Q)`$ in the numerator is itself small, with the result that $`\widehat{h}(Q)`$, and hence also the correlation function $`h(r)`$, is approaching zero. This, in turn, means that the radial distribution function $`g(r)=h(r)+1`$ is very close to unity in this limit and deprived of any significant structure for all values of $`r`$ and it only has some small structure at small $`r`$, which is in fact more pronounced for $`Q^\pm `$-potentials than for $`Q^+`$ ones. In this limit, the hypernetted chain (HNC) closure becomes exact, as the exact relation $`g(r)=\mathrm{exp}[\beta v(r)+h(r)c(r)B(r)]`$, combined with the limits $`g(r)1`$, $`h(r)0`$ and $`c(r)\beta v(r)`$ forces the bridge function $`B(r)`$ to vanish. Moreover, eqs. (7) and (10) reveal that the systems obey a scaling law, namely that the functions $`S(Q)`$ and $`th(r)`$ do not depend on $`\overline{\rho }`$ and $`t`$ separately but only on the ratio $`\overline{\rho }/t`$.

Systems in the $`Q^\pm `$-class freeze before the spinodal is reached. In order to make quantitative predictions, we invoke the empirical Hansen-Verlet freezing criterion hansen:verlet:69 ; hansen:schiff:73 , which states that a system crystallizes when $`S(Q)`$ at its main peak attains, approximately, the value $`S(Q_{})=S_\mathrm{m}3`$. Although this criterion was originally put forward for hard, atomic interactions (HS, Lennard-Jones etc.), recent detailed analyses have demonstrated that it holds for the freezing and the remelting transitions of ultrasoft particles such as star polymers watzlawek:etal:prl:99 ; watzlawek:etal:jpcm:98  and even for the nondiverging Gaussian interaction lang:etal:gcm . Hence, we assume that it is valid for the general class of systems we consider here and combining it with eq. (7), we obtain the equation of the freezing line $`t_\mathrm{f}(\eta )`$ as

$$t_\mathrm{f}(\eta )=\frac{6|\widehat{\varphi }(Q_{})|}{\pi (1S_\mathrm{m}^1)}\eta 2.864|\widehat{\varphi }(Q_{})|\eta .$$
(11)
The value $`|\widehat{\varphi }(Q_{})|`$ determines the slope of the freezing line at the high $`(t,\eta )`$ part of the phase diagram.

## IV Comparison with simulations

We now wish to put these arguments in a strong test, using the concrete family of the FDM, given by eq. (1). First of all, we have calculated the Fourier transform of the potential $`v_\xi (r)`$ of eq. (1) numerically, establishing that members of the FDM with $`\xi <\xi _c`$ belong to the $`Q^\pm `$ class and members with $`\xi >\xi _c=0.49697`$ to the $`Q^+`$ one. The GCM is also a member of the latter class that we will discuss in what follows.

### IV.1 Systems displaying clustering

As examples of systems displaying clustering transitions, we have taken the extreme (and by now well-studied case) case $`\xi =0`$ (the PSM) as well as the case $`\xi =0.1`$ of the Fermi distribution model of eq. (1). We have performed standard Monte Carlo (MC) $`NVT`$ simulations for a large number of values for the temperature and density. We begin with the PSM for which the analytical expression (7) takes the form

$$S(Q)=\left[1+24\eta t^1\left(\frac{\mathrm{sin}(Q\sigma )(Q\sigma )\mathrm{cos}(Q\sigma )}{(Q\sigma )^3}\right)\right]^1.$$
(12)
The high temperature-high density freezing line of eq. (11) takes for this choice of $`\xi `$ the form $`t_\mathrm{f}(\eta )=1.033\eta `$. To test the analytical expression of eq. (12), we move along the ‘diagonal’ $`t=\eta `$, a combination that lies almost on the Hansen-Verlet estimate for the location of the freezing line. In Fig. 1 we show the comparison of the analytical results with those obtained from the MC simulations for $`S(Q)`$ and we also demonstrate that the MC curves for the quantity $`th(r)`$ all collapse onto a single line, amply demonstrating the validity of the mean-field approximation for the PSM. In order to further investigate the validity of the MFA, we have performed MC simulations in a variety of thermodynamic points and we present a selection of the obtained results. We present a selection of these in Figs. 2 and 3 and discuss them below.

In Fig. 2(a) we show a comparison between the MC and MFA results for the radial distribution function $`g(r)`$ along the ‘diagonal’ $`t=\eta `$. It can be seen that the agreement between the two is already very good at $`t=\eta =4.0`$ and thereafter it improves markedly with increasing temperature and density. The results obtained from the present theory are of the same quality as those obtained by Fernaud et al. fernaud:jcp:00 , who used the sophisticated zero-separation (ZSEP) closure to investigate the liquid structure of the system. This closure involves three self-consistency parameters, determined in such a way that the virial-compressiblity, Gibbs-Duhem and zero-separation consistency conditions are fulfilled. At the same time, the present results are of the same quality as the recently obtained results of Rosenfeld et al. rosenfeld:etal:pensph , based on ideas of the universality of the bridge functional.

In Fig. 2(b) we perform the same comparison but now at fixed temperature $`t=5.0`$ and increasing packing fraction $`\eta `$. As can be seen, at this temperature, the MFA, which was originally formulated as a high-density approximation, proves to perform extremely well even at intermediate packings, $`\eta =0.5`$ for instance. This is a direct consequence of the boundedness of the interaction combined with a temperature $`t\epsilon `$. Indeed, for small densities, the direct correlation function tends to the Mayer function, $`c(r)\mathrm{exp}\left[\beta v(r)\right]1`$ hansen:mcdonald . If we are dealing with a bounded interaction at high temperature, we can linearize the exponential, obtaining $`c(r)\beta v(r)`$ at low densities, which matches with the MFA expression, eq. (6), at high densities, thus leading to the conclusion that the MFA is an excellent approximation at all densities. For unbounded interactions the linearization of the exponential is evidently impossible.

In Fig. 3 we present a comparison between MC and MFA at fixed packing fraction $`\eta =3.0`$ and increasing temperature. As can be clearly seen, the validity of the MFA improves with increasing temperature. For bounded interactions, an increasing temperature implies a ‘washing-out’ of the correlation effects caused by the (increasingly weak) interaction effects and a tendency of the system towards the particular ‘high-density ideal gas’ limit characterized by the tendency of the function $`g(r)`$ towards unity. However, it is an interesting peculiarity of these systems that unlike the usual ideal gas, the limit $`g(r)1`$ (or, equivalently, $`h(r)0`$) does not imply a corresponding limit $`S(Q)1`$. Though the Fourier transform of $`h(r)`$, $`\widehat{h}(Q)`$, tends to zero as $`\rho ^1`$, this is compensated by the large density $`\rho `$, so that the structure factor $`S(Q)=1+\rho \widehat{h}(Q)`$ displays the signature of strong ordering through the pronounced peaks seen in Figs. 1(b) and 3(b).

Further, we performed MC simulations at selected points deeply inside the region $`t<t_\mathrm{f}(\eta )`$, finding that the obtained structure factors displayed Bragg peaks and hence confirming the prediction that the system is frozen there. Putting all our results together, we draw in Fig. 4 a semi-quantitative phase diagram of the PSM, accompanied by an assessment of the validity of the MFA at selected thermodynamic points. The MFA appears to be an excellent approximation at all densities above the temperature $`t=3.0`$. Hence, we take as an estimate for the freezing line above $`t=3.0`$ the MFA-Hansen-Verlet line $`t_\mathrm{f}=1.033\eta \eta `$; for lower temperatures, we simply connect the point $`(\eta ,t)=(3.0,3.0)`$ with the point $`(\eta ,t)=(0.5,0)`$, which obtains from the consideration that at $`t=0`$ the PSM reduces to the hard sphere system which is known to freeze at a fluid density $`\eta _{\mathrm{HS}}0.5`$. The monotonic shape of the freezing curve for low temperatures arises from detailed considerations there, which can be found in Ref. likos:pensph:98 .

Next we present in Fig. 5 a comparison for the FDM with $`\xi =0.1`$. For this choice of $`\xi `$, the Hansen-Verlet-based freezing line takes the form $`t_\mathrm{f}=0.712\eta `$. The selected points lie in the fluid region and the comparison indicates once more the excellent accuracy of the MFA both for $`g(r)`$ and for $`S(Q)`$. The radial distribution function $`g(r)`$ of this model is deprived of the jump at $`r=\sigma `$ seen in the PSM; the latter is caused by the discontinuity of the PSM potential there. However, a similarity between the $`g(r)`$’s of the $`\xi =0`$ and $`\xi =0.1`$ models is that they both attain their maximum values at full overlaps between the particles $`r=0`$ and thereafter they decay rapidly, featuring a depletion region around $`r\sigma `$. This is a characteristic pointing to a strong clustering property in the fluid phase, a property thereafter inherited by the incipient thermodynamically stable crystal; the number of particles ‘sitting on top of each other’ and thereby occupying the same crystal site scales linearly with density. In order to corroborate this claim, we can argue in two different ways, using the liquid as a reference point.

First, let us consider the number of particles $`N_\mathrm{c}`$ in the fluid phase whose centers are, on average, within a distance $`\sigma `$ from a given particle. The number $`N_\mathrm{c}`$ is given by the formula:

$$N_\mathrm{c}=1+4\pi \rho _0^\sigma r^2g(r)𝑑r.$$
(13)
In Fig. 6 we show the function $`4\pi r^2g(r)`$ within a particle diameter $`\sigma `$ for a sequence of points along the freezing line of the PSM. As all these curves tend to a common limit with increasing density, the integral $`4\pi _0^\sigma r^2g(r)𝑑r`$ tends to a constant and hence $`N_\mathrm{c}\rho `$ at high densities, where the second term on the rhs of eq. (13) dominates.

Second, we can use the wavevector $`Q_{}`$ at which the fluid structure factor has a maximum in order to estimate for the lattice constant $`a`$ of the incipient crystal through the relation $`aQ_{}^1`$. For the models at hand, this maximum is entirely determined by the pair potential; unlike in usual fluids featuring diverging interactions, for which $`Q_{}`$ scales as $`\rho ^{1/3}`$, in our case $`Q_{}`$ knows nothing about the density. Thus, all post-freezing crystals have the same lattice constant, although their average density is a linear function of the temperature. This clearly shows that clustering must take place in the crystal: by allowing more and more particles to occupy the same lattice sites, a practically constant effective density of clusters is maintained in the crystal, thus leading to a density-independent lattice constant.

### IV.2 Systems displaying reentrant melting

We now turn our attention to the opposite case, namely pair potentials belonging to the $`Q^+`$-class. As an example within the FDM family, we have taken the model with parameter $`\xi =0.6`$ and performed a comparison between MC and MFA results. A characteristic example is shown in Fig. 7. As can be seen in Fig. 7(a), unlike the case of $`Q^\pm `$-class potentials, the radial distribution function is completely deprived of any structure, although the thermodynamic parameters are in the same regime as those presented in Figs. 2, 3 and 5. In fact, in the present case, $`g(r)`$ has a minimum at $`r=0`$, not a maximum. This complete lack of structure is reflected in the shape of $`S(Q)`$, shown in Fig. 7(b).

These characteristic features for the $`Q^+`$-class are not an artifact of the relatively high temperature chosen in the results of Fig. 7. They persist even at extremely low temperatures, provided the density is high enough. This has been amply demonstrated recently for the case of the Gaussian core model, another member of the $`Q^+`$-class lang:etal:gcm . In order to stress this point we present in Fig. 8 the $`g(r)`$ and $`S(Q)`$ of the GCM at $`t=0.01`$ and $`\eta =6.0`$. Though $`g(r)`$ displays some structure up to $`r2\sigma `$, the structure factor $`S(Q)`$ shows no signature of some kind of ordering <sup>1</sup><sup>1</sup>1The small discrepancies between MFA and MC for $`g(r)`$ at small $`r`$ values disappear as the density is increased. Perfect agreement between theory and simulation can be achieved if, instead of the MFA, one employs the HNC or some more refined version of the MFA to calculate the fluid structure, see Refs. lang:etal:gcm  and ard:mft:00 .. At any arbitrarily small but finite temperature, a high enough density can be found for which the MFA is valid and then the assumption that a uniform phase exists leads consistently to a fluid which has ideal-gas behavior, i.e., vanishingly small correlations. These liquids are different from usual ideal gases in that, e.g., their pressure $`P`$ and isothermal compressiblity $`\chi _T`$ scale respectively as $`P\rho ^2`$ and $`\chi _Tt\rho ^2`$. Nevertheless, they are thermodynamically stable. Hence, for potentials in the $`Q^+`$-class, the equilibrium phase for sufficiently high densities at arbitrarily small but finite temperatures is the uniform fluid.

## V Generic phase diagrams

We now turn to the opposite limit of the low temperature-low density part of the phase diagram. There, following the original ideas of Stillinger stillinger:76 , a HS mapping can be performed, as follows. The Boltzmann factor $`\mathrm{exp}[\beta v(r)]`$ of the potential varies monotonically from the value $`\mathrm{exp}(\beta \epsilon )0`$ (since $`\beta \epsilon 1`$ there) at $`r=0`$ to unity at $`r\mathrm{}`$ and has a close resemblance to that of a hard sphere system. We can thus define an effective hard sphere diameter $`\sigma _{\mathrm{HS}}`$ through the relation:

$$\mathrm{exp}[\beta v(\sigma _{\mathrm{HS}})]=1/2.$$
(14)
Writing $`v(r)=\epsilon \varphi (r/\sigma )`$ and using the fact that $`\varphi (x)`$ is a monotonic function in order to establish that the inverse function $`\varphi ^1(x)`$ exists, we can rewrite eq. (14) as:

$$\sigma _{\mathrm{HS}}=\sigma \varphi ^1(t\mathrm{ln}2).$$
(15)
We now use the known fact hard spheres freeze at $`\eta _{\mathrm{HS}}0.5`$ together with eq. (15) above in oder to obtain the low temperature-low density freezing line of the system as:

$$t_\mathrm{f}(\eta )=\frac{1}{\mathrm{ln}2}\varphi \left(\left(2\eta \right)^{1/3}\right).$$
(16)
As the limit $`\varphi (x)0`$ is attained for $`x\mathrm{}`$ only, it follows that the low temperature-low density freezing line of the systems goes to $`\eta =0`$ at $`t_\mathrm{f}=0`$. Eq. (16) is valid for all potentials we consider here; however, for $`Q^+`$-potentials, combining the HS-like freezing at low temperatures and densities with the fact that at high densities the fluid has to be stable, derived in the preceding section, we can draw the conclusion that such systems must display reentrant melting and an upper freezing temperature.

We have now taken eq. (16) for the low-$`t`$ and low-$`\eta `$ freezing line of the FDM and combined it with the analytic expression at the opposite limit, eq. (11), in order to draw schematically the evolution of the phase diagram of the FDM as a function of $`\xi `$, for $`\xi <\xi _c`$. The results are shown in Fig. 9. With increasing $`\xi `$, the slopes of the high-$`t`$ freezing lines decrease; at the limit $`\xi 0`$, corresponding to the PSM, the low-$`t`$ freezing line approaches the horizontal axis vertically, as is dictated by the fact that the PSM becomes equivalent to the HS system there likos:pensph:98 . In the inset of Fig. 9, we show the phase diagram for a system with $`\xi =0.6>\xi _c`$, showing reentrant melting behavior. The evolution of the phase diagram from a clustering to a reentrant melting behavior can be easily visualized from this picture.

Finally, it is important to point out that Stillinger has proven that any system interacting by means of a potential which (i) is differentiable at least four times, (ii) vanishes strongly enough at infinity to be integrable and (iii) is +1 at the origin, will inevitably lead to a reentrant melting phase diagram under the assumption that the competing crystal structures have single lattice site occupancy stillinger:76 . We can therefore now complete the statement and say that if a potential belongs to the $`Q^+`$-class, then it will freeze into crystals of single occupancy and then remelt upon increase of the density. But if it belongs to the $`Q^\pm `$-class, then it will freeze into a clustered solid at any temperature. Clustering appears therefore to be the crucial mechanism for crystal stabilization in these systems.

## VI Summary and concluding remarks

To summarize, we have established a criterion for the topology of the phase diagrams resulting from repulsive, bounded interactions, which is very simple in its formulation and states that if the Fourier transform of the pair potential is positive-definite, then the system shows reentrant melting but if not then it freezes at all temperatures, into clustered crystals with multiply occupied sites. We have also established that at temperatures exceeding the interaction strength the mean-field theory is reliable at all densities and its accuracy improves quickly with increasing temperature. We close with the remark that there is a certain similarity between the ideas put forward here and the considerations on freezing for systems featuring of diverging interactions at infinite spatial dimensions frisch:etal:prl:85 ; klein:frisch:jcp:86 ; frisch:percus:pra:87 ; bagchi:rice:jcp:88 . However, in the latter case, the direct correlation function is given by the Mayer function $`f(r)=\mathrm{exp}[\beta v(r)]1`$ of the interaction potential and not by $`\beta v(r)`$ as in the case at hand.

###### Acknowledgements.

This work was supported by the Österreichische Forschungsfond under Project Nos. P11194-PHY and P13062-TPH. A. L. acknowledges financial support by the Deutsche Forschungsgemeinschaft within the SFB 237.