text:
# UUITP-19/02 ITEP-TH-60/02 Comparing strings in 𝐴𝑑𝑆₅×𝑆⁵ to planar diagrams: an example

## 1 Introduction

The AdS/CFT correspondence establishes an equivalence of the superstring theory in the $`AdS_5\times S^5`$ background and $`𝒩=4`$ supersymmetric Yang-Mills theory (SYM) . This equivalence is a true duality: the interaction strength in the string sigma-model is the inverse of the ’t Hooft coupling in the field theory, which makes any direct comparison of strings with the field theory extremely difficult. On the other hand, planar diagrams of the large-$`N`$ perturbation theory resemble discretized string world sheets , and one may think that a more direct relationship between them exists beyond the fact that they describe one and the same theory in different regions of parameter space. In principle, one has to sum all planar diagrams (solve large-$`N`$ SYM) in order to reach weakly-coupled stringy regime. As became clear recently, it is still possible to make a comparison to string theory without summing all diagrams. This can be done by considering special states with large quantum numbers. In string theory, such states are semiclassical without any reference to the weakness of interactions in the sigma-model. These states, therefore, have a simple string description all the way down to weak coupling where perturbative SYM theory becomes accurate. Though this logic involves a certain stretch in the arguments, it has worked in a number of examples . With the help of the semiclassical string picture, a remarkable progress has been made in understanding operators with large R charge (and also with large spin) in the SYM theory . Correlation functions of these operators have a beautiful world-sheet description, and vice versa, the world-sheet dynamics of the string states with large angular momentum on $`S^5`$ is encoded in certain set of Feynman diagrams which have relatively simple structure .

A local operator in the SYM theory is dual to a closed string state in $`AdS_5\times S^5`$ . The string theory in AdS also has an open string sector which is associated with Wilson loops . By probing Wilson loops with operators that have parametrically large R charge it is possible to reach the semiclassical regime in the open-string sector too . Whether the semiclassical description in the open string sector captures perturbative regime in the field theory or still requires the ’t Hooft coupling to be large is not entirely clear. We will address this question by comparing string-theory calculations with lowest-order Feynman diagrams.

We consider a two-point correlator $`W(C)𝒪_J`$ of the Wilson loop with the operator that carries charge $`J`$ under a $`U(1)`$ subgroup of the $`SO(6)`$ R-symmetry group. The operator of interest is chiral primary

$$𝒪_J=\frac{(2\pi )^J}{\sqrt{J}\lambda ^{J/2}}trZ^J,$$
(1.1)
where $`Z=\mathrm{\Phi }_1+i\mathrm{\Phi }_2`$ is a combination of two of the six adjoint scalars in the $`SU(N)`$ SYM theory. The correlator $`W(C)𝒪_J`$ measures the weight with which operator $`𝒪_J`$ appears in the local operator expansion of the Wilson loop. We should compare this to an overlap of the closed string boundary state created by the Wilson loop with the supergravity state dual to the operator $`𝒪_J`$ . There is one special case in which the exact answer is know: It is possible to compute the correlator exactly for the circular Wilson loop , which is a chiral operator in a certain sense . Because of the supersymmetry cancellations, only diagrams without internal vertices contribute to its expectation value , as well as to its expansion coefficients in chiral primaries . These diagrams can be explicitly resummed which yields the exact results valid at any ’t Hooft coupling. The semiclassical string calculations can also be done to all orders in the sigma-model perturbation theory because of geometric symmetries of the circle . We consider arbitrary contours in this paper, for which the equations of motion in the sigma-model can be solved order by order in $`\lambda /J^2`$, where $`\lambda `$ is the string tension squared***We express the string tension and all other dimensionful quantities in the units of the AdS radius., which according to the AdS/CFT dictionary coincides with the ’t Hooft coupling in the SYM theory: $`\lambda =g_{SYM}^2N`$. We will compare the expansion of the string amplitude in $`\lambda /J^2`$ to the ordinary planar perturbation theory.

The paper is organized as follows. In sec. 2, we set up the notations and review how the correlator of a Wilson loop with the local operator is computed in string theory. In sec. 3, we do a one-loop calculation on the SYM side. The classical solution of the sigma-model, which describes the correlator in string theory, is constructed in sec. 4. We then compare the string amplitude determined by this solution with one-loop SYM perturbation theory. We also show in sec. 5 that diagrams of the SYM perturbation theory exponentiate in accord with predictions of the string theory. We draw the conclusions and discuss the results in sec. 6.

## 2 Wilson loop correlator in string theory

The supersymmetric Wilson loop operator is defined as

$$W(C)=tr\mathrm{P}\mathrm{exp}[ds(iA_\mu x^{}{}_{}{}^{\mu }+\mathrm{\Phi }_1|x^{}|)].$$
(2.1)
This operator is a hybrid of an ordinary non-Abelian phase factor and the unique scalar loop operator which is conformally and Lorentz covariant The most general Wilson loop depends on scalars through an arbitrary linear combination $`\mathrm{\Phi }_i\theta ^i`$, where $`\theta ^i`$ is a unit six-vector, which may depend on $`s`$ and thus parameterizes a closed contour in $`S^5`$. We choose a particular $`\theta ^i`$ for definiteness; results will not be much different for any other constant $`\theta ^i`$. It would be interesting to consider varying $`\theta ^i`$ , but this goes beyond the scope of the present paper..

The Wilson loop operator is dual to a macroscopic string in $`AdS_5\times S^5`$ . The local operator is dual to a supergravity mode. If we want to compute their correlator, we should find vertex operator associated with the supergravity mode, insert it into the world sheet of the string, combine it with the propagator of the supergravity mode, and then integrate over all string world-sheets and all positions of the vertex operator . In the limit when the distance between the loop and the insertion of the local operator goes to infinity the propagator factorizes, and we are left with the one-point correlation function in the sigma-model:

$$W(C)O_J(x)=\frac{1}{|x|^{2J}}d^2\sigma _0D𝒳\mathrm{e}^{S_{\sigma m}[𝒳]}V_J(𝒳(\sigma _0))+o(1/|x|^{2J}),$$
(2.2)
where $`𝒳`$ denotes coordinates, bosonic and fermionic, of the string in the AdS superspace; $`S_{\sigma m}[𝒳]`$ is the Green-Schwarz action of the AdS sigma-model; and $`V_J(𝒳(\sigma ))`$ is the vertex operator which creates string mode dual to the operator $`𝒪_J`$. We assume that the conformal gauge is already fixed and regard corresponding Faddeev-Popov factor as a part of the measure in the path integral.

The metric of $`AdS_5`$ in the Póincare coordinates is

$$ds^2=dp^2+\mathrm{e}^{2p}dx^2.$$
(2.3)
The boundary is at $`p=\mathrm{}`$ and the horizon of $`AdS_5`$ is at $`p=\mathrm{}`$. We adopt the following parameterization of $`S^5`$: $`\theta ^i=(\mathrm{cos}\psi \mathrm{cos}\phi ,\mathrm{cos}\psi \mathrm{sin}\phi ,\mathrm{sin}\psi 𝐧)`$, where $`𝐧`$ is a unit four-vector. The sigma-model action in this parameterization is

$$S_{\sigma m}=\frac{\sqrt{\lambda }}{4\pi }d^2\sigma \left[\left(P\right)^2+\mathrm{e}^{2P}\left(X^\mu \right)^2+\left(\mathrm{\Psi }\right)^2+\mathrm{cos}^2\mathrm{\Psi }\left(\mathrm{\Phi }\right)^2+\mathrm{}\right].$$
(2.4)
The action for the rest of the angles in $`S^5`$ and for the world-sheet fermions will not be important for us. The operator $`𝒪_J`$ is chiral and hence corresponds to the supergravity (ten-dimensional massless) mode. The vertex operator of such string state is a solution of the massless wave equation in $`AdS_5\times S^5`$ . Symmetries and the form of the dual SYM operator $`𝒪_J`$ essentially determine the necessary solution: it should have zero momentum along the boundary, since we consider the large-distance asymptotic of the correlator; it should scale as $`V\mathrm{e}^{J\omega }V`$ under $`PP+\omega `$, since the operator $`𝒪_J`$ has dimension $`J`$; and it should transform as $`V\mathrm{e}^{iJ\phi }V`$ under an $`SO(6)`$ rotation $`\mathrm{\Phi }\mathrm{\Phi }+\phi `$. These simple arguments determine the vertex operator with an exponential accuracy, which will be sufficient for our purposes:

$$V\mathrm{e}^{iJ\mathrm{\Phi }JP}.$$
(2.5)
Fixing normalization is possible , but requires much more work.

We will calculate the correlator (2.2) in the double-scaling limit of large string tension and large R charge with their ratio $`j=J/\sqrt{\lambda }`$ fixed. The string path integral is semiclassical in this limit, so we need just to solve classical equations of motion in the sigma-model, which fortunately leaves behind many hard questions like correct normalization of the measure in the path integral or exact form of the vertex operator. Since the exponent in (2.5) is of the same order as the action, we should treat the action and the vertex operator on the same footing . This effectively adds a source to the action:

$$S_{\mathrm{eff}}=S_{\sigma m}+J\left(P(\sigma _0)i\mathrm{\Phi }(\sigma _0)\right)JS.$$

The action with the source term added is the functional that we should minimize in order to compute the string amplitude in the semiclassical approximation. The action evaluated on the classical solution determines the correlator $`W(C)O_J(x)`$ with an exponential accuracy at large $`\lambda `$ and large $`J`$. Since the rescaled action $`S`$ and, consequently, the equations of motion depend on $`\lambda `$ and $`J`$ only through the ratio $`j=J/\sqrt{\lambda }`$, the correlator will have the form:

$$W(C)O_J(x)\frac{1}{|x|^{2J}}\mathrm{e}^{JS(j;C)}.$$
(2.6)
Solving classical equations of motion of the AdS sigma-model for generic boundary conditions is a hard problem, but it is possible to find an approximate solution when the parameter $`j`$ is large. Then the minimal surface is a cylinder to a first approximation : the shape of the contour is the same in any slice of AdS parallel to the boundary. It was shown in  that this solution predicts the following dependence of the correlator on $`J`$ and $`\lambda `$:

$$W(C)O_J(x)\frac{\lambda ^{J/2}}{J!}.$$
(2.7)
Surprisingly, this is the same as one would get by counting combinatorial factors in the lowest-order diagram (fig. 1) of the SYM perturbation theory. We should stress that the perturbative and the string calculations are valid in different regions of parameters: the semiclassical approximation in string theory requires both $`J`$ and $`\lambda `$ to be large, but the ratio $`\lambda /J^2`$ can be small. The perturbation theory, strictly speaking, is valid only when the coupling is small. Of course, it may turn out that $`\lambda /J^2`$ is the true parameter of perturbative expansion, then $`\lambda `$ should not necessarily be small for the perturbation theory to work.

We will develop a systematic method to solve the sigma-model equations of motion order by order in $`1/j^2`$. This will yield an expansion of the Wilson loop correlator in $`\lambda /J^2`$, which resembles an ordinary perturbative series. To compare the two, we will first compute the one-loop correction on the SYM side.

## 3 One-loop calculation in the SYM theory

The one-loop planar correction to the correlator of a Wilson loop with the chiral primary operator is described by the single diagram in fig. 2. The diagrams with corrections to the external lines exactly cancel, as shown in appendix B. Combining the tree-level diagram 1 with the one-loop correction givesOur notations and conventions for the SYM perturbation theory are collected in appendix A.:

$`W(C)𝒪_J(x)`$ $`=`$ $`{\displaystyle \frac{\lambda ^{J/2}}{\sqrt{J}(4\pi |x|^2)^J}}{\displaystyle 𝑑s_1\mathrm{}𝑑s_J\theta _c(s_1,\mathrm{},s_J)}`$ (3.8)
$`\times \left(1+\lambda {\displaystyle \underset{n=1}{\overset{J}{}}}{\displaystyle _{s_n}^{s_{n+1}}}𝑑s{\displaystyle _s^{s_{n+1}}}𝑑rG(s,r)+O(\lambda ^2)\right)`$
$`=`$ $`{\displaystyle \frac{\sqrt{J}}{J!}}\left({\displaystyle \frac{\sqrt{\lambda }}{4\pi |x|^2}}\right)^J`$
$`\times \left[(2\pi l)^J+\lambda {\displaystyle _0^{2\pi l}}𝑑s{\displaystyle _0^{2\pi l}}𝑑q(2\pi lq)^JG(s,s+q)+O(\lambda ^2)\right],`$
where $`\theta _c`$ equals one, if its arguments are cyclically ordered along the contour, and equals zero otherwise, $`2\pi l=𝑑s|x^{}|`$ is the length of the contour $`C`$, and $`G(s,r)`$ is the sum of the gauge-boson and the scalar propagators inserted between points $`x(s)`$ and $`x(r)`$ on the contour:

$$G(s,r)=\frac{1x^{}(s)x^{}(r)}{8\pi ^2|x(s)x(r)|^2}.$$
(3.9)
From now on, we use the natural parameterization of the contour $`C`$: $`x^{}{}_{}{}^{2}=1`$. In going from the second to the third line in (3.8), we used the formula

$$_0^s𝑑s_1\mathrm{}_0^{s_{J1}}𝑑s_J=\frac{s^J}{J!}$$
(3.10)
to integrate over the end-points of external legs in the diagrams 1, 2.

The answer considerably simplifies in the limit of large $`J`$, because then

$$_0^{2\pi l}𝑑q(2\pi lq)^JF(q)=\frac{(2\pi l)^J}{J}F(0)+O(1/J^2)$$
(3.11)
for any smooth function $`F(q)`$, and the propagator (3.9) can be expanded as

$$G(s,s+q)=\frac{x^{\prime \prime \prime }x^{}}{16\pi ^2}+O(q)=\frac{x^{\prime \prime }^2}{16\pi ^2}+O(q).$$
(3.12)
We get:

$$W(C)𝒪_J(x)=\frac{\sqrt{J}}{J!}\left(\frac{\sqrt{\lambda }l}{2|x|^2}\right)^J[1+\frac{\lambda l}{8\pi J}_0^{2\pi l}dsx^{\prime \prime }{}_{}{}^{2}+O\left(\frac{\lambda }{J^3}\right)+O(\lambda ^2)].$$
(3.13)
The fact that the large-$`J`$ limit is only sensitive to the local limit of the propagator has a simple explanation: since the operator $`𝒪_J`$ contains $`J`$ scalar fields, there are $`J`$ external vertices on the Wilson loop, and the average distance between them is of order $`l/J`$. A distance between the end-points of the propagator $`G(s,r)`$ is, on average, of the same order: $`|rs|l/J`$. The expansion of the propagator $`G(s,r)`$ in $`(rs)`$ therefore generates $`1/J`$ expansion of the Feynman integral.

## 4 String side

In this section, we will develop a systematic method to solve equations of motion in the sigma-model order by order in $`1/j^2`$ and to compute the classical string action as a series in this parameter.

We can choose the coordinates $`\tau `$ and $`s`$ on the world sheet, such that $`s`$ coincides with the natural parameter on $`C`$ when restricted to the boundary at $`\tau =0`$. We can also assume that the world sheet has a topology of the cylinder, and that the vertex operator is placed at $`\tau =\mathrm{}`$. The action then takes the following form§§§The other $`S^5`$ angle $`\mathrm{\Psi }`$ in (2.4) is zero on the classical solution.:

$`S`$ $`=`$ $`{\displaystyle \frac{1}{4\pi j}}{\displaystyle _\epsilon ^{\mathrm{}}}d\tau {\displaystyle _0^{2\pi l}}ds[\dot{P}{}_{}{}^{2}+P^{}{}_{}{}^{2}+\mathrm{e}^{2P}(\dot{X}{}_{}{}^{2}+X^{}{}_{}{}^{2})+\dot{\mathrm{\Phi }}{}_{}{}^{2}+\mathrm{\Phi }^{}{}_{}{}^{2}]{\displaystyle \frac{l}{j\epsilon }}`$ (4.14)
$`+\left(Pi\mathrm{\Phi }\right)|_{\tau =\mathrm{}}`$
The cutoff on the integral over $`\tau `$ is necessary to regularize the divergence of the area at the boundary of AdS space. The area diverges as $`2\pi l/\epsilon `$. With our normalization of the action, we should subtract $`l/j\epsilon `$ to make the action finite. As was shown in , the simple subtraction is the correct way to deal with the divergences in the AdS sigma-model, at least within the semiclassical approximation.

The equations for the minimal surface should be supplemented by boundary conditions. The boundary conditions at $`\tau =0`$ are set by the Wilson loop:

$$X^\mu (s,0)=x^\mu (s),P(s,0)=+\mathrm{},\mathrm{\Phi }(s,0)=0.$$
(4.15)
The source term in the action determines the behavior of the string coordinates at infinity. Let $`z`$ be a local coordinate near $`\sigma _0`$, so that the operator is inserted at $`z=0`$. Then $`X^\mu `$ are regular at $`z=0`$, but $`P`$ and $`\mathrm{\Phi }`$ have logarithmic singularities determined by the source. With the action normalized as in (2.4), we have:

$$P(z)=j\mathrm{ln}|z|+\mathrm{regular},\mathrm{\Phi }(z)=ij\mathrm{ln}|z|+\mathrm{regular}.$$
(4.16)
The exponential parameterization $`z=\mathrm{e}^{(\tau +is)/l}`$ maps the world sheet to a cylinder. The point of operator insertion $`z=0`$ is mapped to infinity with the following boundary conditions on the string coordinates:

$$X^\mu (s,\tau )x_{\mathrm{}}^\mu ,P(s,\tau )\frac{j\tau }{l},\mathrm{\Phi }(s,\tau )\frac{ij\tau }{l}(\tau \mathrm{}).$$
(4.17)
We will solve the equations of motion perturbatively in $`1/j^2`$. In order to do that, it is convenient to introduce new variables in which the solution is non-singular at zeroth order (at infinite $`j`$). We rescale:

$`t=j\tau ,`$ (4.18)
$`Q=P\mathrm{ln}j.`$ (4.19)
In this way, we get:

$`S`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _{j\epsilon }^{\mathrm{}}}dt{\displaystyle _0^{2\pi l}}ds[\dot{Q}{}_{}{}^{2}+{\displaystyle \frac{1}{j^2}}Q^{}{}_{}{}^{2}+e^{2Q}(j^2\dot{X}{}_{}{}^{2}+X^{}{}_{}{}^{2})+\dot{\mathrm{\Phi }}{}_{}{}^{2}+{\displaystyle \frac{1}{j^2}}\mathrm{\Phi }^{}{}_{}{}^{2}]`$ (4.20)
$`{\displaystyle \frac{l}{j\epsilon }}+\left(Qi\mathrm{\Phi }\right)|_{t=\mathrm{}}+\mathrm{ln}j.`$
The dynamics of the angular coordinate $`\mathrm{\Phi }`$ is trivial:

$$\mathrm{\Phi }=\frac{it}{l}.$$
(4.21)
Upon substitution of the solution for $`\mathrm{\Phi }`$, the action combined with the vertex operator (which just imposes the boundary condition for $`Q`$ at $`t\mathrm{}`$) can be written in the following form:

$$S=\frac{1}{4\pi }_{j\epsilon }^{\mathrm{}}dt_0^{2\pi l}ds[(\dot{Q}+\frac{1}{l})^2+\frac{1}{j^2}Q^{}{}_{}{}^{2}+\mathrm{e}^{2Q}(j^2\dot{X}^2+X^{}{}_{}{}^{2})]+Q(j\epsilon )\frac{l}{j\epsilon }+\mathrm{ln}j.$$
(4.22)
Here, we used the equality

$$\left(Q(t)i\mathrm{\Phi }(t)\right)|_{t=\mathrm{}}=\frac{1}{4\pi }_{j\epsilon }^{\mathrm{}}𝑑t_0^{2\pi l}𝑑s\frac{2}{l}\left(\dot{Q}+\frac{1}{l}\right)+Q(j\epsilon ).$$

Note that potentially dangerous boundary terms in (4.20) have been completely absorbed into the bulk part of the action. This would have been impossible without cancellations between $`AdS_5`$ and $`S^5`$ contributions which occur because the vertex operator is marginal. Potential divergences at $`t\mathrm{}`$ are of UV nature, since they come from the vicinity of the operator insertion, though in the coordinates we use they might look as an IR effect. Anyway, they cancel and the action is saturated by $`t1`$. The typical AdS scale $`\mathrm{e}^P`$ is also small: $`Q`$ is finite at large $`j`$ and $`P\mathrm{ln}j`$ according to (4.19). Hence, the largest contribution to the action comes from the region of AdS space close to the boundary.

The equations of motions are

$`\ddot{Q}{\displaystyle \frac{1}{j^2}}Q^{\prime \prime }+\mathrm{e}^{2Q}(j^2\ddot{X}{}_{}{}^{2}+X^{}{}_{}{}^{2})=0,`$ (4.23)
$`j^2\left(\mathrm{e}^{2Q}\dot{X}\right)\dot{}+\left(\mathrm{e}^{2Q}X^{}\right)^{}=0.`$ (4.24)
They are readily solved at large $`j`$. First, we notice that $`X^\mu `$ must be $`t`$-independent to the leading order. Hence, the minimal surface is a cylinder with the contour $`C`$ as the base:

$$X_0^\mu (s,t)=x^\mu (s).$$
(4.25)
The equation (4.23) reduces to

$$\ddot{Q}_0+\mathrm{e}^{2Q_0}=0,$$
(4.26)
which is solved by

$$Q_0=\mathrm{ln}\left(l\mathrm{sinh}\frac{t}{l}\right).$$
(4.27)
To develop a systematic procedure to solve the equations of motion order by order in $`1/j^2`$, we write:

$$Q=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{j^{2n}}Q_n,X^\mu =\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{j^{2n}}X_n^\mu .$$
(4.28)
The structure of equations (4.23), (4.24) is such that we can easily solve them recursively. On each step we will need to solve an ordinary linear differential equation for $`Q_n`$, $`X_n^\mu `$ with an unknown dependence on time only. This is somewhat surprising, since we are dealing with partial differential equations. The reason for the simplification is an enhancement of time derivatives by a factor of $`j^2`$ compared to the derivatives in $`s`$. The iterative solution can be constructed as follows: suppose that we know $`Q_m`$, $`X_m^\mu `$ for $`m<n`$. They must solve the first equation (4.23) with an accuracy $`O(1/j^{2n})`$, and the second equation (4.24) with an accuracy $`O(1/j^{2n2})`$. The next order in the second equation (4.24) allows us to express a linear combination of time derivatives of $`X_n^\mu `$ through the known functions. $`Q_n`$ drops out at this order because $`\dot{X}{}_{0}{}^{\mu }=0`$. The first equation $`(\text{4.23})`$ at order $`1/j^{2n}`$ then reduces to a linear equation for $`Q_n`$ and its second time derivative whose coefficients depend on $`X_n^\mu `$, which we already know, and on $`Q_m`$, $`X_m^\mu `$ with $`m<n`$.

The iterative procedure is best exemplified by the first three steps:

$`\left(\mathrm{e}^{2Q_0}\dot{X}_1^\mu \right)\dot{}`$ $`=`$ $`(\mathrm{e}^{2Q_0}X_0^\mu {}_{}{}^{}){}_{}{}^{}X_1^\mu ,`$
$`\ddot{Q}_12\mathrm{e}^{2Q_0}Q_1`$ $`=`$ $`\mathrm{e}^{2Q_0}\left(\dot{X}_1^2+2X_0^{}X_1^{}\right)Q_1,`$
$`\left[\mathrm{e}^{2Q_0}\left(2Q_1\dot{X}_1^\mu +\dot{X}_2^\mu \right)\right]\dot{}`$ $`=`$ $`\mathrm{e}^{2Q_0}(2Q_1X_0^\mu {}_{}{}^{}+X_1^\mu {}_{}{}^{}){}_{}{}^{}X_2^\mu ,`$
$`\mathrm{}`$
It is straightforward to integrate these equations, though calculations become increasingly difficult with the order of iteration. After rather lengthy algebra we obtain:

$`X_1^\mu `$ $`=`$ $`{\displaystyle \frac{l^2}{4}}(\mathrm{e}^{2\stackrel{~}{t}}+2\stackrel{~}{t}1)x^\mu {}_{}{}^{\prime \prime },`$
$`Q_1`$ $`=`$ $`{\displaystyle \frac{l^2}{4}}{\displaystyle \frac{4\stackrel{~}{t}\mathrm{e}^{2\stackrel{~}{t}}+\mathrm{e}^{4\stackrel{~}{t}}1}{1\mathrm{e}^{2\stackrel{~}{t}}}}x^{\prime \prime }{}_{}{}^{2},`$
$`X_2^\mu `$ $`=`$ $`{\displaystyle \frac{l^4}{16}}[(2\stackrel{~}{t}^2+2\stackrel{~}{t}\mathrm{e}^{2\stackrel{~}{t}}+\mathrm{e}^{2\stackrel{~}{t}}1)x^\mu ^{\prime \prime \prime \prime }`$ (4.30)
$`+4(\stackrel{~}{t}\mathrm{e}^{2\stackrel{~}{t}}+\stackrel{~}{t}+\mathrm{e}^{2\stackrel{~}{t}}1)(x^{\prime \prime }{}_{}{}^{2}x_{}^{\mu }{}_{}{}^{})^{}`$
$`(\mathrm{e}^{4\stackrel{~}{t}}+8\stackrel{~}{t}\mathrm{e}^{2\stackrel{~}{t}}+4\mathrm{e}^{2\stackrel{~}{t}}+4\stackrel{~}{t}5)x^{\prime \prime }{}_{}{}^{2}x_{}^{\mu }{}_{}{}^{\prime \prime }],`$
where $`\stackrel{~}{t}=t/l`$.

The classical action can be written as a power series in $`1/j^2`$:

$$S=\mathrm{ln}j+\underset{n=0}{\overset{\mathrm{}}{}}\frac{S_n}{j^{2n}}$$
(4.31)
withWe can forget about regularization everywhere except in $`S_0`$.

$`S_0`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _0^{2\pi l}}ds{\displaystyle _{j\epsilon }^{\mathrm{}}}dt[(\dot{Q}_0+{\displaystyle \frac{1}{l}})^2+\mathrm{e}^{2Q_0}X_0^{}{}_{}{}^{2}]+Q_0(j\epsilon ){\displaystyle \frac{l}{j\epsilon }},`$
$`S_1`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _0^{2\pi l}}𝑑s{\displaystyle _0^{\mathrm{}}}𝑑t\mathrm{e}^{2Q_0}\left(\dot{X}_1^2+2X_0^{}X_1^{}\right),`$
$`S_2`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _0^{2\pi l}}𝑑s{\displaystyle _0^{\mathrm{}}}𝑑t\mathrm{e}^{2Q_0}\left(Q_1\dot{X}_1^2+\dot{X}_1\dot{X}_2+2X_0^{}X_2^{}\right).`$ (4.32)
Calculating the integrals over $`t`$, we get:

$`S_0`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{2}{l}}1,`$
$`S_1`$ $`=`$ $`{\displaystyle \frac{l}{8\pi }}{\displaystyle _0^{2\pi l}}dsx^{\prime \prime }{}_{}{}^{2},`$
$`S_2`$ $`=`$ $`{\displaystyle \frac{l^3}{64\pi }}{\displaystyle _0^{2\pi l}}ds\left[2x^{\prime \prime \prime }{}_{}{}^{2}(x^{\prime \prime }{}_{}{}^{2})^2\right].`$ (4.33)
The appearance of the dimensionful quantity under the logarithm in $`S_0`$ may seem strange, but in fact is required to reproduce the correct scaling dimension of the correlator, for which we get

$`W(C)𝒪_J(x)`$ $``$ $`{\displaystyle \frac{1}{J!}}\left({\displaystyle \frac{\sqrt{\lambda }l}{2|x|^2}}\right)^J\mathrm{exp}{\displaystyle _0^{2\pi l}}ds\{{\displaystyle \frac{\lambda l}{8\pi J}}x^{\prime \prime }{}_{}{}^{2}{\displaystyle \frac{\lambda ^2l^3}{64\pi J^3}}\left[2x^{\prime \prime \prime }{}_{}{}^{2}(x^{\prime \prime }{}_{}{}^{2})^2\right]`$ (4.34)
$`+O\left({\displaystyle \frac{\lambda ^3}{J^5}}\right)\}.`$
Here we used the Stierling formula to approximate $`\mathrm{e}^J/J^J1/J!`$. Expanding the correlator in $`\lambda `$, we find the complete agreement with the one-loop perturbation theory (3.13)! We also get a prediction for the two-loop contribution. The two-loop calculation on the SYM side is beyond the scope of the present paper. Instead, we will check that a certain set of diagrams of all orders in $`\lambda `$ exponentiate, as required by the general structure of the correlator in string theory.

## 5 Exponentiation

If we believe that the string description is valid at weak coupling for large $`J`$, perturbative series for $`W(C)O_J`$ must exponentiate. This is a rather non-trivial statement, since the usual expansion in Feynman diagrams computes the correlator itself, not its logarithm. Perturbative series exponentiate for Abelian Wilson loops, but in the non-Abelian theory, especially at large $`N`$, exponentiation is not at all obvious. The string-theory prediction for the sum of all planar diagrams is

$$W(C)O_J(x)\frac{1}{|x|^{2J}}\mathrm{e}^{JS(j;C)},$$
(5.35)
where $`S(j;C)`$ has a regular expansion in $`1/j^2=\lambda /J^2`$, eq. (4.31), which can be regarded as the weak-coupling expansion. If we reexpand the correlator in $`\lambda `$:

$$W(C)𝒪_J(x)\frac{1}{J!}\left(\frac{\sqrt{\lambda }l}{2|x|^2}\right)^J\underset{n=0}{\overset{\mathrm{}}{}}a_n\lambda ^n,$$
(5.36)
we find that the contribution which is least suppressed in $`1/J`$ at any given order of perturbation theory is determined by just one number, $`S_1`$:

$$a_n=\frac{S_1^n}{J^nn!}+O\left(\frac{1}{J^{n+1}}\right).$$
(5.37)
This is because the exponent in (5.35) contains an overall factor of $`J`$, and if we take the limit of large $`J`$ at a fixed order of perturbation theory, only the contribution of the leading term in the expansion of $`S(j)`$ survives. Our goal will be to prove this universality by a direct expectation of planar Feynman diagrams in the SYM theory. To do that, we will consider diagrams with largest combinatorial factors at a given order of perturbation theory.

The leading-order diagram contains $`J`$ scalar propagators that connect the local operator to $`J`$ points on the Wilson loop. In view of the standard equivalence of planar Feynman graphs and discretized two-dimensional surfaces , this diagram can be thought of as a prism with $`J`$ facets. The base of the prism is the Wilson loop and its apex is the point where operator $`𝒪_J`$ is inserted. This point is moved to infinity in our approximation. Higher-order diagrams correspond to decorating the facets with extra propagators and vertices. Until the order of perturbation theory becomes comparable to $`J`$, most of the facets will remain empty. It is clear that adding a new element to one of the empty facets produces a combinatorial factor of order $`J`$, while adding a propagator or a vertex to an already decorated facet produces a combinatorial factor of order one. Hence, diagrams with the biggest combinatorial factor have $`n`$ propagators distributed among $`n`$ different facets. Diagrams of this type with internal interaction vertices cancel by a trivial extension of the argument in appendix B. These cancellations leave only the diagrams with $`n`$ loop-to-loop propagators sandwiched between $`J`$ external legs (fig. 3).

Thus, to the leading order in $`1/J`$:

$`a_n`$ $`=`$ $`{\displaystyle \frac{(J1)!}{(2\pi l)^J}}{\displaystyle 𝑑q_1\mathrm{}𝑑q_J\theta _c(q_1,\mathrm{},q_J)}`$ (5.38)
$`\times {\displaystyle \underset{1p_1\mathrm{}p_nJ}{}}{\displaystyle _{q_{p_1}}^{q_{p_1+1}}}ds_1{\displaystyle _{s_1}^{q_{p_1+1}}}dr_1\mathrm{}{\displaystyle _{q_{p_n}}^{q_{p_n+1}}}ds_n{\displaystyle _{s_n}^{q_{p_n+1}}}dr_n`$
$`\times G(s_1,r_1)\mathrm{}G(s_n,r_n)`$
We can change the order of integration and integrate over $`q_i`$ using eq. (3.10). Then:

$`a_n`$ $`=`$ $`{\displaystyle \frac{(J1)!}{(2\pi l)^J}}{\displaystyle 𝑑s_1𝑑r_1\mathrm{}𝑑s_n𝑑r_n\theta _c(s_1,r_1,\mathrm{},s_n,r_n)G(s_1,r_1)\mathrm{}G(s_n,r_n)}`$ (5.39)
$`\times {\displaystyle \frac{J}{n}}{\displaystyle \underset{k_1+\mathrm{}+k_n=J}{}}{\displaystyle \frac{(s_2r_1)^{k_1}\mathrm{}(s_1r_n)^{k_n}}{k_1!\mathrm{}k_n!}}`$
$`=`$ $`{\displaystyle \frac{1}{n(2\pi l)^J}}{\displaystyle 𝑑s_1𝑑r_1\mathrm{}𝑑s_n𝑑r_n\theta _c(s_1,r_1,\mathrm{},s_n,r_n)G(s_1,r_1)\mathrm{}G(s_n,r_n)}`$
$`\times \left[2\pi l(r_1s_1)\mathrm{}(r_ns_n)\right]^J.`$
Applying (3.11) repeatedly $`n`$ times and keeping only the leading order in $`1/J`$, we get

$`a_n`$ $``$ $`{\displaystyle \frac{(2\pi l)^n}{n(J+1)\mathrm{}(J+n)}}{\displaystyle 𝑑s_1\mathrm{}𝑑s_n\theta _c(s_1,\mathrm{},s_n)\frac{x^{\prime \prime }{}_{}{}^{2}(s_1)}{16\pi ^2}\mathrm{}\frac{x^{\prime \prime }{}_{}{}^{2}(s_n)}{16\pi ^2}}`$ (5.40)
$``$ $`{\displaystyle \frac{1}{n!}}({\displaystyle \frac{l}{8\pi J}}{\displaystyle }dsx^{\prime \prime }{}_{}{}^{2})^n,`$
in agreement with the prediction of string theory: the diagrams with largest combinatorial weights indeed exponentiate.

## 6 Discussion

We made a rather detailed comparison between diagrams of perturbation theory and the semiclassical string amplitude which both compute a two-point correlator of an arbitrary Wilson loop with a local operator in $`𝒩=4`$ SYM theory. The results completely agree to the one-loop accuracy, as well as to all orders in perturbation theory for diagrams with the largest combinatorial weights. We should stress once again that there are no apparent reasons for such an agreement. The string theory and the perturbation theory compute different regimes, which can be most easily seen from the general form of the semiclassical string amplitude (2.6). On the string side, we have to assume that the exponent in (2.6) is large, otherwise the semiclassical approximation breaks down. To compare with SYM perturbation theory, we expand the exponential and thus assume that the exponent is small. Perhaps, one can invoke arguments similar to those of  to explain the agreement between the two calculations. The arguments rely on the observation that the ’t Hooft coupling appears only in the combination $`\lambda /J^2`$ on both sides of the correspondence, which for example means that $`\lambda /J^2`$, and not $`\lambda `$ itself, is a parameter of the planar perturbation theory. Still, a direct comparison to the string theory requires a non-trivial resummation of perturbative series.

Incidentally, we found that the string calculations are technically simpler than perturbative SYM calculations. It is relatively easy to compute the string amplitude to the two-loop order and it is definitely possible to push string calculation to higher orders in $`\lambda /J^2`$, while on the SYM side already a two-loop calculation constitutes an enormously hard problem. The simplicity of the string calculation may indicate that various diagrams cancel leaving a simple net result. We indeed observed cancellations at one loop, but those look rather accidental, at least in the way we found them.

### Acknowledgements

We are grateful to A. Tseytlin for discussions. One of us (V.P.) would like to thank the Department of Theoretical Physics at Uppsala University, where this work was partially done, for kind hospitality. This work was supported in part by STINT grant IG 2001-062. V.P. was also supported in part by INTAS grant 00-561, RFBR grant 01-02-17488a, RFBR grant for Support of Young Scientist 02-02-06517, and Russian President’s grant 00-15-99296. K.Z. was supported in part by RFBR grant 01-01-00549, and grant 00-15-96557 for the promotion of scientific schools.

## Appendix A  $`N=4`$ SYM theory

In this appendix we summarize our notations and conventions for the SYM perturbation theory. Feynman rules follow from the SYM action

$$S=\frac{N}{\lambda }d^4xtr\left\{\frac{1}{2}F_{\mu \nu }^2+\left(D_\mu \mathrm{\Phi }_i\right)^2\frac{1}{2}[\mathrm{\Phi }_i,\mathrm{\Phi }_j]^2+\mathrm{fermions}\right\}.$$
(A.1)
We do all calculations in the Feynman gauge in the coordinate representation, where the propagators are

$$\mathrm{\Phi }_i^{AB}(x_1)\mathrm{\Phi }_j^{CD}(x_2)=\frac{g^2}{8\pi ^2}\frac{\delta ^{AC}\delta ^{BD}\delta _{ij}}{|x_1x_2|^2},$$
(A.2)
$$A_\mu ^{AB}(x_1)A_\nu ^{CD}(x_2)=\frac{g^2}{8\pi ^2}\frac{\delta ^{AC}\delta ^{BD}\delta _{\mu \nu }}{|x_1x_2|^2}.$$
(A.3)
The capital letters denote $`U(N)`$ indices.

## Appendix B  Cancellation of diagrams with internal vertices

There are four types of diagrams (fig. 4) that were not taken into account in sections 3 and 5. We show here that they cancel in the large-distance asymptotic of the correlator $`W(C)𝒪_J(x)`$. We can amputate those legs of the propagators in diagrams (a), (b), and (c) which couple to the local operator and replace them by $`1/(8\pi ^2|x|^2)`$. Indeed, the region of integration with intermediate points close to $`x`$ corresponds to renormalization of the operator $`𝒪_J`$. Since the operator is chiral and is not renormalized, these contributions mutually cancel. The intermediate region of integration corresponds to descendants and contributes at higher orders in $`1/|x|`$. Thus we amputate the external legs and also take into account one half of the diagram (d). Another half participates in cancelling the renormalization of the operator $`𝒪_J`$.

We start with the diagram (d). The wave function renormalization does not vanish in the Feynman gauge, and was computed in . We can take the wave-function renormalization into account by multiplying each propagator in the tree-level amplitude by a factor $`1+F_d`$, where

$$F_d=\lambda d^4xD^2(x),$$
(B.1)
and $`D(x)`$ is appropriately regularized scalar propagator, which goes to $`1/(4\pi ^2|x|^2)`$ when regularization is removed.

The diagram (a) with external legs amputated gives a factor of

$$F_a(x_1,x_2)=\frac{\lambda }{2}d^4yd^4zD(yz)\frac{}{y^\mu }D(yx_1)\frac{}{z^\mu }D(zx_2).$$
(B.2)
for each pair of external legs. Here, $`x_1`$ and $`x_2`$ are adjacent vertices on the contour $`C`$. The diagram (b) gives

$$F_b(x_1,x_2)=\frac{\lambda }{2}d^4yD(yx_1)D(yx_2).$$
(B.3)
Finally, the diagram (c) contributes

$$F_c(x_1,x_2)=\frac{\lambda }{2}_{s_1}^{s_2}𝑑sx_{}^{}{}_{}{}^{\mu }(s)d^4yD(yx)\frac{}{y^\mu }\left(D(yx_2)D(yx_1)\right),$$
(B.4)
where $`x(s)`$ is a point on the contour between $`x_1x(s_1)`$ and $`x_2x(s_2)`$. Integration by parts yields:

$$F_c=F_d2F_b.$$
(B.5)
This equation, together with $`F_a=F_b`$, implies that

$$F_a+F_b+F_c+F_d=0.$$
(B.6)


format: text

text:
# Quantum Mechanical Motion of Relativistic Particle in Non-Continuous Spacetime

## 1 Introduction

Continuity and a metric of Lorentz signature are fundamental properties associated with spacetime in the framework of Special and General Relativity. Though intuitive, the assumption of continuity is debatable because of quantum mechanical arguments. Discrete models of spacetime have been investigated thus by various authors .

The modern idea that space could be discrete goes back to Riemann  who reflected about a natural measure of space. More recent attempts to consider spacetime as discrete either aim at resolving divergence problems in Quantum Field Theory at a fundamental level or seek to reconcile General Relativity and Quantum Theory to form a unified theory.

In this paper we investigate possible effects, which the hypothetical discreteness of spacetime could have on Quantum Theory. To this end, Feynman’s path integral approach to Quantum Tehory \[7, fey2\] provides an ideal framework since it naturally accounts for spacetime properties. For continuous spacetime, the path integral formalism has been shown to be completely equivalent to ordinary wave mechanics. However, for a non-continuous spacetime one should not necessarily expect equivalent findings. To gain insight into the possible differences between quantum mechanics in the framework of common continuous spacetime and a non-continuous spacetime model, we seek solutions for the Dirac equation.

As has been pointed out by Feynman and Hibbs , the retarded propagator of the $`1+1`$ dimensional Dirac equation

$$i\mathrm{\Psi }/t=i\sigma _z\mathrm{\Psi }/x\sigma _x\mathrm{\Psi }$$
(1)
(with units $`c=\mathrm{}/m=1`$) can be obtained from a random walk model in which the particle motion is restricted to movements either forward or backward at the speed of light. We consider this model in the framework of a particular non-continuous spacetime model  and demonstrate by an explicit calculation that the model yields common solutions of the Dirac equation.

The paper is organized as follows: Section 2 gives a short review of relevant properties of the non-continuous Minkowskain spacetime model described in . In section 3 we calculate solutions of the Dirac equation. Section 4 summarizes and discusses the results.

## 2 Non-Continuous Minkowskian Spacetime

The two dimensional spacetime model considered here is defined by the subset $`M`$ of $`𝐑^2`$

$`M=\{t,x\}`$ $`=`$ $`\{{\displaystyle \frac{n}{m}}(p^2+q^2),{\displaystyle \frac{n}{m}}(p^2q^2)\}𝐑^2,`$
$`n,m,p,q𝐙\{0\}`$
In the following, $`n,m,p,q`$ may take any value of the indicated range and are not to be considered as fixed. Elements of $`M`$ correspond to spacetime points with temporal and spatial coordinates $`(t,x)`$. $`M`$ is rational, hence countable, i.e. of cardinality $`\mathrm{}_0`$ while $`𝐑^2`$ is of cardinality $`\mathrm{}_1`$. With respect to this difference we denote $`M`$ as non-continuous.

The set $`M`$ is invariant under the transformations

$`\varphi :M`$ $``$ $`M,`$ (3)
$`\varphi (s)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1v^2}}}\left(\begin{array}{cc}1& v\\ v& 1\end{array}\right)s`$ (5)
$`=`$ $`{\displaystyle \frac{1}{2pq}}\left(\begin{array}{cc}p^2+q^2& (p^2q^2)\\ (p^2q^2)& p^2+q^2\end{array}\right)s,`$
$`p,q𝐙\{0\}.`$
The parameter

$$v=v(t,x)=\frac{x}{t}=\frac{p^2q^2}{p^2+q^2}$$
(9)
corresponds to the velocity defined as ratio of space $`x`$ and time $`t`$ displacements measured from the origin of the coordinate system. As illustrated by the second line of (5), the transformation $`\varphi `$ maps rational spacetime coordinates $`(t,x)`$ onto rational ones $`(t^{},x^{})`$.

It is straightforward to verify that (i) the elements $`\phi \varphi `$ map onto $`M`$, that (ii) that $`\varphi `$ with respect to matrix multiplications $``$ possesses group structure and that (iii) $`(\varphi ,)`$ is a subgroup of the common $`1+1`$ dimensional Lorentz group (using natural units). According to these properties and the non-continuity of $`M`$ we denote $`M`$ as non-continuous Minkowskian spacetime. In  it has been shown that in the framework of the d-space formalism  the set $`M`$ corresponds to a spacetime model which is nowhere diffeomorphic to $`𝐑^2`$ but possesses the key properties of the common macrophysical spacetime. In the same context $`M`$ may be considered as a massless solution of the Einstein equations. Note that in contrast to other discrete spacetime models (e.g.  and references therein) which intrinsically violate relativistic covariance and introduce an observationally not supported minimal length, the spacetime model considered here exhibits a generalized form of covariance. Since the set $`M`$ is dense in $`𝐑^2`$, it doesn’t introduce a minimal length as well.

The following sections are based on the spacetime model (2). We will consider light-cone coordinates $`(r,l)=(\frac{1}{2}(t+x),\frac{1}{2}(tx))`$ in order to simplify calculations. In light-cone coordinates $`M`$ becomes

$`M^{}=\{r,l\}=\{{\displaystyle \frac{n}{m}}p^2,{\displaystyle \frac{n}{m}}q^2\}𝐑^2`$ (10)
$`m,n,p,q𝐙\{0\}.`$
To end this section we note that since the set $`M`$ is dense, it is possible to formulate differential equations without referring to a calculus of finite differences.

## 3 Quadratic Checkerboard Model

In  Feynman and Hibbs described a model for the quantum mechanical motion of a relativistic electron. Solutions of the 1+1 dimensional Dirac equation are obtained by summing over all possible particle trajectories consisting of movements either forward or backward at the speed of light. Assuming natural units $`c=\mathrm{}/m=1`$, the motion of the particle corresponds to a sequence of straight path segments of slope $`\pm 45^{}`$ in the x-t plane. The retarded propagator $`\psi _{\delta \gamma }(x,t)`$ of the Dirac equation is obtained from the limiting process (e.g. )

$$\psi _{\delta \gamma }(x,t)=\underset{N\mathrm{}}{lim}A_{\delta \gamma }(ϵ)\underset{R0}{}N_{\delta \gamma }(R)(\text{i}ϵ)^R$$
(11)
$`N`$ is the number of segments with constant length $`ϵ=t/N`$ between the start point (which is assumed to be the origin of the coordinate system) and the end point $`(x,t)`$ of the path. $`R`$ denotes the number of bends while $`N_{\delta \gamma }(R)`$ stands for the total number of paths consisting of $`N`$ segments with $`R`$ bends. The indices $`\gamma `$ and $`\delta `$ correspond to the directions forward or backward at the path’s start and end points, respectively, and refer to the components of $`\psi `$. $`A_{\delta \gamma }(ϵ)`$ accounts for the appropriate normalization.

Since the early 80’ies, the checkerboard model of relativistic particle motion has been subject to anew interest (see e.g. ). In  it has been observed that paths with fixed start and end points have $`R1`$ degrees of freedom, i.e. the last bend of a path is fully specified by the location of the $`R1`$ preceding bends and the end point of the path. As well it has been shown that the normalization constant $`A_{\delta \gamma }(ϵ)1`$ if only the $`R1`$ bends enter the calculation which actually define the path of the particle. The general intuition behind this is the idea that the origin of the divergence problems of Quantum Field Theory could be related to an over-specification of the theory. Here, the corresponding interpretation of the checkerboard model and the calculation scheme described in  are adapted.

In the following we demonstrate by an explicit calculation that expression (11) when generalized to account for the key properties of $`M`$ yields solutions of the Dirac equation. Before starting with the calculation, some remarks about the limiting process (11) and the terminology are in place. First observe the structure of (11). $`\psi _{\delta \gamma }(x,t)`$ is defined in the limit $`N\mathrm{}`$ where $`N`$ is the number of path segments a path has. The length of each of these path segments is $`ϵ=t/N`$. Bends of the particle trajectory occur only at boundaries of path segments. This is equivalent to consider a particle moving on a rectangular spacetime lattice with equal spacing length (or ’resolution’) $`ϵ`$. Consider now the expression

$$\underset{R0}{}N_{\delta \gamma }(R)(\text{i}ϵ)^R.$$
(12)
which is subject to the limit $`N\mathrm{}`$. For a given $`N`$, it is weighted sum of the number $`N_{\delta \gamma }(R)`$ of possible paths with $`N`$ segments and $`R0`$ bends linking the start and the end point under consideration. The complex weight (amplitude) each path is contributing is $`(\text{i}ϵ)^R`$, i.e. each bend of a path contributes $`(\text{i}ϵ)`$. As has been demonstrated (e.g. ), taking the limit $`N\mathrm{}`$ yields solutions of the Dirac equation.

When adapting Feynman’s checkerboard model to the discrete spacetime $`M`$ we will follow as closely as possible the scheme outlined above. Two differences have to be taken into account, the first of which being a generalization. First, it is supposed that path segments can be of non-uniform length $`ϵ_n\text{const.}`$ and that the contribution of each bend to the overall amplitude of a path is proportional to the segment length $`ϵ_n`$ immediately preceding the bend. This leads to the generalized form

$$\psi _{\delta \gamma }(x,t)=\underset{N\mathrm{}}{lim}A_{\delta \gamma }(ϵ)\underset{R0}{}N_{\delta \gamma }(R)(\text{i})^R\underset{n=1}{\overset{R}{}}ϵ_n$$
(13)
of (11). Second, only particle paths on $`M`$ will be considered. These path are characterized by bends occurring at spacetime points $`(t,x)`$ satisfying

$$\frac{x}{t}=\frac{p^2q^2}{p^2+q^2}.$$
(14)
In terms of $`M^{}`$, the light-cone representation of $`M`$, this condition becomes

$$\frac{r}{l}=\frac{p^2}{q^2}$$
(15)
where $`(r,l)=(\frac{1}{2}(t+x),\frac{1}{2}(tx))`$. Points of the spacetime lattice should account for restriction (15). We achieved this by considering quadratic spacetime lattices whose lattice points $`(j_r,j_l)𝐍^2`$ have, up to scaling, spacetime lightcone coordinates $`(r,l)=(j_r^2,j_l^2)`$. The spacetime coordinates of lattice points of this kind naturally account for condition (15). Figure 1 shows an example for a path with $`N=8`$ segments. With respect to the modified properties of the spacetime lattice, the model is denoted as quadratic checkerboard model.

From the definition of the quadratic spacetime lattice (or directly from Figure 1) it is evident that path segments $`ϵ_n`$ are of non-uniform length and are depending on the lattice location $`(j_r,j_l)`$. We only will be interested in the length of paths segment preceding a bend and bends of a path to the right and to the left will be considered separately in the following. Observing that the length of a path segment preceding a bend to the right (left) is uniquely specifying by the lattice coordinate $`j_r`$ ($`j_l`$), the indexes $`r,l`$ can be dropped. Then, the length of a path segment immediately following a bend can be written as

$$ϵ_j=(2j1)ϵ_0$$
(16)
The constant $`ϵ_0`$ accounts for the scaling of the path. It depends on the end point $`(t,x)`$ and on the number of segments $`N`$ a path has. For an explicit expression of $`ϵ_0`$ in terms of the number of segments to the left and right, see equation (29). Note that $`ϵ_0`$ scales not just the length of paths segment followed by a bend but the whole quadratic spacetime lattice. From this it follows immediately that for $`N\mathrm{}`$ the spacetime lattice under consideration becomes a rectangular subset of $`M`$.

In line with the generalization (13) of Feynman’s checkerboard model suppose now that each bend defining a path on $`M`$ contributes an amplitude proportional to the length of the path segment immediately preceding the bend. Again, bends to the right and left are considered separately and the following notation is adopted: The index $`n`$ enumerates the segments of a path followed by a bend to the right (left). $`j_n`$ stands for the lattice coordinate $`j_r`$ ($`j_l`$) of the segment and the set $`\{j_n\}`$ indicates the path segments, after which bends to the right (left) occur. For a schematic example consider Figure 1. The length $`ϵ_{j_n}`$ of the path segment followed by the $`n^{th}`$ bend to the right (left) according to (16) becomes

$$ϵ_{j_n}=(2j_n1)ϵ_0$$
(17)
and the corresponding contribution of the $`n^{th}`$ bend to the overall amplitude is

$$\varphi _{j_n}=\text{i}ϵ_{j_n}=(2j_n1)ϵ_0$$
(18)
The total amplitude of a path contributed by bends to the right (left) is given by the product

$$\varphi =\underset{n}{}(\text{i}ϵ_{j_n})$$
(19)
Having established this notation, we start with the evaluation of (13). As in the ’linear’ case () a path with $`R`$ bends that starts with a positive velocity (i.e. to the right) and ends with a negative velocity (i.e. to the left) consists of exactly $`(R1)/2+1`$ bends to the left and $`(R1)/2`$ bends to the right. The $`(R1)/2`$ bends to the right can occur after an arbitrary path segment to the left. $`(R1)/2`$ of the $`(R1)/2+1`$ bends to the left occur in the same manner after path segments to the right while the additional bend to the left must occur after the last segment to the right. Let $`P`$ denote the total number of path segments to the right $`(+)`$ and $`Q`$ those to the left $`()`$. In total, the path has $`N=(P+Q)`$ segments. The contribution of the $`R^+=(R1)/2`$ bends to the right $`\psi _+`$ is

$$\psi _+(R^+)=\underset{j_1<\mathrm{}<j_{R^+}}{\overset{P1}{}}(2j_11)\mathrm{}(2j_n1)(\text{i}ϵ_0)^{R^+}$$
(20)
Next consider the situation where the path consists of a large number of segments to the right, i.e. $`P1`$ or equivalently $`N\mathrm{}`$. This limiting process corresponds to the spacetime lattice becoming dense. For $`P1`$, $`\psi _+(R^+)`$ is approximated by

$`\psi _+(R^+)`$ $``$ $`{\displaystyle \frac{1}{R^+!}}{\displaystyle \underset{j_1\mathrm{}j_{R^+}}{\overset{P}{}}}(j_1\mathrm{}j_n)2^{R^+}(\text{i}ϵ_0)^{R^+}`$ (21)
$``$ $`{\displaystyle \frac{2^{R^+}(\text{i}ϵ_0)^{R^+}}{R^+!}}\left({\displaystyle \underset{j=1}{\overset{P}{}}}j\right)^{R^+}`$ (22)
$``$ $`{\displaystyle \frac{2^{R^+}(\text{i}ϵ_0)^{R^+}}{R^+!}}\left({\displaystyle \frac{P^2}{2}}\right)^{R^+}`$ (23)
$`=`$ $`{\displaystyle \frac{P^{2R^+}(\text{i}ϵ_0)^{R^+}}{R^+!}}`$ (24)
The contribution of the $`R^{}=(R1)/2+1`$ bends to the left is calculated similarly. The additional bend (occurring after the last segment to the right) does not enter the calculation since a path is fully determined by the $`R1`$ bends to the right and left, respectively. We find

$`\psi _+(R^{})`$ $``$ $`{\displaystyle \frac{2^{(R^{}1)}(\text{i}ϵ_0)^{(R^{}1)}}{(R^{}1)!}}\left({\displaystyle \frac{Q^2}{2}}\right)^{(R^{}1)}`$ (25)
$`=`$ $`{\displaystyle \frac{Q^{2(R^{}1)}(\text{i}ϵ_0)^{(R^{}1)}}{(R^{}1)!}}`$ (26)
To get the total contribution of the $`R^{}`$ bends to the left and the $`R^+`$ bends to the right, expressions $`\psi _+(R^{})`$ and $`\psi _+(R^{})`$ are multiplied yielding

$$\psi _+(R)=(\text{i}ϵ_0)^{R1}\frac{(PQ)^{R1}}{[((R1)/2)!]^2}$$
(27)
which represents the contribution of all paths with exactly $`R`$ bends that start to the right and end to the left. Finally, the summation over all possible paths that start to the right and end to the left is preformed by summing $`\psi _+(R)`$ over the total number of bends $`R`$. This leads to

$$\psi _+=\underset{\text{odd }R}{}(\text{i}ϵ_0)^{R1}\frac{(PQ)^{R1}}{[((R1)/2)!]^2}$$
(28)
where

$$ϵ_0=\frac{t}{P^2+Q^2}$$
(29)
accounts for the proper scaling of the lattice. Notice that for $`N\mathrm{}`$ the scaling factor $`ϵ_00`$. Defining the classical velocity of the particle $`v=\mathrm{\Delta }x/\mathrm{\Delta }t=x/t=(P^2Q^2)/(P^2+Q^2)`$ we find $`PQ=(P^2+Q^2)/2\gamma `$ where $`\gamma =1/\sqrt{1v^2}`$. Finally $`\psi _+`$ becomes

$`\psi _+`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(1)^k{\displaystyle \frac{(t/2\gamma )^{2k}}{[(k)!]^2}}`$ (30)
$`=`$ $`J_0(t/\gamma )`$ (31)
A similar calculation of $`\psi _+`$ leads to the same result. This can be seen from interchanging the roles of $`P`$, $`Q`$, $`R^{}`$ and $`R^+`$.

For $`\psi _{++}`$, the number of bends to the right and to the left is $`R/2`$ for each direction where $`R`$ is even. However, the path again is defined by $`R^+=R/2`$ bends to the right and $`R^{}=R/21`$ to the left. Thus,

$`\psi _{++}`$ $`=`$ $`{\displaystyle \underset{\text{even }R}{}}(\text{i}ϵ_0)^{R1}{\displaystyle \frac{P^RQ^{R2}}{(R/2)!(R/21)!}}`$ (32)
$`=`$ $`{\displaystyle \underset{0,2,4,\mathrm{}}{}}(\text{i}ϵ_0)^{R+1}{\displaystyle \frac{P^2(PQ)^R}{(R/2+1)!(R/2)!}}`$ (33)
$`=`$ $`{\displaystyle \frac{P}{Q}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(1)^k{\displaystyle \frac{(PQϵ_0)^{2k+1}}{(k+1)!(k)!}}`$ (34)
$`=`$ $`\text{i}{\displaystyle \frac{P}{Q}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}(1)^k{\displaystyle \frac{(t/2\gamma )^{2k+1}}{(k+1)!(k)!}}`$ (35)
$`=`$ $`\text{i}{\displaystyle \frac{P}{Q}}J_1(t/\gamma ).`$ (36)
With $`P/Q=(x+t)/(t^2x^2)^{1/2}`$ the component $`\psi _{++}`$ becomes

$$\psi _{++}=\text{i}\frac{(t+x)}{t}\gamma J_1(t/\gamma ).$$
(37)
A similar calculation for $`\psi _{}`$ leads to

$$\psi _{}=\text{i}\frac{(tx)}{t}\gamma J_1(t/\gamma ).$$
(38)
which completes the calculation.

To relate the components $`\psi _{\delta \gamma }`$ to the Dirac equation (1) consider its explicit representation with

$$\sigma _x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _z=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(39)
In this representation, $`\psi _1`$ and $`\psi _2`$ defined as

$$\psi _1=\left(\begin{array}{c}\psi _{++}\\ \psi _+\end{array}\right),\psi _2=\left(\begin{array}{c}\psi _+\\ \psi _{}\end{array}\right)$$
(40)
are two common independent solutions of the Dirac equation (1). Since $`M`$ is dense in $`𝐑^2`$, the differential form of (1) does not pose a problem even though $`M`$ is non-continuous.

The four components $`\psi _{\delta \gamma }`$ are closely related to the matrix elements of the retarded continuum propagator of the Dirac equation (see e.g. ). Despite this there is no full equivalence and a physical interpretation of $`\psi _{\delta \gamma }`$ in terms of propagator components is not straightforward. In the last section we will come back to this and related questions in more detail.

As a side remark note that for the applied calculation scheme it is essential to consider only those bends which actually define the path of the particle. While it is straightforward to ’renormalize’ Feynman’s original model if all bends are considered, this is not the case here. Taking into account all bends, the normalization ‘constant’ $`A_{\delta \gamma }(ϵ_0)`$ for the checkerboard restricted to $`M`$ turns out to not only depend on $`ϵ_0`$ but also on $`\delta `$, $`\gamma `$ and the classical velocity $`v`$.

## 4 Discussion and Conclusions

It is counter-intuitive that Feynman’s spacetime based path integral approach to quantum theory yields in the case of the discussed example intriguing similarities between the common continuous and the discrete spacetime model $`M`$. The solutions $`\psi `$ obtained from the generalized checkerboard model and the common retarded continuum propagator of the Dirac equation are of the same form.

However, despite the same functional form, the components $`\psi _{\delta \gamma }`$ are not fully equivalent to the retarded continuum propagator. Reworking the outlined calculation with physical units $`\mathrm{}/m1`$ shows the dimensionality of $`\psi _{\delta \gamma }`$ to be $`lenght^0`$ which is different from the dimensionality $`lenght^1`$ of the retarded continuum propagator. With respect to the original checkerboard model, the difference is related to two origins: First, taking into account only $`R1`$ bends means accounting for one factor $`mϵ`$ less. Second, Feynman’s original checkerboard model is ’renormalized’ by dividing components like e.g. (28) and (37) by $`2ϵ`$ in order to yield the retarded continuum propagator.

The dimensionality $`lenght^0`$ of $`\psi _{\delta \gamma }`$ may be seen as suggesting a picture consistent with the cardinality $`\mathrm{}_0`$ of $`M`$. Profoundly taking into account the cardinality $`\mathrm{}_0`$, a continuum-like propagator is inconsistent with $`M`$. This is due to the fact that because of its cardinality there exists no measure on $`M`$. As a consequence, the concept of integration is not available. Instead, one can expect (infinite) sums to play the role of integration. This suggests that the components of a discrete counterpart of the continuum propagator (i.e. $`\psi _{\delta \gamma }`$) should have dimensionality $`lenght^0`$. Fully clarifying the physical meaning of $`\psi _{\delta \gamma }`$ in the framework of the discrete spacetime model $`M`$ is not straightforward. One is lead into number theoretical problems that are not easy to tackle.

Other open questions remain. For example: To what extend does the path integral approach capture the relation between spacetime and quantum mechanics? This question is raised not only by the present results but already by Feynman’s original checkerboard model. After all, according to its formulation, bends occur only at spacetime points with rational coordinates. As well, the velocity spectrum of the particle is rational. In this sense, Feynman’s checkerboard model is of cardinality $`\mathrm{}_0`$ and thus discrete in the same manner as $`M`$ is. To our knowledge this is a fact that the literature has not paid attention to.

In  it has been shown that the solution of the Dirac equation can be understood in terms of bit strings detached from any notion of spacetime. The present results seem to point into the same direction, namely, that there is (at least in the special case of the Dirac equation) only a loose connection between spacetime properties and quantum mechanics, except probably for the signature of the metric.

Other open issues are related to the spacetime model $`M`$ itself. The ratio $`x/t`$ interpreted as classical velocity of the particle (starting from the origin) implies a discrete velocity spectrum

$$\{v\}=\frac{p^2q^2}{p^2+q^2},p,q𝐙\{0\}.$$
(41)
What are the implications of the non-continuous spectrum (41)? Progress in this and other directions again is complicated by the fact that $`M`$ is of cardinality $`\mathrm{}_0`$. Finally, a question not addressed here is where the bends of a path occur, that contribute most to the overall amplitude. Because of the structure of $`M`$ it is not evident that the locations are distributed as in the continuous case.

In summary, it has been shown that Feynman’s path integral approach restricted to the subset $`M`$ of $`𝐑^2`$ recovers common solutions of the Dirac equation. The particle path on $`M`$ is limited to bends at positions $`i`$ with coordinate ratios $`x_i/t_i=(p^2q^2)/(p^2+q^2)`$ where $`p,q𝐙\{0\}`$. The bends of paths on $`M`$ thus occur at a subset of the rational space and time coordinates only. It should be noted that both the characteristics of $`M`$ and the discrete velocity spectrum (41) are invariant under the Lorentz subgroup $`\varphi `$. In particular, bends of paths occur only at rational space and time coordinates and the velocity remains rational in all coordinate systems. In this sense the model exhibits a generalized form of Lorentz invariance missing in Feynman’s original checkerboard model and most other lattice models of spacetime. This demonstrates that there is a discrete, two-dimensional spacetime model of cardinality $`\mathrm{}_0`$ accounting for some of the key properties of special relativity and quantum mechanical features of relativistic particle motion.

## Acknowledgements

The author thanks an unknown referee for helpful comments.


text:
# Hierarchical population model with a carrying capacity distribution

## 1 Introduction of the model

As a concrete example, consider bacteria on a line with local population $`N(x,t)`$ defined on $`x[0,1]`$. The population is normalized, so that $`N`$ is to be interpreted as a density rather than the number of bacteria. Generalization to two or three space dimensions is straightforward, as will become clear from the simple structure of the model. Assume, for simplicity, a homogeneous initial condition $`N_0=N(x,0)>0`$, representing inoculation along the line of a thin uniform layer. At $`t>0`$ the growth process is characterized by a probability of birth, $`B`$, which is applied to segments or “patches” along the line. The stochastic character of the growth reflects that the surface is inhomogeneous. Some areas are rich in nutrient (e.g., agar), others not. The division process is effective only in nutrient-rich patches. These patches are assumed to be of characteristic size $`1/\lambda `$, where $`\lambda `$ is a scale factor central to the model. Alternative ways of inducing spatially non-uniform growth or depletion are to apply, for photosynthetic bacteria, favourable or unfavourable (UV) illumination according to a structured or random spatial pattern .

Besides nutrient non-uniformity, the environmental conditions are locally modified in a random way, for example by spraying drops of an antibiotic which subsequently spreads by diffusion or transport resulting from bacterial motility. In places with low chemical concentration the population stagnates, elsewhere it shrinks. The action of the bactericidal product is expressed through a probability of death, $`D`$. An interesting variant on antibiotics, and relevant in particular to inhibition of biofilm formation, is a nutrient remover like, for example, lactoferrin. This protein, which is abundant in human external secretions, mops up traces of iron and thereby deprives bacteria of an essential substance. Applied in concentrations below those that kill or prevent growth, it just slightly increases the division time of the bacteria and enhances their surface motility, causing them to “twitch” or wander around seeking food. As a result, microcolonies do not form and matrix-encased communities specialized for surface persistence (biofilms) do not develop .

The model is composed as follows. The characteristic size of the spatial inhomogeneity of the nutrient is, for simplicity, taken to be the same as that of the applied antimicrobial perturbation. The (unit) line is thus divided in $`\lambda `$ segments, each of which has the size of, say, an antibiotic drop. The probability of birth $`B`$ is applied to each segment, and subsequently the probability of death, $`D`$. Each segment is thus visited twice. We have $`0B,D1`$, and without loss of generality we take $`\lambda =3`$ in this paper.

One generation later this process is repeated, but the maximum local increase or decrease of the population is reduced by a factor of $`\lambda `$. This reduction is repeated in every generation. In this way the growth process eventually comes to a halt, and the population approaches a carrying capacity, as is commonly the case in models for a single species . This rescaling of the growth or the depletion is a crucial ingredient of the model and leads to a hierarchical build-up of the population. While introduced here heuristically, various justifications can be invoked for this feature. It is known that bacterial reproductivity and/or susceptibility to antimicrobial agents can “wear out” in a systematic manner, uniformly in space. For example, the effective birth rate decreases significantly by varying the temperature some 10 degrees away from optimal growth conditions. Also, once a biofilm starts forming, it becomes resistant to nutrient removers like lactoferrin, and it notoriously resists killing by antibiotics.

Further, the rescaling by a factor $`\lambda `$ is applied not only to the increase or decrease of the local population, but also to the typical size of the nutrient-rich or bactericidal patches. The justification for this runs as follows. Since after one generation (of the order of two hours in real time) the nutrient or chemical contamination has been redistributed in space, by diffusion or by bacterial motion, regions which were initially rich (or pure) will show poor (or toxic) spots, while initially barren regions will display small fertile areas. Therefore, we apply the random birth and death rules on a smaller length scale, reduced by a factor of $`\lambda `$ with respect to the previous generation, and repeat this reduction generation after generation.

In sum, every new generation the population fluctuations are reduced by cooling (or heating) and biofilm development, and they occur on a smaller length scale due to diffusive nutrient/antibiotic fragmentation. A single scale factor controls - in this greatly simplified model - all these dynamical complications, so that neither external time-dependent fields nor “interparticle” interactions occur explicitly in what follows. Repeating this process saturation is reached exponentially rapidly, typically in less than 10 generations, as the calculations show. The resulting local population density for long times, $`N(x,\mathrm{})`$, is highly inhomogeneous and largely reflects the initial random distribution of the nutrient distribution and the toxic contamination, brought about by spraying. Considered mathematically, $`N(x,\mathrm{})`$ is a fractal curve with infinite length, but trivial fractal dimension $`d_F=1`$. Figure 1 shows an example of $`N(x,5)`$, the local population after 5 generations, assuming $`B=0.5`$ and $`D=0.2`$.

Note that the natural death of the organisms is not included explicitly in the model, for the simple reason that the induced saturation occurs on a time scale (several hours) significantly shorter than the natural lifetime. In other words, the probability $`D`$ pertains to death caused by bactericidal products.

In Fig.1 it can be seen that what happens in the first generation dominates the future evolution. For example, the green square represents the massive growth on a nutrient-rich segment, $`(2/3,1)`$, while the red square indicates the result of the killing of most of the bacteria on a toxic patch, $`(1/3,2/3)`$. In the first segment, $`(0,1/3)`$, nothing happened in the first generation. After a few generations the population “landscape” consists of towers and deep crests, which get rougher as time progresses, but which hardly change size anymore.

This model is mathematically equivalent to the previously introduced hierarchical random deposition process  with rescaling factor $`\lambda `$, and probabilities $`P=B(1D)`$ for depositing a hill, and $`Q=D(1B)`$ for digging a hole, with $`P+Q1`$. The correspondence between $`(B,D)`$ and $`(P,Q)`$ can be seen by taking into account that in the biological application every segment is visited twice, applying $`B`$ in the first visit followed by $`D`$ in the second, whereas in the deposition model there is only one visit per segment, per generation. Therefore, after a complete visit the population is increased provided offspring was produced and no deaths occurred, whence the product $`B(1D)`$. The population is decreased provided no divisions occurred and the antibiotics killed a fraction of the existing population, whence $`(1B)D`$. Finally, the population remains constant if either offspring was produced but subsequently killed, or no births took place and the antibiotics were absent or had no effect. The probability for this is $`BD+(1B)(1D)=1(P+Q)`$. Note that the total probability for the three outcomes is unity.

Using the deposition model it is easy to derive analytical results in terms of the variables $`P`$ and $`Q`$, and also numerical simulation is straightforward. In generation $`n`$ the process requires just $`\lambda ^n`$ (pseudo-)random numbers. For each segment of width $`\lambda ^n`$ along $`x`$ the population is increased by an amount equal to the segment size (in dimensionless units) with probability $`P`$, or reduced by that amount with probability $`Q`$, or left unchanged with probability $`1(P+Q)`$. This leads to the characteristic rapidly converging population landscape formed by squares of decreasing size (Fig.1).

The initial condition is taken to be $`N_0=1/(\lambda 1)`$, independent of $`x`$ (uniform inoculation). The first quantity of interest is the mean carrying capacity, $`\overline{K}`$, which is the ensemble average of $`N(x,\mathrm{})`$ over all possible random realizations of the process. Since on average the population grows by an amount $`(BD)\lambda ^n`$ in generation $`n`$, we obtain

$$\overline{K}=(1+BD)/(\lambda 1)=(1+PQ)/(\lambda 1),$$
(1.1)
independent of $`x`$. Clearly, the initial condition $`N_0`$ has been chosen large enough so that, for all $`B`$ and $`D`$, the population remains positive for all $`x`$ and $`t`$, and approaches zero only in the extinction limit $`B=0`$, $`D=1`$. Note the simple identity $`BD=PQ`$, making the relations between the growth probabilities and the deposition parameters more transparent.

Generalization of the model to higher dimensions (e.g., substrate dimension 2) is not pursued here, but the results of such extension can easily be anticipated by examining the landscapes shown for the hierarchical deposition model  on planar substrates. In fact, if the bacterial population is compact, the value of $`N(x,t)`$ is simply proportional to the height of the biofilm above the substrate, and the landscapes can be interpreted as three-dimensional in real space. This correspondence holds as long as there are no “overhangs” in the biofilm, and thus works for towering pillar structures but not for mushroom-shaped ones . The similarity between the three-dimensional characteristic landscapes generated in this model, and the structure of its horizontal cross-section (see Fig.6 in that reference), and images of vertical and horizontal sections of a towering pillar biofilm, as obtained by confocal laser scanning microscopy, is remarkable. This suggests that a hierarchical model may be a reasonable first approximation for describing biofilm growth.

The connection between the model and the familiar Verhulst process is elucidated in Section 2. It is made clear that $`B`$ and $`D`$ are biologically relevant parameters determining the carrying capacity $`K`$, whereas $`\lambda `$ simply sets an overall time scale factor. The carrying capacity distribution, the key result of this random growth model, is calculated analytically in Section 3, and concrete examples are discussed in Section 4, together with a check against numerical simulation. The self-similar character of this $`K`$-distribution is demonstrated and the effect of $`B`$ and $`D`$ on the standard deviation of the population is examined in representative cases. Section 5 deals with the most probable carrying capacity and its probability, considered as a function of $`B`$ and $`D`$. Insight is gained in the abundance of large deviations away from the most probable $`K`$. The final section, 6, presents our conclusions and addresses a possible experimental test of this model.

## 2 Connection with a Verhulst process or logistic growth

In this section we are concerned with the uniform, i.e., spatially averaged, population. The time evolution of the mean, or “quenched average”, of $`N(x,t)`$ is given by the change in the mean value, $`\overline{N}(t)`$, from $`t=n`$ to $`t=n+1`$,

$$\overline{N}(n+1)=\overline{N}(n)+\lambda ^{n1}(BD),$$
(2.1)
which is solved by

$$\overline{N}(n)=N_0+(1\lambda ^n)(BD)/(\lambda 1)$$
(2.2)
This is a simple growth, converging exponentially rapidly, for $`n\mathrm{}`$, to the mean carrying capacity. It is useful now to work with the differences

$$\mathrm{\Delta }\overline{N}=\overline{N}N_0,\text{and}\mathrm{\Delta }\overline{K}=\overline{K}N_0,$$
(2.3)
with respect to the initial population $`N_0=1/(\lambda 1)`$. We obtain

$$\mathrm{\Delta }\overline{N}(n)=(1\lambda ^n)\mathrm{\Delta }\overline{K}$$
(2.4)
Using this in (II.1) leads to

$$\mathrm{\Delta }\overline{N}(n+1)\mathrm{\Delta }\overline{N}(n)=\frac{\lambda 1}{\lambda }\frac{\lambda ^n}{1\lambda ^n}\mathrm{\Delta }\overline{N}(n)$$
(2.5)
For large $`n`$, we can neglect $`\lambda ^n`$ compared to 1 in the denominator, and again using (II.4), we get the nonlinear difference equation

$$\mathrm{\Delta }\overline{N}(n+1)\mathrm{\Delta }\overline{N}(n)\frac{\lambda 1}{\lambda }\mathrm{\Delta }\overline{N}(n)(1\frac{\mathrm{\Delta }\overline{N}(n)}{\mathrm{\Delta }\overline{K}})$$
(2.6)
This is clearly a discrete Verhulst process, the continuum time limit of which is of the form of the Verhulst equation

$$\frac{d\mathrm{\Delta }\overline{N}(t)}{dt}=\frac{\lambda 1}{\lambda }\mathrm{\Delta }\overline{N}(t)(1\frac{\mathrm{\Delta }\overline{N}(t)}{\mathrm{\Delta }\overline{K}})$$
(2.7)
This analogy allows us to identify and interpret the model parameters unambiguously. The characteristic time scale of the model is $`\lambda /(\lambda 1)`$, a trivial constant. Therefore, the length rescaling factor $`\lambda `$ in the model has no particular biological relevance, and can be given an arbitrary integer value $`\lambda >1`$. The carrying capacity (difference) given by $`\mathrm{\Delta }\overline{K}=(BD)/(\lambda 1)`$ is an important property, and we find that it is determined simply by the difference of the birth and death rates. This link between a phenomenological parameter of the Verhulst equation, the carrying capacity, and the stochastic “microscopic” model parameters $`B`$ and $`D`$ is a key ingredient of our model. Finally, the effective birth rate $`1/\tau `$ depends on the size of the population, and eventually vanishes as the carrying capacity is approached, as is typical of logistic growth. Asymptotically for long times, it is given by

$$1/\tau \frac{\lambda 1}{\lambda }(1\frac{\mathrm{\Delta }\overline{N}(t)}{\mathrm{\Delta }\overline{K}})=\frac{\lambda 1}{\lambda }\lambda ^n$$
(2.8)
More precisely, the exact expression for all times is

$$1/\tau =\frac{\lambda 1}{\lambda }\frac{\lambda ^n}{1\lambda ^n}$$
(2.9)
Note, and we wish to stress, that for our process the Verhulst equation is nothing but an easily interpretable nonlinear approximation to the actual linear differential equation describing our dynamics,

$$\frac{d\mathrm{\Delta }\overline{N}(t)}{dt}=\frac{\lambda 1}{\lambda }\mathrm{\Delta }\overline{N}(t)(1\frac{\mathrm{\Delta }\overline{N}(t)}{\mathrm{\Delta }\overline{K}})\frac{\mathrm{\Delta }\overline{K}}{\mathrm{\Delta }\overline{N}(t)}=\frac{\lambda 1}{\lambda }(\mathrm{\Delta }\overline{K}\mathrm{\Delta }\overline{N}(t))$$
(2.10)
The precise way (in our case exponential) in which the effective birth rate vanishes is only one possible way in which this rate in a biological system can be driven away from its optimal value. We envisage that this drift is achieved by a temperature variation, upwards or downwards, bringing the bacterial growth to a halt on a time scale of a small number of generations. In the introduction we already alluded to the fact that biofilm development can be a factor leading to resistance to antibiotics, which can lead to a vanishing death rate. We simply assume that both effects (temperature change and biofilm formation) are present uniformly throughout the sample. Our justification for assuming a uniform external perturbation responsible for the saturation is that in our random model the effective birth rate (II.9) does not depend on the local ($`x`$-dependent) population, but only on time. Alternatively, the self-limitation of the population could of course also be due to local overpopulation, but this model is in its present form too simple to allow for an internal local feed-back process.

## 3 The carrying capacity distribution

Even though the average growth of the population is characterized by trivial exponential saturation, the local population $`N(x,t)`$ for a particular realization of the quenched randomness typically shows interesting large fluctuations. These deviations, highly nonuniform in space, can possibly be measured experimentally by probing the local population of bacteria after spraying a few drops of an antibiotic. The fluctuations of the total population $`𝒩(t)=_0^1N(x,t)𝑑x`$ are studied in this section and the next two.

Our aim in this section is to explore the effect of random nutrient inhomogeneity and random antibiotic spray on the carrying capacity $`K`$, which is the asymptotic value for long times of the total population. The simple stochastic nature of the model allows to obtain the full $`K`$-distribution in analytic form in terms of the probability of birth $`B`$ and that of death $`D`$, with modest mathematical effort. Consider a particular random evolution. In generation $`m`$ the total population grows by an amount $`(\lambda ^mk_m)/\lambda ^{2m}`$ where the possible values of the integer $`k_m`$ are contained in the set $`\{0,1,\mathrm{},2\lambda ^m\}`$. For example, in the deterministic limits $`B=1`$ and $`D=0`$, one always has $`k_m=0`$, and for $`B=0`$ and $`D=1`$ one invariably encounters $`k_m=2\lambda ^m`$. The set $`\{k_1,\mathrm{},k_n\}`$ will be referred to as a “growth sequence” of length $`n`$. Of course, growth sequences are in general degenerate in the sense that different population landscapes can possess the same growth sequence. The total population after $`n`$ generations, denoted by $`𝒩(n)`$ is then given by

$$𝒩(n)=\frac{1}{\lambda 1}+\underset{m=1}{\overset{n}{}}\frac{\lambda ^mk_m}{\lambda ^{2m}}$$
(3.1)
At this point it is useful to note that different growth sequences may accidentally result in the same value for the total population. For example, a single local increase of $`N`$ in one generation can be deleted by $`\lambda ^2`$ local decreases in the next generation. Therefore, an alternative sequence in which nothing changes in those two generations would lead to the same total population.

If we work with the deposition probabilities, related to $`B`$ and $`D`$ by $`P=B(1D)`$ and $`Q=D(1B)`$, and define $`S=P+Q`$, the probability for depositing in generation $`m`$ any configuration with, say, $`h`$ holes, $`f`$ flat segments (neither hill nor hole), and $`\lambda ^mfh`$ hills is given by the following expression which takes into account the number of ways in which the holes and flat segments can be put in the standard combinatorial way,

$$𝒫_m(h,f)=\left(\begin{array}{c}\lambda ^m\\ h\end{array}\right)\left(\begin{array}{c}\lambda ^mh\\ f\end{array}\right)P^{\lambda ^mfh}(1S)^fQ^h$$
(3.2)
Note that summing this expression over all possible $`h`$ and $`f`$ gives the total probability,

$$\underset{h=0}{\overset{\lambda ^m}{}}\underset{f=0}{\overset{\lambda ^mh}{}}𝒫_m(h,f)=1$$
(3.3)
which is unity, as it should.

The probability for having a fixed population increment, that is, a particular value $`k_m=k`$, in generation $`m`$, is given by the following partial sum over the $`𝒫_m(h,f)`$,

$$\widehat{𝒫}_m(k)=\underset{h=j(k)}{\overset{[k/2]}{}}\left(\begin{array}{c}\lambda ^m\\ h\end{array}\right)\left(\begin{array}{c}\lambda ^mh\\ k2h\end{array}\right)P^{\lambda ^mk+h}(1S)^{k2h}Q^h$$
(3.4)
Note that $`k`$ equals the number of flat segments plus twice the number of holes, so that $`f=k2h`$. Here, $`[k/2]`$ is equal to $`k/2`$ if $`k`$ is even, and equal to $`(k1)/2`$ if $`k`$ is odd. Further, the integer $`j(k)`$ equals 0 if $`k\lambda ^m`$, and $`j(k)=k\lambda ^m`$ for $`k>\lambda ^m`$. One can verify the normalization

$$\underset{k=0}{\overset{2\lambda ^m}{}}\widehat{𝒫}_m(k)=1$$
(3.5)
Since the $`k_m`$, for different generations $`m`$, are independent variables, the probability of a given growth sequence is the product

$$\widehat{\pi }_n(𝐤)\underset{m=1}{\overset{n}{}}\widehat{𝒫}_m(k_m),$$
(3.6)
where $`𝐤=(k_1,k_2,\mathrm{},k_n)`$. Consequently, the probability $`\mathrm{\Pi }`$ of realizing the total population $`𝒩(n)`$ is the sum of the probabilities of all the growth sequences $`𝐤`$ which produce the same value $`𝒩(n)=𝒩`$,

$$\mathrm{\Pi }_n(𝒩)=\underset{𝐤}{\overset{(𝒩)}{}}\widehat{\pi }_n(𝐤),$$
(3.7)
where the sum takes care of the degeneracy of $`𝒩`$. In practice this sum contains only a few terms, because typically not many growth sequences lead to the same value for the total population. For large $`n`$ the populations $`𝒩(n)`$ converge to the carrying capacities $`K`$, so that the distributions $`\mathrm{\Pi }_n(𝒩)`$ converge to the carrying capacity distribution $`\mathrm{\Pi }_{\mathrm{}}(K)`$, or in short $`\mathrm{\Pi }(K)`$, which we set out to obtain.

Note that the total probability for the entire process up till any given generation $`n`$ is correctly normalized, since

$$1=\underset{m=1}{\overset{n}{}}\underset{k=0}{\overset{2\lambda ^m}{}}\widehat{𝒫}_m(k)=\underset{𝐤}{}\underset{m=1}{\overset{n}{}}\widehat{𝒫}_m(k_m)=\underset{𝐤}{}\widehat{\pi }_n(𝐤)=\underset{𝒩}{}\mathrm{\Pi }_n(𝒩)$$
(3.8)
In particular,

$$\underset{K}{}\mathrm{\Pi }(K)=1$$
(3.9)
In closing this section a few technical remarks are in order. Clearly, in the special cases $`P=0`$, $`Q=0`$, or $`P+Q=1`$ the trinomial formulae we have derived, in particular (III.2) and (III.4), cannot be used in the present form. Instead, simple binomial expressions should be written down in those cases. Further, in view of the relations $`P=B(1D)`$ and $`Q=D(1B)`$ one can verify that the “biological” domain in the $`(P,Q)`$-plane occupies only the kite-shaped region below the line $`\sqrt{P}+\sqrt{Q}=1`$. On this line we have $`B=1D`$. We will discuss later that this is a symmetry line in the $`(B,D)`$-plane. Therefore, the present biological application uses only a subset of the parameters of the physical deposition model. This is illustrated in Figure 2. The relevant kite-shaped region is filled in black. The white area below $`P+Q=1`$ is not accessible in the population model, but it is available in the deposition model.

## 4 Bacterial population distributions versus antibiotic efficiency

In order to illustrate the carrying capacity distribution and its properties we turn to a few concrete examples. Consider a process with $`B=0.5`$, representing equally probable nutrient-rich and nutrient-poor patches, and $`D=0.5`$, corresponding to medium antibiotic strength. It follows that $`P=0.25`$ and $`Q=0.25`$. Figure 3a shows the population distribution after one generation, $`\mathrm{\Pi }_1(𝒩)`$, for $`\lambda =3`$. Since in the first generation the population is uniquely determined by the growth sequence $`k_1`$, through $`𝒩(1)=N_0+(\lambda k_1)/\lambda ^2`$, with $`k_1=0,...,2\lambda `$, we have

$$\mathrm{\Pi }_1(𝒩)=\pi _1(k(𝒩))$$
(4.1)
The probabilities associated with the seven peaks in the figure thus follow directly from (III.4). Note that the distribution is symmetric about the initial population $`N_0=0.5`$, due to the fact that $`B=D`$.

In the language of fractal geometry  Fig.3a can be called the “generator” of the distribution. If this generator were now applied on a smaller scale to split every peak into 7 new peaks, we would obtain a perfect self-similar object, when iterated ad infinitum. However, in the second generation the generator is different and has $`2\lambda ^2+1`$ peaks, 19 in our example. Its width is reduced by a factor of $`\lambda `$, as can be seen from inspection of (III.1). It is the convolution of the two generators, according to the product expression (III.6) and the sum (III.7), which gives the population distribution after two generations. The result is shown in Fig.3b.

In this distribution the 7 probability peaks of the first generation are still clearly visible, but they have been split due to the relatively small population shifts obtained in the second generation. The 7 structured peaks thus have a width and fine structure determined by the generator of the second generation, and all 7 are similar to each other, differing only in their amplitude which is apparent in Fig.3a. Note that the total number of peaks, 73, is significantly smaller than the product $`7\times 19`$. This is due to overlap of peaks, which results, as discussed previously, from different growth sequences leading to the same population.

As the generations progress, the individual probability peaks can no longer be resolved on the normal population scale. This is the case already after three generations, as Fig.3c shows. There are 703 peaks after three, and 6481 peaks after four generations. The global scaling properties of the distribution can be easily understood if one takes into account that the number of peaks increases by approximately a factor of 9, since the inter-peak distance, or minimum population shift, decreases by a factor of $`\lambda ^2`$ in each generation, as implied by (III.1). Consequently, in order to preserve the normalization of the total probability the height of the peaks must decrease roughly by a factor of 9 also.

After four generations the population distribution looks nearly identical to that after three generations, provided the height is scaled by roughly a factor of 9. Therefore, Fig.3c (three generations) already reveals the shape of the population distribution in the limit of a large number of generations, i.e., it gives the carrying capacity distribution $`\mathrm{\Pi }(K)`$ that we are interested in.

Besides the analytic calculation, we have carried out simulations of these random processes and obtained population distributions in the form of histograms. Obviously, for a histogram to present accurately a discrete distribution of peaks, the bin width must be smaller than the inter-peak separation in population. For example, from these simulations, typically involving $`10^4`$ processes for given $`B`$ and $`D`$, a histogram with 1000 bins on $`(0,1)`$ is obtained which reproduces the analytic results accurately. In view of the fact that little or no changes can be perceived in the shape of the (properly scaled) distribution for generation numbers larger than three, it follows that a histogram with fixed bin size (1000) will converge (without scaling) already after three or four generations and reveal the carrying capacity distribution. In actual experiments a bin size is defined by practical resolution limitations on the counting of bacteria. Therefore, the biologically relevant frame-work is that which employs a fixed bin size.

The example we have chosen in Figs.3a-c is a typical case in the category of broad carrying capacity distributions with an average value close to (in this case at) the initial homogeneous population $`N_0`$. For medium nutrient levels and medium antibiotic strength the fluctuations are apparently too large to gain “control” over the outcome of the experiment. The situation changes drastically when bacteria in a nutrient-poor environment are exposed to strong antibiotics. In this case a thorough elimination of bacteria results almost certainly. This is shown in Fig.4.

For birth and death probabilities $`B=0.05`$ and $`D=0.95`$, Fig.4 gives the carrying capacity distribution after three generations. This corresponds to $`P=0.0025`$ and $`Q=0.9025`$. Clearly, the population has been suppressed to far below its initial value $`N_0=0.5`$, and the probability is concentrated on the lowest accessible population values, around $`𝒩0.06`$. This result is to be expected, and in line with common experience with the effect of antibiotics on bacteria in a nutrient-poor medium.

Possibly a more surprising result from this model is what happens when medium nutrient level is combined with strong (or very strong) antibiotics. If we take $`B=0.5`$ and $`D=0.99`$, leading to $`P=0.005`$ and $`Q=0.495`$, we obtain the carrying capacity distribution shown in Fig.5. The following properties are conspicuous: i) the population is entirely shifted to well below its initial value $`N_0=0.5`$, ii) the distribution is broad; iii) the average population is about half the initial population. It appears that the antibiotics has failed in two respects. The population has not been suppressed completely, and it is difficult to predict the eventual asymptotic population due to the large fluctuations. Mathematically, this insufficiency can be traced to be caused, not by the medium abundance of nutrient-rich patches in itself, since only very little growth is observed ($`P=0.005`$ is very small), but by a substantial reduction of the antibiotic efficiency. This reduction, from $`D=0.99`$ in biological parameter space, to $`Q=0.495`$ in physical parameter space, is due to the indirect effect of nutrient-rich spots on antibiotic strength through $`Q=D(1B)`$. Indeed, on $`50\%`$ of the biofilm the antibiotic efficiency is spent on taking out only the newly produced bacteria, since $`B=0.5`$, and only on the remaining $`50\%`$ of the area the initially existing population gets (almost) killed eventually, so that we can understand that $`\overline{N}(n)0.255`$, for $`n\mathrm{}`$.

The few representative cases we have discussed up to now are the most interesting ones which can be studied within this very simple model. We do not need to devote special attention to the limits $`BD0`$ or, equivalently, $`BD1`$, for which the distribution remains sharply concentrated around the initial value, so that the outcome of experiments is simple to predict. In the next section we elaborate on the behaviour of the most probable carrying capacity and its probability as a function of the birth and death rates, i.e., of $`B`$ and $`D`$.

## 5 Most probable carrying capacity and its probability

With regard to an actual experiment it is important to predict the most probable carrying capacity, $`K^{}`$, that will be observed, and its probability of occurrence, $`\mathrm{\Pi }(K^{})`$, given the birth and death probabilities $`B`$ and $`D`$. In general $`K^{}`$ differs from the mean $`\overline{K}`$. We have studied these properties of the model for rescaling factors $`\lambda =2`$ and $`\lambda =3`$. Differences in rescaling factor are seen to lead to quantitative differences only, not essential for our discussion. Here we report on the results for $`\lambda =3`$.

Since the population distributions converge rapidly with increasing generation number, the results for the first generation (the “generator”) determine to a very good approximation the relative importance of the peaks in the final distributions, as well as the population values associated with the peaks. Therefore, we can restrict our attention to the first generation when studying the most probable carrying capacity and its probability qualitatively, and work with the most probable population $`𝒩^{}`$ and its probability $`\mathrm{\Pi }_1(𝒩^)`$. In Appendix A we give the analytic expressions for the probabilities of the populations after one generation.

There are two lines of symmetry in the model, which merit separate attention. The line $`B=D`$ is a symmetry line in the sense that the population distributions for $`B>D`$ are the mirror images, reflected about $`N=0.5`$, of the distributions with $`B`$ and $`D`$ interchanged. On the line $`B=D`$ the average population remains equal to the initial population. There is no net growth (see Section II). One would therefore suspect that, when the antibiotic strength precisely compensates the nutrient abundance, the most probable population after the experiment equals the initial population. This is indeed the case. However, unless $`B`$ is sufficiently close to zero or to 1, the fluctuations are large and the probability $`\mathrm{\Pi }_1(0.5)`$ of observing the most probable population is less than $`50\%`$. This is illustrated in Fig.6 (case $`B=D`$).

The second symmetry concerns the substitution of $`B`$ by $`1D`$, and of $`D`$ by $`1B`$. Clearly, this operation leaves $`P`$ and $`Q`$ unchanged, and therefore the results are invariant. The line $`B+D=1`$ thus acts like a mirror in the $`(B,D)`$-plane. The probability $`\mathrm{\Pi }_1(𝒩^{})`$ along this line shows interesting corner-like singularities where the most probable population makes a jump. In a statistical mechanical context what we are doing is similar to minimizing the “free energy” $`F=\mathrm{\Pi }_1(𝒩)`$ with respect to $`𝒩`$, for every $`B`$, and obtaining the minimum-free-energy (maximum-$`\mathrm{\Pi }_1`$) curve which displays corners at points where two “phases” (values of $`𝒩^{}`$) coexist. Here, the relevant $`𝒩`$ are restricted to the seven population values generated in the first generation (Fig.3a), so that the curve can have (at most) six corners. Fig.6 shows the probability of the most probable population along this line (case $`B+D=1`$). Although this curve has been calculated only for the first generation, it is important to keep in mind that it is already a good approximation to the probability of the most probable carrying capacity $`K^{}`$, defined in the limit of a large number of generations. That curve, $`\mathrm{\Pi }(K^{})`$ must also displays corner-like singularities, representing jumps in $`K^{}`$, at values of $`B`$ (or $`D`$) close to the singularities exhibited in Fig.6. These jumps in $`K^{}`$ occur when two main peaks in the distribution (see Fig.3c for example) exchange maximum height.

A global plot of the most probable population $`𝒩^{}`$ after one generation, in the domain $`0B,D1`$ is presented in Fig.7. The two symmetries we discussed can be seen in this plot. The singularities, in the form of terrace border lines, represent the “coexistence” of two distinct values of $`K^{}`$. These lines run roughly parallel to lines of constant $`BD`$, so that, to a crude first approximation, the most probable population depends mainly on the difference $`BD`$. Recall that the average population, $`\overline{N}`$, depends only on $`BD`$, according to (II.1).

The probability $`\mathrm{\Pi }_1(𝒩^{})`$ as a function of $`B`$ and $`D`$ is shown in Fig.8. It displays lines of corner-like singularities precisely where the terrace borders occur in Fig.7. Note that a sharp population distribution is only found near the four points where $`B`$ and $`D`$ are close to 0 or 1.

In order to quantify this further we have calculated the standard deviation $`\mathrm{\Sigma }_1`$ of the distribution after one generation, as a function of $`B`$ and $`D`$. This quantity is defined through

$$(\mathrm{\Sigma }_1)^2=\frac{1}{2\lambda +1}\underset{k=0}{\overset{2\lambda }{}}(𝒩_k\overline{N}(1))^2\widehat{𝒫}_1(k),$$
(5.1)
with $`𝒩_k1/(\lambda 1)+(\lambda k)/\lambda ^2`$ and, from (II.2), $`\overline{N}(1)=1/(\lambda 1)+(BD)/\lambda `$. Explicit expressions for the $`\widehat{𝒫}_1`$ are given in the Appendix. The standard deviation is shown in Fig.9. Like the average population, the standard deviation is a smooth function of $`B`$ and $`D`$, and illustrates clearly that a broad population distribution is generic.

## 6 Conclusions and a possible experimental test

Keeping in mind the possible relevance of this hierarchical model to towering pillar biofilms, as discussed in the Introduction, we recapitulate here the main model ingredients. The model contains three parameters, birth probability $`B`$, death probability $`D`$ and rescaling factor $`\lambda `$. $`B`$ essentially reflects the fraction (or quality) of nutrient-rich area on the surface, and $`D`$ the density (or quality) of bactericidal agents. The rescaling factor $`\lambda `$ sets the characteristic time scale (Section II) and controls the fragmentation - in space - of the nutrient-rich patches and the antimicrobial activity, mimicking diffusion and transport. Further, the same $`\lambda `$ is used to control the saturation of the population by decreasing uniformly the amount of offspring in each generation and simultaneously increasing uniformly antimicrobial resistance. This mimicks, for example, a gradual temperature drop and biofilm development, respectively, as discussed in Section I. Through the action of $`\lambda `$ the model becomes hierarchical and this leads to a finite carrying capacity (“freezing” or stagnation). Undoubtedly, the use of just one parameter $`\lambda `$ for inducing all these dynamical effects, in the interest of simplicity and transparency of the model, cannot be more than a first crude step towards a more realistic and refined approach.

The growth of the population thus acquires a hierarchical structure, apparent in the fractal population “landscapes” (Fig.1), and expressed by a (nearly) self-similar carrying capacity distributions (Figs.3-5), the mathematics of which has been discussed using elementary notions of fractal geometry. The main conclusion from studying these distributions is that for generic values of $`B`$ and $`D`$ a broad range of carrying capacities can be observed, and the outcome depends largely on the sample used. Ensemble averaging is necessary for predicting statistically relevant properties of this type of growth.

The question now arises whether this model is relevant to an experiment in which bacteria capable of forming a biofilm on an inhomogeneous nutrient field are exposed to antibiotic or protein spray and subsequently put in a refrigerator or an oven (for about a day, with $`\mathrm{\Delta }T1^{}`$C/hr). The applicability of the model relies primarily on the determination of the model parameters starting from experimental system parameters. In our discussion we have taken $`\lambda =3`$, having in mind that the linear size of a nutrient-rich patch or an antibiotic drop is roughly one third of the inoculation line. This length ratio can be determined easily in practice and the rescaling factor can be adjusted, or, of course, two separate rescaling factors can be introduced if necessary. The external cooling (or heating) rate, taking the bacteria away from optimal growth conditions, is assumed to be adjusted so that the processes come to a halt in about 4 to 5 generations. The quality of the nutrient and the strength of the antibiotic, incorporated in $`B`$ and $`D`$, can be tuned experimentally by dilution, for example.

State-of-the-art experiments allow a simultaneous counting of many samples (say, $`10^4`$) on a micro-array in order to obtain carrying capacity distributions, much like is done in ensemble averaging. In sum, a direct experimental check of this model is feasible, and would be worthwhile, before introducing theoretical refinements which undoubtedly are necessary to make the model more realistic, but on the other hand compromise the insight that can be gained using only very few key parameters.

Acknowledgments.
We thank Ralf Blossey for suggesting the possible applicability of the hierarchical population model to biofilms. We are grateful to Reiner Kree and Joachim Krug for constructive comments on the model. This research is supported by the Flemish Programme FWO-G.0222.02 “Physical and interdisciplinary applications of novel fractal structures”. K. S.-W. is supported by the Foundation for Polish Science (FNP Scholarship).

## Appendix A Probabilities of the population after one generation: case $`\lambda =3`$

In this Appendix we give the analytic expressions used in the calculations of the most probable population and its probability after one generation (Section V). The rescaling factor is $`\lambda =3`$. In terms of the auxiliary probabilities

$$P=B(1D)\text{and}Q=D(1B)$$
(A.1)
with $`0B,D1`$, and in terms of $`S=P+Q`$, we obtain, using (III.4),

$`\widehat{𝒫}_1(0)`$ $`=`$ $`P^3`$
$`\widehat{𝒫}_1(1)`$ $`=`$ $`3P^2(1S)`$
$`\widehat{𝒫}_1(2)`$ $`=`$ $`3P(1S)^2+3P^2Q`$
$`\widehat{𝒫}_1(3)`$ $`=`$ $`(1S)^3+6P(1S)Q`$
$`\widehat{𝒫}_1(4)`$ $`=`$ $`3(1S)^2Q+3PQ^2`$
$`\widehat{𝒫}_1(5)`$ $`=`$ $`3(1S)Q^2`$
$`\widehat{𝒫}_1(6)`$ $`=`$ $`Q^3`$ (A.2)