text:
# NUC-MINN-2001/18-T Domain Wall Dynamics of Phase Interfaces

## 1 Introduction

Different phases of matter are separated in space and time by dividing layers called domain walls. The dynamics of the domain walls is the determining mechanism governing phase transitions in an intermediate range between very slow, quasi-static and very rapid, dynamical processes. In this intermediate range the phase transition speed and the speed of external constraints are comparable to each other.

The dynamics of phase transitions is an involved subject even in macroscopic systems. First of all, phase transitions can be different. They may include slow burning or deflagration, detonation, condensation, evaporation, and many other forms of transition. Nevertheless, the basic conditions of all these transitions have some similarities. These arise from the basic conservation laws and from the requirement of local, or at least approximately local, equilibrium.

In a dynamical situation the approach using the equation of state with a first order phase transition is identical both in compression and in expansion. If the compression is supersonic, shock or detonation waves are formed where the final new phase is immediately created. The phase transition speed influences only the width of the shock front. On the other hand, for slow dynamics and rapid phase transition the shock front width is primarily determined by the transport coefficients, viscosity and heat conductivity, and not by the phase transition speed.

At high energies a relativistic treatment is frequently necessary. It is important to mention that a system is also relativistic if the matter is radiation dominated, meaning that the rest mass of the constituent particles is zero or negligible compared to the fourth root of the energy density. These types of systems must be treated as relativistic even if the collective velocities are small. This was one of the important new features recognized in ref. . If conserved charges do not exist the application of the conventional theory of phase transition dynamics  is not possible, and this leads to essential differences in the phase transition dynamics. Most importantly the flow is tied to the energy flow; thus energy or heat current and the fluid flow are identical and heat conduction (with respect to the flow) may not occur. Therefore the coefficient of heat conductivity does not exist, and cannot govern the phase transition speed.

The spatial configuration of phase transition instabilities varies in large homogeneous systems. When the conditions for the occurrence of a new phase are established the most common form of the appearance of the new phase happens via the formation of small critical size bubbles or droplets. (Sub-critical size bubbles or droplets will shrink and vanish but supra-critical size ones will grow.) This configuration is called homogeneous nucleation.

As the level of supercooling increases and the amount of the new phase increases other geometries become energetically more favorable, such as elongated cylindrical objects (spagetti) or layers (lasagna). The phase transition dynamics is then referred to as spinodal decomposition, indicating that systems which supercool (or superheat) and reach the adiabatic or isothermal spinodals on the phase diagram start the formation of the new phase in these configurations immediately.

Finally, when the two phases are about equally abundant and/or the transition is extremely rapid the two phases form a somewhat random occupation of the configuration space called percolation.

Most frequently explicit dynamical calculations are performed for the homogeneous nucleation geometry as this is usually the initial and the slowest of all. Nevertheless, the domain walls and their dynamical properties play an important role in all of the above mentioned configurations. One important feature of the nucleation studies that is not always recognized is that within the dynamical domain wall or droplet wall we do not have thermal equilibrium. Here, quantities like energy, entropy or particle density make sense, while other quantities that rely on thermal equilibrium, such as the equation of state, may not.

If the fluid is perfect and consists of a single phase of matter in full phase equilibrium, the flow is adiabatic . Nevertheless entropy production is still possible if the system is out of phase equilibrium , or out of thermal equilibrium. Thus, in domain walls entropy may be produced even in perfect fluids as here we do not have an equation of state which would satisfy the requirements of equilibrium thermodynamics.

In this paper we would like to discuss connections between the two most frequently encountered approaches to phase transitions using effective field theoretical models and phenomenological thermal and fluid dynamical methods. In section 2 we first demonstrate with a simple effective field theoretical model how one can obtain a Landau-model type of phenomenological description of critical dynamics. In section 3 we elaborate on the assumptions and error of different approaches used recently in the literature. We present a summary and conclusion in section 4.

## 2 The Free Energy Functional

A common approach to dealing with a phase transition in relativistic quantum field theory is to assume a uniform condensate $`\varphi `$ for some scalar field and then to compute quantum and thermal fluctuations about that condensate . This is especially so for systems exhibiting spontaneous symmetry breaking where the condensate field serves as an order parameter to distinquish the two phases. Examples include the Higgs field in electroweak theory and the sigma field or scalar quark condensate in strong interaction physics. The resulting free energy density $`f`$ is obtained from the partition function $`Z`$ in the usual way.

$$f(\varphi ,T)=T\mathrm{ln}Z(\varphi ,T,𝒱)/𝒱$$
(1)
Here $`𝒱`$ is the volume. Since the high temperature symmetric phase usually corresponds to $`\varphi =0`$, it is customary to define the effective potential as the deviation of the free energy from its value in the symmetric phase.

$$V(\varphi ,T)=f(\varphi ,T)f(0,T)$$
(2)
The point $`\varphi =0`$ is either a global minimum or only a local minimum depending on whether $`T`$ is greater than or less than the critical temperature $`T_c`$, respectively. Therefore the results of computations are often found to be, or parametrized as, a fourth order polynomial in $`\varphi `$.

$$V(\varphi ,T)=\underset{n=2}{\overset{4}{}}a_n(T)\varphi ^n$$
(3)
Occasionally one will find additional terms of order $`\varphi ^4\mathrm{ln}\left(\varphi /\mathrm{\Lambda }\right)`$ or higher powers of $`\varphi `$, but these are not common and will not change our analysis qualitatively.

To give a specific example we will use the parametrization of ref. .

$$V(\varphi ,T)=\frac{1}{2}\gamma (T^2T_0^2)\varphi ^2\frac{1}{3}\alpha T\varphi ^3+\frac{1}{4}\lambda \varphi ^4$$
(4)
The $`\gamma `$, $`T_0^2`$, $`\alpha `$ and $`\lambda `$ are temperature independent constants, to be specified shortly. In an equilibrium state the free energy density is the negative of the pressure, $`f(T)=p(T)`$. Normalizing to the high temperature symmetric phase, and denoting the equilibrium pressure of that phase by $`p_h(T)`$, we have

$$f(\varphi ,T)=p_h(T)+V(\varphi ,T).$$
(5)
It is not enough to specify the effective potential; the equilibrium pressure as a function of temperature must be specified too.

If one knows the free energy as a function of $`T`$ then the energy density $`e`$ can be calculated from the basic thermodynamic identities. In particular, if we want to know the energy density for a specific value of the condensate field, including out of equilibrium configurations also, then we must hold it fixed during the temperature differentiation.

$`e(\varphi ,T)`$ $`=`$ $`f(\varphi ,T)T{\displaystyle \frac{}{T}}f(\varphi ,T)`$ (6)
$`=`$ $`e_h(T){\displaystyle \frac{1}{2}}\gamma (T^2+T_0^2)\varphi ^2+{\displaystyle \frac{1}{4}}\lambda \varphi ^4`$
The first line above is the general thermodynamic identity; the second line applies to the specific parametrization under discussion. If one is at a local minimum of the effective potential, either $`\varphi _h(T)=0`$ or $`\varphi _l(T)>0`$ corresponding to the high and low temperature phases, respectively, then it doesn’t matter whether the order parameter is held fixed during the differentiation or not because $`f(\varphi ,T)/\varphi =0`$ at those points. These two points correspond to thermodynamical equilibrium.

If the system undergoes a first order phase transition at a critical temperature $`T_c`$ then at that temperature $`f`$ will have two degenerate minima. For a finite range of temperature above $`T_c`$ there will persist a higher, local minimum at $`\varphi _l(T)>0`$. For a finite range of temperature below $`T_c`$ there will persist a higher, local minimum at $`\varphi _h(T)=0`$. The most important physical quantities are the latent heat $`L=\mathrm{\Delta }e(T_c)=e_h(T_c)e_l(T_c)`$, the correlation lengths $`\xi ^2=^2f(0,T_c)/\varphi ^2=^2f(\varphi _l(T_c),T_c)/\varphi ^2`$, which are equal for the fourth order effective potential at $`T_c`$, and the surface energy $`\sigma `$. At the critical temperature the planar interfacial profile has a nice analytical solution on account of the fact that the effective potential becomes symmetric. In the usual way one finds that the profile field $`\overline{\varphi }(x)`$ is

$$\overline{\varphi }(x)=\frac{\varphi _l(T_c)}{2}\left[1\mathrm{tanh}(x/2\xi )\right].$$
(7)
This interpolates through intermediate non-equilibrium states from $`\varphi _l(T_c)`$ for $`x\xi `$ to $`\varphi _h(T_c)=0`$ for $`x\xi `$. The surface energy can then be expressed as

$$\sigma =_{\mathrm{}}^{\mathrm{}}𝑑x\left(\frac{d\overline{\varphi }}{dx}\right)^2.$$
(8)
The parameters in the effective potential can be expressed in terms of these physical parameters as follows.

$`\gamma `$ $`=`$ $`\left(1+{\displaystyle \frac{L\xi }{6\sigma }}\right){\displaystyle \frac{1}{T_c^2\xi ^2}}`$
$`T_o^2`$ $`=`$ $`{\displaystyle \frac{L\xi }{L\xi +6\sigma }}T_c^2`$
$`\alpha `$ $`=`$ $`{\displaystyle \frac{1}{3\sigma \xi ^3}}`$
$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{T_c\xi }}\sqrt{{\displaystyle \frac{3}{2\sigma \xi ^3}}}`$ (9)
For purposes of illustration, and with the QCD phase transition in mind, we choose the following numerical values.

$`\xi `$ $`=`$ $`{\displaystyle \frac{1}{T_c}}`$
$`\sigma `$ $`=`$ $`{\displaystyle \frac{17}{6}}{\displaystyle \frac{\pi ^2}{45}}T_c^3`$
$`L`$ $`=`$ $`{\displaystyle \frac{68\pi ^2}{45}}T_c^4`$ (10)
For the pressure in the high temperature phase we take 37 effective massless bosonic degrees of freedom corresponding to gluons and two flavors of quark, plus a bag constant $`B=L/4`$ to simulate confinement.

$$p_h(T)=\frac{37\pi ^2}{90}T^4B$$
(11)
The numerical value for the latent heat $`L`$ corresponds to a transition from the 37 degrees of freedom mentioned above to 3 degrees of freedom for massless pions. It is convenient to make everything dimensionless by measuring energies in units of $`T_c`$ and lengths in units of $`1/T_c`$, which practice we adhere to in the rest of this section.

The effective potential (4) as a function of $`\varphi `$ is plotted in figure 1 for three different temperatures: the critical one and 1% above and below it. This displays the traditional behavior of a strong first order phase transition. The plot includes negative values of $`\varphi `$ which may or may not be allowable, depending on the origin of $`\varphi `$ and its physical interpretation.

The free energy density (5) as a function of $`\varphi `$ is plotted in figure 2 for the same temperatures. The only difference between this and figure 1 is the addition of the pressure as evaluated in the high temperature symmetric phase.

The energy density (6) as a function of $`\varphi `$ is plotted in figure 3. It has a maximum at $`\varphi =0`$. This follows directly from the expression for $`e`$ given earlier: For small values of $`\varphi `$ the deviation from $`e_h(T)`$ is quadratic in $`\varphi `$ with negative curvature. Note that the energy density goes negative when $`\varphi `$ is greater than about 2.0.

Rather than using the field $`\varphi `$ as the order parameter one might consider other choices more appropriate for the problem at hand. For example, one could take as the order parameter the deviation of the energy density from its equilibrium value and expand the free energy in a power series in this difference. This is the Landau approach to the description of fluctuations and to departures from equilibrium states . Such an approach to the nuclear liquid-gas phase transition was implemented by Goodman et al. . It is also this same function that was used in working out the nucleation of relativistic first order phase transitions by two of us . With these motivations we plot the free energy versus the energy in figure 4 by using $`\varphi `$ as a parameter. There are a number of points to be made concerning this figure. First, there are two minima at each of the chosen temperatures. The minima occur at the energy densities $`e_l(T)`$ and $`e_h(T)`$. The lower of these two represents the equilibrium state while the other represents a metastable state except at $`T_c`$ when they are degenerate. The value of $`f`$ at a minimum is equal to the negative of the pressure in that particular phase: $`p_l(T)=f(e_l(T),T)`$ and $`p_h(T)=f(e_h(T),T)`$. Second, the value of $`e_l(0.99T_c)<0`$! This is an indication of the inadequacy of the specific parametrization of the effective potential at this temperature, and probably at lower temperatures too. This is an important point to always be aware of when writing down any formula for an effective potential. Third, the free energy has a cusp at the location of the high density symmetric phase. This is an unavoidable consequence of the inference of the free energy from the effective potential. From figure 3 or from the corresponding formula it is clear that the energy density is a maximum at $`\varphi =0`$. This means that the plot of the free energy ends at the equilibrium density $`e_h(T)`$ of the symmetric phase with a cusp. Physically this must be a restriction of fluctuations to those associated with $`\varphi `$ and only $`\varphi `$. This is too restrictive by far; there are certainly other physical processes not taken into account. For example, placing particles in a box in contact with a heat reservoir allows for the exchange of energy between the particles and the reservoir, resulting in fluctuations in energy and pressure. These processes are always present in systems at fixed temperature and are not accounted for, or associated with, fluctuations in an order parameter.

Now let us follow the Landau construction of the free energy away from equilibrium states using the energy density as the order parameter with no reference to $`\varphi `$ whatsoever. This is the construction made in  in the context of constructing spherical surfaces separating two phases at noncritical temperatures. The total free energy of interaction is the sum of a gradient energy and a free energy density integrated over space.

$$F_I\{e(𝐱)\}=d^3x\left[\frac{1}{2}K(e(𝐱))^2+f(e(𝐱),T)\right]$$
(12)
The $`f`$ is most economically written as a fourth order polynomial in $`e`$.

$$f(e,T)=\underset{n=0}{\overset{4}{}}b_n(T)e^n$$
(13)
The coefficients are functions of temperature. They are determined by several requirements. The first is that $`f`$ has minima located at $`e_l(T)`$ and $`e_h(T)`$, such that $`f`$ at those points be equal to the negative of the corresponding equilibrium pressure.

$`f(e_h(T),T)`$ $`=`$ $`p_h(T)`$
$`f(e_l(T),T)`$ $`=`$ $`p_l(T)`$ (14)
This results in four equations involving the five $`b_n`$. Two more equations result from fixing the correlation length (which is the same in the two phases only at $`T_c`$) and the surface energy $`\sigma `$ (which, strictly speaking, is only well defined at $`T_c`$). This set of six equations is not inconsistent because the coefficient of the gradient energy, $`K`$, must also be determined.

$`^2f/e^2|_{e=e_h(T_c)}`$ $`=`$ $`^2f/e^2|_{e=e_l(T_c)}={\displaystyle \frac{6\sigma }{\xi L^2}}`$
$`K`$ $`=`$ $`{\displaystyle \frac{6\sigma \xi }{L^2}}`$ (15)
We insist upon the same equilibrium energy densities and pressures as used in the effective potential approach, and use the same correlation length and surface energy at $`T_c`$ too.

$`\overline{e}(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[e_h(T_c)+e_l(T_c)+\mathrm{\Delta }e(T_c)\mathrm{tanh}(x/2\xi )\right]`$ (16)
$`\sigma `$ $`=`$ $`K{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x\left({\displaystyle \frac{d\overline{e}}{dx}}\right)^2`$ (17)
Then the Landau expansion of the free energy as a function of the energy density is obtained and is plotted in figure 5. The $`f(e,T)`$ is now a smooth and well-behaved function. Its only failing is that the equilibrium energy density at $`0.99T_c`$ is negative, but that is a consequence of insisting that it be the same as for the effective potential, which was negative. (Actually, fluctuations into states with negative energy can be avoided if one uses the Laurent expansion instead of the fourth order polynomial approximation , but only if the lower minimum is at positive energy density.) In fact, with the Landau approach we have a closer connection between the physical observables and the free energy. For example, we can easily specify the equilibrium energy density and pressure in each phase, making it simple to avoid such unwelcome behavior as a negative energy density in equilibrium. Of course the direct connection with the condensate field $`\varphi `$ is lost.

The main difference between the two approaches is that the field theoretical approach parametrizes the out-of-equilibrium configurations in terms of an effective potential as a function of $`\varphi `$ that shows smooth, quadratic minima around the two equilibrium states, whereas the Landau type of thermodynamical approach does the same in terms of the free energy as a function of the energy density $`e`$. These two approaches are not fully equivalent: the highly nonlinear $`f(e)`$ dependence obtained from this particular field theoretical model does not yield a smooth quadratic minimum in terms of $`e`$ at $`e_h`$ (where $`\varphi =0`$). This leads to unphysical estimates for energy density fluctuations, in our opinion. Which approach one takes depends on the physical situation being addressed.

## 3 Dynamics of First Order Phase Transitions

A phase transition occurs because of some change in the global properties of the system. Examples of such changes are: an expansion chamber or a sudden quench in temperature in the laboratory, the expansion of hot matter in a high-energy heavy ion collision, and the expansion of the early universe. Generally we need to compare the phase transition rate to the rate of expansion or quench. The speed at which the phase transition proceeds comes into play when the speed at which the external variables, such as the volume, becomes comparable to or exceeds that of phase conversion.

If the rate of change of an external variable is slow there is sufficient time to maintain phase equilibrium on a coarse-grained time scale. This also means that all other equilibration processes are completed as these invariably require less time and less interaction than the conversion of one phase of matter to another.

Thus, in the case of slow external dynamics and rapid phase equilibration the matter is in complete equilibrium, including phase equilibrium, and the equation of state of the matter undergoing a first order phase transition is given by the Maxwell construction. Then we have a fully developed mixed phase, and the phase abundances are typically given by the conservation of energy and entropy. No information on dynamical processes is needed. This type of transition is adiabatic. Even in moderately fast dynamical expansions there are only small deviations from the ideal and complete phase equilibrium (Maxwell construction). This deviation results in some delay in the creation of the new phase, leading to supercooling or superheating and extra entropy production.

For heavy ion reactions the first attempt to explicitly evaluate the phase transition speed of homogeneous nucleation was described in . The homogeneous nucleation mechanism correctly describes the initial stage of the phase transition when the abundance of the newly created phase is still small and when the rate of phase conversion is the slowest process.

Here a couple of remarks are necessary. There are several requirements to form critical-sized bubbles or droplets of the new phase. Both pressure balance and temperature equality should be established between the phases; this requires the transfer of energy and momentum across the phase boundary. If local equilibrium is assumed both before and after the formation of the new phase then we cannot relax the requirement of pressure and temperature equilibrium.

Langer’s modern theory of nucleation yields the following formula for the rate:

$$I=\frac{\kappa }{2\pi }\mathrm{\Omega }_0e^{\mathrm{\Delta }F/T}$$
(18)
where $`\mathrm{\Delta }F`$ is the change in the free energy of the system due to the formation of the critical droplet. $`\mathrm{\Omega }_0`$ is a statistical prefactor which measures the available phase space volume. $`\kappa `$ is a dynamical prefactor which determines the exponential growth rate of critical droplets which are perturbed from their quasi-equilibrium radius R.

$$\kappa =\frac{d}{dt}\mathrm{ln}[R(t)R_{}]$$
(19)
The dynamical prefactor has been calculated by Langer and Turski  and by Kawasaki  for a liquid-gas phase transition near the critical point, where the gas is not dilute, to be

$$\kappa =\frac{2\lambda \sigma T}{\mathrm{}^2n_{\mathrm{}}^2R_{}^3}.$$
(20)
This involves the thermal conductivity $`\lambda `$, the surface free energy $`\sigma `$, the latent heat per molecule $`\mathrm{}`$ and the density of molecules in the liquid phase $`n_{\mathrm{}}`$. The interesting physics in this expression is the appearance of the thermal conductivity. In order for the droplet to grow beyond the critical size, latent heat must be conducted away from the surface into the gas. For a relativistic system of particles or quantum fields which has no net conserved charge, such as baryon number, the thermal conductivity vanishes. The reason is that there is no rest frame defined by the baryon density to refer to heat transport. Hence this formula obviously cannot be applied to such systems. The dynamical prefactor for such systems was first evaluated a’ la Langer by two of us  to be

$$\kappa =\frac{4\sigma }{(\mathrm{\Delta }w)^2R_{}^3}\left[\frac{4}{3}\eta +\zeta \right].$$
(21)
Here $`\eta `$ and $`\zeta `$ are the shear and bulk viscosities, respectively, and $`\mathrm{\Delta }w`$ is the enthalpy density difference between the two phases (equal to the latent heat since the pressures are equal at the critical temperature). This fully relativistic expression was subsequently generalized by Venugopalan and Vischer  to include a conserved baryon number.

$$\kappa =\frac{2\sigma }{(\mathrm{\Delta }w)^2R_{}^3}\left[\lambda T+2\left(\frac{4}{3}\eta +\zeta \right)\right]$$
(22)
In the non-relativistic limit and when the viscosities are small compared to the heat conductivity this reduces to Langer’s expression. In the relativistic limit with no net baryon number, where effectively $`\lambda 0`$, this reduces to equation (21).

Unfortunately there are several erroneous expressions for $`\kappa `$ which subsequently appeared in the literature. In  one finds that in the relativistic limit with no net baryon number

$$\kappa =\sqrt{\frac{2\sigma w_h}{(\mathrm{\Delta }w)^2R_{}^3}}.$$
(23)
This was ”derived” under the assumption of vanishing shear viscosities. An expression for $`\kappa `$ in the nonrelativistic limit with non-zero baryon number was found which differed from Langer’s result too. By including the viscosities another paper  finds that

$$\kappa =\sqrt{\frac{2\sigma w_h}{(\mathrm{\Delta }w)^2R_{}^3}}+\frac{1}{c_s^2}\frac{\sigma }{(\mathrm{\Delta }w)^2R_{}^3}\left[\frac{4}{3}\eta +\zeta \right].$$
(24)
Here $`c_s`$ is the speed of sound in the low density phase. In the limit that the viscosities vanish this reduces to the result of , and in the limit that they are large this reduces to the result of  albeit with a factor of $`1/c_s^2`$ rather than 4.

The results found by  and  are suspicious because they predict a dynamical growth factor even in the absence of viscous forces. These results are wrong because of a misinterpretation of a derivative. Both papers follow the approach of  relatively closely except for one crucial point. At some point in the analysis one encounters the derivative $`f(\overline{e},T)/\overline{e}`$ evaluated at the radial profile $`\overline{e}=\overline{e}(r)`$ of the energy density for a spherical bubble or droplet connecting the two phase. Thus $`\overline{e}(r)`$ varies between $`e_l`$ and $`e_h`$. The free energy density varies accordingly. As explained carefully and in some detail in , this derivative is taken at fixed temperature. At the equilibrium points it is zero: see the discussion in section 2 above. However, both  and  equate $`f`$ with $`p`$ irrespective of whether we are at one of the minima or not. They then assume an equation of state $`p=c_s^2e`$ to finally obtain $`f(e,T)/e=c_s^2`$. From that point on the analysis diverges, and the erroneous results for $`\kappa `$ inevitably follow.

It is also possible to argue that the erroneous results for $`\kappa `$ are not only mathematically incorrect but physically incorrect too. What is relevant is the slowest required process, not the fastest. Furthermore it is incorrect to add up the rates of all processes to obtain the highest possible phase transition speed. If not all required processes are completed the phase transition is not complete either. For example to establish pressure (momentum) equilibrium  is not sufficient, because then a subsequent step of establishing thermal equilibrium is required, and only when both processes are over has the phase transition been completed.

In ref.  the slowest process, heat transfer via viscosity, was evaluated. In some cases, when we have non-negligible net baryon charge in our system, the heat conductivity may also contribute to the heat transfer to the new phase, thus speeding up somewhat this slowest of all processes . If the transport coefficients are all vanishing so that thermal balance cannot be achieved, we will never reach both phase and thermal equilibrium, so the rate tends to zero as stated in ref. .

In case other degrees of freedom exist which permit energy transport between the phases leading to a common temperature , these processes can be combined with other transport processes like viscosity and/or heat conduction leading to faster temperature balance and a higher rate. However, sound waves do not lead to dissipation as perfect fluid dynamics is adiabatic (in the absence of shock waves). Thus, the sound wave as the sole mechanism cannot characterize dissipative transport processes, as was mistakenly claimed in  and .

A final remark: If the external dynamical evolution is much faster than the processes involved in the phase and kinetic equilibration then the matter involved in the phase transition loses both phase and thermal equilibrium, and we have to abandon the thermal and fluid dynamical approaches altogether. An example in the context of high energy heavy ion collisions is presented in  where an effective field theoretical approach is used. Generally this will be the case for small, rapidly developing systems only.

## 4 Conclusions

In this paper we have examined two subtle issues in the dynamics of domain walls. One of them is the difference between the effective potential as a function of a scalar field and the free energy as a function of energy density, both evaluated near the critical temperature. The other involves a misconception in the literature about the derivation of the dynamical prefactor, or growth rate, in homogeneous nucleation theory. It is important to understand these subtleties because of their importance in such physical environments as cosmology, astrophysics and high energy nuclear collisions.

## Acknowledgements

This work was supported by the US Department of Energy under grant DE-FG02-87ER40328 and by the Research Council of Norway.


format: text

text:
# Theoretical model for magnetooptic imaging

## I Introduction

For centuries humans have tried to ’see’ magnetic fields using various devices. Iron particles, magnetic fluids, inductive coils, Hall probes and magnetic bacterias are only a few of the approaches which have been adopted. During the past 40 years the research on magnetooptic visualization has been highly motivated by the need to ’see’ magnetic fields from superconductors. To that end, magnetooptic visualization using thin film indicators has been particularly useful. In 1957 P.B. Alers visualized for the first time the magnetic field distribution in a superconductor using a mixture of cerous nitrate-glycerol. A highly successful development followed in 1968 when H. Kirchner noticed that chalcogenide-films could be used to visualize fields with improved sensitivity. In the following 20 years magneto-optic studies were published periodically, mainly looking at flux profiles in superconductors. However, there were no major advances until 1991, when it was discovered that bismuth-substituted iron garnet films with in-plane magnetization are excellent indicators with almost no domain activity. This discovery triggered a large number of quantitative investigations of the flux behaviour in high-temperature superconductors. Two of the most recent advances in this field are ultrafast magnetooptical studies (see e.g ref.) and imaging of single vortices.

Although magnetooptical indicators have mainly been used to image magnetic fields from superconductors, a large number of other possible application areas exists. Thus, it has been realised that such indicators are useful for imaging magnetic defects, magnetic phase transitions, microelectronic circuits and magnetic storage media.

Since the field of magnetooptic visualization has matured for more than 40 years, one may expect that most parts of the image formation process is well understood. This is not the case. In fact, only recently a model for magnetooptic imaging has been published. However, this model is mostly concerned with the magnetization redistribution in the indicator. The purpose of the current paper is to develop a more detailed model for the magnetooptic imaging system (see Fig. 1), including the indicator, and to use this to examine the possibility of imaging weak magnetic fields, e.g. single vortices in superconductors.

## II Magnetization in the indicator

The basic sensor element for detection of magnetic fields is the magnetooptic indicator. It is therefore first necessary to know how the indicator respond to an external field. Let $`𝑴`$ be the magnetization vector in the indicator, and $`M_s`$ be the saturation magnetization. Then we have from Fig. 2

$$𝑴=M_s[cos\beta cos\varphi _M,cos\beta sin\varphi _M,sin\beta ]=[M_rcos\varphi _M,M_rsin\varphi _M,M_z].$$
(1)
The external field is given by

$$𝑯=[H_rcos\varphi _H,H_rsin\varphi _H,H_z],$$
(2)
where $`M_r=\sqrt{M_x^2+M_y^2}`$ and $`H_r=\sqrt{H_x^2+H_y^2}`$. We assume that the indicator consists of a single domain, and that its free energy can be expressed as a sum of the uniaxial anisotropy energy ($`E_u`$), the Zeeman energy ($`E_z`$) and the demagnetizing energy ($`E_d`$)

$$E=E_u+E_z+E_d.$$
(3)
The uniaxial anisotropy energy can be written as,

$$E_u=K_ucos^2\beta ,$$
(4)
where $`K_u`$ is the uniaxial anisotropy constant. The Zeeman energy is given as

$$E_z=\mu _0𝑴\mathbf{}𝑯=\mu _0M_s[H_rcos\beta cos(\varphi _M\varphi _H)+H_zsin\beta ].$$
(5)
Since we neglect the cubic anisotropy of the indicator the in-plane magnetization direction must be the same as that of the external field ($`\varphi _M=\varphi _H`$). The demagnetizing energy, which tends to keep the magnetization in the plane of the film, can be given as

$$E_d=\frac{\mu _0M_s^2}{2}sin^2\beta .$$
(6)
Minimizing $`E`$ with respect to $`\beta `$ results in the following relationship

$$H_asin\beta cos\beta +H_rsin\beta H_zcos\beta =0,$$
(7)
where the socalled anisotropy field is given by

$$H_a=M_s\frac{2K_u^{tot}}{\mu _0M_s}.$$
(8)
If we assume that $`H_r=0`$, then the resulting expression for the magnetization is

$$M_r=M_s\sqrt{1\left(\frac{H_z}{H_a}\right)^2},H_zH_a,$$
(9)
and

$$M_z=M_s\frac{H_z}{H_a},H_zH_a.$$
(10)
If $`\beta `$ is very small (corresponding to very small $`H_z`$, or large $`H_a`$ or $`H_r`$), we may put $`sin\beta \beta `$ and $`cos\beta 1`$, and the magnetization becomes

$$M_rM_s,H_zH_a,$$
(11)
and

$$M_zM_s\frac{H_z}{H_a+H_r},H_zH_a,$$
(12)
Thus we see that an in-plane component $`H_r`$ will only reduce the z-component of the magetization. Therefore, it is necessary to keep $`H_r`$ and $`H_a`$ small if large sensitivity to $`H_z`$ is needed.

## III The magnetooptic imaging system

The basic imaging system is shown in Fig.1. In principle, it is a conventional polarizing microscope, although many different designs exists. Basically, the z-component (or any other component) of the magnetization is converted into a small rotation of the plane of polarization, which can then be detected by an array detector (e.g. CCD). The signal on the detector can be found by an analysis based on Jones matrices together with the methods used in refs.. For a general polarization distribution incident on the objective we write

$$𝐄_\mathrm{𝟎}=\left[\begin{array}{ccc}a(\theta ,\varphi )& & \\ b(\theta ,\varphi )& & \end{array}\right],$$

where $`\theta `$ is the incident angle and $`\varphi `$ is the azimuthal angle. In this paper we assume that the incident polarization after the polarizer consists of a x-component only (a=1 and b=0). The effect of the objective lens is to change the incident polarization into s and p-polarized light with the following transformation

$$𝐋=\left[\begin{array}{ccc}cos\varphi & sin\varphi & \\ sin\varphi & cos\varphi & \end{array}\right].$$

Next, the rays passes through the isotropic substrate on which the indicator is deposited

$$𝐓_{\mathrm{𝐢𝐧}}=\left[\begin{array}{ccc}t_p^{in}& 0& \\ 0& t_s^{in}& \end{array}\right].$$

When the rays are reflected from the magnetooptic layer, some of the s-polarization is transformed into p-polarization (and vice versa). In order to assure maximum reflection, one may evaporate a thin mirror onto the indicator. The reflection matrix from the indicator and the mirror can be expressed as

$$𝐑=\left[\begin{array}{ccc}r_{pp}& r_{ps}& \\ r_{sp}& r_{ss}& \end{array}\right],$$

where $`r_{pp}`$ and $`r_{ss}`$ are the usual Fresnel coefficients. In the first order approximation (corresponding to weak magnetooptic effects), we find the following off-axis components for arbitrary magnetization vectors and angles of incidence

$$r_{ps}=\frac{\pi n_1Q}{\lambda _0}t_{01}^st_{10}^p_Df(z)\left(\frac{a(z)M_z+b(z)M_ytan\theta _1}{M_s}\right)𝑑z,$$
(13)
$$r_{sp}=\frac{\pi n_1Q}{\lambda _0}t_{01}^st_{10}^p_Df(z)\left(\frac{a(z)M_zb(z)M_ytan\theta _1}{M_s}\right)𝑑z,$$
(14)
where $`n_1`$ is the refractive index of indicator, Q is the magneto-optical material constant, $`\lambda _0`$ is the wavelength in vacuum, $`t_{ij}^{s,p}`$ is the s or p-component of the transmission coefficient through the interface from medium i to j (similar for $`r_{ij}^{s,p}`$), $`\theta _1`$ is the propagation angle in the indicator, $`f(z)=exp[i4\pi n_1zcos\theta _1/\lambda _0]`$, and

$$a(z)=\frac{[1+r_{12}^sf(Dz)][1+r_{12}^pf(Dz)]}{[1r_{10}^sr_{12}^sf(D)][1r_{10}^pr_{12}^pf(D)]},$$
(15)
$$b(z)=\frac{[1+r_{12}^sf(Dz)][1r_{12}^pf(Dz)]}{[1r_{10}^sr_{12}^sf(D)][1r_{10}^pr_{12}^pf(D)]},$$
(16)
according to the recipy of Hubert and Traeger. The integration is carried out over the thickness (D) of the indicator. Here we have neglected the transverse Kerr-effect, since this is rather small in most indicators.

For strong magnetooptical effects, one should use a similar formalism, but take into account multiple scattering. This complicates the problem considerably, and will not be pursued here.

On its way back to the array detector, the light will again pass through the substrate

$$𝐓_{\mathrm{𝐨𝐮𝐭}}=\left[\begin{array}{ccc}t_p^{out}& 0& \\ 0& t_s^{out}& \end{array}\right].$$

Furthermore, it is collected by the objective lens, and in total we can write the electric field before the analyzer as

$$E=L^1T_{out}RT_{in}LE_0.$$
(17)
After evaluating the equation above we observe that also the Fresnel coefficients between the substrate and air influences the polarization state of the outgoing beam. In many cases it is possible to deposit a dielectric coating which minimizes the depolarization due to the substrate. Here we assume the most ideal case, in which all transmission coefficients become unity. In this case we may write

$$E_x=\frac{1}{2}(r_{pp}+r_{ss})+\frac{1}{2}(r_{pp}r_{ss})cos2\varphi \frac{1}{2}(r_{sp}+r_{ps})sin2\varphi ,$$
(18)
and

$$E_y=\frac{1}{2}(r_{sp}r_{ps})+\frac{1}{2}(r_{pp}r_{ss})sin2\varphi +\frac{1}{2}(r_{sp}+r_{ps})cos2\varphi .$$
(19)
These equations can be used to calculate the expected polarization pattern at the exit pupil seen in the conoscopic image, which is often used to judge the quality of the imaging system. To find the electric field at the detector plane, ($`r^{},\varphi _c`$), we assume that the imaging system obeys the sine condition, and use the method of refs. to obtain

$$E_{xD}=I_0^a+I_2^acos2\varphi _c+I_2^bsin2\varphi _c,$$
(20)
and

$$E_{yD}=I_0^b+I_2^asin2\varphi _cI_2^bcos2\varphi _c,$$
(21)
where

$$I_0^a=_0^\alpha (r_{pp}+r_{ss})sin\theta cos\theta J_0(kr^{}sin\theta )exp(2ikzcos\theta )𝑑\theta ,$$
(22)
$$I_0^b=_0^\alpha (r_{sp}r_{ps})sin\theta cos\theta J_0(kr^{}sin\theta )exp(2ikzcos\theta )𝑑\theta ,$$
(23)
$$I_2^a=_0^\alpha (r_{ss}r_{pp})sin\theta cos\theta J_2(kr^{}sin\theta )exp(2ikzcos\theta )𝑑\theta ,$$
(24)
$$I_2^b=_0^\alpha (r_{sp}+r_{ps})sin\theta cos\theta J_2(kr^{}sin\theta )exp(2ikzcos\theta )𝑑\theta .$$
(25)
Here $`J_n`$ is the Bessel function of the first kind (of order $`n`$), and a constant in front of the integrals is ommited. Note that the integral is taken over the angular aperture of the objective lens, which has a numerical aperture $`NA=sin\alpha `$. In deriving these formulas, we have assumed that the objective is corrected for all spherical aberrations introduced by the film and the substrate. Here the factor $`cos\theta `$ occurs since we have chosen aplanatic apodization at the exit pupil ($`\sqrt{cos\theta }`$, squared for two passes through the objective). The phase factor $`exp(i2kzcos\theta )`$ is included to account for the defocus of the light beam. Here we will assume that z=0, and therefore that the light beam is at focus on the mirror-plane. When an analyzer is placed in front of the detector, the intensity at the detector plane must be written as

$$i_D=|E_{xD}cos\gamma +E_{yD}sin\gamma |^2,$$
(26)
where $`\gamma `$ is the analyzer’s angle from the x-axis. So far, we have not considered diffraction due to the localized magnetization distribution in the indicator. However, this can be included by using the procedure of Kambersky et al.. In that paper a procedure for calculating the diffracted Kerr or Faraday amplitude is presented. The result must be weighted by the polarization vector of the lens system, and summed up at the exit pupil. This is a large numerical task, at least for a reflection-type microscope, and will not be pursued here. That is, we will neglect diffraction from magnetization gradients in the indicator.

## IV Imaging of single vortices in a superconductor

Only recently magnetooptic imaging of single vortices was achieved experimentally. In this section we study magnetooptic imaging of single vortices using approximative solutions of the theory given above. Interestingly, it has been found that the field from a vortex is similar to that from a magnetic monopole located a distance $`z_0=1.27\lambda `$ ($`\lambda `$ is the penetration depth) below the superconductor surface. In this approximation the scalar potential from a vortex can be written as

$$\varphi _V=\frac{\mathrm{\Phi }_0}{2\pi \mu _0}\frac{1}{\sqrt{x^2+y^2+(zz_0)^2}},$$
(27)
where $`\mathrm{\Phi }_0`$ is the flux quantum.

To get some insight into the behaviour of the magnetooptic response, let us neglect the absorption and the multiple reflections, and assume that $`H_aH_r,H_z`$. Then the off-diagonal components of the reflection becomes

$$r_{ps}Qt_{01}^st_{10}^p\left\{\frac{1}{H_a}\left[\varphi _V(d)\varphi _V(d+D)\right]+Dtan\theta _1sin\varphi _M\right\},$$
(28)
and

$$r_{sp}Qt_{01}^st_{10}^p\left\{\frac{1}{H_a}\left[\varphi _V(d)\varphi _V(d+D)\right]Dtan\theta _1sin\varphi _M\right\},$$
(29)
where $`d`$ is the distance between the superconductor and the indicator. Although these expressions are approximations, they still contain some interesting physics. Note that the second term is constant, i.e. independent of the vortex field, since we have assumed that $`H_aH_r,H_z`$. It is also seen that the off-axis reflection coefficients are strongly dependent on the penetration depth of the superconductor. That is, when the penetration depth is large, the magnetic monopole appears to be located far from the surface of the superconductor, whereas when it is small it is located close to the surface. Thus the signal increases with decreasing penetration depth. As an example of this, Fig. 3 shows $`r_{ps}`$ (without the constant term) as a function of the penetration depth when d=125 nm (dashed line) and d=250 nm (solid line). It is noted that the magnetooptical signal falls off nonlinearly with the penetration depth, but that the gradient is less severe when the distance between the superconductor and the indicator increases.

To find the actual signal at the array detector, one must evaluate Eq. 26. When the polarizer and analyzer are in nearly crossed positions ($`\alpha 90^{}`$), and $`\varphi _c=0^{}`$, the signal can be approximated by

$$i_DQ^2\left[\varphi _V^{im}(d)\varphi _V^{im}(d+D)\right]^2|_0^\alpha t_{01}^st_{10}^psin\theta cos\theta J_0(kr^{}sin\theta )𝑑\theta |^2.$$
(30)
Here $`[\varphi _V^{im}(d)\varphi _V^{im}(d+D)]^2`$ may be interpreted as the geometric image of a given position ($`r,\varphi `$) at the indicator. According to the Huygens-Fresnel principle, each such point generates a spherical secondary wavelet. Since these wavelets must pass through the apertures of the optical system, the ’perfect’ geometric image is smeared out at the image plane. That is, each point in the geometric image must be multiplied by a diffraction integral and the point becomes a distribution at the detector. One should also note that in the present approximation the diffraction integral is modified by the off-axis reflection coefficient through the factor $`t_{01}^st_{10}^p`$. In this paper we have not considered the response of the array detector, as this is assumed to be ideal.

It is very important to note that the signal is proportional to the difference in magnetic potential at the two surfaces of the indicator. For this reason the signal increases with thickness, but the signal-distribution also broadens significantly. This issue will be of particular importance when the vortices are placed close together, and the flux-profiles overlap. Thus one may expect a thin indicator to provide better signal contrast than a thick one, in particular since the contrast is often more important than the strength of the signal. To illustrate this point, Fig. 4 shows $`[\varphi _V^{im}(d)\varphi _V^{im}(d+D)]^2`$ at the detector plane when the vortices are placed 2 $`\mu m`$ apart, and the thickness of the indicator is D=5 $`\mu m`$ (dashed line) and D=1 $`\mu m`$ (solid line). Although the signal from the thickest films is 1.7 times stronger than that of the thinnest film, we see that the contrast is better in the latter case.

Finally, it is important to point out that in deriving Eq. (30) we have done several assumptions in order to obtain a simple expression which is readily interpreted. If a detailed analysis is required, one should adopt the full formalism presented in this paper. It is our hope that this formalism could be useful for future studies on magnetooptic visualization.

## V Conclusion

A theoretical model for investigating magnetooptic imaging of weak magnetic fields has been presented. It takes into account the details of the indicator and the magnetooptic system. To demonstrate the applicability of the model, we have tried to understand some issues concerning imaging of single vortices in superconductors. It is found that the useful signal is strongly dependent on the penetration depth of the superconductor, and also the sensitivity and thickness of the indicator.

###### Acknowledgements.

The author is grateful to P.E. Goa, H. Hauglin, M. Baziljevich and T.H. Johansen for many interesting discussions on this topic. This work was supported by the Norwegian Research Council and Tandberg Data ASA.

Figure 1

Figure 2

Figure 3

Figure 4


text:
# A Correlation Between Inclination and Color in the Classical Kuiper Belt

## 1 Introduction

Beyond Neptune, the Solar System is populated by $`10^5`$ Kuiper belt objects (KBOs, also trans-Neptunians) larger than 100 km in diameter, which represent about 0.1 Earth masses of material and are presumed to contain the least thermally processed material in the solar system (Trujillo, Jewitt & Luu 2001). As of Oct 2001, $`480`$ KBOs are known, less than 1% of the inferred population. Dynamically, the KBOs can be divided into at least three classes: (1) the Classical KBOs represent 2/3 of the observed sample and about 1/2 the inferred population; (2) the Resonant KBOs are comprised primarily of Plutinos ($`a39.4`$ AU) and objects in 2:1 mean-motion resonance with Neptune ($`a47.7`$ AU), representing 15% of the observed sample, but only 6% of the total population; and (3) the Scattered KBOs comprise about 10% of the observed sample, but nearly half the inferred population because their extreme dynamical properties ($`50\text{ AU}<a<200\text{ AU}`$ and eccentricities $`0.2<e<0.9`$) place them at unobservably great heliocentric distances for large fractions of their orbits (Trujillo, Jewitt & Luu 2000).

Since the KBOs are thought to be the least thermally modified bodies in the solar system, they are expected to provide insight into the physical and chemical processes that dominated the primitive solar nebula. Given the fact that the KBOs have been observed for nearly a decade, very little is known about their physical properties, primarily because of their extreme faintness (median $`R`$ magnitude $`m_R=23.0`$). The most striking physical feature observed of the KBOs is their substantial color diversity, from nearly neutral reflectance to the reddest objects in the solar system (Jewitt & Luu 1998 and Tegler & Romanishin 2000). Objects are reddened from the visible throughout the near-IR (Davies et al. 2000), presumably from a single coloring agent (Jewitt and Luu 2001), although the physical and chemical mechanisms for the reddening process have not been identified. Recently, Tegler & Romanishin (2000) reported that the Classical KBOs exhibited a correlation between color and inclination. No other observers have reported this trend, and Tegler & Romanishin estimate the significance of their reported trend using their two-color paradigm, which has not been confirmed by independent sources (Barucci et al. 2001 and Jewitt & Luu 2001). Thus, we chose to conduct our own $`BVR`$ observations of this population. Using our new observations, we demonstrate that the inclination and color correlation is indeed statistically significant. In addition, we note the presence of a perihelion and color correlation which appears to be due to sampling bias.

## 2 Observations

We chose to observe the Classical KBOs — those KBOs which are expected to be the least dynamically modified by Neptune. In practice, 19 Classical KBOs and 5 Scattered KBOs ($`a>50`$ AU, 30 AU $`<q<`$ 40 AU) were observed due to the distribution of objects in the sky when observations were scheduled. Results of this work are not changed significantly by the inclusion or exclusion of the Scattered KBOs, as they are a minority population in our sample. All data were taken with the same equipment, the LRIS (Oke et al. 1995) and LRIS-B (McCarthy et al. 1998) cameras at the 10 meter Keck I telescope during the nights of Feb 26–27, 2001. This dichroic-based instrument combination allows the simultaneous use of blue and visible filters (yielding $`BV`$ color) or blue and red filters (yielding $`BR`$ color), allowing a factor two increase in efficiency over the previous single camera system. The efficiency of the LRIS/LRIS-B instrument has permitted us to approximately double the published sample of broadband optical colors for the Classical KBOs in just two nights.

Standard stars (Landolt 1992) were imaged at a variety of airmasses and photometric reductions were performed using a two dimensional PSF-fitting procedure, which allows for arbitrary non-radially symmetric variations of the PSF. Using this method, $`>20`$ bright field stars were sampled in both spatial dimensions of the focal plane to estimate the true point-spread-function (PSF) for a given image. Since exposure times were fixed at 5 minutes, KBO apparent motion (3 arc-sec/hr or 0.25 arc-sec/exposure) was minimized, yielding no significant difference between the KBO PSFs and the stellar PSFs given the $`1`$ arc-sec $`B`$-band seeing (seeing for the $`R`$ and $`V`$ bands was similar, since focus was optimized for $`B`$). The stellar PSF was fit to the KBO image to determine flux and local background levels. In all cases, KBO PSFs were consistent with the model PSFs within the Poisson noise. Additionally, synthetic data tests were performed to demonstrate that the PSF fitting procedure is robust to errors in KBO centroiding. Approximately 4 $`B`$, 2 $`V`$ and 2 $`R`$ measurements were taken of each KBO, although measurements were repeated if contamination by background sources was suspected. The reduced data appear free of instrumental biases — no correlations were observed between color and data collection parameters such as airmass, seeing, time of observations, apparent motion of objects, location in the focal plane or sun-target-observer angle. The most striking feature of these data appears in Figure 1, which shows an apparent correlation between orbital inclination $`i`$ and $`BR`$ color. With few exceptions, this correlation appears to vary continuously from the reddest (low-inclination) KBOs to the bluest (high-inclination) KBOs. Both the $`BR`$ and $`BV`$ color measurements for each object are presented in Table 1. The physical implications of these observations are discussed in § 3, while an estimation of the statistical significance of the color and inclination trend, as well as other trends, appears in the remainder of this section.

### 2.1 The KBO Inclination and Color Correlation

To estimate the significance of the relationship between color and $`i`$, we computed the Spearman (1904) rank correlation statistic $`r`$ for our data set, which takes on non-zero values between $`1<r<1`$ for correlated quantities. The strength of the Spearman method is that it does not presuppose a functional form for the correlation, as do linear least-squares estimates. In addition, the Spearman method is less sensitive to the presence of deviant data points than a chi-squared rank correlation test, as it is the difference in rank between sorted data lists that is compared, not the actual data values themselves. Two statistical tests (Press et al. 1992), the $`t`$-statistic and the variance of the sum squared difference of ranks, $`\text{Var}(D)`$, can be used to estimate the significance of a non-zero value of $`r`$ (here, we report the mean of the two tests). We measure the correlation between $`BR`$ color (a stronger indicator of color than $`BV`$) and orbital inclination as $`r=0.62`$ for the $`N=24`$ data points in our new sample. The probability of this $`r`$ or a more significant one occurring in an uncorrelated sample is $`P=2.0\times 10^3`$ ($`3.1\sigma `$ significance, assuming Gaussian statistics, since $`1P`$ yields the confidence level).

In order to increase sample size, we repeated this analysis combining our data with all previously published data (Table 2). Color measurements of the same object from multiple investigators have been combined into one value, with error bars reflecting the greater of (1) the standard error of the color values or (2) the mean of the reported errors when combined in quadrature. We note that the process that produces the colors appears stochastic in nature, with individual KBOs departing from the color correlation based on their past interactions with the coloring mechanism. Error bars reflect color measurement error, and not the expected departure of an individual KBO from the overall color / inclination trend, which cannot be computed a priori without more information about the physical nature of the coloring process. We note that in cases of duplication between our objects and published objects, our color values are statistically indistinguishable from other published values. By increasing the sample size using previously published data on all KBOs, the significance of the inclination and color correlation increases to $`3.3\sigma `$ ($`N=60`$).

It appears that the correlation between color and inclination has not been widely reported because of the influence of the Plutinos, for which no correlation is apparent. Removing the Plutinos from the combined sample of data results in a large increase in the correlation statistic to $`4.1\sigma `$ ($`N=44`$). However, we must note that the removal of the Plutinos is still statistically tentative, as a two-dimensional Kolmogorov-Smirnov test (Press et al. 1992) shows that the null hypothesis that the Plutino and other KBO ($`i`$, $`BR`$) values are drawn from the same distribution can be rejected with only 93% ($`<2\sigma `$) confidence. This observation in itself has some potential implications for the origin of color diversity, as discussed in the next section. Since the correlation between inclination and color for the non-Plutino KBOs (1) increases in significance with increased sample size, (2) does not appear to be the product of any observational biases, (3) is observed in both $`BV`$ and the $`BR`$, and (4) has $`<4\times 10^5`$ ($`4.1\sigma `$) chance of being produced in an uncorrelated population, we conclude that it is statistically significant.

### 2.2 Other Correlations

The only other color trend that we observe for the Classical KBOs is a perihelion and color correlation. This trend, however, is almost certainly due to sampling bias. To illustrate this, we have considered two subsamples of the $`BR`$ data, which are shown as the greyed bars in Figure 2. The first subsample (the “constant $`q`$” subsample) shows little perihelion variation but large inclination variation and consists of objects between $`37\text{ AU}<q<40\text{ AU}`$. The constant $`q`$ subsample exhibits a $`2.8\sigma `$ correlation between inclination and color, indicating that the inclination and color correlation cannot be attributed to perihelion variations alone (Figure 3). The second subsample (the “constant $`i`$” subsample) shows little inclination variation but large perihelion variation and consists of objects between $`7\mathrm{°}<i<15\mathrm{°}`$. The constant $`i`$ subsample shows no correlation between perihelion and color, indicating that the perihelion and color trend is most likely a spurious correlation induced by sampling bias.

## 3 Interpretation

The link between color and inclination rules out very simple mechanisms of producing the observed KBO color diversity. The presence of a coloring agent (surface frosts, for example) independent of the Kuiper belt dynamical environment cannot be the only cause of color diversity, because in such a scenario, no correlation between color and any orbital parameter would be expected. More complex processes must be invoked such as evolutionary processes or remnant primordial composition effects.

Collisional resurfacing has been suggested as a possible mechanism for producing the wide range of KBO colors (Luu & Jewitt 1996 and Jewitt & Luu 2001). In this model, a close competition between radiative processing, which reddens the surface, and collisional resurfacing, which restores neutral ices to the surface, produces the observed color variations. Assuming that collisions are most likely to take place near the ecliptic where the KBO number densities are highest, a color and inclination correlation seems likely. An object’s speed perpendicular to the ecliptic follows $`\mathrm{\Delta }v\mathrm{sin}i`$. From laboratory simulations (Fujiwara, Kamimoto & Tsukamoto 1977) scaled to solar system events (Fujiwara et al. 1989), the mass ejected during a collision is proportional to impact energy squared, $`mE^2\mathrm{\Delta }v^4`$. Thus, ejecta mass could be a strong function of object inclination, $`m\mathrm{sin}^4i`$, leading to a correlation between color and inclination. If this were the case, one might expect to see a correlation between color and population averaged impact velocity $`v_{\mathrm{imp}}\sqrt{i^2+e^2}`$ (Stern 1995). Unfortunately, because of the paucity of color data, it is not currently possible to separate the effects of an $`i`$ correlation from a $`v_{\mathrm{imp}}`$ correlation because the two parameters are so strongly correlated. In addition, it is not clear why the collisional environment would be different for the Plutinos, which show no correlations between $`i`$ and color. Further modeling is required to fully quantify this effect. Unfortunately, the competition between collisions and radiative processing is quite difficult to model, as we do not have a thorough understanding of all of the astrophysical effects that take place. For instance, no dynamical models have yet estimated collision rates based on population dynamics in a manner that could demonstrate a dynmaics-induced color correlation. In addition, models of radiative processing have not addressed the effect of the complete carbonization of surfaces, which would result in a dark, neutrally colored surface. Furthermore, no collisional resurfacing models have adequately reproduced the magnitude of color variations presently observed for the KBOs (Jewitt & Luu 2001).

The best alternative to evolutionary processes is that the colors are primordial in nature. An $`i`$ correlation might be expected if the KBOs were composed of multiple sub-populations with different primordial colors and inclinations. Such a model is compatible with two recently reported dynamical groupings: (1) the Levison and Stern (2001) large and small size components of the Classical KBOs and (2) the Brown (2001) high-$`i`$ and low-$`i`$ components of the Classical KBOs. However, we find no correlations between color and any orbital parameter other than $`i`$, as might be expected for a cosmogonic origin to the colors. Even if the KBO colors are dominated by primordial processes, the role of collisions in altering KBO surfaces must still be addressed given that (1) collisions are thought to be responsible for reducing the mass of the Kuiper belt by a factor $`100`$ since the epoch of formation (Stern & Colwell 1997 and Kenyon & Luu 1999) and (2) it is expected that a large fraction of KBO surfaces are covered with impact signatures (Durda & Stern 2000).

Antonin Bouchez provided observational assistance at the telescope. David Sprayberry and Gabrelle Saurage provided technical and operational assistance at the telescope.   Table 1: Photometric Data From This Work     Object $`a`$aasemimajor axis $`e`$bbeccentricity $`i`$cc$`i`$ refers to the inclination with respect to the ecliptic. Transforming the inclination to the invariable plane (mean correction 0.6°, Burkhardt 1982) does not have a significant impact on results. $`q`$ddperihelion distance $`R`$ $`BR`$ $`BV`$     \[AU\]  \[deg\] \[AU\] \[magnitudes\] \[magnitudes\]   1999 $`\mathrm{HU}_{11}`$ 43.955 0.061 0.4 41.285 $`23.04\pm 0.06`$ $`2.03\pm 0.16`$ $`1.36\pm 0.16`$  1999 $`\mathrm{CO}_{153}`$ 43.976 0.086 0.8 40.181 $`22.69\pm 0.03`$ $`1.94\pm 0.09`$ $`1.03\pm 0.10`$  1997 $`\mathrm{CU}_{29}`$eeThese objects were also measured by another group, see Table 2. 43.771 0.028 1.5 42.542 $`22.59\pm 0.03`$ $`1.82\pm 0.05`$ $`1.10\pm 0.05`$  1994 $`\mathrm{EV}_3`$eeThese objects were also measured by another group, see Table 2. 42.792 0.045 1.7 40.871 $`23.64\pm 0.06`$ $`1.73\pm 0.15`$ $`1.15\pm 0.23`$  2000 $`\mathrm{FS}_{53}`$ 42.960 0.035 2.1 41.437 $`23.29\pm 0.05`$ $`1.85\pm 0.18`$ $`1.23\pm 0.20`$  1999 $`\mathrm{CM}_{119}`$ 44.623 0.134 2.7 38.640 $`23.58\pm 0.06`$ $`1.78\pm 0.17`$ $`0.90\pm 0.19`$  1999 $`\mathrm{CJ}_{119}`$ 45.647 0.072 3.2 42.375 $`22.92\pm 0.21`$ $`2.07\pm 0.22`$ $`1.38\pm 0.18`$  2000 $`\mathrm{CM}_{105}`$ 42.406 0.068 3.8 39.543 $`22.46\pm 0.03`$ $`1.98\pm 0.23`$ $`1.14\pm 0.23`$  2000 $`\mathrm{CL}_{105}`$ 43.234 0.043 4.2 41.365 $`23.25\pm 0.06`$ $`1.52\pm 0.09`$ $`1.11\pm 0.08`$  1999 $`\mathrm{HT}_{11}`$ 43.602 0.113 5.1 38.670 $`23.05\pm 0.04`$ $`1.83\pm 0.10`$ $`1.14\pm 0.14`$  1999 $`\mathrm{GS}_{46}`$ 44.074 0.009 5.2 43.680 $`22.64\pm 0.02`$ $`1.76\pm 0.07`$ $`1.13\pm 0.09`$  1999 $`\mathrm{CV}_{118}`$ 53.152 0.292 5.5 37.607 $`22.78\pm 0.06`$ $`2.13\pm 0.09`$ $`1.08\pm 0.08`$  1999 $`\mathrm{JD}_{132}`$ 42.996 0.010 10.5 42.546 $`22.26\pm 0.02`$ $`1.59\pm 0.09`$ $`1.11\pm 0.10`$  1999 $`\mathrm{CQ}_{133}`$ 41.517 0.090 13.2 37.784 $`23.07\pm 0.05`$ $`1.35\pm 0.07`$ $`0.75\pm 0.07`$  1999 $`\mathrm{HB}_{12}`$ 55.359 0.411 13.2 32.631 $`22.12\pm 0.02`$ $`1.43\pm 0.04`$ $`0.88\pm 0.06`$  1999 $`\mathrm{CG}_{119}`$ 49.926 0.295 16.6 35.201 $`23.21\pm 0.04`$ $`1.53\pm 0.08`$ $`0.87\pm 0.09`$  1999 $`\mathrm{HW}_{11}`$ 52.422 0.253 17.2 39.168 $`22.93\pm 0.03`$ $`1.33\pm 0.06`$ $`0.80\pm 0.08`$  2000 $`\mathrm{CO}_{105}`$ 47.183 0.132 19.2 40.958 $`22.58\pm 0.18`$ $`1.52\pm 0.18`$ $`0.82\pm 0.06`$  2000 $`\mathrm{CQ}_{105}`$ 57.105 0.389 19.6 34.895 $`23.02\pm 0.05`$ $`1.10\pm 0.09`$ $`0.71\pm 0.09`$  1999 $`\mathrm{CF}_{119}`$ 90.197 0.571 19.7 38.699 $`22.90\pm 0.04`$ $`1.36\pm 0.07`$ $`0.77\pm 0.07`$  1999 $`\mathrm{KR}_{16}`$eeThese objects were also measured by another group, see Table 2. 48.477 0.298 24.9 34.006 $`21.31\pm 0.03`$ $`1.91\pm 0.04`$ $`1.13\pm 0.06`$  2000 $`\mathrm{CG}_{105}`$ 46.595 0.049 27.9 44.334 $`23.40\pm 0.16`$ $`1.17\pm 0.17`$ $`0.78\pm 0.09`$  1998 $`\mathrm{HL}_{151}`$ 40.728 0.092 28.0 36.967 $`23.96\pm 0.08`$ $`1.06\pm 0.13`$ $`0.59\pm 0.19`$  2000 $`\mathrm{AF}_{255}`$ 49.350 0.263 30.9 36.394 $`22.92\pm 0.03`$ $`1.78\pm 0.06`$ $`1.16\pm 0.07`$      Note. — Objects are ordered by increasing inclination. Orbital elements were provided by the Minor Planet Center.       Table 2: Photometric Data From All Other Published Works     Object $`a`$aasemimajor axis $`e`$bbeccentricity $`i`$cc$`i`$ refers to the inclination with respect to the ecliptic. Transforming the inclination to the invariable plane (mean correction 0.6°, Burkhardt 1982) does not have a significant impact on results. $`q`$ddperihelion distance $`BR`$ $`BV`$ Sources   \[AU\]  \[deg\] \[AU\] \[magnitudes\] \[magnitudes\]     1998 $`\mathrm{WX}_{24}`$ 43.612 0.035 0.9 42.097 $`1.79\pm 0.07`$ $`1.09\pm 0.05`$ 6  1994 $`\mathrm{ES}_2`$ 45.883 0.111 1.1 40.779 $`1.65\pm 0.21`$ $`0.71\pm 0.15`$ 1  1994 $`\mathrm{VK}_8`$ 43.133 0.031 1.5 41.800 $`1.68\pm 0.07`$ $`1.01\pm 0.06`$ 6  1997 $`\mathrm{CU}_{29}`$ 43.771 0.028 1.5 42.542 $`2.19\pm 0.48`$ $`1.22\pm 0.10`$ 5, 7  1994 $`\mathrm{EV}_3`$ 42.792 0.045 1.7 40.871 $`2.04\pm 0.20`$ $`1.50\pm 0.19`$ 1  1995 $`\mathrm{WY}_2`$ 46.873 0.128 1.7 40.857 $`1.65\pm 0.21`$ $`1.01\pm 0.16`$ 1, 7  1992 $`\mathrm{QB}_1`$ 44.092 0.072 2.2 40.907 $`1.59\pm 0.14`$ $`0.85\pm 0.10`$ 1, 6, 7  1997 $`\mathrm{CS}_{29}`$ 44.167 0.015 2.2 43.511 $`1.90\pm 0.17`$ $`1.10\pm 0.04`$ 3, 5, 7  1997 $`\mathrm{CQ}_{29}`$ 45.327 0.123 2.9 39.762 $`1.96\pm 0.33`$ $`0.99\pm 0.09`$ 5, 7  1996 $`\mathrm{TK}_{66}`$ 42.848 0.012 3.3 42.313 $`1.70\pm 0.08`$ $`1.04\pm 0.07`$ 6, 7  1996 $`\mathrm{TS}_{66}`$ 44.240 0.134 7.3 38.317 $`1.64\pm 0.14`$ $`1.02\pm 0.05`$ 2, 3, 7  1999 $`\mathrm{DE}_9`$ 55.891 0.423 7.6 32.242 $`1.51\pm 0.04`$ $`0.94\pm 0.03`$ 7  1993 $`\mathrm{FW}`$ 43.650 0.045 7.8 41.665 $`1.59\pm 0.08`$ $`0.95\pm 0.06`$ 1, 5  1997 $`\mathrm{SZ}_{10}`$ 48.190 0.366 11.8 30.546 $`1.79\pm 0.09`$ $`1.14\pm 0.08`$ 6  1998 $`\mathrm{WH}_{24}`$ 46.131 0.110 12.0 41.052 $`1.88\pm 0.31`$ $`0.94\pm 0.04`$ 5, 6  1997 $`\mathrm{QH}_4`$ 42.831 0.084 13.2 39.247 $`1.69\pm 0.12`$ $`1.03\pm 0.10`$ 6, 7  1998 $`\mathrm{SM}_{165}`$ 47.816 0.371 13.5 30.082 $`1.76\pm 0.12`$ $`1.01\pm 0.10`$ 6  1997 $`\mathrm{CR}_{29}`$ 47.310 0.216 19.1 37.082 $`1.36\pm 0.27`$ $`0.67\pm 0.20`$ 7  1996 $`\mathrm{TL}_{66}`$ 84.917 0.587 23.9 35.033 $`1.02\pm 0.11`$ $`0.70\pm 0.04`$ 2, 3, 4, 7  1999 $`\mathrm{OY}_3`$ 43.631 0.166 24.3 36.407 $`1.08\pm 0.02`$ $`0.71\pm 0.01`$ 6  1999 $`\mathrm{KR}_{16}`$ 48.477 0.298 24.9 34.006 $`1.84\pm 0.06`$ $`1.10\pm 0.05`$ 7  1996 $`\mathrm{TO}_{66}`$ 43.364 0.115 27.4 38.372 $`1.07\pm 0.05`$ $`0.69\pm 0.03`$ 2, 3, 4, 7  1996 $`\mathrm{RQ}_{20}`$ 44.036 0.104 31.6 39.467 $`1.54\pm 0.17`$ $`0.96\pm 0.13`$ 7     11footnotetext: Luu & Jewitt 1996 22footnotetext: Jewitt & Luu 1998 33footnotetext: Tegler & Romanishin 1998 44footnotetext: Barucci et al. 1999 55footnotetext: Barucci et al. 2000  66footnotetext: Tegler & Romanishin 2000 77footnotetext: Jewitt & Luu 2001   Note. — Objects are presented in order of increasing inclination. Due to space considerations, Plutinos have been excluded from this table. Orbital elements were provided by the Minor Planet Center.


text:
# 1 Overview

## 1 Overview

I use the word Complexity to refer to the study of complex systems, of which there is no uniformly accepted definition because, well, they are complex. Roughly speaking, one says that a system is complex if it consists of many interacting components (sub-units) and if it exhibits behaviour that is interesting but at the same time not an obvious consequence of the known interaction among the sub-units.

That sounds very vague, especially the use of words like ”interesting” and ”obvious”, but it reflects an evolutionary perspective. For example, a hundred years ago one might have described the study of how a substance changes under heat (phase transitions) as a difficult and interesting problem that required one to deal with systems with a large number of interacting components (atoms). However by now very powerful tools, such as thermodynamics and statistical mechancis, have been developed to deal with such equilibrium systems leading to impressive quantitative agreement between theory and experiment. Though such systems are not commonly referred to as complex, they still provide valuable examples and concepts that have been used in complexity studies.

Current interest has shifted to dynamical systems that are (generally) out-of-equilibrium and thus highly non-linear. Such sytems actually form the bulk of natural phenomena but for which the theoretical tools are as yet poorly developed. Some examples of such complex systems or phenomena are: The economy, the stock-market, the weather, ant colonies, earthquakes, traffic jams, living organisms, ecosystems, turbulence, epidemics, the immune system, river networks, land-slides, zebra stripes, sea-shell patterns, and heartbeats.

There is no single ”Theory of Complexity”, and it is unlikely that there will ever be one. Rather one hopes that apparently different complex systems can be grouped according to some common features that they have, so that intuition and insight gained in studying one can be transferred to another. Thus one of the main aims of complexity studies is to develop concepts, principles and tools that allow one to describe features common to varied complex systems. This leads to exciting interdisciplinary studies because it turns out that ideas developed to handle complex systems in the physical sciences have relevance also for systems in the biological and social sciences, and vice versa!

What are some of the characteristics of complex systems? One often quoted concept is that of emergence, which refers to the appearance of laws, patterns or order through the cooperative effects of the sub-units of a complex system. Thus the emergent phenomena or laws are not an intrinsic property of the sub-units but rather something that is a property of the system as a whole. Simple examples are those of ”temperature” and the ”gas laws”: At the individual microscopic level, none of those make any sense, but they are features of a large system. More sophisticated examples are of ”intelligence” and ”conciousness” – where do they come from ?

Sometimes one sees the phrase ”the whole is more than the sum of its parts”, as a definition of emergence. This again reflects the non-linearity of the system, whereby the output is not proportional to the input, small changes can give rise to large effects, and the non-obvious results that can be produced in a large system.

It is important to realise that the universe consists of many hierarchial levels of complexity linked to each other. Each level has its own emergent patterns and laws: As one goes down from galaxies, solar-systems, planets, ecosystems, organisms, organs, cells, and atoms to quarks, different effective laws emerge. However these laws would not be useful if there was not some degree of universality, that is, one hopes that at each level of complexity the same laws apply to varied systems rather than each following its own tune. It is the apparent universality of the laws of physics, for example, that makes the world comprehensible and gives us faith in its ultimate simplicity. For example, at the atomic level weird quantum mechanics rules, but larger systems are well described according to Newtonian laws, while engineers often use empirical rules, and so do the social scientists.

It appears that nature has chosen to be economical (or is that an illusion on our part?), so that the branching of trees or the air-passages in our lungs, the shape of coastlines or clouds, the form of cauliflower or a mountain range, can be described by fractal geometry: Such shapes are self-similar over a wide range of scales, thus implying scale-invariance, whose hallmark is the appearance of ”power-laws”. In an equilibrium system scale-invariance naturally appears at the critical point of a second-order phase transition, such as that between the liquid and vapour phases of water. However natural systems are out-of-equilibrium and the common appearance of fractals and power-laws in such systems is not as well understaood. Self-organised criticality is the idea that many out-of-equilibrium systems naturally organise themselves, without external tuning or prodding, into a state which is at the threshold between complete disorder and complete order: That is, the system arranges itself into a critical state, and so displays scale-invariance and power-laws.

Living systems are the most complex examples one can think of and it is remarkable how such systems tend in their development towards greater order, organisation and complexity, in contrast to the arrow of time dictated by the Second Law of Thermodynamics. Of course there is no conflict as the increase in disorder and entropy required by the Second Law refers to closed equilibrium systems. Living systems are neither closed nor in equilibrium, but rather use an inflow of energy to drive processes that increase their order (thus decreasing their entropy), and dissipate heat and other waste products that lead to an overall increase in entropy of the universe. One can say that organisms are dissipative structures, and have a tendency towards self-organisation and pattern formation.

Ant-colonies are classic examples of self-organisation. Without a leader (the queen is actually an egg-laying machine) orchestrating everything, and without any of the ants having taken a course in engineering or social science, each ant seems to do its own thing, following a few simple rules that determines its interaction with its environment or its ant-mates. Yet, an incredibly complex and organised society emerges from such an interaction of the many ants. Ant-colonies display remarkable adaptation to changing circumstances, using both feedback mechanisms and parallel analysis of options. In recent years social and computer scientists have taken a keen interest in studying ant colony behaviour in order to help solve problems in their own fields.

Not all systems in nature appear organised or have some pattern to them. Indeed many seem disorderly or ruled by random events. However some of that randomness might only be on the surface. Chaos refers to the property of some non-linear dynamical systems whereby they become extremely sensitive to initial conditions and display long-term aperiodic behaviour that seems unpredictable. Though chaotic behaviour might appear essentially random, there is actually hidden order, apparent only in ”phase-space” rather in ordinary space. Furthermore, many chaotic systems show universality in their approach to chaos, giving one some predictive power. Thus discovering that some random-like events are actually chaotic means one has uncovered a simple determinstic basis for the system and so enabled its understandability.

Often one encounters debates between reductionism and holism. Reductionists like to get right down to the bottom, meaning they are interested in the basic sub-units that make up the whole and believe that that is where all that is of interest lies, the whole itself being just a complicated and uninteresting consequence of the fundamental laws applied to a large system. In short, knowing the microscopic explains all to the extreme reductionist. Particle physicist are such, and in more recent times some molecular biologists involved with genomics are another example. While it is undoubtedly true that knowledge of the microcomponents of a system and the basic interactions among those is essential for us to progress, it is also a fact that such knowledge by itself is insufficient to predict all the diversity and novelty that can arise in a large system. (Take for example the task of predicting superconductivity from Schrodinger’s equation – it is a problem that required much effort after the fact–one knew what to look for. Similarly knowing the whole genome code is not going to predict for us every feature of an organism or a society).

The problem of precisely deducing the whole (of a large system) from its parts is at least two-fold. Firstly it is a computational problem. Problems with a large number of degrees of freedom are too complicated for exact solutions, and for systems far from equilibrium, as complex systems are, they are also not solvable by the probabilistic averaging methods used for equilibrium systems. In recent years the growth of computer power at low cost has produced the first tool that allows large systems to be simulated or solved numerically. However this brings the second problem: Often one does not have full knowledge of the fundamental dynamics, or the initial conditions, or the problem is still too complicated to be handled directly even by computers.

Often what is required is some guesswork or intuition to reduce the actual problem to a simpler model which can then be tested on a computer. Computer simulations of simplified models let one test assumptions quickly, and when the results appear similar to the real world one can take it as plausible validity of the model and the assumptions. Qualitative similarities of course do not constitute a proof, because other models with different assumptions might give similar results, but at least the insight gained helps one to make further guesses and tests in a particular direction rather than being lost in a mess of detail. In fact one of the most important lessons computer simulations have taught us is that a large system with very simple local rules can give rise to collective behaviour of great complexity and variety, showing on the one hand that complex phenomena need not require complicated rules, but at the same time reminding us how difficult it is (without computers) to deduce the emergent behaviour from the sub-units and their interactions.

Thus studying the whole is as interesting as studying its parts, as novel structures and emergent laws arise at each level of complexity. The condensed matter physicist studying superconductivity is not going to be replaced by the string theorist, and neither is the ecologist going to be become obsolete because of the molecular geneticist. Explaining dynamic patterns, order and emergent laws of a complex system by understanding the organising principles among the sub-units is what might be called holism, the counterpoint to reductionism.

## 2 Examples

Let us briefly look at examples to illustrate some of the points above. These are just apetisers, details are in the main course .

### 2.1 Schooling of Fish

Try out the applet at Ref.. Does it not look like a very realistic simulation of fish swimming? The motion of each individual fish is not scripted right from the beginning but rather each individual follows just three simple local rules: cohesion, alignment and separation. Each of the rules is sensible from the biological perspective and so the model is plausible. What is remarkable of course is how the realistic and complicated collective behaviour emerges from the few simple local rules. There is no leader and none of the individuals has a global plan or perspective, (the motion is not orchestrated from the beginning).

This example is an example of self-organisation. Other examples are the herding behaviour of humans, say for example in the stock-market, and the alignment of magnetic spins to form a ferromagnet. Many such examples are studied in .

### 2.2 Bacterial Colonies

Look at the picture of a bacterial colony. It shows a branching structure, which has the property that if one zooms into any region, that part looks similar to the whole. The bacterial colony is an example of a random fractal. Exact fractals appear the same at different magnification scales while random fractals appear only statistically similar at different magnification scales.

Fractals are ubiquitous in Nature. Another example is the branching network of air-passages in the human lung. The advantages of such a structure are an increase in surface to volume ratio which maximises functional efficiency while minimising material and space costs.

The word ”fractal” itself means more than just self-similarity at different scales. It also implies a fractal dimension.

### 2.3 Forest Fires

Look at the figures which show the number of fires as a function of the area burnt, in different regions of the United States of America and Australia. On a log-log plot one sees that the data is well approximated by straight lines, meaning that the number of fires as a function of area is a power law:

$$NA^\alpha $$
(1)
with $`\alpha 1.31.5`$. It is important to note that the straight line fits are for a wide range of the data (one can always fit a straight line to a small range), and for different geographical regions. This suggests a universality in the phenomena that requires a explanation.

Power laws are observed in many other natural phenomena such as earthquakes and solar flare activity. It has been suggested that these phenomena are examples of self-organised criticality, that is, the systems are attracted to a state which is between that of total order and total disorder. The word ”critical” is borrowed from well-studied thermal equilibrium systems that undergo second-order phase transitions at critical points and display power laws. However in the case of forest fires (or earthquakes etc.), the systems are far from equilibrium and the power-law behaviour, that is criticality, does not require fine-tuning — it is self-organised.

Given the complexity fo the actual systems, it is impossible at present to study those systems from first principles. Rather one studies simplified models to check whether power-laws emerge naturally.

Power laws imply a self-similarity at different scales, so it is natural to suppose that self-organised criticality might be the common dynamical mechanism behind the wide occurence of fractal structures in nature.

Power laws also appear in social contexts, for example the frequency of occurence of words in a literary text, and the magnitude of wars.

### 2.4 The Double Pendulum

Most people are familiar with the simple pendulum: A small heavy object suspended at the end of a thin light string and set into oscillation. For small oscillations (and in the absence of air friction) the motion is that of a ”harmonic oscillator”: That is, periodic motion with a period proportional to the sqaure-root of the length of the pendulum.

For larger oscillations, the motion of the simple pendulum is still periodic but no longer given by a simple formula. Indeed, for large oscillations the equations governing the motion fo the pendulum are non-linear in contrast to the linear equations for small oscillations. However while the equations are nonlinear, the motion is still regular and predictable.

A double pendulum consists of two simple pendulums in tandem: One attaches a single pendulum to the end of another! The equations of otion are again non-linear for large oscillations but now the motion becomes quite irregular and very sensitive to the initial conditions. This kind of behaviour is the hallmark of chaos. See the simulation of the double pendulum in the references.

Chaos occurs in many nonlinear systems and it implies that even systems with a few degrees of freedom, and hence naively simple, can show complicated behaviour which is essentially unpredictable on long time scales. However chaos is very different from randomness: The former arises in perfectly deterministic systems while the later is intrinsically nondeterministic, and the distinction between the two at the practical level can be seen by looking at the ”phase space” of the system, as we shall see later.

### 2.5 The Leopard Spots

How did the spots on a leopard, the stripes of the zebra or tiger, or the patterns on sea-shells, form? Is there some simpel general framework which can explain these beautiful patterns? Yes, these patterns are called Turing structures, named after the Bristish mathematician who came up with a model to explain such structures.

Some chemical systems, when maintained far from equilibrium display oscillatory behaviour. When the chemical waves of such a reaction are allowed to diffuse through a medium at different rates andthe resulting pattern stabilised, one obtains Turing structures.

Many other interesting far from equilibrium systems show cyclical behaviour.

## 3 Summary

The aim of the Complexity course in Ref. is two fold:

(1) To provide the student with a relatively gentle introduction to the concepts mentioned above so that they can continue on their own in greater depth.

(2) To broaden the students horizons by introducing to them the interdisciplinary nature of complex systems studies, which exemplifies in a concrete sense the often quoted ”unity of knowledge”.

As in all popular fields, one often finds in descriptions of ”Complexity” misleading or hyped statements, and metaphorical deviations. So some caution and cynicism is required in filtering the raw data from a search, especially over the Internet.

In conclusion, explore the Web-book , the links and exercises therein, and hope for some enlightenment.

## 4 References

1. Complexity: A Web-book at http://staff.science.nus.edu.sg/$``$parwani/complexity.html
2. The dreams of reason : the computer and the rise of the sciences of complexity, by Heinz R. Pagels.
3. The Quark and The Jaguar, by M. Gell-Mann.
4. Web of Life, by F. Capra
5. How Nature Works, by P. Bak
6. The Self-Made Tapestry, by P. Ball
7. The Computational Beauty of Nature, by G. Flake
8. Butterfly Economics, by P. Omerod
9. A Fish schooling applet at
http://www.codepuppies.com/$``$steve/aqua.html
10. The double pendulum applet at
http://www.maths.tcd.ie/$``$plynch/SwingingSpring/doublependulum.html


text:
# Peculiarities of anharmonic effects in the lattice thermodynamics of fcc metals

## Abstract

Explicit expressions for anharmonic contributions to the thermodynamic properties with allowance for higher-order phonon-phonon interactions for closed-packed crystals are given, and the calculations for some fcc metals near the melting (Ir, Rh) and martensite phase transition (Ca, Sr) points are carried out. A detailed comparison of anharmonic and electron contributions to the heat capacity of these metals is carried out. The computational results for high-temperature heat capacity agree well with the available experimental data.

Institute of Metal Physics, 620219 Ekaterinburg, Russia

Russian Research Centre ”Kurchatov Institute”, 123182 Moscow, Russia

Corresponding author: M. I. Katsnelson, e-mail: Mikhail.Katsnelson@usu.ru

PACS: 63.20.-e

Keywords: anharmonic thermodynamics, martensite phase transition, heat capacity, computation

Anharmonic effects (AE) in the lattice dynamics and thermodynamics have been extensively discussed in the literature \[1-5\]. However the quantitative information on AE, both experimental and theoretical, is rather poor so far. The experimental difficulties are due to the fact that at high temperatures ($`TT_m`$, where $`T_m`$ is the melting temperature) when the AE in the thermodynamic properties become noticeable, it is usually difficult to identify their contribution to the heat capacity and thermal expansion from the contributions of vacancies and other thermally excited lattice defects \[2,5-8\]. The theoretical calculations, on the other hand, are based, as a rule, on rough models and, therefore, do not give any definite information on the AE in thermodynamic of real crystals. At the same time this information is needed for understanding the high temperature properties of solids. For example, it is important to clear up what is the role of self-anharmonic effects associated with phonon-phonon interactions, and what is the role of quasiharmonic ones which are only due to the change in the volume because of thermal expansion. Simple estimations  show that the AE contribution to heat capacity at high temperatures of the order of the electron contribution to heat capacity and has the same (linear) temperature dependence; the real ratio of these two contributions essentially depend on the features of the electron structure near the Fermi level $`E_F`$. However this problem has not been studied yet. All these problems are of special importance for refractory metals ($`T_m>2000`$ K) and for the metals with high-temperature martensite transitions which are usually accompanied by essential increase in the AE in the lattice dynamics . In the present work the AE features (in the comparison with the electron contributions) in the thermodynamical properties of refractory fcc metals Ir and Rh and in the bcc and fcc phases of Ca and Sr near the points of structural transformation are studied. The calculations were made in terms of microscopic models used in  which describe a wide range of the lattice properties of these metals.

The initial Hamiltonian of phonon subsystem is

$$H=\underset{\lambda }{}\left(\frac{P_\lambda ^2}{2M}+\frac{M\omega _\lambda ^2Q_\lambda ^2}{2}\right)+\underset{n=3}{\overset{\mathrm{}}{}}H^{(n)}$$
(1)
$$H^{(n)}=\underset{\lambda _1\mathrm{}\lambda _n}{}\frac{\mathrm{\Phi }^{(n)}(\lambda _1,\mathrm{},\lambda _n)}{n!}Q_{\lambda _1}\mathrm{}Q_{\lambda _n}$$
(2)
where $`M`$ is the mass of atoms; $`\omega _\lambda `$ is the phonon frequency, $`\lambda 𝐪\xi `$ where $`𝐪`$ – is the wave vector; $`\xi `$ is the branch number, $`Q_\lambda `$, $`P_\lambda `$ are the normal coordinates in the harmonic approximation and the corresponding pulses, $`\mathrm{\Phi }^{(n)}`$ are amplitudes of intermode interactions. In the leading order of the parameter $`\eta T/E_{at}`$ where $`E_{at}`$ is the energy of the order of binding energy, the basic contribution of AE to free energy $`F`$ is determined by the second order by $`H^{(3)}`$ and the first order by $`H^{(4)}`$

$$F_{an}=F^{(3)}+F^{(4)}$$
(3)
where at $`T\mathrm{\Theta }_D`$ ($`\mathrm{\Theta }_D`$ is the Debye temperature),

$$F^{(3)}=\frac{T^2}{12M^3}\underset{\lambda \mu \nu }{}\frac{|\mathrm{\Phi }^{(3)}(\lambda ,\mu ,\nu )|^2}{\omega _\lambda ^2\omega _\mu ^2\omega _\nu ^2},$$
(4)
$$F^{(4)}=\frac{T^2}{4M^2}\underset{\lambda \mu }{}\frac{\mathrm{\Phi }^{(4)}(\lambda ,\lambda ,\mu ,\mu )}{\omega _\lambda ^2\omega _\mu ^2};$$
(5)
as usually, summing over phonon quasimomenta is fulfilled taking into account the conservation law for $`\lambda (\xi ,𝐤)`$, $`\mu (\eta ,𝐪)`$, $`\nu (\zeta ,𝐤+𝐪)`$. The explicit microscopic expressions for $`\mathrm{\Phi }^{(3)}`$ and $`\mathrm{\Phi }^{(4)}`$ are given in .

The heat capacity of the crystal at $`T\mathrm{\Theta }_D`$, without the standard harmonic contribution $`C_V^{harm}=3R`$ (where $`R`$ is the gaseous constant), has the form:

$$\mathrm{\Delta }C=C_PC_V^{harm}=(C_PC_V)+C_V^{an}+C_V^e+C_V^d$$
(6)
where $`C_V^{an}`$ are the anharmonic contributions to the heat capacity

$$C_V^e(T)=\frac{R}{T}\underset{\mathrm{}}{\overset{\mathrm{}}{}}dE\left[\frac{f(E)}{E}\right](E\mu )^2N(E)$$
(7)
is the electron heat capacity ($`N(E)`$ is the electron density of states, is the chemical potential, $`f(E)`$ is the Fermi distribution function), $`C_V^d`$ is the contribution of defects; the calculation formula for $`C_PC_V`$ in the quasiharmonic approximation is presented in . The density of states $`N(E)`$ was calculated by the FP-LMTO method in the local approximation for the density functional (LDA) . To estimate the contribution of defects to the heat capacity a common approximation of independent monovacancies was used with

$$C_V^d(T)=R\left(\frac{E_v}{T}\right)^2\mathrm{exp}\left(S_v\frac{E_v}{T}\right)$$
(8)
where $`E_v`$, $`S_v`$ are the energy and entropy of monovacancy formation, respectively.

The calculations of AE in the Ir lattice dynamics, made in , show that near the melting temperature $`T=T_m`$ the temperature dependences of the phonon frequency and phonon damping drastically change over the Brillouin zone reaching the values of the order of 15-20%. This value is high enough so that the question on the role of higher orders in the anharmonic perturbation theory could be set. The theory of self-consistent phonons, often used for this purpose (see, e.g., ) is not adequate in this case as it only accounts for a part of $`H^{(2n)}`$ contributions to (1) and neglects the $`H^{(2n+1)}`$ contributions. At the same time, the calculations  show that in Ir the contributions of three-phonon processes to the temperature shifts of the phonon frequencies always dominate over the contribution of four-phonon processes. In the recent work  another approximation is proposed which accounts both the three-phonon and four-phonon processes, and, as stated, reasonably describes the thermodynamics of crystals with the Lennard-Jones type interaction. The applicability of this approximation in the case considered seems somewhat doubtful as it is based on the replacement of anharmonic frequency shift by a certain value average over the Brillouin zone while the calculations in  show that the q-dependence of AE is essential. Moreover, it is known that the approximations adequate for the description of thermodynamics of systems with “hard” potentials of Lennard-Jones type are not usually applicable for metals as they have relatively “soft” potentials .

To calculate the effects of higher orders by the Hamiltonian $`H^{(3)}`$ the approximation similar to that of Gell-Mann and Brueckner (“bubble approximation”) in the theory of electron liquid  was used. It can be shown that the accurate expression for the contribution of three-phonon processes to the free energy at $`T\mathrm{\Theta }_D`$ has the form:

$$F^{(3)}=\frac{T}{3}\underset{0}{\overset{1}{}}\frac{d\alpha }{\alpha }\underset{\lambda }{}\left[G_\lambda (\omega =0)\mathrm{\Sigma }_\lambda (\omega =0)\right]|_{H^{(3)}\alpha H^{(3)}}$$
(9)
where $`G_\lambda `$, $`\mathrm{\Sigma }_\lambda `$ are the Green function and self-energy part for phonon (see ) connected by the Dyson equation

$$G_\lambda ^1(\omega )=\frac{\omega ^2\omega _\lambda ^2}{2\omega _\lambda }\mathrm{\Sigma }_\lambda (\omega ).$$
(10)
The calculations  show that the AE in fcc metals have pronounsed peaks as functions of $`𝐪`$ at points X and L. The largest contribution to the free energy comes from the “bubble” diagrams; this situation is similar to the electron gas of high density with the only difference that in the latter case the “dangerous” region is $`𝐪0`$ . Taking this into account, in (9)-(10) we may restrict ourselves to the expression of the lowest (second) order of the perturbation theory for $`\mathrm{\Sigma }_\lambda `$. Then we have

$$F^{(3)}=\frac{T}{6}\underset{\lambda }{}\mathrm{ln}\left[1\frac{T}{2M^3}\underset{\mu \nu }{}\frac{|\mathrm{\Phi }^{(3)}(\lambda \mu \nu )|^2}{\omega _\lambda ^2\omega _\mu ^2\omega _\nu ^2}\right]$$
(11)
Unlike the bcc metals the fcc metals have no soft modes in the phonon spectrum, which makes the four-phonon processes relatively small comparing with the three-phonon ones . Therefore for $`F^{(4)}`$ it is sufficient to use the simplest corresponding to the Hartry-Fock approximation (the simplest procedure of disengagement in averaging Hamiltonian $`H^{(4)}`$). Then for $`T\mathrm{\Theta }_D`$ we have

$$F^{(4)}=T\underset{\lambda }{}\left[\mathrm{ln}\left(1+\frac{\mathrm{\Delta }_{4\lambda }}{\omega _\lambda }\right)\frac{1}{2}\frac{\mathrm{\Delta }_{4\lambda }}{\omega _\lambda }\right],$$
(12)
where

$$\mathrm{\Delta }_{4\lambda }=\frac{T}{4\omega _\lambda }\underset{\mu }{}\frac{\mathrm{\Phi }^{(4)}(\lambda ,\lambda ,\mu ,\mu )}{M^2\omega _\mu ^2}$$
(13)
is the corresponding correction to the phonon frequency.

The computational results are shown in Figs. 1-5. One can see that the AE higher order contributions to the free energy of Ir and Rh are not large (do not exceed 20%). In Ir and Rh the anharmonic contribution to heat capacity is small comparing with the electron contribution while in the sp-metals (as the hypothetical fcc phase K shown in Fig. 3) these contributions are comparable. Among the metals studied only for Ir the detailed experimental information is available, the theoretical data being in good agreement with these results (Fig. 4). In the estimation of $`C_V^d`$ according to (8) the the first principles calculation result from $`E_v=1.88`$ eV  was used for $`E_v`$, and $`S_v=2.5`$ for $`S_v`$ . It follows from Fig. 5 that in the difference of heat-capacity for bcc and fcc phases of Ca (and similarly Sr) a strong compensation for different contributions – electron, anharmonic ones – and the contribution of $`C_PC_V`$ takes place. Thus, for the discussion of thermodynamics of high temperature structure transformations in metals the microscopic calculations of corresponding values including the anharmonic effects are needed.

The authors are grateful to K. Yu. Khromov for his assistance in performance of calculations and M. N. Khlopkin for helpful discussion of experimental results.

The work is sponsored by the Russian Basic Research Foundation, grant 01-02-16108, and the Netherlands Research Organization, grant NWO 047-008-16.