text:
# Theory of Diffusion Controlled Growth

## I Introduction

Diffusion controlled growth first attracted attention in the literature on solidification, where the advance of a solidification front can be limited by diffusion of either latent heat or compositional excess ahead of the front. Under these conditions a planar front is linearly unstable with respect to long wavelength corrugation, the Mullins-Sekerka instability Mullins and Sekerka (1963), leading to a rich variety of problems in pattern formation. Viscous fingering, arising when a viscous fluid is driven through a porous medium by less viscous one, is recognised as being a problem in the same class. The Diffusion Limited Aggregation model (DLA) of a rigid cluster growing by the irreversible accretion of dilute diffusing particles, introduced by Witten and Sander Witten and Sander (1981) focussed attention on the extreme limit of these problems, where all of the diffusion is ahead of the growth and quasi-static, with the added simplification that the Mullins-Sekerka instability applies on all lengthscales above the size of the accreting particles.

Mathematically these problems share the same general form for the equations governing their growth, with their local interfacial velocity controlled by a conserved gradient flux,

$$v_n\left|_n\varphi \right|^\eta ,^2\varphi =0,\text{ }\varphi _{\text{interface}}0$$
(1)
where \[at least naively\] $`\eta =1`$ Ball (1986a). The generalisation to a range of positive $`\eta `$ was introduced by Niemeyer, Pietronero and Wiesmann Niemeyer et al. (1984) to model dielectric breakdown patterns, and in our recent letter Ball and Somfai (2002) we introduced the idea that this can support equivalences between models where the Mullins-Sekerka instability is controlled locally in mathematically quite different ways.

The DLA model has attracted enormous attention because it contains no limiting lengthscale (save the particle size) and so pattern formation must continue non-trivially on all larger lengthscales, the Mullins-Sekerka instability ruling out simple planar growth. Theoretical interest has been fueled by the fractal and multifractal Amitrano et al. (1991); Halsey et al. (1986) scaling properties of the clusters produced, with controversial claims Plischke and Rácz (1984); Coniglio and Zannetti (1989); Mandelbrot et al. (2002) (and counter-claims Meakin and Sander (1985); Somfai et al. (1999); Ball et al. (2002)) of anomalous scaling, and by the long-standing absence of an overall theoretical framework to understand the problem. A simple mean field theory Ball et al. (1984) does not capture the fractal aspects, which are better understood through various relations between exponents Ball and Witten (1984); Turkevich and Scher (1985); Ball (1986b) and Halsey (1987). The Cone Angle Theory Ball (1986b) gives a plausible argument for the fractal dimension of DLA, whilst the Screened Growth Model Halsey and Leibig (1992) and Makarov’s Theorem Makarov (1985) give insight into how the multifractal spectrum of the growth is generated.

The presence of a cut-off lengthscale $`a`$ below which the physics dictates smooth growth is a crucial ingredient of diffusion controlled growth; it is known that otherwise infinitely sharp cusps develop in the interface within finite time Shraiman and Bensimon (1984). In DLA this cutoff is fixed and set by the size of accreting particles, but in solidification and also in viscous fingering it is set by surface free energy leading to a local offset in the value of $`\varphi _{\text{interface}}`$: this is dependent on the local interfacial curvature and leads to a velocity-dependent selection of approximately constant $`va^2`$ for growing tips Langer (1980). We will gather these different possibilities together with the generalised cut-off law, that advancing tips have radius

$$a\left|_n\varphi \right|^m\text{.}$$
(2)
In terms of $`m`$, simple DLA corresponds to $`m=0`$ and solidification and viscous fingering correspond to $`m=1/2`$; in the theory below in two dimensions we will map onto the case where $`a`$ is such that each growing tip has fixed integrated flux, corresponding to $`m=1`$.

This paper explains and expands the new theory announced in our recent letter Ball and Somfai (2002). In Sec. II we establish mappings between models with different $`\eta `$ and $`m`$. This opens up the opportunity to discuss the full class of models via the section at $`m=1`$, which we will show in Sec. IV is particularly amenable to continuum theoretical description in two dimensions. Section III discusses the other key input to Sec. IV, that we regard diffusion controlled growth as a turbulence problem with self-organising fluctuations. In Sec. V we discuss how to cast the theory of Sec. IV in renormalised form, with divergent factors factored out consistently. Closure approximations are required to obtain explicit theoretical predictions, and in Sec. VI we show how the simplest Gaussian closure leads to a complete theory of the fractal and multifractal scaling. This turns out to be quantitatively quite accurate for the zone of active growth. In Sec. VII we show how one outstanding exponent, the tip scaling of the harmonic measure, can be pinned down through the use of the Electrostatic Scaling Law leading to a very surprising prediction with which numerical data seems compatible. Section VIII gives a more detailed discussion of what happens to the DBM exponents at large $`\eta `$. In Sec. IX we show examples of how the theory can be deployed to tackle deeper quantities such as the relative penetration depth in DLA, which has been the subject of several numerical studies and some controversy. Angular resonances appear leading us to draw parallels in Sec. XI with the earlier Cone Angle Theory.

## II Scaling Properties and Mappings between models at different $`\eta `$ and $`m`$

It is central to fractal (and multifractal) behaviour in DLA that the measure given by the diffusion flux \[density\] $`j_n\varphi `$ onto the interface has singularities Halsey et al. (1986), such that the integrated flux onto the growth within distance $`r`$ of a singular point is given by

$$\mu (r)(r/R)^\alpha ,$$
(3)
where $`R`$ is the overall linear size of the growth. Multifractal scaling of the flux density leads to a whole spectrum of $`\alpha `$ values, with the number of regions of lengthscale $`r`$ having $`\mu (r)(r/R)^\alpha `$ varying as $`(r/R)^{f(\alpha )}`$, but in the following we focus particularly on advancing tips and their associated exponent value $`\alpha _{\text{tip}}.`$

Applying this phenomenology to the scaling around growing tips, we can establish an equivalence between models at different $`\eta `$ and $`m`$ by requiring that the relative advance rates of different growing tips are matched. Consider two growths, growing governed by different parameters $`(\eta ,m)`$ and $`(\eta ^{},m^{})`$ respectively, which at a given moment have the same overall geometry down locally to the level of \[the coarser of\] their cut-off lengthscales. For any given growing tip (labelled $`k`$), the tip radius $`a_k`$ and flux density $`j_k`$ in the unprimed growth will be related to those in the primed growth by $`j_ka_k^{d1}/a_k^\alpha =`$ $`j_k^{}a_k^{d1}/a_k^\alpha `$, where $`\alpha `$ is the local scaling exponent (as per Eq. (3)) of the harmonic measure between lengthscales $`a_k`$ and $`a_k^{}`$. We take this exponent to have value $`\alpha =\alpha _{\text{tip}}`$ on the grounds that this is locally a tip of the growth. Now let us focus on two different growing tips labelled $`1,2`$ whose radii and flux densities are interrelated in the unprimed growth according to Eq. (2) by $`a_1j_1^m=a_2j_2^m`$, and similarly in the primed growth by $`a_1^{}j_1^m^{}=a_2^{}j_2^m^{}`$. If we now insist that the advance velocities are in the same ratio (between tips $`1`$ and $`2`$) in both models, this requires $`(j_1/j_2)^\eta =(j_1^{}/j_2^{})^\eta ^{}`$, which forces the parameter relation

$$\frac{1+m(1+\alpha d)}{\eta }=\frac{1+m^{}(1+\alpha d)}{\eta ^{}}.$$
(4)
For the two models to be equivalent in the relative velocities of all tips requires their parameters be related as above, where $`\alpha =\alpha _{\text{tip}}`$ is the singularity exponent associated with growing tips.

Although we have not strictly proved the equivalence of the models related above, we have shown that any such relationship must follow Eq. (4) and we will assume in the rest of this paper that this equivalence holds. All such models are then classifiable in terms of a convenient reference such as $`\eta _0`$, the equivalent $`\eta `$ when $`m=0`$, corresponding to the original Dielectric Breakdown Model. For example dendritic solidification with $`\eta =1`$ and $`m=1/2`$ corresponds to $`\eta _0=\frac{2}{3+\alpha d}`$: it is thus not equivalent to DLA, but to another member of the DBM class.

Another puzzle resolved by our classification is a recent study showing conflicting scaling between DLA and different limits of a “Laplacian growth” model Barra et al. (2001). The latter model grows bumps of width proportional to flux density, so in the present terminology it corresponds to $`m=1`$. The bumps are also grown with protrusion proportional to flux density. When the coverage of the growing surface (per time step of growth) is low, then as the bumps are also distributed proportional to flux density, this limit corresponds to $`\eta =3`$. By contrast high coverage (with significant suppression of overlapping bumps) corresponds to $`\eta =1`$. Using $`\alpha =0.7`$ (see below) these map through Eq. (4) into $`\eta _0=2.31`$ and $`\eta _0=0.77`$ respectively, so the way their scaling was observed Barra et al. (2001) to bracket that of DLA is quite expected.

## III The Role of Noise

DLA and DBM have been widely regarded as models in statistical physics, in that the local advance rate in Eq. (1) has been implemented as the probability per unit time for the growth locally to make some unit of advance, entailing an inherent shot noise. Here we argue that diffusion controlled growth is a problem of turbulence type, with noise self-organising from minimal input. This was first suggested by Sander et al. Sander et al. (1985) but was only pursued in the case with surface tension cut-off, $`m=1/2`$ in the present terminology, where it has been recognised more recently as chaotic viscous fingering Kessler and Levine (2001).

The renormalisation of noise with lengthscale has hitherto been discussed (at least for DLA) in the context of noise reduction Barker and Ball (1990); Cafiero et al. (1993); Ball et al. (2002), focussing on the idea that as one goes up in lengthscale an equivalent coarse-scale model must have lower level of noise than crude shot noise. The limiting or “fixed-point” level of noise in DLA is small (at least according to references Barker and Ball (1990) and Ball et al. (2002)) but certainly not zero, so it is natural to ask whether it can be approached from below as well as from above. The data in Fig. 3 show clearly that for DLA grown with very low noise by the methods of reference Ball et al. (2002), the relative fluctuations do indeed approach their limiting value from below as well as from above, and the same result was implicit in the earlier Renormalisation Group results of Barker and Ball (1990).

The above discussion leads us to conjecture that for the full range of diffusion controlled growth under a continuum description of the interface, disorder in the initial conditions alone should suffice to feed instability, leading to the same limiting levels of structural fluctuation on larger lengthscales as in the discrete models. The agreement we obtain below from simulations of the continuum model without temporal noise provides direct evidence for this idea, which might also be argued obvious on the grounds that the Mullins-Sekerka instability Mullins and Sekerka (1963) corresponds to a preponderance of positive Lyaponov exponents in the dynamics.

## IV Continuum Theoretical Description and Self-Organisation

The new ideas above, that we can balance changing the cut-off exponent $`m`$ by adjustment of $`\eta `$, and that noise can be left to self organise, are the key to a new theoretical formulation of the problem, at least in two dimensions of space to which we now specialise. In two dimensions the Laplace equation in Eq. (1) can be solved in terms of a conformal transformation between the physical plane of $`z=x+iy`$ and the plane of complex potential $`\omega =\varphi +i\theta ,`$ in which we take the growing interface to be mapped into the periodic interval $`\theta =[0,2\pi ),\varphi =0`$ and the region outside of the growth mapped onto $`\varphi >0`$. Then adapting reference Shraiman and Bensimon (1984), Eq. (1) leads to a closed equation for the dynamics of just the interface $`z(\theta )`$ at $`\varphi =0`$,

$$\frac{z(\theta )}{t}=i\frac{z}{\theta }𝒫\left[\left|\frac{\theta }{z}\right|^{1+\eta }\right].$$
(5)
The linear operator $`𝒫`$ is most simply described in terms of Fourier transforms: $`𝒫\left[_ke^{ik\theta }f_k\right]=_kP(k)e^{ik\theta }f_k=f_0+2_{k=1}^Ke^{ik\theta }f_k`$, where we have introduced here an upper cut-off wavevector $`K`$. It is easily shown that on scales of $`\theta `$ greater than $`K^1`$ a smooth interface is linearly unstable with respect to corrugation for $`\eta >0`$ (the Mullins-Sekerka instability Mullins and Sekerka (1963)), whereas for scales of $`\theta `$ less than $`K^1`$ the equation drives smooth behaviour (corresponding locally to the case $`\eta =1`$). This cutoff on a scale of $`\theta `$, the cumulative integral of flux, corresponds in terms of tip radii and flux densities to $`ajK^1`$, that is an $`m=1`$ cutoff law. Thus the parameter $`\eta `$ in Eq. (5) is more specifically $`\eta _1=\alpha \eta _0`$, using Eq. (4) with $`d=2`$.

We will present a numerical study of Eq. (5) after a variable change in Sec. VII, where disorder was supplied only through the initial conditions. The results clearly confirm that the equation self-organises critical scaling behaviour, without the supply of time dependent noise. The surprising form of the scaling is interpreted below.

## V Renormalisation and the DBM scaling law

We now turn to a theoretical analysis of Eq. (5), and for generality we will consider growth in a wedge of angle $`2\pi c`$ (with periodic angular boundary conditions) so that $`c0`$ corresponds to growth along a channel whilst $`c=1`$ corresponds “radial growth”, that is growth out from a point in the plane. The primary theoretical requirement is that we must obtain results explicitly independent of the cut-off as $`K\mathrm{}`$ : this is hard because we will see that the mean advance rate of the interface diverges as a power of $`K`$, and on fractal scaling grounds one would expect the same divergent factor to appear in the rate of change of other simple variables. One can of course take ratios of rates of change and look to order terms such that divergences cancel. To make this work we have been forced to introduce yet another change of variables,

$`i{\displaystyle \frac{z}{\theta }}`$ $`=`$ $`cR\mathrm{exp}\left(ic\theta \lambda (\theta )\right)`$ (6)
$`=`$ $`cR\mathrm{exp}\left(ic\theta {\displaystyle \underset{k>0}{}}\lambda _ke^{ik\theta }\right),`$
which corresponds to Fourier decomposing the logarithm of the flux density. Here $`R`$ is the effective radius of the growth and the non-analytic factor $`e^{ic\theta }`$ gives the mean winding of the conformal map through the wedge angle $`2\pi c`$, leaving $`\lambda (\theta )`$ as a simple Fourier series (except one-sided; see details in Appendix A). The key to the success of the “logarithmic variables” $`\lambda `$ is that they decompose the flux density itself multiplicatively and, as we shall see, quite naturally capture its multifractal behaviour. In terms of these, time rescaled through $`dt=(cR)^{2y}d\widehat{t}`$ and $`y=(1+\eta )/2`$, the equation of motion (5) becomes

$$\begin{array}{c}\frac{d\lambda _k}{d\widehat{t}}=\underset{j<k}{}(kj)\lambda _{kj}P(j)\left(e^{y(\lambda +\overline{\lambda })}\right)_j\hfill \\ \hfill +2(kc)\left(e^{y(\lambda +\overline{\lambda })}\right)_k\end{array}$$
(7)
where subscripts on bracketed expressions imply the taking of a Fourier component, by analogy with $`\lambda _k`$. The advance rate of the mean interface is correspondingly given by

$$\frac{dR}{d\widehat{t}}=cR\left(e^{y(\lambda +\overline{\lambda })}\right)_0.$$
(8)
Details of the above analysis are given in Appendix A.

At this point we can evaluate the multifractal spectrum in terms of these logarithmic variables. The multifractal spectrum of the harmonic measure follows from computing the general moment Halsey et al. (1986) $`Z(q,\tau )=\left|\delta \theta \right|^q\left|\delta z\right|^\tau `$, in the limit where all the intervals $`\left|\delta z\right|`$ and correspondingly $`\delta \theta `$ approach zero; then the locus $`(q1)D(q)=\tau (q)`$ separates the limiting behaviour $`Z(q,\tau )\mathrm{}`$ from $`Z(q,\tau )0`$. In our case it is convenient to fix $`\delta \theta `$ (admitting wide variations in $`\left|\delta z\right|`$) and we must focus on the restricted range $`K^1<\delta \theta <1`$. For $`\delta \theta 1`$ we have trivially $`Z(cR)^\tau ,`$ whilst for $`\delta \theta K^1`$ the growth begins to look smooth so we can approximate $`Z(q,\tau )\left|\delta \theta \right|^{q\tau }\left|\frac{z}{\theta }\right|^\tau K^{q+\tau +1}𝑑\theta \left|\frac{z}{\theta }\right|^\tau `$. Averaging gives $`\left|\frac{z}{\theta }\right|^\tau =(cR)^\tau e^{(\lambda +\overline{\lambda })\tau /2}`$, yielding

$$Z(q,\tau )K^{q+\tau +1}(cR)^\tau e^{(\lambda +\overline{\lambda })\tau /2}$$
(9)
The multifractal spectrum is then readily obtained from the separator behaviour discussed above, provided we can evaluate the average in Eq. (9). In Sec. VI we show how this can be done quite explicitly in the Gaussian Closure Approximation.

We can now derive an elegant combination of Halsey’s Electrostatic Scaling Law Halsey (1987) combined with the tip scaling law of Ball and Witten Ball and Witten (1984); Turkevich and Scher (1985); Ball (1986b), that for DLA $`\tau (3)=D_f=1+\alpha _{\text{tip}}`$. This new derivation (unlike the earlier results) is not restricted to the $`m=0`$ case. The key idea is that we match the advance rate of the forward tips of the growth (governed by $`\alpha _{\text{tip}}`$) to that of the mean radius governed by a multifractal moment through Eq. (8). For consistency with the rest of this paper it is convenient to present the argument for the $`m=1`$ representation, leading to tip velocity $`\frac{dR}{dt}=j^{\eta _1}\left(K^1/a\right)^{\eta _1}`$ where the tip radius $`a`$ is set by the condition $`K^1aj(a/R)^{\alpha _{\text{tip}}}`$ and hence $`\frac{dR}{dt}R^{\eta _1}K^{\eta _1(1/\alpha _{\text{tip}}1)}`$. The overall advance rate of the growth (in terms of its effective radius) is given from Eq. (8) by $`\frac{dR}{dt}(cR)^{\eta _1}e^{y(\lambda +\overline{\lambda })}`$. For the special value $`\tau /2=y`$ this can be substituted into Eq. (9), giving $`\frac{dR}{dt}R^{\eta _1}K^{q(\tau )\tau 1}`$, when $`Z`$ is kept at its $`K`$-independent value $`(cR)^\tau `$. Comparing the two results leads to the DBM scaling law

$$q(\tau =1+\eta _1)=2+\eta _1/\alpha _{\text{tip}}=2+\eta _0,$$
(10)
or its inverse form

$$\tau (q=2+\eta _0)=1+\eta _1=1+\alpha _{\text{tip}}\eta _0.$$
(11)
The same result can be found, much more tortuously, from growth at general $`m`$.

Now we turn back to the equation of motion of the logarithmic variables. Let us suppose some ignorance of the initial conditions and describe the system in terms of a joint probability distribution over the $`\lambda _k`$, and let us denote averages over this \[unknown\] distribution by $`..`$. We can in principle determine the distribution through its moments, whose evolution we now compute. For simplicity in this paper we assume translational invariance with respect to $`\theta `$, so that only moments of zero total wavevector need be considered, of which the lowest gives:

$$\begin{array}{cc}\hfill \frac{d}{d\widehat{t}}\lambda _k\overline{\lambda }_k=& (\underset{j<k}{}(kj)P(j)\lambda _{kj}\overline{\lambda }_ke_j^{y(\lambda +\overline{\lambda })}\hfill \\ & +2(kc)\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })})+(\text{c. conj.}).\hfill \end{array}$$
(12)
All of the higher moments lead to the same form of averages on the RHS, $`\text{multinomial(}\lambda ,\overline{\lambda }\text{)}e^{y(\lambda +\overline{\lambda })}`$, and all of these terms are conveniently expressed in terms of cumulants Kubo (1962) as detailed in Appendix B. The key helpful feature is that the expressions we require all naturally divide by one factor of $`e^{y(\lambda +\overline{\lambda })}=\frac{1}{c}\frac{d}{d\widehat{t}}\mathrm{ln}R`$, which is what we need in order to remove divergences by eliminating increment of time $`d\widehat{t}`$ in favour of $`\frac{1}{c}d\mathrm{ln}R.`$ The latter quantifies the update of large scale geometry, in terms of the advance in radius relative to the circumference of the wedge (or width of channel). The evolution of the second moments is then given in renormalised form by

$$\begin{array}{cc}\hfill c\frac{d}{d\mathrm{ln}R}\lambda _k\overline{\lambda }_k=& (2(kc)\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })}_c\hfill \\ & \underset{j<k}{}(kj)P(j)[\lambda _{kj}\overline{\lambda }_ke_j^{y(\lambda +\overline{\lambda })}_c\hfill \\ & +\lambda _{kj}e_{jk}^{y(\lambda +\overline{\lambda })}_c\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })}_c])\hfill \\ & +\left(\text{c. conj.}\right)\hfill \end{array}$$
(13)
where $`_c`$ denotes the Kubo cumulant Kubo (1962).

The above result (13) is the key analytical step in this paper, because it removes divergent factors form the equations of motion. It is not dependent on the closure approximation discussed below, and should support other approaches also. Moreover Eq. (13), with the hierarchy of analogous equations for the evolution of higher moments, offers a new entry point towards exact results in the class of DLA and DBM models.

## VI Gaussian Closure Approximation

To obtain simple tractable results we need to introduce some closure approximation(s) and we present here the simplest, neglecting all cumulants higher than the second, equivalent to assuming a joint Gaussian distribution (of zero mean) for $`\lambda `$. This is entirely characterised by its second moments $`S(k)=\lambda _k\overline{\lambda }_k`$ which by Eq. (7) we find evolve according to

$$c\frac{dS(k)}{d\mathrm{ln}R}=2y^2S(k)\left(\frac{\eta _1}{y^2}(kk^{})kS(k)2\underset{j<k}{}jS(j)\right)$$
(14)
where $`k^{}=c(1+1/\eta _1)`$, and again the details are in Appendix B.

Equation (14) evolves to a unique steady state. The key to understanding this is to note that for $`k<k^{}`$ the whole factor in large braces is negative definite, so $`S(k)=0`$ is the unique attractor. Then for $`k=k_{\text{min}}`$, the first integer value above $`k^{}`$, $`S(k_{\text{min}})=0`$ is unstable and the zero of the last factor leads to the global attractor having $`k_{\text{min}}S(k_{\text{min}})=`$ $`\frac{\eta _1}{y^2}(k_{\text{min}}k^{})`$. The attractor values for higher $`k`$ now follow by induction: denote the factor in large braces by $`B(k)`$ and assume that the attractor has $`B(k)=0`$ and $`0kS(k)`$ $`\frac{\eta _1}{y^2}`$, which are true for $`k=k_{\text{min}}`$. Then it follows that $`B(k+1)=\frac{\eta _1}{y^2}kS(k)(k+1)S(k+1)`$ and for $`kS(k)<\frac{\eta _1}{y^2}`$ the attractor must in turn have $`B(k+1)=0`$ and hence a value of $`S(k+1)`$ conforming to $`0(k+1)S(k+1)`$ $`\frac{\eta _1}{y^2}`$. In the case $`kS(k)=\frac{\eta _1}{y^2}`$ the only and stable solution is $`S(k+1)=0`$ which leads to the same conclusions. Thus by induction the only stable attractor of the system has $`B(k)=0`$ for all $`k>k^{}`$ and the corresponding steady state values are

$`kS(k)`$ $`=0,`$ $`k`$ $`<k_{\text{min}}`$
$`kS(k)`$ $`={\displaystyle \frac{\eta _1}{y^2}}(k_{\text{min}}k^{}),`$ $`k`$ $`=k_{\text{min}},k_{\text{min}}+2,\mathrm{}`$
$`kS(k)`$ $`={\displaystyle \frac{\eta _1}{y^2}}(1+k^{}k_{\text{min}}),`$ $`k`$ $`=k_{\text{min}}+1,k_{\text{min}}+3,\mathrm{}`$ (15)
Note that when $`k^{}`$ is integer, as particularly in the case $`c=0`$ corresponding to growth in a channel as discussed in Ball and Somfai (2002), $`k_{\text{min}}=1+k^{}`$ and alternate values of $`kS(k)`$ are zero: in the case of channel growth this absence of even $`k`$ is readily interpreted in terms of the dominance of one major finger and one major fjord.

Within the Gaussian approximation and its predicted variances (15) we can now compute all \[static\] properties of diffusion controlled growth. From Eq. (9) we obtain $`Z(q,\tau )(cR)^\tau \mathrm{exp}\left(\tau ^2/4_k^KS(k)\right)R^\tau K^{\frac{\tau ^2\eta _1}{8y^2}}`$, using the values from Eq. (15), and hence $`Z(q,\tau )R^\tau K^{q+\tau +1+\frac{\tau ^2\eta _1}{8y^2}}`$. Thus the separator of limiting behaviour (now as $`K\mathrm{}`$) is given by

$$q(\tau )=1+\tau +\tau ^2\frac{\eta _1}{2\left(1+\eta _1\right)^2}$$
(16)
It is also easy to see that any closure scheme based on keeping cumulants of $`\lambda `$ up to some finite order leads to a corresponding degree polynomial truncation of $`q(\tau )`$.

From the Legendre Transform of the inverse function $`\tau (q)`$ (as detailed in Appendix C) we obtain the corresponding spectrum of singularities,

$$f(\alpha )=2\frac{1}{\alpha }+\frac{1}{2}\left(\eta _1+\frac{1}{\eta _1}\right)\left(2\alpha \frac{1}{\alpha }\right)$$
(17)
which in Fig. 4 is compared to measured data for DLA Ball and Spivack (1990), which later measurements Jensen et al. (2002) reinforce. For the region of active growth $`\alpha 1`$ ($`q0`$) the theory is quantitatively accurate. At $`\alpha =1`$ it conforms to Makarov’s theorem Makarov (1985), and in contrast to the Screened Growth Model Halsey and Leibig (1992) it does this without adjustment. For $`\alpha >1`$ the spectrum is only qualitatively the right shape, and for such screened regions our equations based on tip scaling may not hold.

## VII Scaling predictions for DBM

To compare with the conventional DBM at $`m=0`$ and parameterised by $`\eta _0`$ , we still need to compute theoretically the value of the tip scaling exponent which enters through $`\eta _1=\alpha _{\text{tip}}\eta _0`$.

We can use the DBM scaling law (10) with the Gaussian Closure Approximation (16) for $`q(\tau )`$ to fix the value of $`\alpha _{\text{tip}}`$, and the resulting prediction is $`\alpha _{\text{tip}}=2/3`$ *independent of* $`\eta _1`$. For DLA in two dimensions this value is respectably close to (but outside) measured values, $`\alpha _{\text{tip}}=D1=0.71\pm 0.01`$ known from large direct simulations of DLA Ossadnik (1991); Somfai et al. , but its suggested independence of $`\eta `$ over a range of DBM is quite shocking. Numerical evidence, however, appears to lend support.

We have investigated numerically what value of $`\alpha _{\text{tip}}`$ is seleced by the dynamics of Eq. (5), with disorder supplied only through the initial condition. Changing variables to $`\psi =\left(i\frac{z}{\theta }(cR)^1e^{ic\theta }\right)^{(1+\eta _1)/2}=e^{y\lambda }`$, we obtain

$$\frac{\psi }{\widehat{t}}=i\frac{\psi }{\theta }𝒫[\psi \overline{\psi }]+iy\psi \frac{}{\theta }𝒫[\psi \overline{\psi }]cy\psi \left(𝒫[\psi \overline{\psi }](\psi \overline{\psi })_0\right)$$
(18)
where the rescaled time $`\widehat{t}`$ is defined in terms of the evolution of cluster radius in Eq. (8). The tri-linear form of the RHS enables us to compute numerically the motion within a purely Fourier representation.

In Sec. V we have seen that the cutoff dependence of the tip velocity is $`vK^{\eta _1(1/\alpha _{\text{tip}}1)}`$. This can be compared with the growth rate of the effective radius, $`\frac{dR}{dt}e^{y(\lambda +\overline{\lambda })}=\psi \overline{\psi }=\psi \overline{\psi }_0`$. So $`\alpha _{\text{tip}}`$ can be obtained from the $`K`$-dependence of $`\psi \overline{\psi }_0`$ (measured at fixed $`R`$). We can obtain it even from simulations with single $`K`$: the truncated sum $`v_{\text{cum}}(k)=_{j<k}\left|\psi _j\right|^2`$ is expected to scale with $`k`$ in the same way as the full sum $`\psi \overline{\psi }_0`$ does with $`K`$, because the Fourier components far below the cutoff should be insensitive to the value of $`K`$. Figure 5 shows the measured variation of $`v_{\text{cum}}(k)=_{j<k}\left|\psi _j\right|^2`$ vs $`k^{\eta _1}`$: this is expected to exhibit a power law with exponent $`(1/\alpha 1)`$ and remarkably we obtain $`\alpha 0.74\pm 0.02`$ with no significant dependence on $`\eta _1`$ in the range studied.

It is a remarkable success for the Gaussian Theory to have predicted the completely unexpected insensitivity of $`\alpha _{\text{tip}}`$ to $`\eta `$ . Whether this result can be truly an exact “superuniversality” is another matter, as certainly the Gaussian value for $`\alpha _{\text{tip}}`$ is only approximate and, as we discuss below, matters become more complicated for large $`\eta `$.

## VIII Breakdown of the DBM model at large $`\eta `$.

Sanchez and Sander Sanchez et al. (1993) first noted that at high enough $`\eta _0`$ the DBM degenerates because all growth is dominated by the most active site, showing by direct simulations (at $`m=0`$) that this happened around $`\eta _0^c4`$, a value reinforced by later discussions Hastings (2001a); Halsey (2002) and new data Hastings (2001b). These discussions are particular to the $`m=0`$ representation and we believe can be associated with the degeneration of the moment governing the rate of gain of cluster mass: this scales with exponent $`\tau (\eta _0)`$ which degenerates to value $`\eta _0\alpha _{\text{min}}`$ when the moment becomes dominated by the (left) end point $`f(\alpha )=0`$ of the multifractal spectrum. The fractal dimension given by $`d_f=1+\eta _0\alpha _{\text{tip}}\tau (\eta _0)`$ then degenerates to $`d_f=1+\eta _0\left(\alpha _{\text{tip}}\alpha _{\text{min}}\right)`$. If the least screened sites are the tips, $`\alpha _{\text{tip}}=\alpha _{\text{min}}`$, then this also leads to $`d_f=1`$ when $`\eta =\eta _0^c.`$

The electrostatic scaling law leads to earlier transition in the behaviour, that is at lower $`\eta `$, which is also more generic in that it does not depend on growth at some particular value of $`m`$. This transition should also limit the applicability of calculating Hastings (2001a); Halsey (2002) exponents perturbatively about $`\eta _0^c`$. The screening transition arises because the moment governing the mean screening of sites has exponent $`\tau (2+\eta _0)`$ which duly appears in the DBM version Eq. (11) of the electrostatic scaling law, and this moment must hit the end of the $`f(\alpha )`$ spectrum before that corresponding to $`\tau (\eta _0)`$ discussed above. Once we have hit this regime, at $`\eta _0\eta _0^s`$, we have $`\tau (2+\eta _0)=(2+\eta _0)\alpha _{\text{min}}`$ and the electrostatic scaling law degenerates to a form which can be rearranged to give

$$2\frac{1}{\alpha _{\text{min}}}=\eta _1\left(\frac{1}{\alpha _{\text{min}}}\frac{1}{\alpha _{\text{tip}}}\right).$$
(19)
This then leads us to choose between two scenarios: either (i) $`\alpha _{\text{tip}}>\alpha _{\text{min}}`$ in which case the behaviour remains non-trivial, or else (ii) $`\alpha _{\text{min}}=1/2`$ which means the arms of the growth are essentially straight and we might suspect self-affine structure.

The Gaussian Closure Approximation, with $`\alpha _{\text{tip}}`$ set by the electrostatic scaling law, has $`\alpha _{\text{tip}}>\alpha _{\text{min}}`$ for almost all $`\eta ,`$ leading to scenario (ii) above. Figure 6 shows the predicted variation of $`\alpha _{\text{tip}}`$, $`\alpha _{\text{min}}`$ and the value $`\alpha _{\text{screening}}`$ corresponding to the exponent $`\tau (2+\eta _0).`$ The screening transition where $`\alpha _{\text{screening}}`$ hits $`\alpha _{\text{min}}`$ occurs (it can be checked exactly) at $`\eta _1^s=2`$ corresponding to $`\eta _0^s=3.`$ Beyond this point $`\alpha _{\text{tip }}`$is no longer quite constant and stays clear of $`\alpha _{\text{min}}`$, because of the changed functional form for $`\tau (2+\eta _0),`$ whilst of course $`\alpha _{\text{screening}}`$ follows $`\alpha _{\text{min}}.`$

The corresponding predicted behaviour of the fractal dimension calculated from $`d_f=1+\eta _0\alpha _{\text{tip}}\tau (\eta _0)`$ is shown in Fig. 7 (upper curve) as a function of $`\eta _0`$ up to $`\eta _0^c5.4,`$ where it has not fallen to unity because $`\alpha _{\text{tip}}>\alpha _{\text{min}}`$ is maintained. There is some change of functional form across $`\eta _0=\eta _0^s`$ but it is scarcely noticeable graphically. Also shown for comparison (lower curve) is the behaviour when we force $`\alpha _{\text{tip}}=\alpha _{\text{min}}`$ instead of obeying the electrostatic scaling law: in this case the fractal dimension does smoothly approach unity as $`\eta _0\eta _0^c,`$ but unfortunately having sacrificed the electrostatic scaling law we cannot see anything relating to the screening transition.

The predicted screening transition at $`\eta _0^s=3`$ seems to mark a break in the match to the simulation data of Hastings Hastings (2001b): below this conforming to the electrostatic scaling law gives the better agreement, whereas beyond this better agreement comes from forcing $`\alpha _{\text{tip}}=\alpha _{\text{min}}`$. The simplest interpretation is that the exact answer conforms to both conditions, and it is just their relative importance which changes around the screening transition.

## IX The Penetration Depth

The multifractal spectrum suggests that the Gaussian approximation is good in the growth zone, so we have computed as a further test the relative penetration depth $`\mathrm{\Xi }`$, defined for DLA as the standard deviation $`\xi `$ of radius of deposition divided by the effective radius $`R`$. For DBM more generally, we have for tractability used as measure the diffusion flux rather than the local growth rate.

The key idea behind the calculation is that we calculate the relative distortion of the conformal map of the interface $`z(\theta )=Re^{ic\theta }\left(1+_{k>0}w_ke^{ik\theta }\right)`$, away from the circular arc $`z_0(\theta )=Re^{ic\theta }`$. Following Ref. Somfai et al. (1999), the squared relative penetration depth is then given by

$$\frac{\xi ^2}{R^2}=_0^{2\pi }\frac{d\theta }{2\pi }\left(Re\frac{z(\theta )z_0(\theta )}{z_0(\theta )}\right)^2$$
(20)
which can be expressed in terms of the coefficients $`w_k`$ as

$$\frac{\xi ^2}{R^2}=\frac{1}{2}\underset{k>0}{}|w_k|^2.$$
(21)
From the definition (6) of $`\lambda (\theta )`$, the coefficients can be identified as $`w_k=\frac{c}{ck}\left(e^{\lambda (\theta )}\right)_k`$, giving

$$\begin{array}{cc}\hfill \frac{\xi ^2}{R^2}=& \frac{1}{2}\underset{k>0}{}\frac{c^2}{(ck)^2}_0^{2\pi }\frac{d\theta }{2\pi }e^{ik\theta }_0^{2\pi }\frac{d\varphi }{2\pi }\hfill \\ & \times e^{ik\varphi }\left(e^{_p\lambda _pe^{ip\theta }_p\overline{\lambda }_pe^{ip\varphi }}\right).\hfill \end{array}$$
(22)
This expression remains to be averaged over the distribution of cluster geometries. Appendix D details the averaging of this under the Gaussian Closure Approximation, leading to the results shown in Fig. 8.

For radial DLA the theory predicts $`\mathrm{\Xi }=(\xi /R)_{\text{r.m.s.}}=0.20`$ in rather modest agreement with $`0.13`$ extrapolated from simulations Somfai et al. (1999), and the prediction over a range of DBM parameter is shown in Fig.8. The predicted variation of $`\mathrm{\Xi }`$ for DLA grown in a wedge is also shown, and in the limit of zero wedge angle the value for the penetration depth relative to the width of the channel is theoretically $`0.13`$ compared with $`0.14`$ measured Somfai et al. . What is perhaps more interesting is the prediction of resonant features which arise when $`c`$, the wedge angle relative to $`2\pi `$, is an integer multiple of $`\frac{\eta _1}{1+\eta _1},`$ corresponding to integer $`k^{},`$ because the *idea* of resonant angles in DLA has been much discussed Kessler et al. (1998); Turkevich and Scher (1985); Ball (1986b): this is the first direct and explicit prediction, which we look forward to seeing tested by simulations.

## X Growth Fluctuations

Numerical confidence in the scaling properties of DLA was greatly bolstered by the idea of an intrinsic but low level of self-organised noise Ball et al. (2002), so it is natural to ask if the present theory can address this. The simulation studies of noise have rested on tracking the extremal radius, which is hard to extract from our analytic formulation, so we have had to compromise on something more accessible theoretically. The relative penetration depth has fluctuations which reflect the differing geometry of the growth, and for the case of growth in channel we have measured these fluctuations to be $`(\delta \mathrm{\Xi }/\mathrm{\Xi })_{\text{simulation}}=0.18`$.

The corresponding theoretical calculation is a fairly straightforward generalisation of the penetration depth calculation itself: we simply calculate the average square of the expression in Eq. (22), minus the square of Eq. (28) to obtain the variance of $`\mathrm{\Xi }^2`$. The evaluation of this under the Gaussian Closure Approximation is detailed in Appendix D.

The most useful comparison is for the channel, $`c0`$, for which the penetration depth itself happened to be given rather accurately by the theory. In this case we obtain the relative variance of the square of the penetration depth as $`(\xi /R)^4/(\xi /R)^2^21=1.26`$, leading to $`(\delta \mathrm{\Xi }/\mathrm{\Xi })_{\text{theory}}=0.56`$ which is substantially higher than the simulation value. It is clear physically that the penetration depth comes predominantly from the lowest index modes of $`\lambda `$, and this is apparent from our expressions above if we linearise Eq. (22). Keeping only $`\lambda _1`$for the channel would then make $`\xi /R`$ be the magnitude of a single Gaussian distributed complex scalar, leading to very similar $`\delta \mathrm{\Xi }/\mathrm{\Xi }`$. Thus it seems to be quite fundamentally the Gaussian form of our closure approximation which leads to overstatement of the penetration depth fluctuations.

## XI Relation to Cone Angle Theory

The angular resonances predicted in the penetration depth turn out to be in interesting correspondence with part of the earlier Cone Angle Theory (CAT) of DLA Ball (1986b). In that theory a growing cluster was viewed as having an identifiable number of major arms $`n`$, and it was then further supposed that the growth should be marginally stable with respect to the loss of major arms through competition for growth. The strongest mode of such competition is where alternate fingers gain and lose, and the condition for this mode to marginally stable is in the present notation

$$\left(\frac{n/c}{2}1\right)\eta _0\alpha =1,$$
(23)
as calculated in ref Ball (1986b) for the case $`c=1`$ and $`\eta _0=1`$. In the CAT fractional $`n`$ was presumed an acceptable approximation and condition (23) was combined with geometrical approximations to predict $`\alpha `$, but here let us focus on the values of $`n`$ implied. Using $`\eta _1=\alpha \eta _0`$ this gives

$$\frac{n}{2}=c\frac{1+\eta _1}{\eta _1},$$
(24)
so our resonance condition corresponds directly to the case where the number of marginally stable major arms is an even number - which is of course required for the alternating mode stability calculation to be strictly applicable.

The CAT was closed in Ref. Ball (1986b) by approximating the cluster as a solid polygon of $`n`$ sides, for which $`\alpha =\frac{1}{1+2\nu /n}`$, leading to $`\alpha _{\text{CAT}}=\frac{1+\eta _0+\sqrt{1+6\eta _0+\eta _{0}^{}{}_{}{}^{2}}}{4\eta _0}`$ which is clearly quite different in principle from the Gaussian Closure prediction of constant $`\alpha _{\text{tip}}`$. However it is not easy to distinguish between them on the basis of previously published DBM data, as shown in Fig. 9, and unlike GCA the CAT does not predict any other exponents.

## XII Discussion

For DLA and its associated Dielectric Breakdown Models we have shown a theoretical framework which is complete in the sense that essentially all measurable quantities can be calculated. This extends to amplitude factors such as the relative penetration depth for which there is no theoretical precedent. For the full spectrum of exponents the practical advance over the Screened Growth Model is the elimination of fitting parameters, and it remains an open challenge to extend our theory to give quantitatively credible predictions for the large $`\alpha `$ part of the spectrum. For the exponent $`\alpha _{\text{tip}}`$ we have in the Gaussian approximation a striking new result that this is predicted constant over a wide range of $`\eta `$, which begs direct confirmation by (expensive) particle-based simulations. We look forward to addressing this in a following paper.

For the DBM at high $`\eta `$ we find structure more rich than discussed hitherto, with a Screening Transition intervening before the upper critical value $`\eta _0^c`$ is reached. Beyond the screening transition the scenario where $`\alpha _{\text{min}}=\alpha _{\text{tip}}`$ looks prospectively solvable (at least in terms of exponents) given the degenerate form of the Electrostatic Scaling Law which applies. The Gaussian Closure Approximation leads to the richer but quite possibly misleading scenario $`\alpha _{\text{min}}<\alpha _{\text{tip}}`$, so sorting out the truth of this inequality would be very interesting.

We have shown that the GCA naturally exhibits angular resonances which are in interesting correspondence with the ideas of the earlier Cone Angle Theory. Notably these resonances now have clear predicted consequences such as we demonstrated for the relative penetration depth, and they can be explored either by growing in a wedge of variable angle or by varying $`\eta `$ \- so once again behaviour vs $`\eta `$ is a key probe of our understanding of the problem. For DLA in particular the best theoretical value of $`\alpha _{\text{tip}}`$ remains $`1/\sqrt{2}0.71`$ from the CAT Ball (1986b), but the greater test now lies in the $`\eta `$ dependence on which CAT and GCA differ qualitatively.

Within DLA and DBM we look forward to calculating more properties such as the response to anisotropy, which is fairly readily incorporated into our equations of motion. The hard part is that in breaking angular symmetry we can no longer exclude non-zero first cumulants $`\lambda _k`$, and a full matrix of second cumulants, but the calculation is in principle straightforward. A conceptually more challenging avenue is to improve on the Gaussian approximation itself which we have used to obtain explicit theoretical results. Truncating at a cumulant of higher order than the second is hard, and more seriously it does not correspond to a positive (semi-)definite probability distribution. An alternative route of improvement which we are exploring is closure at the level of the full multifractal spectrum.

Whilst our main use of the equivalences within the class of $`(\eta ,m)`$ models has been to facilitate calculation, through mapping onto $`m=1`$, the particular associated claim that surface tension control is included through $`m=1`$ may prove controversial. This would imply that the scaling properties of the chaotic viscous fingering regime can be predicted from suitable DBM simulations. The DBM simulations required are relatively accessible and the greater difficulty in pursuing this agenda lies in obtaining suitably calibrated experimental data or accurate direct simulations of fingering out to high degrees of ramification.

There are possibilities for wider application of ideas in this paper, where we have formulated DLA and DBM as a turbulent dynamics governed by a complex scalar field in 1+1 dimensions. Decomposing this field multiplicatively (through Fourier representation of its logarithm) was the crucial step to obtain renormalisable equations and theoretical access to the multifractal behaviour, even though other representations offered equations of motion (18) with weaker non-linearity. It is natural to speculate whether the same strategy might apply to turbulent problems more widely, where the key issue appears to be identifying suitable fields to decompose multiplicatively which are of local physical significance, and subject to closed equations of motion.

###### Acknowledgements.

This research has been supported by the EC under Contract No. HPMF-CT-2000-00800.

## Appendix A Logarithmic variables

In this Appendix we present some details about the logarithmic variables in Eq. (6) and their equation of motion (7) and (8).

We begin by explaining the choice of the analytic form in Eq. (6). In terms of $`\omega =\varphi +i\theta ,`$ in which complex plane the region exterior to the growth is mapped to a \[half\] strip, Eq. (6) analytically continues to

$$\frac{dz}{d\omega }=cRe^{c\omega }e^{_{k>0;\text{ }kc}\lambda _ke^{k\omega }}.$$
(25)
Expanding the second exponential factor above to all orders gives

$$\begin{array}{c}\frac{dz}{d\omega }=cRe^{c\omega }(1\lambda _1e^\omega (\lambda _2\lambda _1^2/2)e^{2\omega }\hfill \\ \hfill (\lambda _3\mathrm{})e^{3\omega }\mathrm{})\end{array}$$
(26)
and integrating with respect to $`\omega `$ gives a conformal map from the right hand part of the strip to the exterior region of the wedge as required. When $`c`$ is positive integer, periodic boundary conditions put a restriction on the $`\lambda `$’s: $`\lambda _1=0`$ when $`c=1`$, and for the less interesting cases $`c=2,3,\mathrm{}`$, the prefactor in front of $`e^{c\omega }`$ has to vanish. The Gaussian Closure solutions (15) automatically satisfy this constraint, because they have $`\lambda _k=0`$ for $`kc`$.

In the region far from the growth the leading term dominates, giving $`z(\omega )=Re^{c\omega }`$ which shows the significance of $`R`$: it is the apparent radius of the growth (corresponding to $`\varphi =0`$) as seen from far away.

To obtain the transformed equation of motion (7), first take the logarithm of Eq. (6), giving $`\lambda (\theta )+\mathrm{ln}R=\mathrm{ln}\left(\frac{z}{\theta }\right)+\mathrm{ln}\left(i/ce^{ic\theta }\right)`$. Then differentiating both sides with respect to time (at constant $`\theta `$) gives

$$\begin{array}{cc}\hfill \frac{\lambda }{t}+\frac{R}{Rt}& =\left(\frac{z}{\theta }\right)^1\frac{}{\theta }\frac{z}{t}\hfill \\ & =\left(c+i\frac{\lambda }{\theta }\right)𝒫\left[(cR)^{2y}e^{y\left(\lambda +\overline{\lambda }\right)}\right]\hfill \\ & i\frac{}{\theta }𝒫\left[(cR)^{2y}e^{y\left(\lambda +\overline{\lambda }\right)}\right]\hfill \end{array}$$
(27)
where we have used $`\left|\theta /z\right|=(cR)^1e^{Re\lambda }`$ and subsequently the powers of $`R`$ can be taken outside $`𝒫`$. It is then trivial to take Fourier components of both sides to obtain Eqs. (7) and (8), the latter coming from the zeroth component which we chose to be absent from $`\lambda `$.

## Appendix B Cumulant expansion

The Kubo Cumulants for independent variables $`X_i`$ are given in terms of their moments by $`\mathrm{ln}e^{_i\beta _iX_i}=e^{_i\beta _iX_i}1_c,`$ where $`\beta _i`$ are arbitrary (scalar) parameters Kubo (1962). From this it is well known that by differentiation with respect to parameters one obtains $`X_1e^W=X_1e^W_ce^W`$ and $`X_1X_2e^W=\left(X_1X_2e^W_c+X_1e^W_cX_2e^W_c\right)e^W`$, where $`W_i\beta _iX_i`$, and similar results for higher moments such as $`X_1X_2X_3e^W`$.

To obtain the renormalised equation of motion (13) from the unrenormalised equation (12) we need first to apply the above to $`\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })}`$. Writing this as $`\frac{d\theta }{2\pi }e^{ik\theta }\overline{\lambda }_ke^{y(\lambda (\theta )+\overline{\lambda }(\theta ))}`$ allows us to apply the cumulant identities directly, giving $`\frac{d\theta }{2\pi }e^{ik\theta }\overline{\lambda }_ke^{y(\lambda (\theta )+\overline{\lambda }(\theta ))}_ce^{y(\lambda (\theta )+\overline{\lambda }(\theta ))}`$ and hence $`_j\overline{\lambda }_ke_{kj}^{y(\lambda +\overline{\lambda })}_ce_j^{y(\lambda +\overline{\lambda })}`$. In the translationally invariant case considered here only the $`j=0`$ term survives in the latter summation, leading to $`\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })}=\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })}_ce_0^{y(\lambda +\overline{\lambda })}`$ as required. The calculation of $`\lambda _{kj}\overline{\lambda }_ke_j^{y(\lambda +\overline{\lambda })}=\left(\lambda _{kj}\overline{\lambda }_ke_j^{y(\lambda +\overline{\lambda })}_c+\lambda _{kj}e_{jk}^{y(\lambda +\overline{\lambda })}_c\overline{\lambda }_ke_k^{y(\lambda +\overline{\lambda })}_c\right)\times e_0^{y(\lambda +\overline{\lambda })}`$ proceeds along precisely analogous lines.

## Appendix C Legendre transform of an inverse function

In Eq. (16) we are confronted with a slightly unusual situation: we wish to find the Legendre Transform $`f(\alpha )`$ of the function $`\tau (q)`$, that is $`f(\alpha )=q\alpha \tau (q)`$ where $`\alpha =d\tau /dq`$, given a simple form for the inverse function $`q(\tau ).`$ Let $`g(x)`$ be the Legendre Transform of $`q(\tau )`$, that is $`g(x)=x\tau q(\tau )`$ where $`x=dq/d\tau `$. Thus $`x=1/\alpha `$ and we have $`g(1/\alpha )=\tau /\alpha q=f(\alpha )/\alpha `$, so the two Legendre Transforms are very simply related.

The example needed from Eq. (16) has the form $`q(\tau )=1+\tau +b\tau ^2`$ leading to $`x=1+2b\tau `$ and hence $`g(x)=1+b\left(\frac{x1}{2b}\right)^2.`$ Then we have $`f(\alpha )=\alpha g(1/\alpha )`$ leading directly to Eq. (17).

## Appendix D Calculation of the penetration depth

We start from the expression given in Eq. (22) for the square of the relative penetration depth without any averaging over clusters. As the expression is simple exponential in the $`\lambda _k`$ it is straightforward to average over a Gaussian distribution leading to

$$\mathrm{\Xi }^2=\frac{\xi ^2}{R^2}=\frac{1}{2}\underset{k0}{}\frac{c^2}{(ck)^2}_0^{2\pi }\frac{d\sigma }{2\pi }e^{ik\sigma }\left(e^{_pS(p)e^{ip\sigma }}1\right).$$
(28)
Note that the term $`1`$ in the integrand makes no difference for $`k0`$, but we can include $`k=0`$ in the outer summation. This enables us to rearrange the combination of summation over $`k`$ and integration, giving

$$\mathrm{\Xi }^2=\frac{c^2}{2}_0^{\mathrm{}}𝑑xxe^{cx}\left(e^{_pS(p)e^{px}}1\right),$$
(29)
as is readily verified upon expanding the exponential $`e^{_pS(p)e^{px}}`$to all orders. The summation in this exponential can be evaluated in closed form using the variances from Eq. (15) and the standard forms $`u+u^3/3+u^5/5+\mathrm{}..=\mathrm{ln}\sqrt{\frac{1+u}{1u}}`$ and $`u^2/2+u^4/4+u^6/6+\mathrm{}.=\mathrm{ln}\sqrt{\frac{1}{1u^2}}`$ . This leaves one numerical quadrature to obtain the results shown in Fig. 8.

To compute the fluctuations in relative penetration depth we again start from Eq. (22) and now average its square to give

$$\begin{array}{cc}\hfill (\xi /R)^4=& \frac{1}{4}\underset{k>0}{}\frac{c^2}{(ck)^2}_0^{2\pi }\frac{d\theta }{2\pi }\frac{d\varphi }{2\pi }e^{ik(\theta \varphi )}\underset{k^{}>0}{}\frac{c^2}{(ck^{})^2}_0^{2\pi }\frac{d\theta ^{}}{2\pi }\frac{d\varphi ^{}}{2\pi }e^{ik^{}(\theta ^{}\varphi ^{})}\hfill \\ & \times \mathrm{exp}\left(\underset{p}{}\lambda _pe^{ip\theta }\underset{p}{}\overline{\lambda }_pe^{ip\varphi }\underset{p}{}\lambda _pe^{ip\theta ^{}}\underset{p}{}\overline{\lambda }_pe^{ip\varphi ^{}}\right)\hfill \end{array}$$
(30)
where averaging the last factor gives $`\mathrm{exp}\left(_pS(p)[e^{ip(\varphi \theta )}+e^{ip(\varphi \theta +\theta \theta ^{})}+e^{ip(\varphi ^{}\theta ^{}+\theta ^{}\theta )}+e^{ip(\varphi ^{}\theta ^{})}]\right)`$. One integration with respect to an absolute angle is now redundant, and the integrations with respect to $`\varphi \theta `$ and $`\varphi ^{}\theta ^{}`$ can be rearranged using the same trick as in calculating $`(\xi /R)^2`$ previously, leading after some cancellations to

$$\begin{array}{cc}\hfill \delta \left(\mathrm{\Xi }^2\right)^2& =\left(\frac{\xi }{R}\right)^4\left(\frac{\xi }{R}\right)^2^2\hfill \\ & =\frac{c^4}{4}_0^{\mathrm{}}𝑑xxe^{cx}_0^{\mathrm{}}𝑑x^{}x^{}e^{cx^{}}_0^{2\pi }\frac{d\psi }{2\pi }\left(e^{_pS(p)[e^{px}+e^{p(xi\psi )}+e^{p(x^{}+i\psi )}+e^{px^{}}]}e^{_pS(p)[e^{px}+e^{px^{}}]}\right).\hfill \end{array}$$
(31)
Above we used the simplification $`_0^{2\pi }\frac{d\psi }{2\pi }e^{_pS(p)e^{p(xi\psi )}}=1`$, which follows upon expanding the exponential where only the leading term survives. This can be further applied leading to the more compact form $`\delta \left(\mathrm{\Xi }^2\right)^2=_0^{2\pi }\frac{d\psi }{2\pi }\left|\frac{c^2}{2}_0^{\mathrm{}}𝑑xxe^{cx}e^{_pS(p)e^{px}}\left(e^{_pS(p)e^{p(xi\psi )}}1\right)\right|^2`$. Particularly in order to address the limit of channel growth, $`c0,`$ it is convenient to bypass normalisation conventions by looking at the relative fluctuations in $`\mathrm{\Xi }^2`$ which are now given by

$$\left(\frac{\delta \left(\mathrm{\Xi }^2\right)}{\mathrm{\Xi }^2}\right)^2=\frac{_0^{2\pi }\frac{d\psi }{2\pi }\left|_0^{\mathrm{}}𝑑xxe^{cx}e^{_pS(p)e^{px}}\left(e^{_pS(p)e^{p(xi\psi )}}1\right)\right|^2}{\left[_0^{\mathrm{}}𝑑xxe^{cx}\left(e^{_pS(p)e^{px}}1\right)\right]^2}.$$
(32)
From this we obtained the result cited in Sec. IX by numerical quadrature.