text:
# Hidden Sector Supergravity Breakdown

## 1 Introduction

In the early 1980’s considerations of the hierarchy problem in the standard model of strong and electroweak interactions (as well as grand unified theories) strengthened the interest in applications of supersymmetry to particle physics. I was at SLAC at that time and I remember being intrigued by the idea that the scale of dynamically broken supersymmetry  would explain the weak mass scale of the standard model and induce the breakdown of $`SU(2)\times U(1)`$. Still this goal was not so easy to achieve. After the first excitement of model building in the framework of global supersymmetry, problems started to appear. One of them was the problem of the cosmological constant that was strictly positive in spontaneously broken supersymmetry. In addition, the tree level supertrace formulae  made it difficult to consistently push up the superpartners to a high mass scale and there was no known sensible mechanism to generate gaugino masses. A new symmetry (R-symmetry) was needed to assure proton stability, but the necessary breakdown of that symmetry led to problems with unacceptable axions. It soon became clear that realistic models had to be quite complicated and that supersymmetry breakdown should be somewhat remote from the standard model itself.

With this in mind I decided that the right thing to do was to learn supergravity. After all gravity existed in nature and should be included sooner or later in the framework of the standard model. It was particularly imppressive for me that the local version of supersymmetry included gravity automatically (with the gravitino as gauge particle). In fall 1981 I was moving to CERN as a fellow and this was the right place to learn supergravity. The papers to read were  and two of the authors, Sergio Ferarra and Luciano Girardello were at CERN.

In this talk I would like to present the origin and status of hidden sector supergravity in three steps. The suggestion of the scheme  came in spring 1982 and led to hectic activities, so that in the beginning of 1983 the picture of what is today called “gravity mediation” was pretty complete. A second important step of the development came in 1985 after the “string revolution” initiated by Green and Schwarz . Here the “hidden $`E_8`$-sector” of the heterotic string played a crucial role and the moduli fields of the $`d=10`$ dimensional string theory entered the stage. This will be discussed in chapter 3. More recently, with the advances of string dualities and the formulation of M-theory  some of the remaining questions could be clarified. The hidden sector now became a hidden wall and this version of the scheme will be presented in chapter 4.

## 2 1982: Hidden Sector Supergravity

In a supersymmetric framework we wanted to understand the breakdown scale of $`SU(2)\times U(1)`$ (thus the W- and Z-mass) as a consequence of supersymmetry breakdown ($`M_\mathrm{S}`$). The hierarchy problem was rephrased: why is $`M_\mathrm{S}`$ so small compared to the grand unified scale $`M_{\mathrm{GUT}}`$ or the Planck scale $`M_{\mathrm{Planck}}`$? A satisfactory scheme would have been a dynamical of breakdown supersymmetry, such as fermion-antifermion condensation or the condensation of gauge fermion bilinears: gaugino condensation. Unfortunately, developments in 1981 led to the conclusion, that in the simplest model, globally supersymmetric pure SU(N) gauge theory, supersymmetry was not broken by the formation of gaugino condensates .

Meanwhile my study of supergravity made progress and I could handle some simple models. As an exercise I considered the effective Langrangian for pure supersymmetric Yang-Mills theory as formulated by Veneziano and Yankielowicz  in its local supersymmetric version. Some modifications had to be made to fine-tune the cosmological constant. The result was surprising: supergravity was broken in the presence of a gaugino condensate . Apart from the Polonyi superpotential  we had a second example of local supersymmetry breakdown with a vanishing cosmological constant.

But how could that be consistent with the apperent no-go theorem in the globally supersymmetric theory? Let us assume that $`\mathrm{\Lambda }`$, the renormalization group invariant scale of the gauge interaction sets the scale for gaugino condensation, thus

$$<\chi \chi >\mathrm{\Lambda }^3$$
(1)
where $`\chi `$ denotes the gaugino. The supersymmetry breakdown scale $`M_{\mathrm{SUSY}}`$ (as a vacuum expectation value of an auxiliary field) was here given by

$$M_{\mathrm{SUSY}}^2\frac{\mathrm{\Lambda }^3}{M_\mathrm{P}}$$
(2)
where $`M_\mathrm{P}`$ denotes the Planck scale.

Thus everything was consistent, in the global limit $`M_\mathrm{P}\mathrm{}`$ one obtains $`M_{\mathrm{SUSY}}0`$ and unbroken supersymmetry. In spontaneously broken supergravity the gravitino becomes massive:

$$m_{3/2}\frac{\mathrm{\Lambda }^3}{M_\mathrm{P}^2}$$
(3)
This was very exciting since we had the possibility of dynamically broken supergravity in a “hidden sector” with an additional Planck mass suppression for the gravitino mass. What would that mean in case of coupling to the “observable sector” (the supersymmetric extension of the standard model). In the Polonyi case this had not been studied yet: in any case one expected (at that time!) Str$`M^2=0`$ for the chiral multiplets and (in addition) no convincing way to generate gaugino masses (i.e. for the gauginos in the observable sector).

In the present case, however, there was a way to generate gaugino masses in the observable sector. Supergravity contains 4-fermion interactions

$$\frac{h(\varphi )}{M_\mathrm{P}^2}(\chi \chi )(\chi \chi ),$$
(4)
where $`h(\varphi )`$ is a function of scalar fields $`\varphi `$.

Imagine now a mixing term with one pair of observable sector gauginos (like the gluinos, winos, etc.) and a second pair of hidden sector gauginos with a nonvanishing vacuum expectation value. With this one obtains soft supersymmetry breaking gaugino masses in the observable sector and a completely satisfactory model . Masses for squarks and sleptons appear via radiative corrections from the soft gaugino masses. The scale of soft supersymmetry breaking parameters is determined by the gravitino mass (with potential suppression factors). A detailed discussion can be found in ref. . In this way the general picture of hidden sector supergravity breakdown (or gravity mediated supersymmetry breakdown) emerged. One needs two sectors:

* the observable sector including quarks, leptons, gauge particles and Higgs bosons as well as their supersymmetric partners and
* a hidden sector coupled to the observable sector through interactions of gravitational strength.

The hidden sector provides supersymmetry breakdown at an intermediate scale and leads to a gravitino mass at a rather low scale. Soft breaking terms of order of the gravitino mass (or a few orders of magnitude smaller) appear in the observable sector and lead to a phenomenological satisfactory version of the supersymmetric extension of the standard model of strong and electroweak interactions.

Up to this point the supergravity action had only been constructed for one chiral multiplet. A more detailed analysis of the low energy effective action would require a generalization of this work. Such a generalization was subsequently presented by Cremmer, Ferrara, Girardello and van Proyen  stressing in particular that the supertrace formulae of global supersymmetry was no longer valid in this case. Chamseddine, Arnowit and Nath  quoted similar results in the framework of a grand unified theory, while Barbieri, Ferrara and Savoy  presented a supersymmetric extension of the standard model with soft scalar mass terms, exploiting the new Str$`M^2`$ formulae in local supersymmetry. Srednicki, Wyler and myself  completed this analysis of soft SUSY breaking terms with among others the introduction of the A-parameters. The relevance of the gauge kinetic function ($`f`$) for the generation of soft gaugino masses was pointed out in ref. . A general comprehensive discussion was given in . So within nine months the full picture had been completed. We now had a consistent scheme with (dynamically) broken supersymmetry and a vanishing cosmological constant but without tachyonic instabilities or harmful axions. Gaugino masses and A-parameters were generated naturally: a consequence of the breakdown of continuous R-symmetry to discrete R-parity. Problems with flavour changing neutral currents (the so-called flavour problem) could be avoided by assuming certain degeneracies of soft scalar mass terms, which appeared automatically in case of gaugino masses as the dominant source of supersymmetry breakdown.

Still today “Hidden Sector Supergravity” appears as the most promising d=4 set-up for a low energy effective theory that extends the standard model and avoids the hierarchy problem.

Remains the question of how this fits into “Grand Scenario”, i.e. Strings and/or extra dimensions.

## 3 1984: Towards Higher Dimensions and String Theory

In their classic paper  on anomaly cancellation in $`S0(32)`$ superstring theory, Green and Schwarz initiated a vast research effort in that field. Very soon new string constructions, such as the heterotic string theories  appeared and gave hope for an application in particle physics. Especially the heterotic $`E_8\times E_8`$ theory seemed very promising in this respect. Compactification of six dimensions on a Ricci flat manifold was shown to lead to supersymmetric models in $`d=4`$ dimensions that could serve as candidates for a supersymmetric extension of the standard model . In its simplest form (the so-called standard embedding) this led to a theory with gauge group $`E_6\times E_8`$ and quarks and lepton superfields in the 27-dimensional representation to $`E_6`$. Thus all the observable fields were contained in the $`E_6`$ sector.

With Jean Pierre Derendinger and Luis Ibanez at CERN, we were excited about this fact. We had now a candidate theory which

* naturally contained an observable ($`E_6`$) sector and a hidden ($`E_8`$) sector, that were
* connected only via interactions of gravitational strength, and
* in the simplest model the hidden sector consisted of a pure supersymmetric gauge theory.

This was exactly the picture of hidden sector supergravity we have discussed in the previous section: supersymmetry is potentially broken by gaugino condensation in the hidden sector. It naturally appeared in the framework of superstring theories. Still, it was necessary to clarify the details. For this one needed a simple way to arrive at an at least qualitative description of the low energy effective action. Such a method was supplied in the scheme of reduction and truncation as defined in ref.  which we will describe shortly. With this machinery at hand we could explicitely discuss the question of supersymmetry breakdown in this set-up, as given in ref. . Similar results were found independently in ref. . What we see is that the messenger fields that are responsible for the mediation of supersymmetry breakdown are a combination of the string theory dilaton and moduli fields that describe the size and shape of the extra compactified dimensions. In the simple case considered here they are given by the moduli fields $`S`$ and $`T`$ that will be defined below. We shall be quite brief in our discussion here and refer the reader to ref.  for more details. For a more recent review see ref. .

.

### 3.1 Supersymmetry breakdown in heterotic string theory

In our review of supersymmetry breakdown in the $`d=10`$ weakly coupled $`E_8\times E_8`$ theory we shall start from the $`d=10`$ effective field theory and go to $`d=4`$ dimensions via the method of reduction and truncation explained in ref. . In string theory compactified on an orbifold this would describe the dynamics of the untwisted sector. We retain the usual moduli fields $`S`$ (the dilaton) and $`T`$ (the overal size of compactified extra dimensions) as well as matter fields $`C_i`$ that transform nontrivially under the observable sector gauge group. In this approximation, the Kähler potential is given by

$$G=\mathrm{log}(S+\overline{S})3\mathrm{log}(T+\overline{T}2C_i\overline{C}_i)+\mathrm{log}\left|W\right|^2$$
(5)
with superpotential

$$W(C)=d_{ijk}C_iC_jC_k$$
(6)
and the gauge kinetic function is given by the dilaton field

$$f=S.$$
(7)
We assume the formation of a gaugino condensate $`<\chi \chi >=\mathrm{\Lambda }^3`$ in the hidden sector where $`\mathrm{\Lambda }`$ is the renormalization group invariant scale of the confining hidden sector gauge group. The gaugino condensate appears in the expression for the auxiliary components of the chiral superfields

$$F_j=\left(G^1\right)_j^k\left(\mathrm{exp}(G/2)G_k+\frac{1}{4}f_k(\chi \chi )\right)+\mathrm{}$$
(8)
which are order parameters for supersymmetry breakdown. Minimizing the scalar potential we find $`F_S=0`$, $`F_T0`$ and a vanishing cosmological constant. Supersymmetry is broken and the gravitino mass is given by

$$m_{3/2}=\frac{F_T}{T+\overline{T}}\frac{\mathrm{\Lambda }^3}{M_P^2}$$
(9)
and $`\mathrm{\Lambda }=10^{13}`$ GeV would lead to a gravitino mass in the TeV – range. A first inspection of the soft breaking terms in the observable sector gives a surprising result. They vanish in this approximation. Scalar masses are zero because of the no–scale structure  in (5) (coming from the fact that we have only included fields of modular weight $`1`$ under T–duality in this case ). In a more general situation we might expect scalar masses $`m_0`$ comparable to the gravitino mass $`m_{3/2}`$ and the above result $`m_0=0`$ would just be an artifact of the chosen approximation at the classical level. Thus soft scalar masses are strongly model dependent. Gaugino masses $`m_{1/2}`$ are easier to discuss. They are given by

$$m_{1/2}=\frac{\frac{f}{S}F_S+\frac{f}{T}F_T}{2\mathrm{R}\mathrm{e}f}$$
(10)
and with $`f=S`$ and $`F_S=0`$ we obtain $`m_{1/2}=0`$. One loop corrections will change this picture as can be seen already by an inspection of the Green–Schwarz anomaly cancellation counter terms, as they modify $`f`$ at one loop. In the simple example of the so–called standard embedding with gauge group $`E_6\times E_8`$ we obtain

$$f_6=S+ϵT;f_8=SϵT.$$
(11)
This dependence of $`f`$ on $`T`$ will via (10) lead to nonvanishing gaugino masses which, however, might be small compared to $`m_{3/2}`$ and $`m_0`$ since $`ϵT`$ is considered a small correction to the classical result. This might be problematic when applied to the class of supersymmetric extension of the standard model, where $`m_0`$ is large compared to $`m_{1/2}`$. With the large difference of the soft scalar and the gaugino masses a sizeable fine tuning is needed to induce the breakdown of electroweak symmetry at the correct scale . The smallness of the gaugino masses might also lead to a problem in the context of relic abundances of the lightest superparticles (LSP) . Ideally one would like to consider again a model where the main source of supersymmetry breaking is given by the gaugino masses.

### 3.2 The role of the antisymmetric tensor field

There are a few things that should be retained from this discussion that are not so transparent in the short presentation given above. In particular it is the role of the two index antisymmetric tensor field $`B_{MN}`$ that is important. In $`d=10`$ supergravity the (3 index) field strength $`H_{MNP}`$ is not just the curl of $`B`$ but has to be supplied with so-called Chern Simons terms .

We now want to stress the following points that will also be important in the discussion in the next section.

* The superpotential of the $`d=4`$ effective action originates in the co-called Yang-Mills-Chern-Simons term $`\omega ^{\mathrm{YM}}`$.
* The Green Schwarz counterterms of anomaly cancellation, will provide (among others) a one-loop corection to the gauge kinetic function, that is classicaly just given by $`S`$
* In $`d=4`$ the real part of $`S`$ is given by $`g^2`$, thus the gauge coupling constant $`g`$ is a field dependent quantity.
* In the $`d=10`$ dimensional theory the field strength $`H`$ and the gaugino bilinears $`(\chi \chi )`$ appear in a complete square, such that at the classical level the potential is positive definite.
* As a result of this a nonvanishing vacuum expectation value of $`H`$ will trigger supersymmetry breakdown and the dilaton will adjust its vacuum expectation value such that the cosmological constant vanishes.
* A vanishing value of $`<H>`$ will then lead to a situation where $`S`$ runs to infinity with a vanishing gauge coupling constant.

### 3.3 Open questions

Certainly the potential runaway of the dilaton $`S`$ is an open problem in this scenario (as in any theory based on perturbative string theory ), and up to today no really satisfactory solution has been proposed. Other problems concern the exact form of the scalar potential. While at the classical level the vacuum energy vanishes, this is certainly not obvious in the quantum theory. As a result of this, it is quite difficult to obtain a statement about even the qualitative form of the potential in a model independent way. This, of course, implies that the soft supersymmetric breaking parameters (like squark and slepton masses as well as A-parameters) are strongly model dependent. In contrast, the situation for the gaugino masses is much simpler (as usual) and maybe they could even be dominant in a realistic framework. Their origin in the simple set-up is a direct consequence of the Green-Schwarz counterterms in the low-energy effective action. In perturbative string theory they can be obtained through a one loop calculation  that agrees with the former result in the limit of large $`T`$.

In a next step, one would like to use more explicit input from string theory. One approach would be to exploit the presence of stringy symmetries , like $`T`$-duality, that should be valid in perturbative string theory. In fact, this picture of $`T`$-duality can be described successfully for the classical effective action . There appear, however, problems when one incorporates loop effets, as now the $`S`$ and $`T`$-moduli mix. A successful implemenation of $`T`$-duality then seems to require an understanding of $`S`$-duality as well. $`S`$-duality, however, is not an explicit symmetry of the perturbative low energy effective action, and this is probably the reason, why this approach can not be completed . The problem of the runaway dilaton, and thus the determination of the vacucum expextation value of $`S`$ seems to require a nonperturbative treatment and no really convincing scheme has emerged. As we have seen above this problem might be intimately related to the origin of a vacuum expectation value of the three-index tensor field strengh $`<H>`$ and the apearence of Chern-simons terms. Still, a better understanding of string symmetries would be the road to a better control of the low energy effective action.

## 4 1995: Towards $`d=11`$ and (Heterotic) M–theory

A new step towards the understanding of nonperturbative dualities in string theory came in 1995 with the formulation of the concept of M-theory  and its embedding in $`d=11`$ dimensional space time. For our purposes we are especially interested in the $`E_8\times E_8`$ version as constructed by Horava and Witten . The new aspect of this picture is the fact that now hidden and observable sector become separated geometrically. Thus hidden sector supergravity breakdown now becomes the breakdown of supergravity at a hidden brane, and again gaugino condensation is the preferred mechanism as suggested by Horava  and Olechowski, Yamaguchi and myself . Early related work has been done in

Let us now reconsider these questions in the strongly coupled $`E_8\times E_8`$$`M`$–theory in somewhat more detail. The effective action is given by (for details see )

$`L={\displaystyle \frac{1}{\kappa ^2}}{\displaystyle d^{11}x\sqrt{g}\left[\frac{1}{2}R+\mathrm{}\right]}+\mathrm{}`$ (12)
$`+{\displaystyle \frac{1}{2\pi (4\pi \kappa ^2)^{2/3}}}{\displaystyle d^{10}x\sqrt{g}\left[\frac{1}{4}F_{AB}^aF^{aAB}+\mathrm{}\right]}.`$
Compactifying to $`d=4`$ we obtain  (with the correction pointed out in ref. )

$`G_N`$ $`=`$ $`8\pi \kappa _4^2={\displaystyle \frac{\kappa ^2}{8\pi ^2V\rho }},`$ (13)
$`\alpha _{GUT}`$ $`=`$ $`{\displaystyle \frac{(4\pi \kappa ^2)^{2/3}}{V}}`$ (14)
with $`V=R^6`$ and $`\pi \rho =R_{11}`$. Fitting $`G_N`$ and $`\alpha _{GUT}=1/25`$ then gives $`R_{11}M_{11}`$ of order 10 and $`M_{11}R2.3`$. The rather large value of the $`d=4`$ reduced Planck Mass $`M_P=\kappa _4^1`$ is obtained as a result of the fact that $`R_{11}`$ is large compared to $`R`$.

We now perform a compactification using the method of reduction and truncation as above. For the metric we write

$$g_{MN}=\left(\begin{array}{ccc}e^\gamma e^{2\sigma }g_{\mu \nu }& & \\ & e^\sigma \delta _{mn}& \\ & & e^{2\gamma }e^{2\sigma }\end{array}\right)$$
(15)
with $`M,N=1\mathrm{}11`$; $`\mu ,\nu =1\mathrm{}4`$; $`m,n=5\mathrm{}10`$; $`2R_{11}=2\pi \rho =M_{11}^1e^\gamma e^\sigma `$ and $`V=e^{3\sigma }M_{11}^6`$. At the classical level this leads to a Kähler potential as in (5)

$$K=\mathrm{log}(𝒮+\overline{𝒮})3\mathrm{log}(𝒯+\overline{𝒯}2C_i\overline{C}_i)$$
(16)
with

$`𝒮`$ $`=`$ $`{\displaystyle \frac{2}{\left(4\pi \right)^{2/3}}}\left(e^{3\sigma }\pm i24\sqrt{2}D\right),`$ (17)
$`𝒯`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{\left(4\pi \right)^{4/3}}}\left(e^\gamma \pm i6\sqrt{2}C_{11}\right)`$ (18)
where $`D`$ and $`C_{11}`$ fields are defined by

$`{\displaystyle \frac{1}{4!}}e^{6\sigma }G_{11\lambda \mu \nu }`$ $`=`$ $`ϵ_{\lambda \mu \nu \rho }\left(^\rho D\right),`$ (19)
$`C_{11i\overline{j}}`$ $`=`$ $`C_{11}\delta _{i\overline{j}}`$ (20)
and $`x^i`$ ($`x^{\overline{j}}`$) is the holomorphic (antiholomorphic) coordinate of the Calabi–Yau manifold. The imaginary part of $`𝒮`$ (Im$`𝒮`$) corresponds to the model independent axion, and the gauge kinetic function is $`f=𝒮`$. This is very similar to the weakly coupled case. Before drawing any conclusion from these formulae, however, we have to discuss a possible obstruction at the one loop level. It can be understood from the mechanism of anomaly cancellation . For the 3–index tensor field $`H`$ in $`d=10`$ supergravity to be well defined one has to satisfy $`dH=\mathrm{tr}F_1^2+\mathrm{tr}F_2^2\mathrm{tr}R^2=0`$ cohomologically. In the simplest case of the standard embedding one assumes $`\mathrm{tr}F_1^2=\mathrm{tr}R^2`$ locally and the gauge group is broken to $`E_6\times E_8`$. Since in the M–theory case the two different gauge groups live on the two different boundaries of space–time such a cancellation point by point is no longer possible. We expect nontrivial vacuum expectation values (vevs) of

$$(dG)\underset{i}{}\delta (x^{11}x_i^{11})\left(\mathrm{tr}F_i^2\frac{1}{2}\mathrm{tr}R^2\right)$$
(21)
at least on one boundary ($`x_i^{11}`$ is the position of $`i`$–th boundary). In the case of the standard embedding we would have $`\mathrm{tr}F_1^2\frac{1}{2}\mathrm{tr}R^2=\frac{1}{2}\mathrm{tr}R^2`$ on one and $`\mathrm{tr}F_2^2\frac{1}{2}\mathrm{tr}R^2=\frac{1}{2}\mathrm{tr}R^2`$ on the other boundary. This might pose a severe problem since a nontrivial vev of $`G`$ might be in conflict with supersymmetry ($`G_{11ABC}=H_{ABC}`$). The supersymmetry transformation law in $`d=11`$ reads

$$\delta \psi _M=D_M\eta +\frac{\sqrt{2}}{288}G_{IJKL}\left(\mathrm{\Gamma }_M^{IJKL}+\mathrm{}\right)\eta +\mathrm{}$$
(22)
Supersymmetry will be broken unless e.g. the derivative term $`D_M\eta `$ compensates the nontrivial vev of $`G`$. Witten has shown  that such a cancellation can occur and constructed the solution in the linearized approximation (linear in the expansion parameter $`\kappa ^{2/3}`$) which corresponds to the large $`T`$–limit in the weakly coupled theory<sup>*</sup><sup>*</sup>*For a discussion beyond this approximation in the weakly coupled case see refs. ..

The supersymmetric solution can be shown to lead  to a nontrivial dependence of the $`\sigma `$ and $`\gamma `$ fields with respect to $`x^{11}`$:

$$\frac{\gamma }{x^{11}}=\frac{\sigma }{x^{11}}=\frac{\sqrt{2}}{24}\frac{d^6x\sqrt{g}\omega ^{AB}\omega ^{CD}G_{ABCD}}{d^6x\sqrt{g}}$$
(23)
where the integrals are over the Calabi–Yau manifold and $`\omega `$ is the corresponding Kähler form. Formula (23) contains all the information to deduce the effective $`d=4`$ supergravity theory, including the Kähler potential and the gauge kinetic function . A definition of our $`𝒮`$ and $`𝒯`$ fields in the four–dimensional theory now requires an average over the 11–dimensional interval. We therefore write

$`𝒮`$ $`=`$ $`{\displaystyle \frac{2}{\left(4\pi \right)^{2/3}}}\left(e^{3\overline{\sigma }}\pm i24\sqrt{2}\overline{D}\right),`$ (24)
$`𝒯`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{\left(4\pi \right)^{4/3}}}\left(e^{\overline{\gamma }}\pm i6\sqrt{2}\overline{C}_{11}\right)`$ (25)
where bars denote averaging over the 11th dimension. It might be of some interest to note that the combination $`𝒮𝒯^3`$ is independent of $`x^{11}`$ even before this averaging procedure took place.

$`\mathrm{exp}(3\sigma )`$ represents the volume of the six–dimensional compact space in units of $`M_{11}^6`$. The $`x^{11}`$ dependence of $`\sigma `$ then leads to the geometrical picture that the volume of this space varies with $`x^{11}`$ and differs at the two boundaries. In the given approximation, this variation is linear, and for growing $`R_{11}`$ the volume on the $`E_8`$ side becomes smaller and smaller. At a critical value of $`R_{11}`$ the volume will thus vanish and this will provide us with an upper limit on $`R_{11}`$. For the phenomenological applications we then have to check whether our preferred choice of $`R_{11}`$ that fits the correct value of the $`d=4`$ Planck massWith $`V`$ depending on $`x^{11}`$ we have to specify which values should be used in eq. (14). The appropriate choice in the expression for $`G_N`$ is the average value of $`V`$ while in the expression for $`\alpha _{GUT}`$ we have to use $`V`$ evaluated at the $`E_6`$ border. satisfies this bound. Although the coefficients are model dependent we find in general that the bound can be satisfied, but that $`R_{11}`$ is quite close to its critical value. A choice of $`R_{11}`$ much larger than $`(\mathrm{few}\times 10^{15}\mathrm{GeV})^1`$ is therefore not permitted.

This variation of the volume is the analogue of the one loop correction of the gauge kinetic function (11) in the weakly coupled case and has the same origin, namely a Green–Schwarz anomaly cancellation counterterm. In fact, also in the strongly coupled case we find , as a consequence of (23), corrections for the gauge coupling constants at the $`E_6`$ and $`E_8`$ side.

Gauge couplings will no longer be given by the (averaged) $`𝒮`$–field, but by that combination of the (averaged) $`𝒮`$ and $`𝒯`$ fields which corresponds to the $`𝒮`$–field before averaging at the given boundary:

$$f_{6,8}=𝒮\pm \alpha 𝒯$$
(26)
at the $`E_6`$ ($`E_8`$) side respectivelyWith the normalization of the $`𝒯`$ field as in (25), $`\alpha `$ is a quantity of order unity.. The critical value of $`R_{11}`$ will correspond to infinitely strong coupling at the $`E_8`$ side $`𝒮\alpha 𝒯=0`$ (Notice the similarity to (11) in the weakly coupled limit). Since we are here close to criticality a correct phenomenological fit of $`\alpha _{\mathrm{GUT}}=1/25`$ should include this correction $`\alpha _{\mathrm{GUT}}^1=𝒮+\alpha 𝒯`$ where $`𝒮`$ and $`\alpha 𝒯`$ give comparable contributions. This is a difference to the weakly coupled case, where in $`f=S+ϵT`$ the latter contribution was small compared to $`S`$. Observe that this picture of a loop correction $`\alpha 𝒯`$ to be comparable to the tree level result still makes sense in the perturbative expansion, since $`f`$ does not receive further perturbative corrections beyond one loop .

### 4.1 Supersymmetry breakdown in M–theory

In a next step we are now ready to discuss the dynamical breakdown of supersymmetry via gaugino condensation in the strongly coupled M–theory picture. In analogy to the previous discussion we start investigating supersymmetry transformation laws in the higher–dimensional (now $`d=11`$) field theory :

$`\delta \psi _A`$ $``$ $`D_A\eta +G_{IJKL}\left(\mathrm{\Gamma }_A^{IJKL}8\delta _A^I\mathrm{\Gamma }^{JKL}\right)\eta `$ (27)
$``$ $`\delta (x^{11})\left(\overline{\chi }^a\mathrm{\Gamma }_{BCD}\chi ^a\right)\mathrm{\Gamma }_A^{BCD}\eta +\mathrm{}`$
$`\delta \psi _{11}`$ $``$ $`D_{11}\eta +G_{IJKL}\left(\mathrm{\Gamma }_{11}^{IJKL}8\delta _{11}^I\mathrm{\Gamma }^{JKL}\right)\eta `$ (28)
$`+`$ $`\delta (x^{11})\left(\overline{\chi }^a\mathrm{\Gamma }_{ABC}\chi ^a\right)\mathrm{\Gamma }^{ABC}\eta +\mathrm{}`$
where gaugino bilinears appear in the right hand side of both expressions. It can therefore be expected that gaugino condensation breaks supersymmetry. Still the details have to be worked out. In the $`d=10`$ example, the gaugino condensate and the three–index tensor field $`H`$ contributed to the scalar potential in a full square. This has led to a vanishing cosmological constant as well as the fact that $`F_S=0`$ at the classical level. Hořava has observed  that a similar mechanism might be in operation in the $`d=11`$ theory After a careful calculation this leads to a vanishing variation $`\delta \psi _A=0`$. In our model (based on reduction and truncation) we can now compute these quantities explicitly. We assume gaugino condensation to occur at the $`E_8`$ boundary

$$\overline{\chi }^a\mathrm{\Gamma }_{ijk}\chi ^a=\mathrm{\Lambda }^3ϵ_{ijk}$$
(29)
where $`\mathrm{\Lambda }<M_{\mathrm{GUT}}`$ and $`ϵ_{ijk}`$ is the covariantly constant holomorphic 3–form. This leads to a nontrivial vev of $`G_{11ABC}`$ at this boundary and supersymmetry is broken<sup>§</sup><sup>§</sup>§One might speculate that a nontrivial vev of $`D_A\eta `$ might be operative here as in the case without gaugino condensation (see discussion after eq. (22)). However, the special values of $`H_{ijk}ϵ_{ijk}`$ necessary to cancel the contribution of the gaugino condensate do not permit such a mechanism (see footnote 6 in ref. ).. At that boundary we obtain $`F_𝒮=0`$ and $`F_𝒯0`$ as expected from the fact that the component $`\psi _{11}`$ of the 11–dimensional gravitino plays the role of the goldstino.

In the effective $`d=4`$ theory we now have to average over the 11th dimension leading to

$$F_𝒯\frac{1}{2}𝒯\frac{𝑑x^{11}\delta \psi _{11}}{𝑑x^{11}}$$
(30)
as the source of SUSY breakdown.

This will then allow us to compute the size of supersymmetry breakdown on the observable $`E_6`$ side. Gravitational interactions play the role of messengers that communicate between the two boundaries. This effect can be seen from (14): large $`R_{11}`$ corresponds to large $`M_P`$ and $`F_𝒯`$ gives the effective size of SUSY breaking on the $`E_6`$ side ($`R_{11}\mathrm{}`$ implies $`M_P\mathrm{}`$). The gravitino mass is given by

$$m_{3/2}=\frac{F_𝒯}{𝒯+\overline{𝒯}}\frac{\mathrm{\Lambda }^3}{M_P^2}$$
(31)
(similar to (9) in the weakly coupled case) and we expect this to represent the scale of soft supersymmetry breaking parameters in the observable sector . These soft masses are determined by the coupling of the corresponding fields to the goldstino multiplet. As we have seen before, we cannot compute the scalar masses reliably. In the naive approximation used here we obtain $`m_0=0`$ because of the no–scale structure . This might be an artifact of the approximation. Fields of different modular weight will receive contribution to $`m_0`$ of order $`m_{3/2}`$. For the mass of a field $`C`$ we have

$$m_0^2=m_{3/2}^2F^i\overline{F}^{\overline{j}}\frac{Z_{i\overline{j}}Z_iZ^1\overline{Z}_{\overline{j}}}{Z}$$
(32)
where $`i,j=𝒮,𝒯`$ and $`Z`$ is the moduli dependent coefficient of $`C\overline{C}`$ term appearing in the Kähler potential. Scalars of modular weight $`1`$ will become massive through radiative corrections. This then leads to the expectation that $`m_{3/2}`$ should be in the TeV–region and $`\mathrm{\Lambda }10^{13}`$ GeV In realistic models $`E_8`$ is broken and $`\mathrm{\Lambda }`$ is adjusted by model building.. So far this is all similar to the weakly coupled case.

### 4.2 Gaugino masses and WIMPs

An important difference appears, however, when we turn to the discussion of observable sector gaugino masses (10). In the weakly coupled case they were zero at tree level and appeared only because of the radiative corrections at one loop (11). As a result of this small correction, gaugino masses were expected to be much smaller than $`m_{3/2}`$. In the strongly coupled case the analog of (10) is still valid

$$m_{1/2}=\frac{\frac{f_6}{𝒮}F_𝒮+\frac{f_6}{𝒯}F_𝒯}{2\mathrm{R}\mathrm{e}f_6}$$
(33)
and the 1–loop effect is encoded in the variation of the $`\sigma `$ and $`\gamma `$ fields from one boundary to the other. Here, however, the loop corrections are sizable compared to the classical result because of the fact that $`R_{11}`$ is close to its critical value. As a result we expect observable gaugino masses of the order of the gravitino mass. The problem of the small gaugino masses does therefore not occur in this situation. Independent of the question whether $`F_𝒮`$ or $`F_𝒯`$ are the dominant sources of supersymmetry breakdown, the gauginos will be heavy of the order of the gravitino mass. Thus it could very well be that gaugino masses are the dominant sources of supersymmetry breakdown in the observable sector. The exact relation between the soft breaking parameters $`m_0`$ and $`m_{1/2}`$ will be a question of model building. If in some models $`m_0m_{1/2}`$ this might give a solution to the flavor problem. The no–scale structure found above might be a reason for such a suppression of $`m_0`$. As we have discussed above, this structure, however, is an artifact of our simplified approximation and does not necessarily survive in perturbation theory. At best it could be kept exact (but only for the fields with modular weight $`1`$) in the $`R_{11}\mathrm{}`$ limit. The upper bound on $`R_{11}`$ precludes such a situation. With observable gaugino masses of order $`m_{3/2}`$ we also see that $`m_{3/2}`$ cannot be arbitrarily large and should stay in the TeV – range .

The question of the soft gaugino masses is, of course, very important for the determination of the mass of the lightest supersymmetric particle (LSP). This might be stable and provide a candidate for the dark matter found in the universe (a so-called WIMP, for weakly interacting massive particle). Given the present situation of the supersymmetric extension of the standard model, the LSP is most likely a gaugino of the wino- or bino-type. The relic abundance of these neutralinos in the early universe depends crucially on its mass and its coupling to ordinary matter. The mass is predominantly determined through the soft gaugino masses, while the interactions are mediated by the scalar partners of leptons. Too high a mass of the scalar partners leads then to a weak coupling of the LSP to ordinary matter, reducing the annihilation of WIMPs in the early universe; with the potential danger of an overcritical WIMP density. For this relic density of WIMPs to be both acceptable and cosmologically interesting, only certain ranges of soft scalar masses and gaugino masses are allowed. It turns out, that the heterotic M-theory could lead to such an acceptable scenario (in contrast to the the situation in the weakly coupled case). A detailed discussion can be found in ref. .

### 4.3 The Super–Higgs mechanism

In the previous sections the estimate of the gravitino mass was obtained  using a the simplified approximation according to which the higher dimensional bulk fields were integrated out via an averaging proceedureA corresponding analysis in global supersymmetry has been performed in ref. .. In this picture, the goldstino mode was represented by the lowest Kaluza–Klein $`\mathrm{\Psi }_0`$ mode of a higher dimensional field $`\mathrm{\Psi }`$. In the super–Higgs mechanism this mode supplies the additional degrees of freedom to render the gravitino massive. Qualitatively this simplified approximation does give a consistent picture, but there remain some open questions and potential problems when one looks into details of the super–Higgs mechanism. In this section we would like to point out these potential problems and show how they can be resolved. For a detailed explanation we refer the reader to . We shall concentrate here on a qualitative discussion of the mechanism involved. This will include a discussion of the possible nature of the goldstino (is it a bulk or a wall field), the relation to the Scherk-Schwarz mechanism  in that context  and an upper limit for the gravitino mass in the present picture. We shall argue that a meaningful realization of the super–Higgs mechanism seems to require some modes in the bulk other than the graviton and the gravitino. Finally we shall comment on the phenomenological consequences of this findings, including a discussion of the nature of the soft breaking terms on both walls.

Specifically we want to address the following two questions:

* the nature of the massless gravitino in the presence of several $`F`$terms<sup>\**</sup><sup>\**</sup>\**We generically use the notation $`F`$term for the source of supersymmetry breakdown. Depending on the specific situation this could represent a $`D`$term or a gaugino condensate as well. on different walls that cancel and lead to unbroken supersymmetry,
* the identification of the goldstino in the case of broken supersymmetry.

The first question (i) arises because of a particalur nonlocal effect of supersymmetry breakdown first observed by Horava . A given source of supersymmetry breakdown (parametrized by a vacuum expectation value (vev) of an auxiliary field $`F`$) on one wall could be compensated by a similar but opposite value $`(F)`$ on another (separated) wall. Any calculation and approximation of the system thus has to reproduce this behaviour. The previously mentioned averaging proceedure over the bulk distance does this in a trivial way, leading to unbroken supersymmetry as expected. A detailed inspection of the gravitino, however, reveals a problem. If we start with the situation $`F=0`$ it is easy to define the massless gravitino $`\mathrm{\Psi }_0`$ in the $`d=4`$ theory. Switching on a nontrivial $`F`$ on one brane and $`(F)`$ on the other still should give a massless gravitino, but $`\mathrm{\Psi }_0`$ turns out to be no longer a mass eigenstate. The resolution of this problem and the correct identification of the gravitino can be found in . It is a particular combination of the possible gravitini that appear when one, for example, reduces a 5-dimensional theory to a theory in $`d=4`$ on a finite $`d=5`$ interval. The theory on a $`d=5`$ circle would lead to $`N=2`$ supersymmetry in $`d=4`$ and two massless gravitini (zero modes on the circle). The $`Z_2`$ projection on the interval removes one of the gravitini and is $`N=1`$ supersymmetric. A nonvanishing vev of $`F`$ now interferes with the boundary conditions and the massless gravitino will be a linear combination of the zero mode and all the excited Kaluza-Klein modes whose coefficients will depend on $`F`$ (assuming, of course, unbroken supersymmetry due to a compensating vev $`F`$ on another wall).

The second question (ii) deals with the nature of the goldstino (i.e. the longitudinal components of the gravitino) in the case of broken supersymmetry. Remember that the simplified averaging proceedure leads to a goldstino that corresponds to the lowest Kaluza–Klein mode $`\mathrm{\Psi }_0`$ of a higher-dimensional bulk field $`\mathrm{\Psi }`$. Inspecting the gravitino mass matrix in this case reveals the fact that this field $`\mathrm{\Psi }_0`$ is not a mass eigenstate, but mixes with infinitely many higher Kaluza–Klein modes $`\mathrm{\Psi }_n`$. A consistent manifestation of a super–Higgs mechanism would require a diagonalization of this mass matrix and an identification of the goldstino. This problem has been solved in , by a suitable redefinition of the Kaluza-Klein modes. Thus it is shown that a consistent Super-Higgs mechanism is at work.

## 5 The nature of hidden wall supersymmetry breakdown

The resolution of these puzzles clarifies some of the other questions of the approach. It also allows us to generalize this to string theories of type I with D-branes and supersymmetry breakdown on distant branes. The qualitative picture, of course, is very similar.

* The nonlocality of the breakdown shows some resemblance to the breakdown of supersymmetry via the Scherk–Schwarz  mechanism. Here, however, the real goldstino of the spontaneous breakdown of supersymmetry can be unambiguously identified.
* The possibility to cancel the supersymmetry breakdown on a distant wall by a vev on the local wall tells us, that the mass splittings of broken supersymmetry have to be of order of the gravitino mass $`m_{3/2}`$ on both walls.
* In terms of the physical quantities there is no real extra suppression, once we separate the walls by a large distance $`R`$. In the limit $`R\mathrm{}`$ we will have $`M_{\mathrm{Planck}}\mathrm{}`$ as well. The suppression of the soft breaking parameters will always be gravitational.
* In general, when we have a system of many separated branes with potential sources of supersymmetry breakdown, the actual breakdown will be obtained by the sum of these contributions. The averaging proceedure will be very useful to decide whether supersymmetry is broken or not. The identification of the goldstino, however, is more difficult and requires a careful calculation.
* A successful implementation of the super–Higgs mechanism will require some fields other than gravitino and graviton in the bulk<sup>††</sup><sup>††</sup>††Usually they arise as modes of the higher dimensional supergravity multiplet.. This implies that in the absence of such fields a consistent spontaneous breakdown of supergravity might not be achieved.

We thus have a consistent scenario of supersymmetry breakdown in the framework of heterotic M-theory. It provides us with explicit formulae of soft supersymmetry breaking terms in the low energy effective theory. Interestingly enough these formulae differ from those obtained in the weakly coupled heterotic $`E_8\times E_8`$ theory, most notably with respect to the size of soft gaugino masses. The phenomenological properties of the scenario, therefore have to be reconsidered in various circumstances. This might include a discussion of the mass and nature of the LSP, the question of the universality of the soft scalar masses as well as possible new solutions to the strong CP problem. More work needs to be done to work out these questions in detail. In summary we can say that the picture of supersymmetry breakdown in the M–theoretic limit looks very promising. It is very similar to the weakly coupled case, but avoids the problem of the small gaugino masses. This has important consequences for the phenomenological and cosmological properties of the effective models in four dimensions.

## 6 Outlook

Since its original formulation in 1982, gravity mediation appeared as a very attractive scheme to introduce spontaneous (dynamical) supersymmetry breakdown into the supersymmetric extension of the $`SU(3)\times SU(2)\times U(1)`$ standard model of strong, weak and electromagnetic interactions. Its natural implementation in (heterotic) string theory and M-theory stengthens the confidence in this mechanism of hidden sector supergravity breakdown. It gives us useful hints about the possibility of a local cancellation of the classical vacuum energy, the role of anomaly cancellation via the Green-Schwarz counterterms, the importance of underlying string symmetries and the specific role of the antisymmetric tensor fields in string- and M-theory.

Still there remain some open questions about the possible incorporation of $`T`$\- and $`S`$-duality, the microscopic origin of the vacucum expectation value of the $`<H>`$ and $`<G>`$ fields and their potential to avoid the run-away problem of the dilaton. The control of the cosmological constant beyond the classical level remains a problem and this makes it difficult to determine the soft scalar mass terms in a model independent way. Soft gaugino masses are easier to handle since their size seems to be controlled by the mechanism of anomaly cancellation. The most attractive scheme emerges in the case where gaugino masses are the dominant source of supersymmetry breakdown.


format: text

text:
# Dimension-splitting for simplifying diffusion in lattice-gas models

## I Introduction

In computer modeling of physical systems, there is not always a direct correspondence between simple models and computationally efficient ones. What seems like a direct and intuitive description to a human does not necessarily correspond to an algorithm that performs well on a computer. Simplicity and directness can, however, be valuable assets regardless of whether they result in better performance. Simple models are easier to analyze and to translate into correct computer programs. In addition, complex models composed of well-understood simple component models are easy to design and test. Of course what we would really like is to combine simplicity and performance.

Molecular dynamics—numerically simulating the classical or semiclassical behavior of a collection of atoms or molecules—is an example of a direct technique with wide applicability that is conceptually simple. The technique is often used in situations where a computational solution based on differential equations is infeasible due to the complexity of the physical dynamics or the complexity of the boundary conditions. Often, in molecular dynamics simulations, the microscopic dynamics of the particles is simplified in order to enable more computationally efficient and hence larger-scale simulations, as in the dissipative particle dynamics and lattice-gas approachesCov-diss-prtcl ; BU-lga-meeting . In a lattice-gas automaton (LGA), for example, only a small number of particle velocities are allowed, and particles are confined to the sites of a finite spatial lattice. Steps in which particles concurrently hop between lattice sites alternate with steps in which the particles present at each lattice site interact. This simplification of the microscopic details results in models that combine simplicity and performance, while still capturing the full essence of phenomena of interest.

LGA-like models have long been studied in physics as theoretical models which capture some essence of a physical phenomenon. The classic example is the random walk on a lattice, which has been used as a microscopic model of diffusionWeiss-rwalks . The first true LGA models were introduced for the theoretical study of fluid flowHPP-pra76 . It was realized only much later that the locality, uniformity and spatial regularity of models based on LGA’s make them ideal candidates for large-scale simulation on parallel hardwareFHP-prl86 ; MTV-prl86 ; cam8-waterloo . Since then, LGA simulations have been used to study a variety of physical phenomena, including fluid dynamics, chemical reactions, and changes in phaseRothZales-book ; Rivet-Boon-lgas-book ; Dieter-lgas-book ; Kapral-prl91 ; LGA-reactive-phys-report96 . More generally, lattice dynamics which emulate the locality and uniformity of physical dynamics occupy a borderland between abstract theoretical models and practical algorithms: they are often bothcambook ; Smith-phd ; Dsouza-phd ; Margolus-feynlec . The computer model is the same as the theoretical model, unlike models based on partial differential equations (PDE’s), for example, where an intervening stage of numerical analysis needed. Thus lattice models and related modeling techniques have both a theoretical and a practical aspect, providing simple models that we can directly experiment with, visualize and analyze mathematically.

In this paper, we introduce a simple technique for incorporating diffusive phenomena into LGA and other lattice models. This technique is closely related to the classical random walk model, but differs in important respects. First, it incorporates exclusion: only a limited number of particles can occupy each lattice site. The technique we use in this manuscript to implement diffusion with exclusion was first introduced by Toffolicambook . A second difference involves randomness: it is explicitly provided by concurrently simulating a different invertible LGA as a source of stochasticity. This approach has also been used before in other parallel LGA modelscam8-waterloo ; Smith-phd . Finally, perhaps the most important distinction from a classical random walk model is that this technique splits each $`n`$-dimensional simulation time-step up into $`n`$ fractional time-steps, updating only one dimension during each fractional step. This is similar to some dimension-splitting and fractional-time-step techniques used in solving PDE’sAmes ; prng-ran2 , but here we apply it in the context of a lattice-gas with exact particle conservation. Earlier Cellular Automata (CA) models have used dimension-splittingcam6-rotation ; Fredkin-note , but not for simulating diffusion, and not as a general modeling technique. For a review of other techniques for simulating diffusion with discrete lattice systems see Ref. Band-cadiff .

LGA models of complex materials and phenomena can be constructed in a modular fashion by combining simpler models—something which is difficult to do with differential equation based models. Qualitative models are particularly easy to construct in this manner, but realistic quantitative models of physical phenomena can also be constructed. For example, to reproduce known chemical behavior, several subsystems can be combined. Each subsystem would consist of independently diffusing chemical species, and interaction between subsystems (e.g, chemical reaction steps) would occur in between independent species-diffusing stepsKapral-prl91 ; LGA-reactive-phys-report96 . The diffusing steps can be handled with a dimension-splitting approach, resulting in a very simple treatment of the diffusion component of the model, which requires only a small amount of state information at each lattice site. A smaller amount of state not only allows larger simulations to be achieved with given memory resources, but also makes it easier to couple component subsystems together. The uniform treatment of one dimension at a time provides an additional benefit. Consider, for example, a chemical reaction model of the sort described above. Such a model can be developed and tested as a two-dimensional model and then converted into a three-dimensional model by simply applying the same dynamics to an additional dimension.

In general, LGA models have a different domain of applicability than PDE-based models. Since they are based on a kind of molecular dynamics optimized for exact digital computation on massively parallel hardware, LGA’s can be constructed to handle complex interactions and complex boundaries—circumstances poorly suited to PDE-based simulations. For example, running a fluid LGA simulation through an MRI reconstruction of a real porous rock involves no more computational work than running it with any other boundary condition. An example illustrating the power of the LGA approach for complex systems with complicated boundaries is a lattice-gas model of amphiphilic fluid flow in porous mediaCov-etal-porous .111Note that in Sec. (II.1.1) we introduce the lattice-Boltzmann equation (LBE) for our model and discuss some of the distinctions between the LGA and the LBE approaches.

In this paper we first review a version of Toffoli’s LGA model of diffusion with exclusion. We follow an analysis due to Chopard and DrozChopDroz-diff ; ChopDroz-book  to show explicitly how to start with the microscopic model and recover the appropriate macroscopic diffusion equation in the continuum limit. We then introduce a dimension-splitting technique for use in LGA models that involve diffusion. We show explicitly how to map this dimension-reduced model onto Toffoli’s LGA model. We then discuss the origins of the conservation laws present in the model. Finally we illustrate the use of this technique as a component for building two composite LGA models: a microscopically reversible model of diffusion-limited aggregation, and a model of growth of biological films. These two models were developed in two dimensions but have also been run and studied in three dimensions.

## II Diffusion

Modeling diffusion is a cornerstone for building general models of a wide range of physical phenomena. Examples include models for reaction-diffusion processes, pattern formation, and heat flowcambook ; BLS-transportbook ; Walgraef-patformbook . Models of diffusive transport alone encompass such diverse phenomena as flow through porous mediaRothZales-book , flow of fluids in biological systemsBerg-rwalks , and flow in geological processesRothZales-book ; BLS-transportbook .

Two standard approaches to modeling diffusion come from opposite extremes: the continuum limit, where diffusion is modeled by a partial differential equation; and the microscopic limit where diffusion is modeled as a collection of particles undergoing random walks (which is meant to capture the phenomenon of Brownian motion). One-dimensional diffusive phenomena were first described from a macroscopic perspective by FickCrank . The generalization to multiple dimensions is the well known diffusion equation:

$$\frac{\rho (\stackrel{}{x},t)}{t}=D^2\rho (\stackrel{}{x},t).$$
(1)
Here $`\rho (\stackrel{}{x},t)`$ is density as a function of space and time, and $`D`$ (the coefficient of self-diffusion) describes the diffusion rate.

We are interested in modeling diffusive phenomena beginning at the microscopic scale and recovering Eq. (1) in the appropriate limit. We first consider a lattice-gas automata (LGA) model of particles hopping at random along the sites of a periodic lattice. The particles execute simultaneous random walks while obeying an exclusion principle. This model was introduced by Toffoli, and it was extended and studied extensively by Chopard and Droz. They show that in the limit of the infinite lattice, the microscopic particle model maps directly onto the macroscopic continuum description, Eq. (1). Additionally, they show that the Green-Kubo relation for expressing transport coefficients (such as the diffusion coefficient $`D`$) in terms of autocorrelation functions is recovered exactly.

Below, we present the model and review the basic analysis of Chopard and Droz. Then we present a new model which introduces a dimension-splitting technique. Then we show, by a remapping, that both models are equivalent.

### II.1 An LGA model of diffusion with exclusion

Consider a two-dimensional square lattice.222A similar argument would apply to models defined on the two-dimensional triangular lattice, with the transport split into a sequence of three operations performed one principle lattice direction at a time. We can think of particles hopping between adjacent sites of the lattice with unit velocity. The particles are initialized on lattice sites so at every integer time, the particles occupy a lattice site. The particles can move in one of four directions: north, west, south, or east—corresponding to the four lattice directions. In the language of lattice-gases, we consider a “transport channel”, $`N_i(\stackrel{}{x},t)`$, along each lattice direction, $`\widehat{c}_i`$, at each site $`\stackrel{}{x}`$ at time $`t`$, where $`i\{0,1,2,3\},`$ and $`\widehat{c}_i`$ is the unit vector along that lattice direction. The transport channel can either be occupied (with a particle to be moved) or empty: $`N_i(\stackrel{}{x},t)=1`$ or $`0`$ respectively. As there are four channels at each site, there can be up to four particles at each site.

The dynamics can be decomposed into two portions: the “interaction” phase and the “streaming” phase. We denote the state after interaction but before streaming as $`N_i^{}(\stackrel{}{x},t)`$. Note, the interval for the completion of both phases is considered the unit interval of time. During the streaming phase the particles are transported along the lattice. If there were no interaction phase, the particles would stream along the lattice ballistically, with their velocities never changing.

To implement simultaneous random walks, during the interaction phase the states of the transport channels at each lattice site are randomly permuted.333We should note that this approach can be considered a field-theoretic molecular dynamics: the fields, not the particles, are permuted. See Ref. cambook  for an early discussion of this distinction. The random permutation is applied to all lattice sites simultaneously, with each permuted independently of all other sites. Since the system is initialized with at most one particle in each channel, this restriction is preserved by the dynamics. Thus all the particles execute simultaneous random walks with exclusion. In the model considered by Chopard and Droz, each channel at site $`\stackrel{}{x}`$ at time $`t`$ was simultaneously permuted by a rotation of $`r(\stackrel{}{x},t)\pi /2`$ radians, where $`r(\stackrel{}{x},t)\{0,1,2,3\}`$ is a random number with values that occur with a frequency corresponding to a specified probability distribution. By choosing this probability distribution appropriately, we can control the mean-free path (i.e., the average number of sites a particle advances before having its velocity permuted). Thus after the permutation,

$$N_i^{}(\stackrel{}{x},t)=N_{i+r(\stackrel{}{x},t)}(\stackrel{}{x},t),$$
(2)
where here and in what follows a sum in a subscript, such as $`(i+r(\stackrel{}{x},t))`$, is taken modulo four. During the subsequent streaming phase the particles are transported along the lattice:

$$N_i(\stackrel{}{x},t+1)=N_i^{}(\stackrel{}{x}\widehat{c}_i,t)=N_{i+r(\stackrel{}{x}\widehat{c}_i,t)}(\stackrel{}{x}\widehat{c}_i,t).$$
(3)
Though Chopard and Droz considered a more general model, we restrict ourselves to the case where $`r(\stackrel{}{x},t)`$ is uniformly and independently distributed amongst the four possible values for each $`(\stackrel{}{x},t)`$. Figure 1 shows the qualitative diffusive behavior of the model. The system is initialized with a uniform block of particles occupying an otherwise vacant two-dimensional lattice, as shown in Fig. 1(a). After $`t=360`$ updates of the full lattice we observe the configuration shown in Fig. 1(b).

#### II.1.1 Establishing the lattice Boltzmann equation

Up to here, we have described the dynamics for a single realization of the lattice, as it is actually implemented and run on a computer in practice. However, for purposes of understanding how to the use this technique as a realistic generator of diffusion, we would like to know “what happens on average in the limit” (where these terms have not yet been defined). Figure 1 shows what happens for a particular realization; in this section, we take the average; and in the next section, we take the limit.

We first consider the mathematical idealization of an ensemble of all such simulations. We consider that each member of the ensemble is initialized in a consistent way (as discussed in detail at the end of this section), and evolves independently. There are many quantities one could average over this ensemble to test the model, but the main quantity of interest is the state of each channel: $`N_i(\stackrel{}{x},t)`$. Taking the expectation of Eq. 3 gives:

$$N_i(\stackrel{}{x},t+1)=N_{i+r(\stackrel{}{x}\widehat{c}_i,t)}(\stackrel{}{x}\widehat{c}_i,t)=\frac{1}{4}\underset{j=0}{\overset{3}{}}N_j(\stackrel{}{x}\widehat{c}_i,t)=\frac{1}{4}\underset{j=0}{\overset{3}{}}N_j(\stackrel{}{x}\widehat{c}_i,t).$$
(4)
Note that we used the fact that the expectation of a sum is the sum of expectations, regardless of the correlations between distinct channels. Finally, defining the occupation number as $`n_i(\stackrel{}{x},t)N_i(\stackrel{}{x},t)`$ gives the lattice-Boltzmann equation444A lattice-Boltzmann equation (LBE) describes a dynamical system spatially arrayed on a lattice but with continuous variables at each site (as opposed to discrete variables, which is the case for CA’s). LBE’s have many properties which make them better suited for modeling hydrodynamic behavior. Yet in situations where microscopic noise, kinetic fluctuations and interparticle correlations play an essential role, one must use LGA models. For a discussion of the striking difference in phenomenology accessible to LGA as opposed to LBE see Ref. bmb-LBE-amphifluids . for our system,

$$n_i(\stackrel{}{x},t+1)=\frac{1}{4}\underset{j=0}{\overset{3}{}}n_j(\stackrel{}{x}\widehat{c}_i,t).$$
(5)
Averaging over the ensemble in this way transforms our discrete Boolean variables into continuous probabilities.

Equation (5) is one of the few instances where an interesting physical dynamics is exactly described by a linear lattice-Boltzmann equation (as discussed above, correlations between channels do not matter). Typically, there are interactions between channels, so in order to perform the ensemble average (as in Eq. (4)), one must make the approximation that the channels are uncorrelated (i.e., the assumption of molecular chaos), and the resulting equations are nonlinear. In the present case there are no interactions between the channels because the streaming and permutations are data blind, meaning they are not conditional on the state of the channels. As there is still confusion about the correspondence between microscopic dynamics and the macroscopic limit, especially in the context of the increase in macroscopic entropy which may accompany a microscopically reversible dynamicsLebowitz ; stauffer-q2r , models described by linear equations may be illuminating. Exact solutions can be obtained and discrepancies with the continuous limit discussed rigorouslyChopDroz-book .

#### II.1.2 Mapping onto the diffusion equation

If we now view each $`n_i(\stackrel{}{x},t)`$ as a smooth function of a continuous spacetime, we can derive a partial differential equation that governs the dynamics of the total density, $`\rho (\stackrel{}{x},t)`$, of particles at a given site:

$$\rho (\stackrel{}{x},t)\underset{i=0}{\overset{3}{}}n_i(\stackrel{}{x},t).$$
(6)
Combining this with Eq. (5) we can write out the evolution of the density,

$`\rho (\stackrel{}{x},t+1)`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{3}{}}}n_i(\stackrel{}{x},t+1)`$
$`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}\left[{\displaystyle \underset{j=0}{\overset{3}{}}}n_j(\stackrel{}{x}\widehat{c}_i,t)\right]={\displaystyle \frac{1}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}\rho (\stackrel{}{x}\widehat{c}_i,t).`$ (7)
The limit of interest is the diffusive regime where the time step, $`\mathrm{\Delta }t`$, and lattice spacing, $`\mathrm{\Delta }x=\left|c_i\right|`$, approach zero: $`\mathrm{\Delta }x0`$, and $`\mathrm{\Delta }t0`$, while $`\left(\mathrm{\Delta }x\right)^2/\mathrm{\Delta }t`$ constant. Note in this regime (of diffusion on an infinite lattice), the propagation velocity, proportional to $`\mathrm{\Delta }x/\mathrm{\Delta }t`$, becomes infinite. This unphysical propagation velocity is actually a well known pathology of the diffusion equation, Eq. (1).

Taylor expanding the spatial terms in Eq. (II.1.2), keeping those terms to order $`\left(\mathrm{\Delta }x\right)^2`$, we obtain

$`\rho (\stackrel{}{x},t+1)`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}{\displaystyle \underset{j=0}{\overset{3}{}}}\left[n_j(\stackrel{}{x},t)+\mathrm{\Delta }x\left(\widehat{c}_i\right)n_j(\stackrel{}{x},t)+{\displaystyle \frac{(\mathrm{\Delta }x)^2}{2}}\left(\widehat{c}_i\right)^2n_j(\stackrel{}{x},t)\right]`$
$`=`$ $`\rho (\stackrel{}{x},t)+{\displaystyle \frac{\left(\mathrm{\Delta }x\right)}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}{\displaystyle \underset{j=0}{\overset{3}{}}}(\widehat{c}_i)n_j(\stackrel{}{x},t)+{\displaystyle \frac{\left(\mathrm{\Delta }x\right)^2}{8}}{\displaystyle \underset{i=0}{\overset{3}{}}}{\displaystyle \underset{j=0}{\overset{3}{}}}(\widehat{c}_i)^2n_j(\stackrel{}{x},t).`$ (8)
After noting that $`\left(\widehat{c}_i=\widehat{c}_{i+2}\right)`$ and $`\left[(\widehat{c}_i)^2+(\widehat{c}_{i+1})^2\right]=^2`$, we can simplify Eq. (II.1.2):

$$\rho (\stackrel{}{x},t+1)=\rho (\stackrel{}{x},t)+\frac{\left(\mathrm{\Delta }x\right)^2}{4}^2\rho (\stackrel{}{x},t).$$
(9)
Taylor expanding to order $`\left(\mathrm{\Delta }t\right)`$ we obtain the diffusion equation:

$$\frac{\rho (\stackrel{}{x},t)}{t}=D^2\rho (\stackrel{}{x},t),$$
(10)
where the diffusion constant, $`D=(\mathrm{\Delta }x)^2/4\mathrm{\Delta }t.`$

We have shown above that the ensemble as a whole reproduces the diffusion equation in the limit of an infinitely fine lattice and continuous time. However, for purposes of an actual LGA simulation, all we have is a single member of the ensemble in a discrete spacetime. In particular, the state of each individual channel is a $`0`$ or $`1`$ and not a real number between $`0`$ and $`1`$ as results from a lattice-Boltzmann treatment. So we want to argue that a typical member of the ensemble—as generated by a single run of a simulation—provides a realistic rendition of the diffusion of a collection of particles.

The main thing we want in this respect is that the expected number of particles in a patch of lattice sites is in some sense the same as the expected number in the corresponding region of continuous space. More specifically, we demand equivalence at the mesoscopic scale (the scale at which patches are large compared to the lattice spacing but small compared to the characteristic scale of the system as a whole). Note that it is also at this mesoscopic scale that we demand equivalence between LGA and lattice-Boltzmann equations. We could sum $`\rho (\stackrel{}{x},t)`$, as obtained from the lattice-Boltzmann formulation, over a patch of lattice sites and obtain a quantity representing the number of particles in that patch. On the other hand, we could generate individual members of the ensemble from simulations according to the LGA formulation, count the number of particles in the patch, and then take the ensemble average of this count. Both approaches give exactly the same quantity. Much as is the case for the separate channels feeding into a single site, as discussed following Eq. (4), the expectation over the ensemble for this count in the LGA formulation is the same as the sum of the expectations giving $`n_i(\stackrel{}{x},t)`$ in the lattice-Boltzmann formulation.

This is true as long as each member of the ensemble is chosen to have the appropriate average mesoscopic initial conditions. Equation (5) provides an exact description of any ensemble that starts with the correct average at each lattice site, no matter how strong the correlations between sites. However, not every member of such an initial ensemble is necessarily “typical” in that it could be consistent with the mesoscopic conditions, yet have a very unusual pattern of correlations between the cells; for example, the particles might align to form an image of a face. To prevent this in practice, one initializes each channel at a relevant site $`(\stackrel{}{x},t)`$ independently at random to be a $`1`$ with a probability given by $`n_i(\stackrel{}{x},t)`$ and to be $`0`$ with a probability given by $`1n_i(\stackrel{}{x},t)`$. This generates a sample of the maximum entropy ensemble that is consistent with the initial specified averages. Moreover, the sample is self-averaging in the sense that the average density of particles in any patch around a point will closely match the initial average at that point. Finally, at any point in the simulation, the configuration is likely to be a typical member of the ensemble because the dynamics only serves to increase coarse-grained entropy, not generate highly correlated patterns. It is these typical ensemble members that provide the most realistic rendition of the diffusion of a collection of particles.

While the above analysis is very promising, we should mention some caveats of the model. In particular, the simulation cannot match more complicated statistics of a real physical system. For example, the fluctuations in the particle number are bounded above by the number of channels in any patch whereas in principle, there is no upper bound in a continuum. Also, the specific microscopic details will never be the same as in a real physical system or from one run to the next. While the exact configuration of all the particles do not matter to the average density of the diffusing species alone, they do matter a great deal when they are coupled to a system that is very sensitive to initial conditions, such as the reversible aggregation model that will be described in Sec. III.1, where microscopic fluctuations dictate the dynamics.

### II.2 The dimension-splitting technique

Our approach is to consider the multi-dimensional diffusion process as a sequence of one-dimensional diffusion events. Each one-dimensional diffusion step consists of an interaction (i.e., mixing) phase and a transport phase. We need consider only a single transport channel and its opposite channel: $`M_0(\stackrel{}{x},t)`$ and $`M_1(\stackrel{}{x},t)`$. First consider the mixing phase. An unbiased, random binary variable, $`\eta (\stackrel{}{x},t)\{0,1\}`$, $`\eta (\stackrel{}{x},t)=1/2`$, is sampled to determine whether the states of the two channels at site $`(\stackrel{}{x},t)`$ are interchanged:

$`M_0^{}(\stackrel{}{x},t)`$ $`=`$ $`[1\eta (\stackrel{}{x},t)]M_0(\stackrel{}{x},t)+\eta (\stackrel{}{x},t)M_1(\stackrel{}{x},t)`$
$`M_1^{}(\stackrel{}{x},t)`$ $`=`$ $`[1\eta (\stackrel{}{x},t)]M_1(\stackrel{}{x},t)+\eta (\stackrel{}{x},t)M_0(\stackrel{}{x},t).`$ (11)
Then the updated channels are transported along the $`\widehat{x}`$ direction, as shown in Fig. 2(a):

$`M_0(\stackrel{}{x},t)`$ $`=`$ $`M_0^{}(\stackrel{}{x}\widehat{x},t1)`$
$`M_1(\stackrel{}{x},t)`$ $`=`$ $`M_1^{}(\stackrel{}{x}+\widehat{x},t1).`$ (12)
Once these two substeps are complete, we consider diffusion along the $`\widehat{y}`$ direction which is also composed of the analogous mixing and transport steps (i.e., substitute $`\widehat{y}`$ for $`\widehat{x}`$ in Eq. (II.2)). As shown in Fig. 2, we first use the two channels to transport data along the horizontal direction, then we use these same two channels to transport data along the vertical direction. We thus reinterpret the function of the same two degrees of freedom.

In this decomposition of a two-dimensional dynamics into two one-dimensional dynamics, instead of four transport channels at each site, we need only two. Thus at each lattice site, three bits of state are required: two bits, $`M_\gamma (\stackrel{}{x},t)`$ where $`\gamma \{0,1\}`$, denote the presence or absence of a particle in each of the two channels respectively; the third bit, $`\eta (\stackrel{}{x},t)`$, is a binary pseudorandom variable. Note, the conventional approach requires six bits of state at each site.

The number of degrees of freedom interacting at any one time is simplified (and remains so regardless of overall dimensionality). Yet we end up executing a sequence of four substeps (mix, then transport along $`\widehat{x}`$, mix then transport along $`\widehat{y}`$). In contrast, the full two-dimensional problem requires only two substeps (one mixing, one transport). Note also that extending the dynamics to one more dimension requires only including an additional sequence of mixing followed by transport (i.e., an additional fractional timestep, with the transport along the new dimension).555Note that when dealing with diffusion on a two-dimensional triangular lattice we would apply a sequence of three substeps, using the same two channels to transport data subsequently along each of the three principle lattice directions. Generalization of transport on a more complex higher-dimensional lattice would require addition of a substep along each new principle lattice direction.

In the two-dimensional diffusion model considered, we find that at the end of one complete update (i.e., one update along $`\widehat{x}`$ and one update along $`\widehat{y}`$) a particle has advanced one unit along each dimension—one step along an oblique lattice with lattice spacing $`\left|\widehat{s_i}\right|=\sqrt{2}\left|\widehat{c_i}\right|`$, where $`\left|\widehat{c_i}\right|`$ is the original lattice spacing. The mapping onto this oblique lattice is shown in Fig. 3.

#### II.2.1 Mapping onto the diffusion equation

We show that this dimension-reduced approach to diffusion maps onto the Chopard-Droz analysis discussed in Sec. II.1, and thus also onto the diffusion equation, Eq. (1), in the appropriate regime. For simplicity we again consider the two-dimensional case, and again deal with the average over an ensemble of similar systems, $`m_i(\stackrel{}{x},t)M_i(\stackrel{}{x},t)`$.

One complete update requires the four substeps discussed above. Thus we consider the unit interval of time to be the time required to complete these four substeps. The result at $`t^{}=t+1/2`$ (of transport and mixing along the horizontal direction) is

$`m_0(\stackrel{}{x},t+{\displaystyle \frac{1}{2}})`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[m_0(\stackrel{}{x}\widehat{x},t)+m_1(\stackrel{}{x}\widehat{x},t)\right]`$
$`m_1(\stackrel{}{x},t+{\displaystyle \frac{1}{2}})`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[m_0(\stackrel{}{x}+\widehat{x},t)+m_1(\stackrel{}{x}+\widehat{x},t)\right].`$ (13)
Now considering transport along the $`\widehat{y}`$ direction we find that

$`m_0(\stackrel{}{x},t+1)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[m_0(\stackrel{}{x}\widehat{y},t+{\displaystyle \frac{1}{2}})+m_1(\stackrel{}{x}\widehat{y},t+{\displaystyle \frac{1}{2}})\right]`$ (14)
$`=`$ $`{\displaystyle \frac{1}{4}}\left[m_0(\stackrel{}{x}\widehat{x}\widehat{y},t)+m_1(\stackrel{}{x}\widehat{x}\widehat{y},t)+m_0(\stackrel{}{x}+\widehat{x}\widehat{y},t)+m_1(\stackrel{}{x}+\widehat{x}\widehat{y},t)\right].`$
Analogously,

$$m_1(\stackrel{}{x},t+1)=\frac{1}{4}\left[m_0(\stackrel{}{x}+\widehat{x}+\widehat{y},t)+m_1(\stackrel{}{x}+\widehat{x}+\widehat{y},t)+m_0(\stackrel{}{x}\widehat{x}+\widehat{y},t)+m_1(\stackrel{}{x}\widehat{x}+\widehat{y},t)\right].$$
(15)
Let us redefine the lattice vectors such that $`\widehat{s}_0=\widehat{x}+\widehat{y}`$, $`\widehat{s}_1=\widehat{x}+\widehat{y}`$, $`\widehat{s}_2=\widehat{s}_0`$, and $`\widehat{s}_3=\widehat{s}_1.`$ Note that the unit vectors of the oblique lattice are longer than those on the original lattice: $`\left|\widehat{s}_i\right|=\sqrt{2}\left|\widehat{c}_i\right|.`$ The evolution of the density

$`\rho (\stackrel{}{x},t+1)`$ $`=`$ $`m_0(\stackrel{}{x},t+1)+m_1(\stackrel{}{x},t+1)`$
$`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}\left[{\displaystyle \underset{j=0}{\overset{1}{}}}m_j(\stackrel{}{x}\widehat{s}_i,t)\right]={\displaystyle \frac{1}{4}}{\displaystyle \underset{i=0}{\overset{3}{}}}\rho (\stackrel{}{x}\widehat{s}_i,t).`$ (16)
Comparing this equation to Eq. (II.1.2), we see that the dimension-splitting model maps onto the conventional model since the oblique lattice vectors, $`\widehat{s}_i`$, map onto the original lattice vectors, $`\widehat{c}_i`$. Reference ChopDroz-diff  and Ref. ChopDroz-book  contain many quantitative plots of the resulting diffusion. Instead of reproducing those plots here, we refer the reader to the original works.

#### II.2.2 Diffusion coefficients

We can easily control the coefficient of self-diffusion ($`D`$ in Eq. (10)), a property which will prove useful in the models discussed in Sec. III. In particular, we can increase the diffusion constant by a factor of $`k^2`$ in one of two ways. First we can transport the particles by $`k`$ lattice unit vectors each streaming step. Secondly, recalling that LGA models have transport steps alternating with interaction steps, we can run $`k^2`$ diffusion steps for each particle interaction step. With more effort, we can use a non-uniform $`r(\stackrel{}{x},t)`$ as shown by Chopard and DrozChopDroz-diff ; ChopDroz-book  to increase or decrease the probability of reversing direction at each interaction step and thereby decrease or increase the mean free path, and hence the diffusion constant.

### II.3 Constraints and conserved quantities

As discussed above, the particles execute simultaneous random walks with exclusion. If we start with at most one particle per channel, that constraint is preserved by the dynamics. Moreover, the total number of particles—alternatively viewed as mass, kinetic energy, or total energy—is conserved. Momentum is not conserved because particles can change their states of motion without a compensating change elsewhere.

Our model also has spurious conserved quantities (quantities that we do not necessarily wish to conserve but are conserved nonetheless). The dynamics partitions the lattice into four interleaved sublattices, and the number of particles on each sublattice is independently conserved. Consider a particle initialized at the center of the lattice (marked by the shaded circle) as shown in Fig. 3(b). At the next time step we know it will occupy one of the sites marked with an open circle. Likewise, only a particle initialized on a site marked by an open circle will be able to occupy the site marked with the shaded circle at the next time step. If we think of superimposing a checkerboard on the lattice, where we color the squares on one checkerboard sublattice “red” and those on the other “black”, it is clear that all particles on the red sites will occupy red sites at the end of any complete update of the system. Yet we can further divide each checkerboard into two. Again consider the situation illustrated in Fig. 3(b). A particle initialized on the shaded site cannot occupy that shaded site at the next time step (even though that shaded site is on the appropriate checkerboard). Thus only half of the sites on that checkerboard are accessible in one time step, and so we are simulating diffusion independently on each of four interleaved sublattices. In contrast, the Chopard-Droz model partitions the lattice into two checkerboard sublattices, simulating diffusion independently on each.666The original lattice-diffusion model introduced by Toffoli uses a partitioning scheme which doesn’t split the space into uncoupled sublattices. If we increase the dimensionality of our model, we double the number of independent checkerboards with each additional dimension.

In Fig. 1, we started all four sublattices in the identical state, so the separate particle conservation on interleaved sublattices is not apparent. If we had started with one of the sublattices occupied, at any point in the dynamics only one of the sublattices would be occupied. When diffusion is used as a component of a larger simulation, these conservation laws are often broken by interactions between subsystems. For example, as discussed below, we consider interactions of diffusing vapor particles and a stationary solid, where the interactions mix the sublattices. A vapor particle will aggregate regardless of the sublattice ordering, and subsequent evaporation from the solid returns that vapor particle placed on an arbitrary sublattice. We can also eliminate the spurious conservation laws in a more intrinsic way: we simply simulate a single sublattice and embed it into the original lattice by shrinking it by a factor of two in each direction. Odd and even timesteps must be embedded with slightly different offsets (amounting to one-half the lattice spacing), but this is just a slight complication in the implementation. For a detailed discussion of the embedding technique see Ref. nhm-spaceram .

Often, when developing microscopic models, we wish to create models which are microscopically invertible (i.e., models which conserve microscopic information). If we use an invertible random number generator for our source of random bits, we can build an invertible model of diffusion, since we can then undo the mixing permutation. (In the first example below, we use an independent invertible LGA as the generator of a dynamical random number field). We can always undo the transport steps (by transporting in the opposite direction). Thus both the Toffoli/Chopard-Droz model and our model can be made invertible.

## III Applications

### III.1 The reversible aggregation model

We developed this technique for simulating diffusion while building a model of cluster growth via the aggregation of particles which were previously diffusing along the lattice. This model is an extension of the athermal diffusion-limited aggregation (DLA) modelWittSand-dla . We place a DLA lattice system in contact with a simulated heat bath and allow only microscopically reversible heat exchanges between the two subsystems.

During each aggregation event, a diffusing gas particle becomes a stationary cluster member while latent heat is released. The heat is explicitly modeled by introducing a diffusing heat particle onto a heat bath lattice. The inverse process is also allowed: a cluster member can absorb a heat particle and evaporate, becoming a gas particle. Initially the heat bath is empty, and we observe rapid, nonequilibrium growth. During this initial phase the cluster structures resemble those typical to DLA systems, such as the bushy clusters formed as frost collects on a window pane. During the subsequent slow approach to thermodynamic equilibrium, the cluster structures anneal until reaching the highest entropy macrostate allowed for a connected cluster in a finite volume: a branched polymer. For more details of the reversible aggregation (RA) model see Ref. dsouza-nhm-RA .

This local, microscopic, reversible approach to modeling allows us to explicitly track how information flow at the microscopic scale gives rise to dissipation and properties such as increasing temperature at the macroscopic scale. We demonstrate that even far from thermodynamic equilibrium, as the simulated clusters are growing and annealing, we have a model with a well defined temperature, and that these models can serve as a numerical laboratory for investigating nonequilibrium thermodynamics. Moreover, discrete reversible systems such as this—which include local heat flow, creation of entropy, and macroscopic dissipation—seem to capture the full essence of the phenomena of pattern formation more completely than do irreversible models. For more details see Refs. dsouza-nhm-RA ; Dsouza-phd .

Simulations of the RA model were implemented and run on the CAM-8 lattice-gas computercam8-waterloo . Figure 4 shows a typical RA cluster, first for a two-dimensional implementation, then for a three-dimensional one. Implementing the three-dimensional model was a simple matter of including one more line of code in our program (for transport along the third dimension). Note that the diffusion technique forms the foundation of our model, which involves two interacting diffusion fields (a field of gas particles and a field of heat particles). We can control the diffusion coefficient independently for each field, as discussed in Sec. II.2.2, and observe the change in the cluster structures and time to reach equilibrium as a function of this parameter.

### III.2 A model for sedimentation

Another application which makes use of the dimension-splitting technique in a reaction-diffusion context is the Biofilm CA model. This is a model developed by Pizarro, Griffeath and NogueraPizarro-sludge ; Pizarro-phd , and programmed by Shalizi and Margolus on the CAM-8cam8-waterloo .

The Biofilm CA simulates the growth of biological films (e.g. sludge, dental plaque) by having food and bacteria react and diffuse in two or more dimensions. The diffusion rate for the food depends on whether there is bacteria present, and the model includes birth, food sources, eating, death by starvation, random death, boundaries, attachment to a substrate, attachment to an aggregate, death of unattached cells, and the periodic removal of entire bacterial clusters which have become disconnected from the substrate (this simulates the action of water currents).

The diffusion in this model is handled one dimension at a time. Each diffusing species can occupy one or both (or neither) of two channels at each lattice site, and the contents of the channels are alternately swapped (or not swapped) at random, and then transported oppositely along one dimension. For the food particles, each channel contains a single bit (present or not). For the bacteria, each channel contains two bits (present or not, and if present how hungry). In regions which should have a longer mean free path, the mixing probability is reduced. Food diffuses in from the top of the space, bacterial growth starts on the substrate at the bottom. To reduce the time taken for food to initially reach the bacteria by diffusion, food starts off uniformly distributed across the space. For more details on the model see Ref. Pizarro-phd .

The model was initially developed as a two-dimensional system, which could be run quickly and easily, and conveniently visualized. The model was tested and debugged in two-dimensions until it behaved as expected. A typical state in the growth of a two-dimensional biofilm is shown in Fig. 5(a). The model was then generalized to three dimensions without adding any bits to each lattice site or defining any new rules to apply to each site. The sequence of operations needed to update a single dimension was simply iterated three times instead of two, with the third iteration involving transport of diffusing particles along the added third dimension. An appropriate three-dimensional initial state was prepared, and generic three-dimensional visualization routines were activated. The model immediately behaved as expected in three dimensions—no further debugging was needed. Only fine tuning of initial densities, diffusion rates, etc., was required. A typical state in the growth of a three-dimensional biofilm is shown in Fig. 5(b).

## IV Discussion

We have introduced a dimension-splitting technique for use in LGA models with diffusive particle transport. For an LGA of arbitrary dimensionality, interactions involving only one dimension at a time are sufficient for iterating the dynamics. Thus, a model can be developed in one dimension, and extended to more dimensions simply by adding additional fractional timesteps. The same dynamics are applied in each fractional timestep, with the same diffusing particles moving along each dimension in turn. By reinterpreting and reusing the same degrees of freedom and the same interactions at each fractional timestep, diffusion is performed with a small amount of state information.

Using this scheme, simple qualitative models involving diffusion can easily be constructed using a kind of stylized molecular dynamics. We can appeal to our intuition in constructing such models, simply including all relevant species and making them all interact in a physically reasonable manner. Constructing more realistic models involves more work, however, since the diffusion behavior of our model has only been demonstrated in the mesoscopic limit. To use the model only in this limit would involve running a very large number of diffusion steps for each particle interaction step. This might be resolved by requiring a finer lattice for the diffusers than for other entities in the simulation. For many models, however, we don’t expect that a quantitatively realistic interaction requires explicit mesoscopic averaging of the microscopic constituent simulations. Demonstrating realism for a given kind of simulation requires both theoretical work and experimental verification. Even models which are not completely realistic, however, may be view as “toy” physical systems which can be simulated exactly and analyzed theoreticallydsouza-nhm-RA .

It seems that there is an important cost for our model’s simplicity: the dynamic range in particle density at each lattice site is small. $`N`$ bits of state at each site yield a maximum density of only $`N`$ diffusers per site. Numerical simulations of diffusion, in contrast, allow up to $`2^N`$ diffusers to be represented at each site using the same amount of state, with only a linear increase in the computational complexity of the dynamics. However, this contrast is only apparent. The different species present at a lattice site can of course be given different power-of-two weights: each “1” of species $`A`$ might represent a single particle, each “1” of species $`B`$, two particles, etc. Independent dimension-split diffusion of each such species provides a large dynamic range while retaining exact particle conservation.

LGA models with as few as two bits of diffusing state at each lattice site can be constructed using the dimension-splitting technique. The benefits of similar economies in simulation state have recently been discussed in other contexts. For example, in the context of reversible compilers see Ref. Perumala-rev-comp ; in the context of systems relying on pseudorandom numbers see Ref. Nisan92 .

It might seem that the operation on one dimension at a time would result in slower simulations than a more complex LGA model, which mixes data in several dimensions at once. However, from a computational complexity standpoint, this is not at all obvious, since the more complex LGA model would have more data to interact. If, for example, diffusion is being handled by circuitry capable of complex operations, then using only simple operations may waste resources. If, on the other hand, only very simple operations are available in the circuitry, and only one operation can be performed at a time, then more data would require more operations and mean a slower simulation. If it is ever possible to perform simulations with molecular computers, the latter of these two situations seems much more likely than the former.

Our dimension-splitting diffusion technique is particularly well suited for use in concert with other dimension-splittable primitives. Composite models based only on such primitives share the properties of reduced state and the possibility of developing models in one or two dimensions, and then running them without change in higher dimensions.

A useful primitive of this type is the LGA-fluid attraction and repulsion technique used by Yepez and Appert in their momentum-conserving lattice-gas aggregation modelsyepez-blobs ; apert-blobs . In these models, pairs of particles lying along principle lattice directions and separated by a given distance are brought together to interact. By using a set of distances and specified types of interactions (repulsive interactions turn both particles away from their midpoint, attractive interactions turn them towards it), prescribed potentials can be synthesized operating between particle aggregates.777Their original models didn’t use dimension-splitting but instead used “virtual particles” moving in all directions at once to carry the forces. In translating their models to the CAM-8 machine, we employed dimension-splitting.

Since these models aren’t hydrodynamic, momentum conserving gas interactions used to mix particles going in different directions can also be made to interact one dimension at a time, forming a rest particle from each head-on collision which then decays into two new particles moving out in opposite directions—when the particle decays can be determined stochastically. Using the technique of the present paper, diffusing components can also be included in such models. Note that in all of these models, the fractional timesteps are applied along each of the principle lattice directions, so that in a two-dimensional triangular lattice model, for example, there are three transport steps rather than just two. Thus a better name for the technique might be “lattice-splitting” rather than “dimension-splitting”. It would be interesting to find additional examples of such splitting techniques which can be used as components of LGA models.


text:
# Age Dating of Globular Clusters Using UBV(RI) Main-sequence Photometry in the Two-color Diagram: Age of NGC 6397

## 1. Age Indicator, Isochrones and Observational Data

The main sequence of metal-poor GCs exhibits a bend in the two–color diagram(s). Analysis of the multicolor Yale isochrones of Demarque et al. (1996) shows that color width of such a bend, i.e., a color difference between the turnoff point and main sequence at the same level of the $`UB`$, is a function of age at a given metallicity. Hence age of a GC can be derived by fitting isochrones of appropriate metallicity to the cluster two–color diagram(s). The best fit isochrone must show the most adequate reproduction of fiducial line of the sequence (mainly in its upper part) and location of its three basic points: points of the turnoff and main sequence at the same level of the $`UB`$, and point of the extremum.

Multicolor photometry by Alcaino et al. (1997) of NGC 6397 has been used to apply this approach. The densest part of the cluster field and stars above the turnoff level were excluded from consideration. Since the most probable cluster metallicity (see Gratton et al. 2000; Castilho et al. 2000) is close to $`[Fe/H]=2.0`$, the set of isochrones of appropriate metallicity was used.

## 2. Result

Panels of Fig. 1 show that the fitting criteria described above are best met by isochrone for 16 Gyr age in each diagram presented. A fitting accuracy is of order $`\pm 1`$ Gyr, taking into account that result of the best fit is approximately the same in all the diagrams. However, real error in age may be as large as a few Gyr (VandenBerg 1999, private communication) even if a photometry is precise, with no systematic effects. That may be due to incertaincy in the parameters of theoretical calculations such as $`T_{eff}`$, or to the problems of transformation of isochrones to the colors used, etc. (e.g., see, for details, Grundahl et al. 2000). Note that change of 0.1 dex in $`[Fe/H]`$ causes change in age estimate $`1`$ Gyr so that say for $`[Fe/H]=1.9`$ the age derived for NGC 6397 is $`15`$ Gyr.

### Acknowledgments.

I am deeply thankful for the IAU travel grant. Also, I am very grateful to Dr. D. VandenBerg for valuable discussions and comments.

## References

Alcaino, G., Liller, W., Alvarado, F., et al. 1997, AJ, 114, 1067

Castilho, B. V., Pasquini, L., Allen, D. M., et al. 2000, A&A, 361, 92

Demarque, P., Chaboyer, B., Guenther, D., et al. 1996, Yale Isochrones

Gratton, R. G., Bonifacio, P., Bragaglia, A., et al. 2000, astro-ph/0012457

Grundahl, F., VandenBerg, D. A., Bell, R. A, et al. 2000, AJ, 120, 1884