Laplacian patterns

Francisco Vera
Pontificia Universidad Católica de Valparaíso, Av. Brasil 2950, Valparaíso, Chile


We propose a simple universal mechanism for explaining the spontaneous appearance of complex patterns in many natural systems.

Complex patterns are produced spontaneously in a diversity of different natural systems, some typical examples are: snow crystals, spirals, the tortuous path of lightning, biological systems, etc. It is easy to find many similarities between patterns formed in completely different systems, that suggest the existence of an universal explanation for their appearance. In the last thirty years we have seen many advances in the theoretical description of each particular system, but a general model that explains the appearance of complex universal patterns is still missing. The common belief in the scientific community is that no simple general principle can explain this diversity of complex patterns.

Two years ago we started the search for an universal energy functional from which most of these patterns should appear. After simplifying the typical dynamical equations for several systems, elimination of microscopic details, and neglecting time dependent variables, we were left with a simple laplace equation. Our starting point was then to study a single scalar field that obeys laplace equation.

The most general form of laplace equation includes a conductivity that depends on spatial coordinates, it turns out that in our model this conductivity will have different values inside the pattern and in the region outside the pattern. Laplace equation is a good starting point, because this equation is a consequence of a conservation law, for the flux of a vectorial field derived from the scalar field. It is well known that laplace equation can be obtained from a variational principle, and the energy functional can be found in any textbook of electromagnetism. We found that this energy functional can produce a great variety of complex patterns, when used in a quasi-static model that let the system to probe different configurations for the conductivity and evolve towards the one where the energy functional is maximized.

Using our model and simple numerical algorithms, we have obtained non-trivial results for patterns in dielectric breakdown and snow crystals, our model has no free parameter to adjust, and all the parameters have a direct relationship with measurable quantities.

The model: To fix ideas and understand the underlying physics, let us to explain our model using an electric scalar field $\phi$ inside a parallel plate capacitor. After imposing boundary conditions and given the dielectric permittivity $\epsilon$, which in general depends on coordinates, one must solve laplace equation

\end{displaymath} (1)

to obtain $\phi$ in the region of interest. It is well known that laplace equation can be derived using a variational principle from the energy functional

\end{displaymath} (2)

$U$ is also the total energy in the capacitor when the charge $Q$ is maintained constant. In the case of a parallel plate capacitor this energy will be proportional to $d/\epsilon$ where $d$ is the separation between the plates. Static physical systems evolve trying to reduce the total energy, this implies a force between the plates trying to reduce the distance $d$, and if you insert a slab of a material having a dielectric permitivity $\epsilon$' greater than $\epsilon$, the slab will be pulled into the capacitor. If instead of maintaining $Q$ constant, the potential difference between the plates $V$ is maintained constant, the energy $U$ for this system will be proportional to $\epsilon/d$. In experiments at constant $Q$ or $V$ there are forces trying to reduce d and to increase $\epsilon$, then in experiments at constant $V$, $U$ is not the total energy of the system, because the system evolves trying to increase $U$. There is a missing energy term that comes from a rearrangement of charges in the wires to maintain $V$ constant. When this term is introduced, the total energy is again proportional to $d/\epsilon$. The previous discussion justifies that we study systems at constant $V$ by letting them to evolve towards regions of higher $U$.

If V is large enough the dielectric material will break causing a short circuit, we propose that for large V there is a new possibility for the system to lower the total energy, a short circuit. Assuming that $\epsilon$ can evolve locally towards bigger values, forming a discharge channel, we could obtain complex patterns of dielectric permitivity by letting the system to evolve towards higher values of U. Our model is quasi-static (see below), but there are two implicit time scales in the problem: a slow evolution for the local change in $\epsilon$, and a fast evolution for the electric potential $\phi$ after a local change in $\epsilon$. In general a higher value for the dielectric permitivity corresponds to a higher value for the system conductivity. In the electric case, complex patterns arise because the system is trying to lower the total energy by increasing locally the conductivity. It is not clear at this point why this evolution would produce complex patterns, and not just a straight path of bigger conductivity. We have no simple answer to this question, because the evolution for the system is highly nonlinear, and we have to rely on numerical methods to study this evolution.

We study these systems using a quasi-static treatment as follows: First, set the boundary conditions which are maintained trough all the steps in our simulation. Second, assign a fixed value $\epsilon$ to the dielectric permitivity inside the boundaries. Third, using laplace equation, obtain $\phi$ after changing the dielectric permitivity locally to a greater value $\epsilon'$, near one electrode. Fourth, find the energy U using the new values for the dielectric permitivity and the scalar field. Fifth, repeat steps 3 and 4 to obtain the energy values U for each of the neighbors of one electrode. These energies are compared and the neighbor providing the biggest energy value is added to the channel. These steps are repeated including the new neighbors for the channel, until a pattern develops. We used a square lattice, pattern evolution through the diagonals is not permitted, and the pattern grows adding only one site to the evolving pattern at each time step.

Our first example: The tortuous behavior of lightning. When a charged conductor has a sharp end or tip, the electric field just outside this region is bigger than the field outside other points that are far from the tip, and when the electric field exceeds a certain value, the dielectric medium (air for example) will break and a discharge or spark will be initiated at the tip. A similar scenario occurs in thunderstorms [1], where air currents separate negative and positive charges producing regions of high electric fields, when this field exceeds the value for breaking the air a lightning stroke is produced. As the tip of this leader nears ground, the electric field at sharp objects on the ground increases until it exceeds the breakdown strength of air. At that time, one or more upward-moving discharges are initiated from those points, and an attachment with the downward-moving leader occurs some tens of meters above ground. The lightning discharge in thunderstorms and sparks between charged conductors, span scales from millimeters to kilometers and evolve following a tortuous path, forming a branched structure that closely resembles a fractal [2].

When one is faced with the problem of simulating lightning the first thing to do is to simplify the geometry of the problem. The most used geometries are: a circular conductor (outer electrode) plus a central electrode (see Fig. 1), and an elongated rectangular region with a electrode in one end and the other end acting as the other electrode. Then an electric potential is applied to each electrode, the difference of potential must be big enough to reach the necessary electric field to break the dielectric medium. All conclusions obtained in these two-dimensional geometries can be readily applied to the real 3-dimensional case.

In the theoretical study of the evolution of the discharge channel, the simplest possibility is to expand the inner electrode towards the regions of maximum electric fields. The presently accepted model of lightning [3,4] was developed by Niemeyer, Pietronero and Wiesmann in 1984, and follow these steps but includes a stochastic term that weights a probability, that is a function of the value of the local electric field. This is known as the Dielectric Breakdown Model (DBM) and produces a branched structure whose fractal dimension is similar to the ones obtained experimentally for the same geometry, Pietronero model is based on the Diffusion Limited Aggregation (DLA) model, developed by Witten and Sander in 1981[5].

Our results show that it is possible to obtain a branched structure of lightning that follow from a deterministic treatment[6], that only relies in minimizing the total energy in the system and changing locally the dielectric permitivity of the medium (not the geometry of the inner electrode) at each step of iteration. We note that our model is almost all of the time deterministic, but for certain configurations there is degeneracy for the extreme value of the energy $U$, and the numerical noise will be responsible for selecting the next step in the evolution of the pattern.

We will study dielectric breakdown in the circular geometry of Fig. 1, considering a $20\times20$ two-dimensional square lattice, where the central point is the inner electrode and the outer electrode is modeled as a circle. The boundary conditions, $\phi=0$ in the inner electrode and $\phi=1$ in the outer electrode, are maintained trough all the steps in our simulation. For each different configuration, the numerical solution of laplace equation was obtained using a Gauss-Seidel algorithm and accepted when the numerical residual was less than $10^{-2}$, the values for the dielectric permitivity outside the channel was $\epsilon=1$ and inside the channel was $\epsilon'=5$. The filled boxes represent the discharge channel (sites where the permitivity is $\epsilon'$). The circles show all possible sites where the channel can evolve, the diameter of each of these circles represent the value for the energy U of the system in case this site is added to the channel.

Figure: Discharge channel, for a $20\times20$ square lattice, showing the attachment between the return branch (left black box) and the main branch (inner gray boxes).

Upward-moving discharges initiated from earth attach with the downward-moving leader in real lightning. Our calculation considered two initial branches: the main branch coming from the central lattice site and the return branch coming from the outer electrode. The return branch begin to evolve, after 36 steps of evolution of the main branch, at the site having the black box inside (at the extreme left of the figure). The site between this point and the main branch is the next step in the evolution, completing the path for this discharge channel from the inner electrode towards the outer electrode. Obtaining the return branch in our numerical simulations is a really satisfactory finding, because this is a highly non trivial result and we do not know of any other work that can obtain this attachment.

Our second example: Snow crystals. We all have admired snowflakes and other beautiful examples of solidification patterns, these crystals often come with a six-fold symmetrical shape as a consequence that the most typical basic form of an ice crystal is a hexagonal prism. Other basic forms occur at very low temperatures or very high pressures. The observation of snowflakes and these sub-millimeter snow crystals can be traced to 1611 with the work of Johannes Kepler entitled "The Six-Cornered Snowflake"[7]. After that, several authors have continued the experimental and theoretical study of snow crystals[8]. The preliminary observations together with the latest systematic observations and experiments have provided the basis for understanding the basic mechanism responsible for the creation of these fascinating structures. There are many good reviews on the subject of pattern formation[9,10,11,12,13] and one of the most used methods for studying solidification is the phase field method[14]. There are also two specific models that generate patterns for snow crystals similar to the patterns found in nature[16,15]. The origin of structures in these systems is known to be related to the Mullins-Sekerka instability[17], Mullins and Sekerka studied why a protrusion in the interface causes the interface to grow faster there, but the theory for the evolution of a growing protrusion towards complex structures has remained a mystery.

By analogy with our model of dielectric breakdown, where the scalar fields correspond to the temperature $T$ and the thermal conductivity $k$, we can study systems with boundary conditions for the temperature $T$, by letting them to evolve towards regions of higher $U$[18]. This turn out to be equivalent to the evolution towards a pattern that minimize the total energy in the scalar field $T$, after a local change in the thermal conductivity. Laplace equation for this case reads

\nabla\cdot(k\nabla T)=0,
\end{displaymath} (3)

We will study the evolution of a central seed in the circular geometry of Fig. 2. We consider a two-dimensional square lattice, where the central point is the seed and the outer interface is modeled as a discrete version of the big circle. The boundary conditions, $T=0$ in the seed and $T=1$ in the outer interface, are maintained trough all the steps in our simulation. The first step is to assign a fixed value for the thermal conductivity of each lattice point between the seed and the outer interface, and set the initial pattern as the central point. The second step is to obtain the energies $U$ of the system after changing the thermal conductivity in one of each neighbor of the pattern to a greater value $k'$, these energies are compared and the neighbor providing the biggest energy value is added to the pattern.

Fig. 2 shows the structure of the pattern for a $70\times70$ square lattice after 750 iterations. For each different configuration, the numerical solution of Eq. 2 was obtained using a Successive Over Relaxation algorithm and accepted when the numerical residual was less than $10^{-1}$, the values for the thermal conductivity outside the channel was $k=2$ and inside the channel was $k'=6$. The evolution of the pattern shows that opposite branches are not exactly aligned, we consider this as a prediction of our model and it is possible to find this effect in the experimental results shown in Fig. 3 of ref [9]. Our example also shows that the system develops, forming initially a central compact core, this core supports the evolution of the main venous branches and diagonal branches. We expect that secondary branches emerging from the main branches would appear after some additional iterations. Our simulated structure lacks the six-fold symmetry of real snow crystals, this is a consequence of using a square lattice in our numerical model. We expect that realistic six-fold patterns will be formed in a lattice that imposes the hexagonal symmetry. We have not tested our model changing the supporting lattice, because several months of computing time was needed for completing this example. Repetitive calculations, or calculations using bigger lattices, would require implementing sophisticated numerical methods[14].

Figure: Pattern developed for a $70\times70$ square lattice after 750 iterations for thermal conductivity outside the channel $k=2$ and inside the channel $k'=6$.

The large computing time needed for these simulations is a consequence of the necessity for the system to probe different possible configurations and select the one minimizing the total energy. This is the famous principle of least action that is nicely explained in the textbook of Richard Feynman[19]. In normal systems, like a mass falling towards the earth surface as a consequence of the force of gravity, this principle leads to simple differential equations for the motion of particles. In the case studied in this article there are not trajectories emerging from this principle, instead there are complex configurations of the inner interface as a result of minimizing this energy.

Our third example: Viscous fingers. Two-dimensional experiments, where a less viscous fluid is injected into a more viscous fluid, show the appearance of patterns[20,21]. These patterns were discovered for the first time when water was injected to extract the remaining oil deposited in the bottom of natural wells; the expectation was that the high density fluid (water) would be deposited below the lower density fluid (oil) allowing extracting the oil easily. Contrary to what was expected, water invaded the oil forming patterns and frustrated the extraction. It is very easy to produce these patterns in the laboratory, and many experiments have been reported for many different geometries and fluids. These patterns are believed to form as a consequence of a dynamical instability in the governing equations; here we will offer a completely different explanation for their appearance.

After elimination of all the dynamics in the governing equations, we are left again with laplace equation for the scalar field P

\nabla\cdot(k\nabla P)=0,
\end{displaymath} (4)

where $k=b^{2}/(12\eta)$, $P$ is the pressure, $\eta$ the viscosity, and $b$ is the separation between the plates in a two-dimensional Hele-Shaw experiment.

Figure: Solution of laplace equation showing curves of constant pressure, the conductivity $k=1$ everywhere and $k'=20$ in the central black region at the bottom.

We are in the process of implementing adaptive grid algorithms to use a cluster of computers for solving laplace equation. In fig. 3 we see the solution using $P=0$ at the bottom, $P=1$ at the top, and linear boundary conditions at the left and right boundaries. We divided the domain in a $5\times5$ grid setting $k=1$ everywhere, and $k'=20$ in the central black region at the bottom. Using 70000 nodes for the adaptive grid, we can see smooth curves for regions of constant pressure. Fig. 4 show the grid for the previous solution when the calculation is made using 3 processors in the cluster. Each color represents the grid that is solved by each processor in our cluster of computers. We are trying to get control over the numerical error before attempting to generate patterns using these algorithms.

Figure: Adaptive grid after solving laplace equation using three processors.

General remarks: Our model needs much more computing time than other models in the literature, because the system has to probe different possible configurations and select the one minimizing the total energy. We have been very careful for not making more than one point of the interface to evolve at every time step, because evolving more points in a given time step would be completely arbitrary. We have seen the appearance of a return branch, which is typically found in real lightning. We obtained the breaking of chiral symmetry in a global pattern, as a consequence of time evolution. Almost all the patterns obtained using our model, show branches that are roughly identical in length, this is one of the most important unsolved questions in the area, and using our model we obtained this result straightforward by letting the system to evolve locally towards a configuration of minimal energy. For the same set of parameters we have seen a pattern to grow first in a compact form, and after reaching a critical size the system begin to form branches.

But, why this simple model can produce the main structures seen in such diversity of systems? Can we answer for example why water injected into oil form patterns, but oil injected into water show only a compact grow?

In all the systems studied or mentioned in this work, there is conservation of flux and consequently they obey laplace equation. These systems are in a static configuration for the scalar field, but they are not in a state of minimal energy. All these systems have a conductivity that has the possibility of changing its value locally. When the system is able to increase this conductivity, a channel of bigger conductivity can be formed between boundaries (having different fixed values for the scalar field); this would allow the system to equalize the values for the scalar field at the boundaries, releasing the stored energy. An increase in dielectric permitivity, the freezing of water, substitution of oil by water, etc., have in common that these systems are increasing locally their conductivities for electric field, temperature and pressure respectively.

We hope that this mechanism of reducing energy by increasing conductivities could be applied to many others different pattern forming systems (biological systems, spiral galaxies, etc.), and we believe that this mechanism will be the key for beginning a global quantitative understanding of complexity.

Francisco Vera 2003-08-20