Free boundary problem and hydrodynamic limit

Jeremy Quastel University of Toronto

The following problem comes up in nuclear waste management^{1}. Radioactive waste is placed in a container and buried. Although the container may be very well designed, there is still some leakage, and so radioactive particles leave the container (very slowly we hope!) and then, as a rough approximation, perform independent Brownian motions. In order to contain them, chemical traps are distributed around the area. When a radioactive particle hits a trap it is destroyed. However, the traps also become saturated after a certain number of particles hit them. To keep matters simple we will assume here a homogeneous distribution of traps. The problem is, given the rate of production of particles and the density of traps, what is the shape and size of the set of saturated traps after a long time?

If you believe the shape is spherical, then the radius r(t) as a function of time can be computed. If c(t) is the rate of creation of particles at time t then the expected number of particles which hit the boundary of the sphere up to time T is

F (T) = ∫_{0}^{T} P ( sup_{s≤t≤T} β_{t}/r(t) ≥ 1)c(s)ds  (1)

where β_{t} is a Brownian motion starting at 0 at time s. On the other hand the traps in a sphere of radius r(T) will have eaten about αω_{d}r^{d}(T) particles. Here α is the density of traps which should really be counted according to how many particles they can destroy before they are saturated. The constant ω_{d} is the volume of the sphere of radius one in dimension d. So given λ(t) and ρ, the correct r(t) is the one for which F (T) = αω_{d}r^{d}(T). It looks a little daunting until you realize that (1) can be written as

F (T) = ∫_{0}^{T} c(t)dt − ∫ ρ(x, t)dx

where ρ is the solution of the Stefan problem

\begin{matrix}\frac{∂ρ}{∂t} = Δρ + c(t)δ_{0},  in \{ρ > 0\}, \\ ρ = 0  on s(x) ≥ t, \\ ∇ρ · ∇s = −α  on s(x) = t. \end{matrix}  (2)

If c(t) = t^{(d−2)/2} then in any dimension there is an explicit solution and the desired radius turns out to be K\sqrt{t} where exp\{−K^{2}/4\} = αK^{d}π^{d/2}/Γ(d/2). If c(t) ≡ 1 then r(t) ≈ t^{1/d} if d ≥ 3 and r(t) ≈ \sqrt{t log t} if d = 1. In d = 2 it is just as in the previous line.

On the other hand we have no particular reason to believe that the shape is a sphere, or for that matter that there is any sort of shape at all. So that is the real problem.

We could put the system on a lattice. In other words we replace the physical space R^{d} by Z^{d} and for each site x in the lattice let ζ_{x} and ξ_{x} denote the number of particles and the number of traps in a box of side length one around x. Then we approximate the original model by having the particles produced at the origin at rate c(s) and then perform independent continuous time random walks on the lattice until they reach a site x where ξ_{x} ≥ 1. Then that particle is destroyed and ξ_{x} is reduced to ξ_{x} − 1.

Such models have been studied by Lawler, Bramson and Griffeath under the name Internal Diffusion Limited Aggregation^{3}. In particular, they verified the spherical shape in the case c(s) ≡ 1 and ξ_{x} ≡ 1 in dimensions d ≥ 3. This was extended to d ≤ 2 by Gravner in joint work with the author^{4} and is presently being extended to the case where ξ_{x} is a stationary random field by Ramirez and the author. From now on we will stick to the case ξ_{x} ≡ 1 in order to illustrate the main ideas. First we want to explain the dimensional dependence of the results. In fact the dimension is not the key point, but the asymptotic behaviour of

t^{−d/2} ∫_{0}^{t} c(s)ds.

If this goes to zero as t gets large, then we say the model is subcritical. If it goes to infinity, we say the model is supercritical. In between is the critical case. The density inside the occupied set, in either the Stefan problem, or the Internal DLA, is asymptotically zero in the subcritical case, infinite in the supercritical case, and non-trivial in the critical case. The method of Lawler, Bramson and Griffeath, which is by comparison with a discrete time version of the model, works only in the subcritical regime. The method of Gravner and the author is by comparison with the Stefan problem. Although this comparison is only correct in the critical regime, the Lawler-Bramson-Griffeath method can be turned around and therefore it can be extended to the subcritical regime and recovers even the shape theorems for the discrete time versions. Not much is known in the supercritical regime except in one dimension where there is not much of a shape theorem to prove.

Critical case

Let us consider the critical case of Internal DLA. Particles are produced at the origin of Z^{d} at rate c(t) = t^{(d−2)/2} and then perform independent random walks until they reach a trap at which point the two are mutually destroyed. The number of particles at x at time t is called ζ_{x}(t) and the number of traps at x at time t is called ξ_{x}(t). Initially ξ_{x}(0) = 1 for each site. What we want to show is that in an appropriate weak sense

\begin{matrix}ζ_{⌊ϵ^{−1}x⌋}(ϵ^{−2}t) → ρ(x, t) \\ \{(x, t) : ξ_{⌊ϵ^{−1}x⌋}(ϵ^{−2}t) = 0\} → \{(x, t) : s(x) ≤ t\} \end{matrix}

as ϵ → 0 where ρ(x, t), s(x) is the solution of Stefan problem (2).

There are now several general methods to prove such hydrodynamic scaling limits. The most popular are based on the evolution of the relative entropy of the system with respect to global or local equilibrium measures. Ignoring for a moment the creation of particles at the origin, our model has two mutually singular equilibrium classes. The first has a trap at every site and no particle in sight, and the second class of invariant measures has no trap anywhere, while the particles are distributed according to a product of Poisson distributions with a fixed density ρ which parametrizes the class. Unfortunately the entropy methods are not well adapted to systems such as ours, where the key object of study is an unknown boundary separating regions in which the system is locally in one of two such equilibrium classes which have infinite entropy with respect to each other.

The proof is based on the following clever transformation of Stefan problem (2) which is well known to researchers in that field. If we define

u = ρ + 1(ρ > 0)  and  λ(u) = (u − 1) ∨ 0,

so that ρ = λ(u) then if ρ(x, t), s(x) is a solution of (2) an application of Green’s theorem shows that u(x, t) is a weak solution of

\frac{∂u}{∂t} = Δλ(u) + c(t)δ_{0}.  (2′)

In fact (2’) only has one solution, and the number of unknown functions has magically been reduced from two to one, so it is definitely an improvement over (2). To see that there is only one solution, we notice that it is in the form of a degenerate porous medium equation and the porous medium equation has the nice property that it contracts the H_{−1} norm,

∥f∥_{−1}^{2} = ∫ f(−Δ^{−1})fdx.

In fact if u and u′ are two solutions to (2’) then one can easily check that

∥u(T) − u′(T)∥_{−1}^{2} = − ∫_{0}^{T} ∫ (u − u′)(λ(u) − λ(u′))dxdt.  (3)

The integral on the right is non-negative whenever λ is a non-decreasing function of u. We can try a similar trick for the Internal DLA. Let

η_{x} = ζ_{x} + (1 − ξ_{x})

be the number of live or dead particles at site x in the lattice. Then η by itself is a Markov process and its generator is

Lf(η) = ∑_{x,|e|=1} λ(η_{x}) (f(η^{x,x+e}) − f(η)) + t^{(d−2)/2} (f(η_{+}^{0}) − f(η)) .

Here η^{x,y} means the new configuration after one particle has jumped from x to y, i.e., η_{x}^{x,y} = η_{x} − 1, η_{y}^{x,y} = η_{y} + 1 and η_{z}^{x,y} = η_{z} otherwise. Similarly η_{+}^{0} means the new configuration after one new particle has been created at the origin. Notice that the process takes the form of a degenerate zero-range model.

The idea, which first appeared in work of Chang and Yau^{5} in a different context, is to compare the microscopic system to the macroscopic one using the H_{−1} norm. We work on the lattice ϵZ^{d} on which we have the Internal DLA model η^{ϵ} running at rate ϵ^{−2}. We compare it to the lattice version of the Stefan problem (2’),

\frac{∂u^{ϵ}}{∂t} = Δ_{ϵ}λ(u^{ϵ}) + dP^{ϵ},  (2′′)

where the Poisson process P^{ϵ}(t) is just the number of particles created in the Internal DLA model up to time t, and the lattice Laplacian is given by

Δ_{ϵ}f_{x} = ϵ^{−2} ∑_{|e|=ϵ}(f_{x+e} − f_{x}).

We also define the H_{−1} norm of a function on ϵZ^{d} by

∥f∥_{−1,ϵ}^{2} = ϵ^{d} ∑_{x∈ϵZ^{d}} f_{x}(−Δ_{ϵ}^{−1}f)_{x}.

Now you apply the generator to the H_{−1} of the difference of u^{ϵ} and η^{ϵ} and after some computation you obtain the following remarkable formula,

E[∥η^{ϵ}(T) − u^{ϵ}(T)∥_{−1,ϵ}^{2}] = −2 ∫_{0}^{T} ϵ^{d} ∑_{x} E[V (η_{x}^{ϵ}(t), u_{x}^{ϵ}(t))]dt  (3′)

where

V (n, u) = (n − u)(λ(n) − λ(u)) − λ(n).

It looks as if something has gone wrong. On the right hand side the integrand is the same sort of non-negative object as in (3) except that now it has an extra term −λ(n) appended at the end which seems to spoil everything. Of course this term is really very important. If it was not there then we would have shown that for fixed ϵ the H_{−1} norm of the difference of η^{ϵ} and u^{ϵ} is zero which is patently false because we expect that inside the occupied set the distribution of the particles is more or less a product of Poissons whose parameter is approximately the solution of the Stefan problem. Therefore the convergence can only be weak and correct in the limit as ϵ is sent to zero.

What happens is that because of the time integral, which is really over an enormous time ϵ^{−2} and the fact that typical sites are far from the origin we can pretend that the expectation of V is being taken locally with respect to an invariant measure of the process without creation. Furthermore we can transfer whatever smoothness can be obtained for u^{ϵ} into some spatial averaging, so the invariant measure can also be assumed to be translation invariant. One can use some coupling arguments^{6} to show that all such invariant measures fit into the two classes described a few paragraphs ago and the punch line is that for any measure μ in those classes and any u whatsoever,

E_{μ}[V (η_{0}, u)] ≤ 0.

The reason it happens is that the extra term −λ(n) is exactly the right thing to compensate for the fact that if N is a Poisson random variable with expectation m then E[N^{2}] = m^{2} + m.

The details of the proof can be found in 4. We also learned at the meeting that Funaki^{8} has recently studied the same model without creation and shown convergence to the Stefan problem. We note that our method can be adapted to allow the particles to interact to some extent before they are killed. One replaces the function λ(n) in the generator L by g(n − 1) ∨ 0 where g is some general increasing function with linear growth. Then λ(u) in the Stefan problem (2’) is replaced by γ(u − 1) ∨ 0 where γ is an appropriate function depending on the choice of g and the result is a non-linear version of the Stefan problem. One can also have several creation sites and there are analogous results in that case as well.

Supercritical case

In the supercritical case we do not know how to prove any shape theorem so we are restricted to d = 1. We let c(t) = t^{α} with α > −1/2 and compute the size of the occupied set by studying the asymptotics of (2). Note first of all that

u(x, t) = ∫_{0}^{t} |t − s|^{α} \frac{e^{−x^{2}/4s}}{\sqrt{4πs}} ds

satisfies the equation u_{t} = u_{xx}+δ_{0}. We look for a function r(t) such that ρ̂(x, t) = u(x, t)− u(r(t), t) for |x| ≤ r(t) is an approximate solution of (2). We need −u_{x}(r(t), t) ∼ r′(t) and ∂_{t}u(r(t), t) ∼ 0. Then we can construct out of ρ̂ appropriate sub and super-solutions to show that r(t) is the real size of the boundary. Differentiating and changing variables we get

−u_{x}(x, t) = \frac{t^{α}}{\sqrt{π}} ∫_{x^{2}/4t}^{∞} |u − \frac{x^{2}}{4t} |^{α}u^{−(α+1/2)}e^{−u}du.

Using the fact that ∫_{a}^{∞} |u − a|^{α}u^{−(α+1/2)}e^{−u}du ∼ a^{−(α+1/2)}e^{−a} as a → ∞ one can check that

r(t) = \sqrt{4t log \{t^{−d/2} ∫_{0}^{t} c(s)ds\}}

is the radius of the occupied set. The case c(t) = 1 was previously obtained^{4,7} by much more complicated methods.

Subcritical case

The discrete time version of Internal DLA is a Markov chain A(0), A(1), … whose state space is the set of connected sets in the lattice Z^{d} containing 0 and whose transition probabilities are given by p(A, B) = 0 unless B can be obtained from A by adding on one site x in its boundary, in which case p(A, B) is the probability that a random walk starting at the origin hits the boundary of A at x.

Let A_{t} denote the occupied set (set of saturated obstacles) in the continuous time version of Internal DLA with some creation rate c(t). We now give a rough sketch of the argument of Lawler, Bramson and Griffeath^{3} comparing A_{t} with the discrete time version A(n). One numbers the particles according to their time of birth, and lets A_{t}(n) be the occupied set of the first n particles, with the rule that if an older particle arrives at a site where an younger particle has died, the older particle replaces the younger one as dead at that site and the younger particle goes off on its merry way. A_{t}(n) is always a subset of both A_{t} and A(n) and, for fixed n, A_{t}(n) increases to A(n) as t increases to infinity since each of the first n particles will after some time have found its resting place. If we let P (t) be the number of particles created up to time t in the continuous model, then P (t) is about n(t) ≐ ∫_{0}^{t} c(s)ds and A_{t}(P (t)) = A_{t} so we have

A_{t}(n(t)) ≃ A_{t}

the only error coming from the surplus or surfeit of particles produced up to time t. By ≃ we mean that the measure of the symmetric difference normalized by the measure of the union, goes to zero as t → ∞. In order to compare A_{t} with A(n) the key point is therefore whether

A_{t}(n(t)) ≃ A(n(t))

which is a way of saying that the density of live particles is going to zero, or, in other words, that we are in the subcritical case. The precise argument can be found in reference 6 and the conclusion is that the shape theorems for A(n) and A_{t} can be deduced from each other. Lawler, Bramson, and Griffeath then used probabilistic potential theory methods to show that A(n) is asymptotically a sphere. Actually, their argument works in the much stronger Hausdorff metric, i.e. they show that for sufficiently large n, n^{−1/d}A(n) contains any sphere centered at the origin of volume less than one, and is contained in any sphere centered at the origin of volume greater than one.

What we want to do here is deduce the shape theorems for both A(n) and A_{t} from a hydrodynamic limit.

We introduce a ‘soft’ version of Internal DLA with creation rate c(t). Particles are produced at the origin in Z^{d} at rate ϵ^{−d}c(t) and perform independent random walks at rate ϵ^{−2} until they reach a trap, at which point the two are mutually annihilated as usual. In fact the hydrodynamic methods apply to this critical model without modification and the density and occupied set converge to those of the subcritical Stefan problem (2). If we let A_{t}^{ϵ} denote the occupied set of the soft model and A_{t} denote the occupied set of Stefan problem (2) then what we have is that as ϵ goes to zero,

ϵA_{t}^{ϵ} ≃ A_{t}.

Although A_{t} is spherical, its radius is not quite explicit. But from the subcriticality it is clear that as t goes to infinity, n(t)^{−1/d}A_{t} converges to a ball of volume one. Now let P^{ϵ}(t) be the number of particles created in the soft model up to time t. We have P^{ϵ}(t) ≃ n(ϵ, t) ≐ ϵ^{−d}n(t) and we have A_{t}^{ϵ}(n(ϵ, t)) ≃ A_{t}^{ϵ} and therefore, as ϵ goes to zero then t goes to infinity,

n(ϵ, t)^{−1/d}A_{t}^{ϵ}(n(ϵ, t)) ≃ A ball of volume one.

But also from the subcriticality we have that

A_{t}^{ϵ}(n(ϵ, t)) ≃ A(n(ϵ, t))

which completes the argument.

The reason that all subcritical models have the same asymptotics is that if the density of live particles is going to zero the growth of the boundary is essentially by harmonic measure from the origin. In the critical case where there is a non-trivial density of particles in the region, the growth of the boundary is essentially by the average of the harmonic measure from all points inside the region, weighted by the density.

Notes and References.

1. I learned about this application from Gerard BenArous.

2. We have normalized things so that our Brownian motions have generator Δ.

3. Lawler, G., Bramson, M., Griffeath, D., Internal diffusion limited aggregation, Ann. Prob.20 (1992), no. 4, 2117–2140.

4. Gravner, J., Quastel, J., Internal DLA and the Stefan problem, preprint.

5. Chang, C.C., Yau, H.-T., Fluctuations of one-dimensional Ginzburg-Landau models in nonequilibrium. Comm. Math.. Phys. 145 (1992), no. 2, 209–234., Yau, H.-T., Metastability of Ginzburg-Landau model with a conservation law. J. Statist. Phys. 74 (1994), no. 3-4, 705–742.

6. Liggett, T., Coupling the simple exclusion process, Ann. Prob., 4 (1976), 339–356., Andjel E., Invariant measures for the zero-range processes, Ann. Prob. 10 525–547 (1982).

7. BenArous, G., Ramirez, A., Diffusion and saturation processes in random media, preprint.

8. Funaki, T., preprint.