text:
# Stable propagation of an ordered array of cracks during directional drying

## I Introduction

Cracks appearing when a material shrinks are common in everyday life. The most popular examples are probably the cracks appearing in paint layers, and those in the surface of muds. In the case of paints, there is a superficial layer of material that shrinks on top of a substrate to which it is attached. For muds, the superficial layer and the substrate are the same material, but difference in the humidity concentration produces the cracking at the surface. In these two cases, the shrinking is due to changes in the humidity concentration within the sample, but it can also be due to the existence of non-uniform temperature distributions, which produces stresses and generates cracking. The problem of surface fragmentation has been studied in the last years, both theoretically and experimentally.

There are situations in which cracks appear at the surface of the material, and penetrate into the sample later on. A well known example corresponds to the columnar fracturing of basaltic lava, seen in many different geographical locations. A detailed description of this problem has only recently been foreshadowed, and still some points remain obscure. The two dimensional equivalent of the columnar structure of basalts corresponds to the case of a two dimensional material that dries (or cools down) starting from a free edge. As in the three dimensional case, drying creates internal stresses in the material that generate cracks, first appearing at the free edge. As the drying front propagates to the interior the cracks also do so in turn. There have been different experimental realizations of this phenomenon. In one of them the material was a very thin layer of a colloidal dispersion placed between two glasses, with one free edge, from which humidity can escape. In other experiments, a slurry of Al2O3-water mixture was deposited onto a substrate, and a glass was placed a few millimeters above it. A low, ultra dry $`N_2`$ gas breeze was injected in the slot above the sample. The $`N_2`$ became saturated with water as it passed over the material, and a rather sharp drying front propagates with time in the same direction than the air flow. A third realization corresponds to the drying of thin films of aqueous silica sol-gel, where very nice patterns have been observed. In all these cases the propagation of a set of parallel cracks has been observed.

A detailed theoretical description of this phenomenon is lacking, and we will make an attempt in this direction for the two dimensional case. In the next section we will make a general description of the phenomenon, emphasizing the difference with a more standard fracture mechanics problem. In Section III we present the calculations of elastic and crack energy for the case in which the drying (or cooling) front is sharp. In Section IV we discuss the conditions under which the crack front (formed by the tips of all cracks) is stable, and progresses into the material when the drying front advances. In Section V we show that the array of cracks generated at the surface of the material becomes evenly spaced when penetrates the sample, and that the typical width of the stripes is determined mainly during the first stage of the process, near the surface. In Section VI we briefly comment upon the effects of other drying (or cooling) conditions than the one studied previously. Section VII contains some implications of our work for the study of the three dimensional case. Finally, in Section VIII we summarize and conclude.

## II Quasistatic fracture mechanics. Energetic and stress analysis

One remarkable thing about the experiments cited in the previous section  is the fact that a large number of cracks penetrate the sample in a coordinated and quasistatic manner, as the external conditions (the humidity profile) change. Then this problem is qualitatively different of a typical problem in fracture mechanics, where the propagation of a crack is usually an abrupt phenomenon and typically leads to the failure of the material. In standard fracture mechanics, stress and energetic analysis are two different ways of predicting the evolution of the fracturing process. However, the equivalence of both approaches is not clear, particularly in cases in which the cracks propagate at large speeds, which is almost always the case when failure occurs.

The situation is different in our case. Cracks propagate only because the external humidity profile changes. If the external conditions were stationary, crack advance would be arrested. If we think on the configuration of the system as a point in configuration space (the space spanned by all coordinates of all particles of the material), at each moment the system is at one minimum of the energy landscape. As the humidity profiles changes the landscape changes itself, the minimum on which the system is located shifts, and the configuration of the system adapts so as to remain in the shifted local minimum. This is what we understand as a quasistatic propagation of cracks.

Under these conditions the propagation of cracks in a two dimensional geometry can be studied by two different but equivalent procedures: Stress analysis consists on the calculation of the stress intensity factors $`K_I`$ and $`K_{II}`$ of the opening and shearing modes at the tips of the cracks present in the system. The propagation will occur in the direction along which $`K_{II}=0`$, which coincides with that for which $`K_I`$ is maximum. Propagation actually occurs only if the energy relieved by the advance of the cracks is enough to overcome the fracture energy needed to elongate the cracks. In the quasistatic case, these two energies will differ only infinitesimally.

In the energetic procedure to calculate crack advance, the total energy (including elastic and crack energy) after virtual advances of the cracks is calculated. Cracks will actually propagate only if this propagation reduces the total energy. Under quasistatic advance the energy reduction during propagation is infinitesimally small, and typically, there is only one possible direction of propagation for each crack, all other propagation directions would produce an increase of the total energy of the system.

Both approaches are equivalent in the case of quasistatic advance of the cracks. We will use stress analysis or energetic analysis according to convenience in each case.

## III Estimations of elastic and crack energy

We will consider an isotropic and homogeneous material, and assume that linear elasticity can be applied. It is well known that in this case the material possesses only two independent elastic constants. We will take as these two parameters the bulk modulus $`C`$ and the Poisson ratio $`\nu `$. The expansion or contraction properties of the material are described by a humidity expansion coefficient $`\alpha `$, which is formally equivalent to the thermal expansion coefficient, i.e., the relative change of linear size $`\delta L/L`$ of a piece of material after changing the humidity concentration $`h`$ by some quantity $`\mathrm{\Delta }h`$ is given by

$$\delta L/L=\alpha \mathrm{\Delta }h.$$
(1)
The ideal situation we will address consists of a semi-infinite two-dimensional sample that is being dried from its surface. The drying will be considered to be non-homogeneous, and we will model it by a humidity profile that is given by some function $`h(z,t)`$ depending on depth $`z`$ within the material, and time $`t`$. The precise form of this function will be specified later for different experimental situations, but it is important to point out that we do not consider the case in which the appearance of cracks in the system modifies itself the drying process. Then our analysis applies more to cases as that described in  and possibly , but much less to that in  where evaporation of humidity through the cracks seems to be relevant(see also ). Although in the experimental situation the material usually lies on top of a substrate, we will study the case in which there is no interaction with this substrate. In the experiments of Ref.  this is achieved by introducing a layer of some slippery material between the sample and the substrate. Then, all the stresses on the sample are originated internally, and are due to the existence of a non-uniform humidity concentration.

It is known that in an isolated (namely, not clamped) piece of homogeneous material placed under a constant thermal (or humidity) gradient all stresses vanish in linear elastic approximation . Under these conditions cracks cannot appear at all, or if already present from the beginning its propagation is completely halted. It is in fact crucial for the propagation of cracks that the humidity gradient it is not constant within the sample.

We will take the edge of the semi-infinite sample as the $`z=0`$ line, and $`z>0`$ in the interior. For convenience, we will also refer to the edge as the ‘surface’, and consider it to be horizontally placed, in such a way that the cracks propagate down the material. At $`t=0`$ the whole material is supposed to have a constant humidity concentration $`h_0`$. In this situation the material is unstressed. Humidity concentration is assumed to decrease with time $`h(z,t>0)h_0`$. For any reasonable experimental realization of the drying process occurring from the $`z=0`$ free surface, it is clear that at any time, well inside the material we should reach the original humidity concentration, i.e., $`h(z\mathrm{},t)=h_0`$. Then the majority of the material is always at constant humidity $`h_0`$. In this region the sample must be unstressed, otherwise it will store an infinite amount of energy. Then a boundary condition for our problem will be that stresses go to zero as $`z`$ go to infinity. As a simplification of a possible experimental situation, we will consider the case of an abrupt drying profile (see Fig. 1), namely $`h(z,t)=h_0`$ for $`z>z_0(t)`$, and $`h(z,t)=h_1`$ for $`z<z_0(t)`$, with $`h_1<h_0`$, and $`z_0(t)`$ being an increasing function of time.

A given set of cracks will always correspond to a local minimum of the total energy of the system. The total energy is the sum of two well different parts. One is the fracture energy, namely, the energy spent in the creation of all fractures present in the sample. This is typically proportional to the total length of cracks in the material. The second part is the elastic energy stored in the sample. Let us suppose we have an evenly spaced array of cracks, defining stripes of width $`l`$, which have penetrated down to some distance $`\overline{z}`$, with the drying front being located at some position $`z_0`$.

The elastic energy $`e_{\mathrm{el}}`$ stored in the material per unit of horizontal length must be proportional to the bulk modulus of the material $`C`$, and to the second power of the typical change in linear density caused by the humidity gradient (this change being $`\alpha \mathrm{\Delta }h`$, with $`\mathrm{\Delta }h(h_0h_1)`$). This is at the basis of linear elasticity theory. Using the fact that $`l`$ and $`\overline{z}z_0`$ are the only two relevant lengths for this geometrical configuration, dimensional analysis allows to write down the following expression for $`e_{\mathrm{el}}`$

$$e_{\mathrm{el}}=(\alpha \mathrm{\Delta }h)^2Clg\left((\overline{z}z_0)/l\right),$$
(2)
where $`g(x)`$ is a dimensionless geometrical function. We have determined $`g`$ numerically and the result is shown in Fig. 2 as a continuous line. We can rationalize the general form of the function $`g(x)`$, considering how the elastic energy is distributed along the $`z`$ direction. The total energy $`e_{\mathrm{el}}`$ is the integral over $`z`$ of the density of elastic energy $`\delta e_{\mathrm{el}}(z)`$. In Fig. 3 we show qualitatively the form of $`\delta e_{\mathrm{el}}(z)`$ for the cases $`\overline{z}z_0`$, and $`\overline{z}z_0`$. Let us consider first the case $`\overline{z}z_0`$ (Fig. 3(a)). $`\delta e_{\mathrm{el}}(z)`$ is zero for $`z>z_0`$, since as we already discussed, the material has to be unstressed for $`z\mathrm{}`$. In the region $`z_0>z>\overline{z}`$ there is a rather constant energy density, associated to the change in humidity concentration, which cannot be compensated by a change in density of the material since the material here is attached to the part below $`z_0`$. Around the position of the crack front there is an increase of the stored elastic energy, which is associated to the elastic energy around the tips of the cracks. For $`z<\overline{z}`$ the elastic energy density goes to zero, since here the existence of the cracks has allowed to relieve the elastic energy accumulated previously to the cracks formation. The linear dependence of $`g(x)`$ as $`x\mathrm{}`$ (i.e. for $`\overline{z}z_0`$) comes from the energy stored between $`\overline{z}`$ and $`z_0`$. Let us consider now the case when $`\overline{z}z_0`$ (Fig. 3(b)). The elastic energy becomes independent of $`\overline{z}`$, as the crack tips are in a region of material that is unstressed. The constant value of $`g(x)`$ as $`x+\mathrm{}`$ comes from the energy stored in the independent stripes around the position of the humidity front, which in this limit is well behind the crack front. When $`\overline{z}z_0`$ there is a smooth crossover between the two limiting regimes.

To determine the actual position of the crack front in a realistic situation we will rely on the energetic argument. At any time during the drying process, the crack front will be located at the position that minimizes the total energy of the system. Equation (2) gives the elastic energy of a set of cracks that have penetrated down to the position $`\overline{z}`$. In order to get the total energy, $`e_{\mathrm{el}}`$ has to be added to the energy cost of creating the cracks. This part, when measured per unit of horizontal length, will be called the crack energy $`e_{\mathrm{ck}}`$, and it is simply given in terms of the specific energy fracture of the material $`\eta `$ in the form

$$e_{ck}=\eta \overline{z}/l$$
(3)
In order to determine the value of $`\overline{z}`$ at which the fracture front prefers to be located, we have to minimize the total energy $`e_{\mathrm{tot}}=e_{\mathrm{el}}+e_{\mathrm{ck}}`$. The result we obtain is shown in Fig. 4, where we plot the most convenient position of the crack front as a function of the parameter $`u\eta /[C(\alpha \mathrm{\Delta }h)^2l]`$ . As we see, for very small value of $`u`$, $`x(\overline{z}z_0)/l`$ takes large and positive values, i.e., the crack front is located well below the humidity front. This is due to the negligible contribution of the crack energy compared to the elastic energy in this case. As $`u`$ increases $`x`$ decreases, crossing zero (namely, the crack front coinciding with the humidity front) at $`u0.66`$. $`x`$ tends to $`\mathrm{}`$ as $`u`$ approaches the limiting value $`u_01.14`$. In fact, for $`u>u_0`$ the crack energy is so high that cracks do not penetrate the sample at all.

Under the conditions analyzed in this section, the crack front will be located at a position such that $`(\overline{z}z_0)/l`$ is given by the function plotted in Fig. 4. The crack front advances only due to changes in the position of the humidity front $`z_0`$, keeping always $`\overline{z}z_0`$ as constant.

## IV Stability of the flat crack front

The previous analysis has assumed that all cracks penetrate down to a uniform depth $`\overline{z}`$, and has focused on what the value of $`\overline{z}`$ is, on energetic grounds. It has to be complemented however with a stability analysis of the crack front. The necessity of this is clear from the following example. If a material with an array of vertical cracks is loaded with a uniform horizontal stress (this situation can be thought to be realized in our case if $`\overline{z}z_0`$), there will be typically a single crack that propagates and fractures the material. This is a consequence of the fact that as soon as a single crack moves forward a small distance, the stress at its tip increases, and those at the tips of the other cracks decrease. This generates an unstable, rapid propagation of a single crack. We will see that in our case, this can be compensated by the fact that stresses decrease ahead of the humidity front, and this can stabilize a flat crack front.

The same kind of energetic arguments used in the previous section will be used to determine the stability conditions of the crack front. Consider an evenly spaced set of cracks, labelled sequentially by an index $`j`$, where now the tips of the cracks are at vertical positions $`z^j`$, which can be slightly displaced from the mean position $`\overline{z}`$, i.e., $`z^j=\overline{z}+\delta _ij`$, with $`_j\delta _j=0`$. The horizontal positions of the cracks are given by $`x^j=jl`$. The elastic energy of this configuration contains a term of the form (2), plus a correction that can be expanded in powers of $`\delta _j/l`$. The first order term in this expansion vanishes, as $`_j\delta _j=0`$. The second order term can be written in the form

$$\mathrm{\Delta }e_{\mathrm{el}}=l^2(a\underset{j}{}\delta _j^2+b\underset{j}{}\delta _j\delta _{j+1}+c\underset{j}{}\delta _j\delta _{j+2}+\mathrm{})$$
(4)
Successive term contain ‘interactions’ between more distant cracks. As the elastic energy is a local quantity, we expect $`|a|<|b|<|c|<\mathrm{}`$. Any small displacement of the crack front, defined by the quantities $`\delta _j`$ can be decomposed in a sum of ‘normal modes’, by going to the Fourier representation

$$\delta _j𝑑k\stackrel{~}{\delta }_k\mathrm{exp}(ikx^j),$$
(5)
and the energy decomposes into independent term, in the form

$$\mathrm{\Delta }e_{\mathrm{el}}\underset{k}{}|\stackrel{~}{\delta }_k|^2(a+b\mathrm{cos}(2\pi k/l)+c\mathrm{cos}(4\pi k/l)+\mathrm{}).$$
(6)
In order for the flat crack front to be stable, $`\mathrm{\Delta }e_{\mathrm{el}}`$ must be positive for any choice of the $`\delta _i`$, and thus of the $`\stackrel{~}{\delta }_k`$. This implies that

$$(a+b\mathrm{cos}(2\pi k/l)+c\mathrm{cos}(4\pi k/l)+\mathrm{})$$
(7)
must be positive for all $`k`$. In the limit of very small $`k`$, $`\mathrm{\Delta }e_{\mathrm{el}}`$ is equivalent to a uniform advance of the crack front, and then is has to be a positive quantity, as the curvature of $`e_{\mathrm{el}}(\overline{z})`$ is always positive (see (2) and Fig. 2). On the other hand, It will be proportional to the sum of all coefficients in expression (7), i.e., $`\mathrm{\Delta }e_{\mathrm{el}}^{k0}a+b+c+\mathrm{}`$. Then it is clear that if an instability exists for some value of $`k`$, it will occur at $`k=\pi /l`$, where $`\mathrm{\Delta }e_{\mathrm{el}}^{k=\pi /l}ab+c\mathrm{}`$. Then we will analyze the stability of the flat crack front against a perturbation with $`k=\pi /l`$. We took the equilibrium position of a crack front obtained in the previous section (plotted in Fig. 4) and calculated numerically the quadratic change in energy $`\mathrm{\Delta }e_{\mathrm{el}}=\chi (x)(\alpha \mathrm{\Delta }h)^2C\delta ^2/l`$ when a perturbation of $`k=\pi /l`$ and amplitude $`\delta `$ is introduced. The dimensionless function $`\chi (x)`$ can be seen in Fig. 2 (dashed line). We see that $`\chi `$ is positive (negative) for $`x`$ greater (lower) than $`x_{\mathrm{cr}}0.038`$. In this way we obtain that there are two regions in the plot of Fig. 4. The one at the left of $`u_{\mathrm{cr}}1.02`$ corresponds to a stable situation. If the parameters of the system make the pattern of cracks to lie in that region, then the time evolution of the drying process (i.e., the increase with time of $`z_0`$) will produce a smooth advance of the crack front, keeping always the same value of the distance $`\overline{z}z_0`$ to the humidity front. At the right of $`u_{\mathrm{cr}}`$, the pattern is unstable. Should we have one of those patterns at a given time, it will immediately propagate forward some of its cracks (ideally, one of each two cracks), in order to reach a stable situation. This will imply in particular that some cracks will remain halted. The further evolution of the crack front will correspond to a new crack pattern with less cracks (i.e., with larger $`l`$) being propagated.

## V Appearance and ordering of cracks

Up to now we have assumed that a set of evenly spaced cracks exists, and we focused on its stability conditions. We will study now how this pattern can appear, starting with the process at the surface. Cracks are not expected to appear evenly spaced at the surface. In fact, it is known in one-dimensional models of surface fragmentation  (that can be used to represent the first stage of cracking of our two-dimensional problem) that the distribution of fragment length is strongly dependent on the presence of small inhomogeneities in the material, making the fragment length distribution broad. But we will see that as the superficial cracks penetrate the sample they become evenly spaced. We will first discuss how the process of nucleating new cracks at the surface is, and then argue that the cracks become evenly spaced as they penetrate the sample.

### A The appearance of cracks at the surface

When the drying front penetrates the material from its surface, and before the material gets cracked, the stresses can be very simply calculated. The state of the system corresponds to the material being completely unstrained horizontally (considering the unstrained state as that corresponding to the humidity value $`h_0`$ that occurs well inside the material). Under this condition a uniform horizontal stress $`T_0`$ appears for all $`z<z_0`$, which is simply calculated as $`T_0=C\alpha \mathrm{\Delta }h`$. In order for the first crack to appear, this value has to overcome the uniform traction resistance of the material $`T_r`$ :

$$C\alpha \mathrm{\Delta }h>T_r.$$
(8)
If this relation is not satisfied, no cracks will appear whatsoever. Assuming that the relation (8) is satisfied, the first crack will nucleate at the surface. This crack penetrates the sample as long as this penetration reduces the total energy of the system. In our case, the crack penetrates only down to a distance $`d_0`$ of the order of $`z_0`$, where the humidity front is located. The horizontal stress at the surface is now zero at the position of the crack, and increases as we move away from it, reaching the value $`F_0`$ at large distances of the crack. This means that new cracks will nucleate away of the first one, in regions where relation (8) is still satisfied. The number of cracks nucleated at the surface will be typically the minimum number that makes the horizontal stress at every point of the surface to be lower than $`T_r`$. The typical distance between cracks $`l`$ can then be estimated to be

$$ld_0f(T_0/T_r)$$
(9)
where the geometrical function $`f(T_0/T_r)`$ goes to infinity when $`T_0T_r^+`$, and is of order 1 for $`T_0T_r`$. Typically, $`l`$ is a few times $`d_0`$.

If we stick to the ideal sharp drying front we have been studying, the first crack appears for an infinitesimally small value of $`z_0`$, the only restriction is that the relation (8) is satisfied. Therefore, since the depth of the first crack $`d_0`$ is of the order of $`z_0`$, we should expect a very dense set of superficial cracks. However, this unrealistic situation is removed if we note that any physical drying front will have some typical (finite) humidity gradient $`h`$. At the moment when relation (8) is first satisfied (with $`\mathrm{\Delta }h`$ being now the difference between humidity concentrations right at the surface, and well inside the material), the first crack will penetrate down to a distance

$$d_0\mathrm{\Delta }h/h=T_r/(C\alpha h),$$
(10)
and the typical distance between cracks at the surface will be a few times this distance.

### B How the distribution of cracks gets uniform

The cracks that appear at the surface of the sample nucleate at positions that are strongly influenced by the presence of small inhomogeneities in the sample . But a superficially uneven set of cracks has a tendency to become evenly spaced as it penetrates the sample. The reason for this tendency can be understood on an energetic basis. For a fixed number of cracks the evenly spaced configuration corresponds to the minimum of elastic energy. It is then clear, using the energetic arguments, that this configuration will be approached during the propagation process.

It is important to clarify the different effects that inhomogeneities have on the creation and propagation of cracks. For the superficial layer, inhomogeneities are relevant, and responsible for the broad distribution of fragment lengths, since at the beginning all the surface is uniformly stressed, and tiny differences in the properties of the material will dictate which point of the surface will fail first. However, once the array of cracks has been defined at the surface, inhomogeneities play a secondary role, since new cracks do not appear during penetration, and the elastic energy is only very weakly dependent on the precise distribution of defects. Then the evolution of the pattern is basically the same as if inhomogeneities were absent.

The previous energetic arguments predict the trend towards evenly spaced cracks as they penetrate the sample. Stress analysis allows to re-obtain and make this result quantitative. We have calculated the stresses that are present in the semi-infinite sample with a sharp drying front, in the presence of a particular uneven set of parallel cracks. We took an infinite set of cracks separated sequentially by distances $`l_1`$ and $`l_2`$.

For this geometrical configuration we have determined numerically the direction of maximum opening stress (mode I), at the tips of the fractures, and characterized it by its angle $`\theta `$ with respect to the direction of forward advance. The first result is that $`\theta `$ always points in the direction of making the pattern more evenly spaced. The actual value of $`\theta `$ is absolutely independent of the values of $`C`$ and $`\eta `$ of the material, and also on the value $`\mathrm{\Delta }h`$ of the step of the humidity front. It is in fact a quantity that only depends on the Poisson ratio $`\nu `$ of the material and the values of $`x(\overline{z}z_0)/l`$ and $`(l_1l_2)/(l_1+l_2)`$. The results for the angle $`\theta `$, for a material with $`\nu =1/3`$ can be summarized as follows: $`\theta `$ is rather independent of $`x`$, and within a maximum 10 % error, it can be written as

$$\theta (\mathrm{deg})29(l_1l_2)/(l_1+l_2)$$
(11)
The value of $`\theta `$ vanishes for $`l_1=l_2`$, as in this case the pattern is actually evenly spaced, and the cracks advance straightforwardly by symmetry. The limiting case $`l_1l_2`$ gives $`\theta 29`$ deg. This is the maximum bending we can expect from a given crack, in the process of uniformising the widths of the stripes. Note that as a consequence of the fact that there are no typical lengths in (11), the length needed for the pattern to become evenly spaced will be proportional to the typical width of the stripes $`l`$.

Once the cracks have started to deviate in order to reduce the elastic energy as much as possible, the further detailed prediction of its evolution becomes more difficult, since now we should calculate stresses ahead of a set of non parallel fractures (or alternatively, the elastic energy of this nontrivial configurations of curved cracks). However, the main conclusion that the pattern becomes evenly spaced in a distance of the order of the stripes widths remains.

## VI Dependences on the characteristics of the drying front

The sharp drying front is an idealization that is never exactly realized in practice. For the experiments of Ref. , for instance, a typical distance over which humidity changes is expected. In other cases in which the drying is through surface diffusion, the humidity profiles are still smoother. As we have already discussed the stability reason of a flat crack front in the case of a sharp humidity step is due to the rapid reduction of stresses ahead of the crack front. Smoother humidity profiles will produce a weaker tendency to generate a flat crack front, and it can even occur that a flat crack front is never stable if the humidity profile is smooth enough. In this case a pattern of bifurcations of the crack front has been predicted. Note that an additional stabilizing factor of the crack front exists when the drying is favored by the presence of the cracks themselves, a case we have not addressed here.

## VII Implications for the three dimensional case

The same phenomenon we have been discussing acquires novel characteristics in a three dimensional geometry. It has a beautiful realization in the geological formations named columnar basalts, which can be reproduced in a kitchen experiment using corn-starch. When an originally hot volcanic lava flow starts to cool down (after solidification), the thermal stresses generate cracks at the surface, which propagate progressively towards the interior of the igneous body. There are many coincidences with the two dimensional case, and also some new features . Cracks are known to appear at the surface and at later times propagate to the interior. At the surface the pattern of cracks is rather disordered, but becomes progressively ordered as it penetrates the material. In three dimensions, the tendency to order manifest in progressive tendency to form a polygonal pattern of cracks, in planes parallel to the surface. Qualitatively this tendency to order can be understood on the same basis as the two dimensional case, as a tendency of the system to reduce its total energy as much as possible during the evolution. It poses however some interesting problems, because typically, the expected perfect hexagonal pattern is not reached, but a collection of polygons with different number of sides and areas. It has been recently demonstrated that this is a consequence of the fact that the patterns starts being disordered at the surface, and in the process of minimizing the energy it is not able to reach the absolute minimum (namely, the honeycomb pattern) but is trapped in a metastable minimum. Based in the arguments for the two dimensional case, we can expect also in three dimensions that the typical width of the columns is set by the typical distance between cracks at the surface, and it is not modified with the further penetration. In fact this is what is experimentally observed, since columns are seen to keep its horizontal size over distances that reach hundred times the width of the columns. For the three dimensional case, the cooling through diffusion plus convection (and/or radiation) at the surface is the most realistic model of cooling to consider. Although this type of cooling can lead in two dimensions to the instability of the flat crack front, this is not necessarily so in three dimensions, since in this case all fractures form a connected structure, and this generates an additional tendency to stabilize the flat crack front. These issues and some additional ones particular of the 3D case will be discussed in a forthcoming publication.

## VIII Summary and conclusions

In this paper we have focused on a two dimensional material that cools down or dessicates from one edge. We studied the conditions under which a set of cracks penetrate the sample during the process. We showed that these cracks have a tendency to become evenly spaced, an analyze the conditions for which all the tips of the cracks form a plane front. Our results provide a quantitative framework to analyze recent experiments. Extensions to the three dimensional case are expected to provide both a qualitatively similar scenario, and also some new features that will be fully developed elsewhere.

## IX Acknowledgments

This work was financially supported by CONICET (Argentina). Partial support from Fundación Antorchas is also acknowledged.


format: text

text:
# Optimal design, robustness, and risk aversion

(February 19, 2002)

## Abstract

Highly optimized tolerance is a model of optimization in engineered systems, which gives rise to power-law distributions of failure events in such systems. The archetypal example is the highly optimized forest fire model. Here we give an analytic solution for this model which explains the origin of the power laws. We also generalize the model to incorporate risk aversion, which results in truncation of the tails of the power law so that the probability of disastrously large events is dramatically lowered, giving the system more robustness.

optimization, power laws, designed systems, self-organization, robustness

In a series of recent papers, Carlson and Doyle CD99 ; CD00 ; DC00  have proposed a model for designed systems which they call “highly optimized tolerance” or HOT. The fundamental idea behind HOT is that systems designed for high performance naturally organize into highly structured, statistically unlikely states that are robust to perturbations they were designed to handle, yet fragile to rare perturbations and design flaws. As an example they consider an idealized model of forest fires CD00 . In this model a forester is charged with finding the optimal distribution of the trees on a grid so as to maximize tree harvest in the face of occasional fires that burn complete connected clusters of trees and are started by sparks that arrive with a given spatial distribution. They find that optimizing the harvest, or yield, for the model gives rise to a segmented forest consisting of contiguous patches of trees separated by firebreaks, and that the resulting distribution of fire sizes usually follows a power law. While this type of configuration does typically achieve very good yields, the system is also fragile in the sense that perturbations to the firebreaks or changes in the spark distribution can lead to substantially sub-optimal performance. They argue that these are pervasive phenomena: high-performance engineering leads to systems that are robust to stresses for which they were designed but fragile to errors or unforeseen events.

In this paper we argue that simple yield maximization is problematic even if there are no errors in firebreaks or changes in the spark distribution. Because the power-law distributions generated by yield maximization have fat tails, disastrously large forest fires occur with non-negligible frequency—far greater frequency than one would expect from intuition based on normal distributions. This idea, that yield optimization can lead to ruinous outcomes, is not new. For the classic problem of gambler’s ruin, for example, it is well known that optimizing total return leads to ruin with probability one. By contrast, if one is willing to accept suboptimal returns it is possible to construct gambling strategies that are immune to ruin CT91 . Applying similar ideas in the present context, we show that a risk-averse engineer who is willing to accept some loss in average system performance can effectively limit the large deviations in the event size distribution so that disasters are rare. We call this variation on the HOT theme “constrained optimization with limited deviations”, or COLD. By avoiding total ruin, a COLD design is more robust than a HOT one, even in a world of perfect error-free optimization.

To demonstrate the difference between HOT and COLD we first revisit the HOT forest fire model. We give an analytic solution for the model which shows that the distribution of fire sizes does indeed follow a power law, the cumulative distribution having a exponent $`(1+1/d)`$, where $`d`$ is the dimensionality of the system. Using a fast percolation algorithm NZ00  we perform numerical simulations that confirm the value of this exponent. We then generalize our solution to include risk aversion in the design, thereby breaking the power law scaling and dramatically reducing the frequency of disastrously large events.

Following Refs. CD99,  and CD00,  then, we consider a forest divided into a large number of regions or patches, with firebreaks between them that prevent the spread of fire from one patch to another. Although the original forest fire model was based on a lattice, the model we consider is a continuum one, since this makes the mathematical treatment more tractable. For large system sizes, we expect the behavior of this continuum model to converge to that of the lattice model.

It is assumed that during the lifetime of the forest a single spark lands at a random position $`𝐫`$ and starts a fire that burns the surrounding patch. Let us denote the area of this patch by $`s(𝐫)`$. The forest is then harvested, giving a yield equal to the area of the remaining forest. In units where the total area of the forest is one, the yield is $`1s(𝐫)F`$, where $`F`$ is the cost in terms of yield of constructing the firebreaks.

Because dimensionality is an important property of HOT systems, we consider the model for general dimension $`d`$. If the cost of constructing firebreaks is $`a`$ per unit length (or per unit surface area for $`d>2`$), then the cost of the firebreak surrounding a patch $`m`$ is $`agds_m^{(d1)/d}`$, where $`s_m`$ is the value of $`s(𝐫)`$ in patch $`m`$ and $`g`$ is a geometric factor of order 1 that depends on the geometry of the lattice and the shape of the patches. In the lattice version of the forest fire model, $`a`$ is simply equal to the lattice parameter (i.e., the nearest-neighbor spacing), but in the continuum model we are at liberty to give $`a`$ any value we feel to be appropriate. As we will see, as long as $`a`$ is finite its value does not affect the shape of the distribution of fire sizes.

Because $`s(𝐫)`$ is constant inside each patch, the integral of $`1/s(𝐫)`$ over any patch is identically 1, and hence, summing over all patches, the total area occupied by firebreaks is

$$F=agd\underset{m}{}s_m^{(d1)/d}=agd\underset{m}{}s_m^{(d1)/d}_m\frac{\mathrm{d}^dr}{s(𝐫)}=agds(𝐫)^{1/d}\mathrm{d}^dr.$$
(1)
Letting the normalized probability distribution of sparks be $`p(𝐫)`$, the mean yield is then

$$Y=1p(𝐫)s(𝐫)\mathrm{d}^dragds(𝐫)^{1/d}\mathrm{d}^dr,$$
(2)
where the integrals run over the entire area of the forest.

To find the maximum yield with respect to the patch sizes $`s(𝐫)`$, we set the functional derivative $`\delta Y/\delta s(𝐫)=0`$, giving

$$s(𝐫)=\left[\frac{ag}{p(𝐫)}\right]^{d/(d+1)},$$
(3)
for all $`𝐫`$. The optimal yield is then given by substituting back into Eq. (2) to get

$$Y_{\mathrm{opt}}=1(d+1)(ag)^{d/(d+1)}p(𝐫)^{1/(d+1)}\mathrm{d}^dr,$$
(4)
and the optimal number of patches is

$$n=\frac{\mathrm{d}^dr}{s(𝐫)}=(ag)^{d/(d+1)}p(𝐫)^{d/(d+1)}\mathrm{d}^dr.$$
(5)
For the lattice version of the model, $`a`$ goes as $`L^1`$, where $`L`$ is the (linear) system size, and hence the number of patches should scale as $`L^{d/(d+1)}`$, i.e., as $`L^{1/2}`$ in one dimension or $`L^{2/3}`$ in two. Numerical experiments on a one-dimensional system confirm this, giving a measured exponent of $`0.47\pm 0.03`$.

Now we wish to calculate the distribution $`\rho (s)`$ of fire sizes that arises if we make this choice of patch sizes. We have

$$\rho (s)=p(𝐫)\frac{\mathrm{d}^dr}{\mathrm{d}s}=p(𝐫)\frac{\mathrm{d}^dr}{\mathrm{d}p}\frac{\mathrm{d}p}{\mathrm{d}s}=ag\frac{d+1}{d}p(𝐫)\frac{\mathrm{d}^dr}{\mathrm{d}p}s^{(2+1/d)},$$
(6)
where $`\mathrm{d}^dr`$ here represents the volume of the space between the contours $`p`$ and $`p+\mathrm{d}p`$ on the $`p(𝐫)`$ surface. As we will show, the term $`p(𝐫)\mathrm{d}^dr/\mathrm{d}p`$ is constant or contributes logarithmic corrections for a wide selection of possible distributions $`p(𝐫)`$, while the principal power-law behavior in the event size distribution comes from the factor $`s^{(2+1/d)}`$111As noted in Ref. DC00, , the exponent of this power law *decreases* with increasing dimension $`d`$, which is the opposite of the behavior seen in other models such as the self-organized critical forest fire model of B. Drossel and F. Schwabl, Phys. Rev. Lett. 69, 1629–1632 (1992)..

An alternative method for deriving Eq. (6) is to maximize the simple yield functional $`Y=1p(𝐫)s(𝐫)\mathrm{d}^dr`$, subject to a constraint that fixes the volume $`F`$ occupied by the firebreaks (Eq. (1)). This method, which is similar to the approach taken in Ref. CD99, , is equivalent to the method above, via a Lagrange transform, provided $`F`$ is chosen so as to make the corresponding Lagrange multiplier equal to $`agd`$.

Consider then the case (which covers all the examples in Refs. CD99,  and CD00, ) of a distribution of sparks with a single maximum at the origin, so that the volume $`\mathrm{d}^dr`$ takes the form of an annulus enclosing the origin. If we denote by $`\mathrm{\Omega }_d`$ a $`d`$-dimensional solid angle centered on the origin, then a volume element of the annulus is $`r^{d1}\mathrm{d}r\mathrm{d}\mathrm{\Omega }_d`$. In terms of the thickness $`\mathrm{d}p`$ of the annulus, we can write $`\mathrm{d}r=(r\mathrm{d}p)/(𝐫p)`$, and integrating our volume element over the contour $`p=\text{constant}`$ we find

$$\frac{\mathrm{d}^dr}{\mathrm{d}p}=_p\frac{r^d\mathrm{d}\mathrm{\Omega }_d}{𝐫p}.$$
(7)
For example Carlson and Doyle CD99  studied the case of a spark distribution in two dimensions having the form of the product of two Gaussians with different widths:

$$p(𝐫)=N\mathrm{exp}\left(\left[\frac{x^2}{2\sigma _x^2}+\frac{y^2}{2\sigma _y^2}\right]\right),$$
(8)
where $`N`$ is a normalization constant. For this distribution the denominator of the integrand in Eq. (7) is

$`𝐫p`$ $`=`$ $`N\left[{\displaystyle \frac{x^2}{\sigma _x^2}}+{\displaystyle \frac{y^2}{\sigma _y^2}}\right]\mathrm{exp}\left(\left[{\displaystyle \frac{x^2}{2\sigma _x^2}}+{\displaystyle \frac{y^2}{2\sigma _y^2}}\right]\right)`$ (9)
$`=`$ $`2p\mathrm{log}{\displaystyle \frac{p}{N}},`$
which is constant over our contour of constant $`p`$. The element of solid angle in two dimensions is simply the element of polar angle $`\mathrm{d}\theta `$, and hence Eq. (7) simplifies in this case to

$$\frac{\mathrm{d}^2r}{\mathrm{d}p}=\frac{1}{2p\mathrm{log}(p/N)}_pr^2d\theta =\frac{A(p)}{p\mathrm{log}(p/N)},$$
(10)
where $`A(p)`$ is the area enclosed by the contour. This contour is a line of constant $`(x/\sigma _x)^2+(y/\sigma _y)^2`$, i.e., an ellipse, which has major and minor axes $`a=\sqrt{2\sigma _x^2\mathrm{log}(N/p)}`$ and $`b=\sqrt{2\sigma _y^2\mathrm{log}(N/p)}`$. Thus the area enclosed by the contour is $`A(p)=\pi ab=2\pi \sigma _x\sigma _y\mathrm{log}(N/p)`$. Combining Eqs. (6) and (10) we then find that the distribution of event sizes is

$$\rho (s)=3\pi \sigma _x\sigma _yags^{5/2}.$$
(11)
Thus, for the Gaussian case in two dimensions the model generates a perfect power-law with slope $`\frac{5}{2}`$.

This argument is easily generalized to other spark distributions and other dimensions. We find that the HOT forest fire model generates a perfect power law with slope $`(2+1/d)`$ for all dimensions when a spark distribution of the form $`p(𝐫)=N\mathrm{exp}(_{i=1}^d[x_i/\sigma _i]^d)`$ is used. As a test of this prediction, we show in Fig. 1 numerical results from direct simulations of the forest fire model in one and two dimensions with distributions of this type 222The simulations were carried out using the greedy algorithm of Ref. CD99, , in which trees are placed on the lattice one by one, each in the position that maximizes the yield at that step. This is not a true global optimization, but appears to give results close to the global optimum and in good agreement with theory.. For better visualization and analysis the distributions pictured are cumulative, so the expected slope is $`(1+1/d)`$, rather than $`(2+1/d)`$. As the figure shows, the slopes of the observed distributions are in good agreement with this prediction.

We note in passing that the slope of $`(1+1/d)`$ for the cumulative distribution of fire sizes seen in both our exact solution and our numerical results is different from the slope of approximately $`1`$ found numerically by Carlson and Doyle CD00  in two dimensions. The source of this discrepancy is unclear, although it may be that the simulations of Ref. CD00,  provided too few data points to make an accurate evaluation of the exponent possible. We note also that the value $`\frac{3}{2}`$ for the two-dimensional case is quite different from the slope of $`\frac{1}{2}`$ measured for the cumulative size distribution of real forest fires DC00 ; MMT98 .

Other functions with exponential tails also generate power laws, but give logarithmic corrections as well. For instance, if $`p(𝐫)=N\mathrm{exp}(_{i=1}^d[x_i/\sigma _i]^\gamma )`$ with $`\gamma d`$ then Eq. (10) still applies, but now $`A(p)[\mathrm{log}(N/p)]^{d/\gamma }`$ and hence

$$p(𝐫)\frac{\mathrm{d}^dr}{\mathrm{d}p}[\mathrm{log}(N/p)]^{d/\gamma 1}.$$
(12)
Thus the distribution of event sizes fundamentally still follows a power law with slope $`(2+1/d)`$, but there is a logarithmic correction. Similar logarithmic corrections are noted in Ref. CD99, .

Spark distributions with power-law tails also (unsurprisingly) give power-law event size distributions, but in this case the exponent of the distribution is non-universal, varying with the exponent of the spark distribution. For example, if $`p(𝐫)`$ takes the generalized Lorentzian form

$$p(𝐫)=\frac{N}{_i(x_i/\sigma _i)^\nu +\mathrm{\Gamma }^\nu },$$
(13)
then we find

$$\frac{\mathrm{d}^dr}{\mathrm{d}p}\frac{1}{p(𝐫)}\left(\frac{N}{p}\mathrm{\Gamma }^\nu \right)^{d/\nu 1},$$
(14)
which goes asymptotically as $`p^{d/\nu }`$ in the power-law tail where $`p(𝐫)`$ becomes small. Thus the tail of the distribution of event sizes goes as $`\rho (s)s^{(2+1/dd/\nu )}`$.

We now turn to the COLD variant of the forest fire model, which incorporates risk aversion. In constructing this model we are guided by theories of risk aversion in economics, where the subjective benefit of outcomes is typically a nonlinear function of the loss $`s`$, which is captured by a utility function $`u(s)`$ Savage72 . Sensible utility functions are decreasing with increasing loss: $`u^{}<0`$ 333In traditional economics utility functions are usually defined to be increasing functions of their arguments; this could arranged in the present case, if necessary, by substituting $`ss`$.. Risk aversion also implies that $`u^{\prime \prime }<0`$, so that the negative utility of bad outcomes is weighted more strongly than the positive utility of good outcomes. One standard family of utility curves that achieves this is the one-parameter family

$$u(s)=\frac{(1s)^\alpha }{\alpha }.$$
(15)
Note that, since we will be concerned only with maximizing utility, $`u(s)`$ is arbitrary to within both additive and multiplicative constants. For $`\alpha =1`$, Eq. (15) gives $`u=1s`$ and maximizing utility is precisely equivalent to minimizing loss. For $`\alpha <1`$, we have risk-averse utility functions and for $`\alpha <0`$ we are infinitely averse to losing our entire investment.

Our goal now is to maximize the average utility functional

$$U=p(𝐫)u(s(𝐫))\mathrm{d}^dr,$$
(16)
subject to the constraint of fixed $`F`$, Eq. (1) (or equivalently maximize a combined utility functional similar to Eq. (2)). Carrying out the functional derivatives and using the utility function of Eq. (15), we find that the optimum $`U`$ corresponds to

$$p(𝐫)s(𝐫)^{(d+1)/d}[1s(𝐫)]^{\alpha 1}=\lambda ,$$
(17)
where $`\lambda `$ is a Lagrange multiplier whose value can be calculated from Eq. (1). The distribution of event sizes is given by Eq. (6) as before, and using Eq. (17) we find that the derivative $`\mathrm{d}p/\mathrm{d}s`$, which gives the principal variation in $`\rho (s)`$, is

$$\frac{\mathrm{d}p}{\mathrm{d}s}=\lambda \frac{(\alpha +1/d)s(1+1/d)}{(1s)^\alpha s^{2+1/d}}.$$
(18)
For $`\alpha =1`$ our utility maximization is equivalent to simple yield maximization (HOT), so it is not surprising to observe that when we set $`\alpha =1`$ in the above expression we recover our previous $`s^{(2+1/d)}`$ power-law. For $`\alpha <1`$, we have risk-averse utility functions (COLD), which give rise to event distributions following the $`s^{(2+1/d)}`$ form for small event sizes, but having lower probability of large event sizes. When $`\alpha <0`$, event probability tends to zero as $`s1`$, as we would expect.

In Fig. 2 we compare the distribution of event sizes in HOT and COLD regimes for a variety of values of the risk-aversion parameter $`\alpha `$. The figure shows that the COLD distribution approaches the HOT one as $`\alpha `$ approaches 1. For $`\alpha `$ large and negative the HOT power law is followed for only a small portion of the range of event sizes—about 20% in the case of $`\alpha =5`$.

It is worth noting that while risk aversion truncates the power-law behavior in the event size distribution, the distribution of the *utilities* of events still follows a power law: we find that the tail of the distribution of utilities goes as $`\rho (u)u^\beta `$ with $`\beta =(2\alpha 1)/\alpha `$. Note that this exponent is independent of the system dimension $`d`$ 444A similar result can also be derived within the framework of the “PLR” model of Ref. DC00, , which is a general model of engineered systems subject to stresses that occur with specified probability. With this model a risk-averse utility function of the form (15) again produces a sharp cutoff in the event distribution, but gives a distribution of the utilities of events that has a dimension-independent power-law tail..

While the introduction of the utility function has reduced the risk of large losses, the optimal utility solution does not normally coincide with the optimal yield solution, and hence we pay a cost for risk aversion in terms of yield. For the lattice forest fire model however we find that the cost paid is small 555For other models or parameter regimes, such as the continuum version of the forest fire model with large $`a`$ or the models of Ref. DC00, , the cost paid could certainly be larger.. For example, in a $`128\times 128`$ two-dimensional system with a Gaussian spark distribution, we find numerically that the mean yield at the $`\alpha =1`$ optimum (HOT) is $`0.904`$, dropping to $`0.900`$ for $`\alpha =3`$ and $`0.888`$ for $`\alpha =5`$. It appears therefore that the introduction of risk aversion garners substantial benefits in terms of the reduction of large losses—and the complete elimination of 100% losses—while at the same time costing us only a few percent at most in terms of average system yield.

We conjecture that the suppression of power law tails in the COLD event size distribution will also make the system more robust against the other problems mentioned in the introduction, namely errors in the design and changes in the spark distribution. The truncation of the power law means that the largest patches in the COLD solution are considerably smaller than those in the HOT solution. Thus, if a design flaw, such as a gap in one of the firebreaks, causes two patches to merge, the resulting combined patch is smaller too. Similarly, if the spark distribution is changed, the size of the resulting fires is smaller, and hence the effect on the average yield is not as catastrophic as in the HOT case. (A similar conjecture is also made in Ref. CD99, .)

In this paper we have examined in detail only the forest fire model, but similar principles should apply to other problems as well. We have shown that in order to produce power-law event size distribution, the HOT model requires the auxilliary assumption of risk-neutrality. If humans are risk-averse they will tend to prefer COLD designs, although this does not necessarily mean that HOT designs never occur. It might be, for instance, that blind evolutionary processes of the type found in natural systems would simply optimize yield, without risk aversion. On the other hand, COLD designs are more robust to rare events than HOT designs, and therefore might be selected for on long time-scales. Of course, in the real world, imperfect designs that fail to optimize either yield or utility are always a possibility too.

The authors thank Jean Carlson and John Doyle for many useful comments and suggestions. This work was supported in part by the National Science Foundation under grant number DMS–0109086 (MEJN and MG) and by McKinsey Corporation, Credit Suisse First Boston, Bob Maxfield, and Bill Miller (JDF).


text:
# Microcanonical mean-field thermodynamics of self-gravitating and rotating systems

## Abstract

We derive the global phase diagram of a self-gravitating $`N`$-body system enclosed in a finite three-dimensional spherical volume $`V`$ as a function of total energy and angular momentum, employing a microcanonical mean-field approach. At low angular momenta (i.e. for slowly rotating systems) the known collapse from a gas cloud to a single dense cluster is recovered. At high angular momenta, instead, rotational symmetry can be spontaneously broken and rotationally asymmetric structures (double clusters) appear.

The statistical equilibrium properties of systems of particles interacting via long-range forces (the so-called non-extensive systems) are currently the subject of intense research, both for their highly non-trivial thermodynamics (displaying such features as negative heat capacities thirring70 ; gross174 ) and for the considerable conceptual and technical difficulties they present. It is known that the long-range nature of the potential makes the canonical ensemble inadequate for describing their statics landau96 ; gallavotti99 ; padmanabhan90 , because the usual thermodynamic limit, where (number of particles) $`N\mathrm{}`$ and (volume) $`V\mathrm{}`$ while intensive variables are kept fixed, does not exist. A central issue is hence whether phase transitions and other conventional statistical phenomena are possible in non-extensive systems gross186 .

Among non-extensive systems, self-gravitating gases, i.e. systems of classical particles subject to mutual gravitation, have deserved the most attention. Their usual static description is based on the microcanonical ensemble padmanabhan90 . In this framework, the key problem is finding the most probable equilibrium configuration of a self-gravitating gas enclosed in a finite $`3`$-dimensional box of volume $`V`$ as a function of the conserved quantities (integrals of motion), the simplest (but possibly not the only relevant ones) being the total energy $`E`$ and the total angular momentum $`𝑳`$.

Dynamical methods chandrasekhar69  based on fluid-mechanics techniques suggest (see e.g. lyndenbell66 ) that upon increasing the ratio between rotational and gravitational energy, the stationary distribution can change from a single dense cluster to a double cluster, and that other structures such as disks and rings might appear.

On the other hand, so far static theories could not recover the richness of the dynamical picture. Taking the total energy as the only control parameter (see padmanabhan90  for a review and chavanis01 ; pettini01 ; chavanis02  for more recent work and references) after removing the rotational symmetry artificially e.g. by constraining the system into a non-spherical box, a “collapse” transition has been found thirring70 , where, as the energy (temperature) is lowered, the equilibrium configuration changes from a homogeneous cloud to a dense cluster lying in an almost void background, with an intermediate “transition” regime characterized by negative specific heat. Despite some attempts gross181 , a detailed static theory embodying angular momentum is lacking.

In this work, building essentially on lyndenbell66 ; laliena98 , we aim at bridging the gap between the static and the dynamical approaches by analyzing the statics of a self-gravitating and rotating gas in a microcanonical setting including angular momentum. Within a mean-field approximation, we derive an algebraic integral equation for the density profiles maximizing the microcanonical entropy and solve it numerically as a function of $`E`$ and $`L`$. Along with the usual collapse, occurring at low $`L`$, we find that for sufficiently high $`L`$ the rotational symmetry of the Hamiltonian can be spontaneously broken, giving rise to more complex equilibrium distributions, including double clusters, rings and disks. The global phase diagram of the system is presented. We shall concentrate here on the equilibrium density profiles, deferring a detailed discussion of the related (highly non-trivial) thermodynamic picture to a more extensive report.

We consider the $`N`$-body system with Hamiltonian

$$H_NH_N(\{𝒓_i\},\{𝒑_i\})=\frac{1}{2m}\underset{i=1}{\overset{N}{}}p_i^2+\mathrm{\Phi }(\{𝒓_i\})$$
(1)
with $`\mathrm{\Phi }(\{𝒓_i\})=Gm^2_{i<j}|𝒓_i𝒓_j|^1`$. $`𝒓_i`$, $`𝒑_i`$ and $`m`$ denote, respectively, the position, momentum and mass of the $`i`$-th particle. The system is assumed to be enclosed in a spherical volume $`V`$ (to preserve rotational symmetry and ensure angular momentum conservation). The crucial quantity to be evaluated is the microcanonical “partition sum”

$$W_N(E,𝑳)=\frac{ϵ}{N!}\delta (H_NE)\delta (𝑳\underset{i=1}{\overset{N}{}}𝒓_i\times 𝒑_i)D𝒓D𝒑$$
(2)
where $`D𝒓=_{i=1}^Nd𝒓_i`$, $`D𝒑=_{i=1}^N(d𝒑_i/h^3)`$, and $`ϵ`$ is a constant that makes $`W_N`$ dimensionless. According to Boltzmann, the entropy is given by

$$S_N(E,𝑳)=\mathrm{ln}W_N(E,𝑳)$$
(3)
Our aim is to find the density profiles that maximize $`S_N`$.

Following Laliena laliena98 , the integration over momenta in (2) can be carried out using a Laplace transform. This yields

$$W_N(E,𝑳)=\frac{A}{N!}[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }(\{𝒓_i\})]^{\frac{3N5}{2}}D𝒓$$
(4)
where $`A`$ is a constant and $`𝕀𝕀(\{𝒓_i\})`$ is the inertia tensor, with elements $`I_{ab}(\{𝒓_i\})=m_{i=1}^N(r_i^2\delta _{ab}r_{i,a}r_{i,b})`$ ($`a,b=1,2,3`$).

In order to evaluate the integral over $`V^N`$, we use the following mean field approximation. Letting $`\rho (𝒓)`$ denote the particles’ density inside $`V`$ ($`\rho (𝒓)𝑑𝒓=N`$), we set

$$\mathrm{\Phi }(\{𝒓_i\})\mathrm{\Phi }[\rho ]=\frac{Gm^2}{2}\frac{\rho (𝒓)\rho (𝒓^{})}{|𝒓𝒓^{}|}𝑑𝒓𝑑𝒓^{\mathbf{}}$$
(5)
$$I_{ab}(\{𝒓_i\})I_{ab}[\rho ]=m\rho (𝒓)\left(r^2\delta _{ab}r_ar_b\right)𝑑𝒓$$
(6)
so that (4) can be re-cast in the form of the functional-integral

$$W_N(E,𝑳)=\frac{A}{N!}[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[\rho ]]^{\frac{3N5}{2}}P[\rho ]𝑑\rho (𝒓)$$
(7)
where $`P[\rho ]`$ is the probability to observe a density profile $`\rho \rho (𝒓)`$. To estimate $`P[\rho ]`$, we adopt the logic of Lynden-Bell lyndenbell66 . We subdivide the spherical volume $`V`$ into $`K`$ identical cells labeled by the position of their centers. The idea is to replace the integral over $`V`$ with a sum over the cells. In order to avoid overlapping we assume that each cell may host up to $`n_0`$ particles ($`1n_0N`$). This condition is essentially equivalent to introducing a hard-core for each particle, and projects out all the physics that is expected to play a role at short distances. $`P[\rho ]`$ is now proportional to the number of ways in which our $`N`$ particles can be distributed inside the $`K`$ cells with maximal capacity $`n_0`$. Denoting by $`n(𝒓_k)`$ the number of particles located inside the $`k`$-th cell, a simple combinatorial reasoning lyndenbell66  leads to

$$P[\rho ]N!\underset{\text{cells }k}{}\left(\genfrac{}{}{0pt}{}{n_0}{n(𝒓_k)}\right)$$
(8)
Defining the relative cell occupancy $`c(𝒓)=n(𝒓)/n_0`$, so that $`\rho (𝒓)=Kn_0c(𝒓)/V`$, and applying Stirling’s formula to approximate the factorials (assuming $`n(𝒓_k)1`$ and $`n_0n(𝒓_k)1`$), one obtains

$$\begin{array}{c}P[c]N!e^{\frac{n_0K}{V}{\scriptscriptstyle \left[c(𝒓)\mathrm{log}c(𝒓)+(1c(𝒓))\mathrm{log}(1c(𝒓))\right]𝑑𝒓}}=\hfill \\ \hfill =N!e^{\frac{N}{\mathrm{\Theta }}{\scriptscriptstyle \left[c(𝒙)\mathrm{log}c(𝒙)+(1c(𝒙))\mathrm{log}(1c(𝒙))\right]𝑑𝒙}}\end{array}$$
(9)
where we introduced the dimensionless variable $`x=r/R`$ and defined $`\mathrm{\Theta }=\frac{NV}{n_0KR^3}=c(𝒙)𝑑𝒙`$.

Plugging this expression into (7) one easily arrives at

$$W_N(E,𝑳)e^{S_N(E,𝑳)}e^{N\left(s_1[c]+s_2[c]\right)}𝑑c(𝒙)$$
(10)
where $`s_1`$ and $`s_2`$ are given by

$$s_1[c]=\frac{3}{2}\mathrm{log}[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[c]]$$
(11)
$$s_2[c]=\frac{1}{\mathrm{\Theta }}[c(𝒙)\mathrm{log}c(𝒙)+(1c(𝒙))\mathrm{log}(1c(𝒙))]𝑑𝒙$$

Notice that $`𝕀𝕀[c]`$. We have neglected terms appearing in the exponent which do not scale with $`N`$.

For large $`N`$ one can resort to the saddle-point method to evaluate (10). Variation of the entropy $`S_N`$ with respect to the relative cell occupancy $`c`$, with the constraint on $`\mathrm{\Theta }`$ enforced by a Lagrange multiplier $`\mu `$, leads to the stationarity condition

$$\mathrm{log}\frac{c(𝒙)}{1c(𝒙)}=\frac{\beta }{\mathrm{\Theta }}U(𝒙)+\frac{1}{2}\beta (𝝎\times 𝒙)^2\mu $$
(12)
or, equivalently,

$$c(𝒙)=(1+e^{\frac{\beta }{\mathrm{\Theta }}U(𝒙)\frac{1}{2}\beta (𝝎\times 𝒙)^2+\mu })^1$$
(13)
where $`𝝎`$ is the angular velocity (related to the total angular momentum by the relation $`𝑳=𝕀𝝎`$), and the shorthands $`\beta `$ and $`U(𝒙)`$ are respectively defined as

$$\beta =\frac{3/2}{[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[c]]}$$
(14)
$$U(𝒙)=\frac{c(𝒙^{\mathbf{}})}{|𝒙𝒙^{\mathbf{}}|}𝑑𝒙^{\mathbf{}}$$
(15)
The essence of our mean-field approach becomes clear if we notice that $`\mathrm{\Phi }[c]c(𝒙)U(𝒙)𝑑𝒙`$.

Eqs (12) or (13) are our central result. They hold for any long-range potential $`\mathrm{\Phi }`$ for which our mean-field approximation can be justified. Functions $`c(𝒙)`$ solving them and corresponding to entropy maxima represent our desired equilibrium distribution of particles. Of course, for fixed energy and angular momentum there may exist different solutions, each having its entropy. In such cases, the higher the entropy the more probable the solution.

Technically, one can only hope to solve e.g. (13) by numerical integration. However, the implicit dependence of $`U(𝒙)`$ on $`c(𝒙)`$ via a 3-dim. integral makes this a formidable task. To simplify things, we pass to spherical coordinates, $`𝒙=(x,\theta ,\varphi )`$, and expand our potential and relative occupancy in series of real spherical harmonics:

$$\frac{1}{|𝒙𝒙^{\mathbf{}}|}=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}\frac{4\pi }{2l+1}\frac{(xx^{})^l}{(xx^{})^{l+1}}Y_{lm}(\theta ,\varphi )Y_{lm}(\theta ^{},\varphi ^{})$$

$$c(𝒙)=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}b_{lm}(x)Y_{lm}(\theta ,\varphi )$$
(16)
with $`xx^{}=\mathrm{min}\{x,x^{}\}`$ and $`xx^{}=\mathrm{max}\{x,x^{}\}`$. $`b_{lm}(x)`$ is a radial function whose precise form we shall soon derive. Using the series (16), together with the completeness relation for our basis set $`\{Y_{lm}\}`$, one can easily show that $`U(𝒙)=_{l,m}u_{lm}(x)Y_{lm}(\theta ,\varphi )`$, with

$$u_{lm}(x)=\frac{4\pi }{2l+1}\frac{(xx^{})^l}{(xx^{})^{l+1}}b_{lm}(x^{})(x^{})^2𝑑x^{}$$
(17)
Multiplying both sides of (13) by $`Y_{lm}`$ and integrating over angular variables one finds

$$b_{lm}(x)=Y_{lm}(\theta ,\varphi )\left[1+e^{\frac{\beta }{\mathrm{\Theta }}_{l=0}^{\mathrm{}}_{m=l}^lu_{lm}(x)Y_{lm}(\theta ,\varphi )\frac{1}{2}\beta \omega ^2x^2\mathrm{sin}^2\theta +\mu }\right]^1\mathrm{sin}\theta d\theta d\varphi $$
(18)
where $`l=0,1,\mathrm{}`$ and $`m=l,l+1,\mathrm{},l`$.

The integral-algebraic system (18) can be solved numerically as follows. Fixing $`E`$ and $`L`$ and starting from a reasonable initial guess for $`b_{lm}(x)`$, one can compute $`u_{lm}(x)`$ from (17) ($`1`$-dim. integral). Using this, (18) can be calculated ($`2`$-dim. integral) to obtain a better form for $`b_{lm}(x)`$. This scheme can be iterated until convergence. Clearly, numerical calculations must be performed with a finite number of harmonics. A first reduction is obtained by excluding odd harmonics from the calculation. Exclusion of $`l=1`$ fixes the center of mass in the origin, while absence of higher-order odd harmonics prevents the formation of asymmetric structures (e.g. two clusters of different sizes lying at different distances from the origin). Besides this, from the typical behaviour of $`b_{lm}(x)`$ obtained from the above procedure, shown in Fig. 1, one clearly sees that $`b_{lm}`$ dies out fast as $`l`$ increases. Therefore, all calculations were performed with even harmonics up to and including $`l=16`$. We set the particles’ masses to one and measured energy and angular momentum in units of $`GN^2/R`$ and $`(RGN^3)^{1/2}`$, respectively. Solutions of (13) at fixed $`E`$ and $`L`$ were obtained, choosing $`\mathrm{\Theta }=0.02`$ always. For the sake of simplicity, we took the angular momentum parallel to the $`3`$-axis.

The resulting phase diagram is reported in Fig. 2. In order to discern pure phases from phase-coexistence regions (mixed phases), the Hessian of $`S_N`$, namely

$$[S_N]=\text{det}\left(\begin{array}{cc}_E^2S_N& _L_ES_N\\ _E_LS_N& _L^2S_N\end{array}\right)$$
(19)
must be analyzed gross174 ; gross186 . Indeed, in pure phases one has $`[S_N]>0`$, while in mixed phases $`[S_N]<0`$. In the latter, the specific heat is negative and statistical ensembles (microcanonical and canonical) are inequivalent. For low angular momenta, the system is more likely to be found in a single dense cluster (SC) at low energies, while for higher energies the most probable state is a gaseous cloud (G). The two phases, both pure, are separated by a mixed phase with negative specific heat where different equally-probable equilibrium configurations coexist. One thus recovers the usual collapse scenario thirring70  that is found in theories without angular momentum. For higher angular momenta, instead, the most probable equilibrium configuration is a double cluster (DC, pure phase), although the gas remains the most probable at sufficiently high energies. A sample of equilibrium density profiles is shown in Fig. 3. A central question is clearly that of stability of such structures. A detailed analysis will be given elsewhere, however rings turn out to be unstable, at odds with single and double clusters.

The appearance of the previously-unobserved double-cluster solutions at high angular momenta is particularly remarkable (and stresses again the importance of angular momentum in self-gravitating systems), because in such a state rotational symmetry is spontaneously broken. This is shown explicitly in Fig. 4, where the entropy is plotted as a function of the order parameter $`I_{11}I_{22}`$ (measured in units of $`NR^2`$), $`I_{ab}`$ being the components of the inertia tensor $`𝕀`$ ($`a,b=1,2,3`$). (We remind the reader that we chose the angular momentum to be parallel to the $`3`$-axis.)

If $`L=0`$ (or in absence of angular momentum) the system is isotropic ($`I_{11}=I_{22}=I_{33}`$) and rotational symmetry cannot be broken. When $`L0`$, anisotropies may occur ($`I_{33}I_{11},I_{22}`$) and one can have either rotationally-homogeneous ($`I_{11}=I_{22}`$) or rotationally-heterogeneous ($`I_{11}I_{22}`$) solutions. The latter correspond to double clusters. Indeed, the entropy profile develops two peaks at non-zero values of $`I_{11}I_{22}`$, corresponding to binary-star systems, with the two stars either aligned on the $`1`$-axis or on the $`2`$-axis.

Summarizing, we calculated the static equilibrium density profiles of a self-gravitating and rotating gas using a microcanonical mean-field approach, showing that the formation of double-cluster structures (previously observed only through dynamical approaches) can be obtained from the spontaneous breakdown of a fundamental symmetry of the Hamiltonian (1). We have presented for the first time the global phase diagram as a function of the conserved quantities $`E`$ and $`L`$. The inclusion of angular momentum in this analysis is absolutely crucial for these results, the formation of a double-cluster structure being possible only through the spontaneous breaking of a fundamental symmetry of the Hamiltonian (i.e. the rotational symmetry). However, it goes by itself that the formation and stability of stars, double stars, etc. involves other forces (e.g. nuclear and sub-nuclear) than Newtonian gravity, and hence could require other ingredients than just energy and angular momentum contopoulos60 . The results presented here provide nevertheless the most detailed static analysis of the problem to date, and bridge the existing gap between static and dynamical theories. A more exact picture of the situation could be obtained by introducing correlations, which are ignored in the present mean-field approach.

We wish to thank P.-H. Chavanis and O. Fliegans for useful discussions, comments and suggestions.