text:
# A model metal potential exhibiting polytetrahedral clusters

## I Introduction

The understanding of the structure of metal clusters has seen many developments in recent years,Alonso (2000); Johnston (2002) however there is still much to be learnt. From the theoretical perspective, for example, only relatively recently have global optimization techniques become sufficiently powerful to find the most stable structures of metal clusters with up to 100 atoms,Wales and Scheraga (1999) even when described by relatively simple many-body potentials. This has led to many interesting new structures being revealedGarzón et al. (1998); Doye and Wales (1998); Wilson and Johnston (2000); Soler et al. (2001); Wang et al. (2001); Michaelian et al. (2002); Lai et al. (2002); Doye and Hendy (2003) that go beyond those classified for pair potentials, but there are still probably many classes of structure that remain undiscovered. Furthermore, the many-body nature of the bonding makes it difficult to unravel the relationship between the observed structures and the form of the interactions, although there has been some interesting recent progress.Soler et al. (2000); Baletto et al. (2002); Doye and Hendy (2003) This situation contrasts with that for pair potentials, where, for example, the effects of the width of the potential wellDoye et al. (1995); Doye and Wales (1997) and oscillations in the potentialDoye et al. (2001); Doye and Wales (2001); Doye et al. (cond-mat/0205374) on the cluster structure have been mapped out and rationalized.

One of the more interesting recent discoveries has been that some monatomic metal clustersDassenoy et al. (2000); Wang et al. (2001) can exhibit structures with polytetrahedral order,Nelson and Spaepen (1989); Sadoc and Mosseri (1999) i.e. the whole of the cluster can be naturally divided up into tetrahedra with atoms at their vertices.clo  In the experiments of Dassenoy et al. on cobalt clusters, the electron microscopy images and diffraction patterns could not be explained by any of the usual structural forms for a cluster, but instead were consistent with a polytetrahedral structure.Dassenoy et al. (2000) Furthermore, the magic number at $`N=61`$ in the mass spectrum of strontium clustersBréchignac et al. (2000); Wang et al. (2001) was assigned to a polytetrahedral structure,Wang et al. (2001) and the next magic number at $`N=82`$ also coincides with a polytetrahedral magic number.Doye and Wales (2001) Similar magic number patterns have also been found for some rare-earth metals.Bréchignac et al. (1993)

These results seem somewhat surprising since ordered polytetrahedral structures are restricted to alloys in bulk, e.g the crystalline Frank-Kasper phasesFrank and Kasper (1958, 1959); Shoemaker and Shoemaker (1988); bin  and quasicrystals,Shechtman et al. (1984) structural models of which often invoke substantial polytetrahedral order.Elser and Henley (1985); Sachdev and Nelson (1985); Audier and Guyot (1986) However, as we will see, there are good reasons why polytetrahedral structures are more likely to be seen for monatomic systems in clusters.

Figure 1(a) illustrates that all space cannot be packed with regular tetrahedra. (For this reason polytetrahedral structures are sometimes said to be frustrated.Sadoc and Mosseri (1999)) Therefore, the nearest-neighbour pair distances in a polytetrahedral structure cannot all take the same value. For example to close the gaps associated with the icosahedral packing of regular tetrahedra, the tetrahedra must be distorted such that the edges of the tetrahedra on the surface are 5% longer than the radial edges. For a monatomic system this strain usually leads to an energetic penalty, although for a binary system, if a smaller atom of the right size is placed in the centre of the icosahedron, all pair distances could still take their ideal values.

For clusters of sufficiently small size the energetic penalties associated with the strain inherent to polytetrahedral structures can often be more than compensated by a favourable surface energy. The surface of a polytetrahedral cluster is triangulated leading to high coordination numbers, whereas clusters based on a close-packing scheme may have to expose $`\{100\}`$ faces, which involve atoms with lower coordination number, and hence typically have a higher surface energy.

The 13-atom icosahedral structure, which is very common for systems interacting with both pair potentials and many-body metal potentials, illustrates this principle. The polytetrahedral packing of the icosahedron can be continued by adding capping atoms above the faces and vertices.Northby (1987) This “anti-Mackay” overlayer leads to polyicosahedral structures, where each atom in the interior of the cluster has local icosahedral coordination. So, the 19-atom double icosahedron is also common, but as the size increases these structures become less likely to be observed. For example, for Lennard-Jones clusters they are lowest in energy up to $`N=30`$.Northby (1987) The reason that these clusters become disfavoured is clear from Fig. 1(a). As further regular tetrahedra are added to these packings the gaps involved become larger and larger, and hence the strain energy associated with polyicosahedral packings increases very rapidly with size.

For the above polytetrahedral clusters, there are five tetrahedra around each edge in the interior of the cluster. To balance out the increasing gaps in this type of packing, the bulk Frank-Kasper phases also involve sites where six tetrahedra share a common edge—a negative disclination is said to run along this edge.dis  Although the local distortion required to remove the overlap when six regular tetrahedra are placed around a common edge is larger, the introduction of disclination lines leads to a decrease in the overall strain.

The coordination polyhedra involved in the Frank-Kasper phasesFrank and Kasper (1958, 1959) are depicted in Figure 1(b). Local icosahedral coordination involves no disclinations, but as the coordination number increases negative disclinations must be introduced. Note a disclination line can only end by forming a loop or exiting at a surface. So for the Z14 coordination polyhedron the disclination passes through the central atom, and for the Z15 and Z16 coordination polyhedra the central atom acts as a node for three and four disclinations, respectively. Polytetrahedral packings can therefore be described by a network of disclination lines threading an icosahedrally-coordinated medium. In Frank-Kasper crystals this disclination network is ordered and periodic,maj  and is mediated by the Z14, Z15 and Z16 Frank-Kasper coordination polyhedra.

Polytetrahedral clusters that involve disclinations have previously been found for two very different pair potentials. Firstly, for clusters interacting with a sufficiently long-ranged Morse potential the strain involved in a polytetrahedral structure can be tolerated, and ordered polytetrahedral clusters are found up to $`N70`$.Doye and Wales (1997) Secondly, a modified Dzugutov potential favours clusters that are based on the C14 and C15 Frank-Kasper phases up to at least $`N=250`$.Doye and Wales (2001) This potential has a maximum in the potential that has been chosen to coincide with the next-nearest neighbour distance in close-packed structures, and somewhat resembles the Friedel oscillations that can occur in effective pair potentials derived for metallic systems.Hafner (1998); Moriarty and Widom (1997)

These two examples are for somewhat unusual potentials. However, polytetrahedral structures are perhaps more likely for metals, because the nature of the metallic potentials makes them potentially less sensitive to the strain in the individual pair distances,Soler et al. (2000) for reasons that will be discussed in more detail in Section II.1. For example, similar structures to those of the Morse clusters have also been found for clusters described by a generalized set of metallic interactions.Cune and Apostol (2000) Furthermore, in this paper I report that the lowest-energy structures of clusters interacting with an aluminium glue potentialErcolessi and Adams (1994) are also polytetrahedral up to at least $`N=190`$.

The paper is organized as follows. In Section II the main focus is on the form of the aluminium potential and the implications of this form for cluster structure. In Section III I describe in detail the structures of the putative global minima, and in Section IV I take Al61 as a case study for understanding the reasons for this system’s preference for polytetrahedral structures. Then, in Section V I briefly examine the size evolution of the cluster structure at larger sizes. The emphasis in the paper is less on the clusters as a realistic model of aluminium clusters that will hold up to detailed scrutiny—this is very difficult for any empirical potential—and more on using these clusters as a model system for exploring the types of polytetrahedral structure that metal clusters might exhibit, and the reasons for their formation. However, in Section V I do compare the results with the available experimental information on the structure of aluminium clusters.

## II Methods

### II.1 Potential

To model the aluminium clusters I use a glue potentialErcolessi et al. (1988) that has been constructed using the force-matching method.Ercolessi and Adams (1994) The potential energy is given by

$`E`$ $`=`$ $`E_{\mathrm{pair}}+E_{\mathrm{glue}}`$ (1)
$`=`$ $`{\displaystyle \underset{i<j}{}}\varphi \left(r_{ij}\right)+{\displaystyle \underset{i}{}}U\left(\overline{\rho }_i\right),`$ (2)
where $`\varphi (r)`$ is a short-ranged pair potential, $`U(\overline{\rho })`$ is a many-body glue function and $`\overline{\rho }_i`$ is defined as

$$\overline{\rho }_i=\underset{j}{}\rho \left(r_{ij}\right),$$
(3)
where $`\rho (r)`$ is an “atomic density” function. These three functions have been fitted to match the forces produced by first-principles electronic structure calculations for a large set of configurations, including those that correspond to surfaces, clusters, liquids and crystals.Ercolessi and Adams (1994) That the potential has been designed to model a wide variety of atomic environments, such as occur in clusters, is particularly important for the current application, as opposed to a potential that has just been fitted to bulk data. This flexibility has also led the potential to be used in a wide variety of applications, often with considerable success. Examples include clusters,Sun and Gong (1993) wires,Gülseren et al. (1998) surface thermodynamics,Di Tolla et al. (1995); Raphuthi et al. (1995) grain boundaries,Rittner et al. (1996) bulkSandberg et al. (2002) and surfaceTrushin et al. (2001) diffusion and melt solidification.Liu et al. (2001)

The functions $`\varphi (r)`$, $`\rho (r)`$ and $`U(\overline{\rho })`$ are displayed in Figure 2. The forms of these functions are not assumed, but are determined by the fitting process. This contrasts with other empirical potentials, such as the Gupta or Sutton-Chen potentials. With these latter types of potentials it would be impossible, for example, to reproduce the double well structure of $`\varphi (r)`$, which also occurs for a sophisticated ‘bulk’ potential.Mishin et al. (1999)

The shape of these functions have some fairly straightforward consequences for the cluster structure. $`U(\overline{\rho })`$ has a single minimum of depth $`2.6571`$eV at $`\overline{\rho }=1.0`$. This term provides the dominant contribution to the energy, so it is important that each atom has a value of $`\overline{\rho }`$ as close as possible to one. Of course, the values associated with the surface atoms will be significantly less than this. There is therefore a tendency for surface distances to contract so that the atoms achieve as high a value of $`\overline{\rho }`$ as is possible. This effect will be balanced by the energetic penalty associated with the resulting compression of the interior and the reduction in the pair energy of the surface atoms. Furthermore, the minimum of the glue function is rather broad, and so the energy loss associated with $`\overline{\rho }`$ values that deviate from one is not so severe. Thus, surface contraction will not be as pronounced as for some other metallic systems. For example, the lead glue potential,Lim et al. (1992) whose cluster structures have recently been analysed,Hendy and Doye (2002); Doye and Hendy (2003) has a stronger dependence of $`U`$ on $`\overline{\rho }`$.

$`\rho (r)`$ decreases roughly linearly with distance until $`r4`$Å, after which it decays more slowly to zero at the cutoff distance. As the value of $`\rho (r)`$ in the latter regime is always small, the major contribution to the $`\overline{\rho }`$ values of each atom comes from nearest neighbours, and hence the average value of $`\overline{\rho }_i`$ will be maximized by clusters with a large number of nearest neighbours.

$`\varphi (r)`$ has an interesting double-well structure. The first minimum occurs at $`r_{\mathrm{min1}}=3.095\mathrm{\AA }`$, and has a depth of $`0.0587`$eV. This is separated by a maximum at $`3.669\mathrm{\AA }`$ from a second minimum at $`r_{\mathrm{min2}}=4.303\mathrm{\AA }=1.390r_{\mathrm{min1}}`$ of slightly greater depth, $`0.0633`$eV. As the number of next-nearest neighbours will be significantly larger than nearest neighbours (except for the smallest clusters), the largest contribution to $`E_{\mathrm{pair}}`$ comes from the second neighbour shell.

The ratio of the positions of the two minima would initially suggest that they are almost ideally placed to coincide with the first and second shells for closed-packed structures, which have a ratio of $`\sqrt{2}`$. However, the nearest-neighbour distance is determined more by the glue part of the potential. For example, for an atom with twelve nearest neighbours, these neighbours would need to be at a distance for 2.611 Å to achieve a $`\overline{\rho }`$ value of 1. If the effects of next neighbours are also included, this distance does not not need to be as short, and so the nearest neighbour distance for the equilibrium face-centred-cubic (fcc) crystal, $`r_{\mathrm{nn}}^{\mathrm{xtal}}`$, is 2.851 Å.Ercolessi and Adams (1994) Compared to this distance the position of the second minimum is not so advantageous for closed-packed structures; $`r_{\mathrm{min2}}=1.509r_{\mathrm{nn}}^{\mathrm{xtal}}`$.

One of the important features of the glue energy for the current study is that it does not directly depend on the distribution of pair distances or the coordination number of an atom, but only on the $`\overline{\rho }_i`$ values.Soler et al. (2000) This situation contrasts with that for a pair potential, and makes metallic systems potentially less sensitive to internal strains. This feature has already been cited as an important factor in the stability of disordered clusters for a variety of metallic systems.Soler et al. (2000) The same principle applies for ordered polytetrahedral clusters, which have both strain and high coordination numbers. For example, the ideal glue energy can be achieved for an atom with coordination number higher than twelve simply by expanding the nearest-neighbour distances for this atom with respect to $`r_{\mathrm{nn}}^{\mathrm{xtal}}`$ until $`\overline{\rho }=1`$ for that atom. Similarly, there is no reason why an atom with a wide distribution of nearest-neighbour distances cannot also achieve a $`\overline{\rho }`$ value of one.

### II.2 Global Optimization

The global optimization of the aluminium clusters was performed using the basin-hopping Wales and Doye (1997); Wales and Scheraga (1999) (or Monte Carlo minimization Li and Scheraga (1987)) approach. This method has proved particularly successful in locating putative global minima for a wide variety of cluster systems. As the results are extended up to 190 atoms, a considerable computational effort was required. It proved particularly important to supplement the standard unbiased runs from random starting points for each size, with runs started from low-energy minima of nearby sizes with the appropriate number of atoms added or removed. These latter runs were applied iteratively until no further new global minima were located.

It should be noted that, of course, there is no guarantee that I have been able to locate the true global minima, and the probability that a global minimum has been missed will increase with cluster size, as the size of the search space, and hence the number of minima,Tsai and Jordan (1993); Stillinger (1999); Doye and Wales (2002a) increases exponentially with $`N`$. However, from examination of the statistics of how often independent runs locate the same lowest-energy minimum, I am confident that virtually all of the putative global minima up to $`N=100`$ cannot be bettered, and that beyond this size the putative global minima, if not truly global, are very good estimates. In this latter size range the runs seeded from low-energy clusters of other sizes become more important. Therefore, global minima are most likely to be missed when they are structurally different from those of nearby sizes. For example, the lowest-energy structure found in a basin-hopping run for Al147 was at $`436.6892`$eV, even though the 147-atom Mackay icosahedron has an energy that is $`0.0133`$eV lower. However, this is the only example I found where reoptimization of a structure that is known to often be particularly stable bettered the global optimization results. I also tried to construct large polytetrahedral structures based upon some of the structural motifs present for $`N<100`$, but again this did not lead to any improvements.

## III Global minima for $`N190`$

The energies and point groups for the putative global minima are given in Table 1. Point files are available online at the Cambridge Cluster Database.Web  The energies of the global minima are represented in Figure 3 in such a way that makes particularly stable clusters stand out. A selection of clusters are depicted in Figure 4 that are either particularly stable or have some interesting structural feature. Most of the intermediate sizes involve structures similar to those in this figure, but normally with a somewhat incomplete surface structure.

The only other study that has examined the structures of aluminium clusters interacting with the present potential studied clusters with 13, 43, 55 and 147 atoms.Sun and Gong (1993) For $`N`$=13, 55 and 147, Sun and Gong considered high symmetry structures and found the Mackay icosahedra to be lowest in energy, in agreement with the results here. For $`N`$=43, they performed simulated annealing, generating a lowest-energy structure that is $`0.0562`$eV above the global minimum, and which is in fact the second lowest-energy isomer. They also claimed that this cluster had a ‘glassy’ structure, but in common with the global minimum it in fact has an ordered polytetrahedral structure.

The structures with from 3 to 13 atoms are those typically seen for isotropic potentials. Most are on a polytetrahedral growth sequence leading to the 13-atom icosahedron, except for the six-atom octahedron and the 8-atom dodecehedron. Beyond $`N=13`$ growth occurs around the icosahedron with atoms being added above the faces and vertices of this structure in the so-called anti-Mackay overlayer.Doye et al. (1995) This growth maintains the polytetrahedral character of the clusters and leads to a series of interpenetrating 13-atom icosahedra at $`N=19`$, 23, 26, 29, 32 and 34. For example, the 34-atom structure involves seven interpenetrating icosahedra with their centres in the shape of a pentagonal bipyramid. This growth sequence is very similar to that for Lennard-Jones clustersNorthby (1987) although in a few instances slightly different sites in the anti-Mackay overlayer are occupied.

Beyond this size is where things become different from any previous system. For example, for Lennard-Jones clusters there is a crossover at $`N=31`$ to the Mackay overlayer,Northby (1987) which leads to the 55-atom Mackay icosahedron. For model metal potentials, even those that favour Mackay icosahedra, one often sees decahedral and close-packed structures at sizes away from the complete icosahedra.Doye and Wales (1998) Instead, for these aluminium clusters the polytetrahedral character of the global minima continues. However, it is not by continuation of the anti-Mackay overlayer as for long-ranged Morse clusters,Doye and Wales (1997) but through structures involving ordered arrays of disclinations. Initially some of these structures are similar to that for a modified form of the Dzugutov potential that was designed to produce compact polytetrahedral clusters,Doye and Wales (2001) but they then quickly diverge.

The 40-atom global minimum is a sixfold-symmetric equivalent of Al34. It involves a ring of six interpenetrating icosahedra with two interpenetrating Z14 polyhedra along the symmetry axis. These latter two polyhedra give rise to the single disclination line that passes through the structure. Similar structures have previously been found for Dzugutov clusters,Doye et al. (2001); Doye and Wales (2001) and it is a small fragment of the Frank-Kasper Z phase.

Similar clusters are seen for $`N=44`$ and 48, but with a Z15 and Z16 polyhedron, respectively, lying slightly off-centre; these act as three-fold and four-fold nodes in the disclination networks. Only for Al48 does further growth occur on this cluster leading to the particularly stable 61-atom global minimum with tetrahedral symmetry. This cluster was first identified for a long-ranged Morse cluster,Doye and Wales (1997) and has since been found for clusters with a generalized set of metallic interactions,Cune and Apostol (2000) and in experiments on strontium clusters.Wang et al. (2001)

Two exceptions to this dominance of polytetrahedral structures occur at $`N=54`$ and 55, where a complete Mackay icosahedron and one with a vertex missing atom are most stable. However, away from this magic number structures based on Mackay icosahedra quickly become significantly higher in energy than the polytetrahedral structures.

The 64-atom global minimum is the first to involve not just a single network of connected disclinations. It is based on Al61 but it also has a single linear disclination passing through the clusters. Such disclination lines are a common motif for the larger clusters. At $`N=69`$ the global minimum has an unusual structure with a 19-coordinate central atom that acts as a node for six disclinations going out from it.

Starting at $`N`$=72 a new class of polytetrahedral clusters becomes most stable. They are all based in some way on a set of Frank-Kasper phases, which are layer structuresShoemaker and Shoemaker (1988) and can be considered as a tiling of squares and triangles that are suitably decorated with atoms. The four simplest crystalline tilings are represented in Figure 5(a). The four Frank-Kasper phases correspond to the four possible ways in which a point can be surrounded by squares and triangles. The Z phase has six triangles around each vertex, the A15 phase four squares and the H and $`\sigma `$ phases correspond to the two ways that two squares and three triangles can share a common vertex.

Overlaid on the tilings in Figure 5(a) are the associated disclination networks. All have planar arrays of disclinations that are mediated by Z14 and Z15 polyhedra and that are threaded by linear disclinations that are perpendicular to the plane of the tiling. A number of the global minima can be understood simply in terms of finite square-triangle tilings, as illustrated in Figure 5(b). Al78 has a disclination network that is part of the Z phase; Al72 and Al85 correspond to the H phase, and $`N=138`$, 142, 160 and 172 to a $`\sigma `$-like tiling. The latter slightly deviates from the $`\sigma `$ phase, because, although the local arrangement of squares and triangles around all the vertices in the tiling is the same as for the $`\sigma `$ phase, the central triangle only shares edges with squares. In fact this motif is common in models of dedecagonal quasicrystals.Stampfli (1986); Dzugutov (1993); Roth and Gähler (1998)

The other larger global minima can be basically thought of in terms of these tilings, but with some deviations in the disclination networks at the cluster surface. So $`N=107`$, 120, 132, 134 and 148 have H-like disclination networks, where the Al85 tiling has been extended downward for $`N=134`$ and across for $`N=120`$, 132 and 148. Similarly, $`N=94`$, 116, 151 and 186 also have $`\sigma `$-like tilings. All the $`\sigma `$-like structures are heading towards or based upon Al160.

The deviations from the square-triangle structures involve breaking the planarity of the disclination networks to connect disclinations in different layers of the structure. The most common is to link two parallel disclinations that would otherwise have simply exited the cluster. This connection gives a disclination pattern similar to that for a three triangle tiling, but in a plane perpendicular to the usual plane of the disclination network. For example, the only difference between the structures at $`N=94`$ and 116, and 138 and 151 is the addition of one of these linkages. Similarly, by the addition of two of these linkages one can generate Al186 from Al160. This motif is also present in the structures at $`N=120`$, 132, 134 and 148. The other linkage is similar and involves linking what would otherwise have been three-fold nodes in the disclination network, thus generating Z16 nodes. This type of connection occurs for $`N`$=94 and 116.

It is noticeable that the variations in the energy with respect to the average for clusters with more than 100 atoms are relatively small (Fig. 3). This is because there are lots of ways of creating polytetrahedral packings. This behaviours is unlike what occurs, say, for structures based on the Mackay icosahedra, where, between the magic numbers, structures with incomplete outer layers give rise to a substantial energy variation. Even the most stable large polytetrahedral packings, the $`\sigma `$-like structures centred around $`N=166`$, only give rise to a broad minimum in Fig. 3. The structure at the minimum, Al166, is in fact intermediate between the 160- and 172-atom structures depicted in Fig. 4 with three of the six disclinations extended.

## IV Al61: a case study

To try to unravel why the aluminium glue potential favours polytetrahedral structures, in this section I take Al61 as a case study and compare the energetics of some competing structural forms (Table 2). The structures I consider are the polytetrahedral global minimum (Fig. 4), an icosahedral structure that can be formed by adding a six-atom overlayer to the 55-atom Mackay icosahedron and an fcc cluster. The latter two are illustrated in the insets to Fig. 6.

The decomposition of the potential energy into glue and pair contributions shows (Table 2) that the global minimum is stabilized by a significantly lower $`E_{\mathrm{pair}}`$. This is partially offset by a higher $`E_{\mathrm{glue}}`$ than the icosahedral structure. Figures 6(a) and (b) enable us to understand why the pair energy is lowest for the polytetrahedral clusters, where (a) depicts $`n_{\mathrm{pairs}}(r)`$, the number of pairs of atoms that are separated by less than $`r`$, and (b) $`E_{\mathrm{pair}}(r_{ij}<r)`$, the contribution to $`E_{\mathrm{pair}}`$ from pairs of atoms that are separated by less than $`r`$.

The contribution to $`E_{\mathrm{pair}}`$ from the nearest neighbour shell (defined as those pairs with a separation less than 3.5Å) can be decomposed into two terms

$$E_{\mathrm{pair}}^{\mathrm{nn}}=n_{\mathrm{nn}}\varphi (r_{\mathrm{nn}})+\mathrm{\Delta }E_{\mathrm{pair}}^{\mathrm{dist}},$$
(4)
where $`n_{\mathrm{nn}}`$ is the number of nearest neighbours, $`r_{\mathrm{nn}}`$ is the average separation of pairs in the nearest-neighbour shell. The first term is the $`E_{\mathrm{pair}}^{\mathrm{nn}}`$ value that would result if all the $`n_{\mathrm{nn}}`$ nearest neighbours had the same pair separation, and the second term, $`\mathrm{\Delta }E_{\mathrm{pair}}^{\mathrm{dist}}`$, is the energy cost associated with the nearest-neighbour pair distances having a distribution of values.str

As is well-known icosahedral clusters have larger values of $`n_{\mathrm{nn}}`$ than fcc clusters, because Mackay icosahedra are quasispherical and only have high-coordinate $`\{111\}`$ faces, whereas the most spherical fcc clusters must have a significant proportion of lower-coordinated $`\{100\}`$ faces. This more favourable surface energy is the reason why icosahedral clusters are often lowest in energy for small clusters interacting with pair potentials. The polytetrahedral global minimum has an even higher value of $`n_{\mathrm{nn}}`$ than the icosahedral structure because of the higher coordination number for the interior atoms that lie on disclinations, and the densely packed surface.

As anticipated in Section II.1, $`r_{\mathrm{nn}}`$ is significantly shorter than the distance corresponding to the first minimum in the pair potential—some nearest-neighbour pairs in the polytetrahedral and icosahedral structures even have positive pair energies. $`r_{\mathrm{nn}}`$ also has a significant structural dependence. For the icosahedral and fcc structures, $`r_{\mathrm{nn}}`$ is contracted with respect to $`r_{\mathrm{nn}}^{\mathrm{xtal}}`$, as is usual for metal clusters,Soler et al. (2000) in order to try to increase the $`\overline{\rho }`$ values of the surface atoms. By contrast, $`r_{\mathrm{nn}}`$ for the global minimum is actually slightly larger than $`r_{\mathrm{nn}}^{\mathrm{xtal}}`$. This is not that there is no contraction at the surface, but rather that this contraction is with respect to the appropriate bulk Frank-Kasper phase.

Given the large value of $`n_{\mathrm{nn}}`$ for the global minimum and that $`r_{\mathrm{nn}}`$ lies closest to the potential minimum, it is not so surprising that the global minimum has the lowest value of $`E_{\mathrm{pair}}^{\mathrm{nn}}`$. However, it is only marginally lower than that for the fcc cluster, which in turn is lower than the icosahedral cluster (Figure 6(b)) even though it has fewer nearest neighbours. These effects can be understood by taking account of the distribution of $`r_{ij}`$ values. As $`\varphi (r)`$ has positive curvature in this region, $`\mathrm{\Delta }E_{\mathrm{pair}}^{\mathrm{dist}}`$ must be positive.

As is evident from Figure 6(a) and the $`\sigma _r`$ values in Table 2, the global minimum has the widest distribution of nearest-neighbour distances, and hence by far the largest $`\mathrm{\Delta }E_{\mathrm{pair}}^{\mathrm{dist}}`$. By contrast the fcc cluster has the narrowest distribution and smallest $`\mathrm{\Delta }E_{\mathrm{pair}}^{\mathrm{dist}}`$. This explains why $`E_{\mathrm{pair}}^{\mathrm{nn}}`$ for the fcc cluster is actually lower than that for the icosahedral cluster and only just above that for the global minimum. The large value of $`\sigma _r`$ for the global minimum reflects the distortions of the tetrahedra away from regularity that are necessary in any polytetrahedral packing.

$`E_{\mathrm{pair}}^{\mathrm{nn}}`$, however, only represents roughly a quarter of the total pair energy; the dominant contribution energy comes from next neighbours. There are a number of possible configurations for such neighbours. The nearest next neighbours correspond to those at opposite vertices of an octahedron, which if regular corresponds to a separation of $`\sqrt{2}r_{\mathrm{nn}}=1.4142r_{\mathrm{nn}}`$. These configurations occur most frequently in the fcc clusters, also occur in the twenty fcc tetrahedra that make up the Mackay icosahedron, but do not occur at all in a polytetrahedral structure.

The next set of pair distances correspond to atoms at opposite ends of a trigonal bipyramid, which if it is made of two regular tetrahedra occurs at $`2\sqrt{2/3}r_{\mathrm{nn}}=1.6330r_{\mathrm{nn}}`$. These configurations are absent from fcc clusters, occur around the fivefold axes of the Mackay icosahedra and also across the twin planes between the strained fcc tetrahedra that make up the icosahedron and are most prevalent in polytetrahedral clusters.

Then, there are pairs at the extreme ends of a planar rhombus, which if it consists of two equilateral triangles would be at a separation of $`\sqrt{3}r_{\mathrm{nn}}=1.7321r_{\mathrm{nn}}`$. However, these pairs, which are most common in the fcc structures due to the close-packed planes, are less important to the energetics because, of the three sets of configurations, they are the furthest from the second minimum in the pair potential.

Given that $`r_{\mathrm{min2}}=1.509r_{\mathrm{nn}}^{\mathrm{xtal}}`$, one might expect from the above ideal configurations that the octahedral pairs would lie slightly closer to the second potential minimum than those pairs at opposite ends of a trigonal bipyramid. However, the contraction of $`r_{\mathrm{nn}}`$ with respect to the crystal takes the octahedral distances further away from this minimum. Furthermore, the wide distribution of nearest neighbour distances in the polytetrahedral structures also leads to a wide variation in the dimensions of the trigonal bipyramids. As can be seen from Figure 6(a) a significant fraction of these polytetrahedral pairs for the global minimum lie closer to the pair potential minimum than the octahedral distances. These effects are responsible for the substantially lower pair energy of the global minimum (Figure 6(b)). The trigonal bipyramids present in the icosahedral structure also cause the difference in the total pair energy between it and the fcc structure to narrow, compared to that for $`E_{\mathrm{pair}}^{\mathrm{nn}}`$.

Similar to the analysis for $`E_{\mathrm{pair}}^{\mathrm{nn}}`$ the glue energy can be written as

$$E_{\mathrm{glue}}=NU(\overline{\rho })+\mathrm{\Delta }E_{\mathrm{glue}}^{\mathrm{dist}}$$
(5)
where the first term is the glue energy that would result if all atoms had the same value of $`\overline{\rho }`$ ($`\overline{\rho }`$ is the average value of $`\overline{\rho }`$), and the second term is the correction due to the $`\{\overline{\rho }_i\}`$ having a distribution of values.

We can examine the role of different pair distances on the glue energy by following $`\overline{\rho }(r_{ij}<r)`$, the contribution to $`\overline{\rho }`$ from pairs of atoms that are separated by less than $`r`$. It is clear from Figure 6(c) that, as anticipated, the major contribution to $`\overline{\rho }`$ comes from nearest neighbours, and the values of $`\overline{\rho }(r_{ij}<r)`$ after the nearest-neighbour shell reflect the relative values of $`n_{\mathrm{nn}}`$. $`r_{\mathrm{nn}}`$ is also important. The average contribution to $`\overline{\rho }`$ from each nearest neighbour in the global minimum is less because of the longer $`r_{\mathrm{nn}}`$. Indeed, if the three structures had the same $`r_{\mathrm{nn}}`$ the polytetrahedral structure would have comfortably the lowest glue energy. The differences in $`\sigma _r`$ have little effect, because the function $`\rho (r)`$ is approximately linear for the relevant range of $`r`$ (Figure 2(a)).

As $`\rho (r)`$ is a monotonically decreasing function, the next largest contribution to $`\overline{\rho }`$ comes from octahedral next neighbours. All other pairs have little effect on the relative values of $`\overline{\rho }`$ for the different structures because by these distances $`\rho (r)`$ is small in magnitude.

The effect of next-nearest neighbours is to cause $`\overline{\rho }`$ for the icosahedral structure to just overtake that for the global minimum, and similarly to diminish the difference between $`\overline{\rho }`$ for the global minimum and the fcc structure, compared to the relative values of $`\overline{\rho }(r_{ij}<r)`$ after the nearest-neighbour shell.

The ordering of the $`\overline{\rho }`$ values agrees with that for $`E_{\mathrm{glue}}`$. However, although the global minimum has a similar value of $`\overline{\rho }`$ to that for the icosahedral structure, its glue energy is only just greater than that for the fcc structure. This effect can be explained by taking into account the distribution of $`\overline{\rho }_i`$ values. As $`U(\overline{\rho })`$ has positive curvature, $`\mathrm{\Delta }E_{\mathrm{glue}}^{\mathrm{dist}}`$ must be positive, and is in fact largest for the global minimum because it has the greatest variance in $`\overline{\rho }_i`$ values (Table 2). To a first approximation the distribution of $`\overline{\rho }_i`$ values is bimodal with similar values for surface atoms, and similar values for interior atoms. As the polytetrahedral cluster has the most atoms in the interior of the cluster due to the high coordination numbers, it has the largest $`\sigma _{\overline{\rho }}`$.

As we have pointed out above the value of $`r_{\mathrm{nn}}`$ is important for the values of both $`E_{\mathrm{pair}}`$ and $`E_{\mathrm{glue}}`$, so it is important to understand what determines $`r_{\mathrm{nn}}`$. At $`r_{\mathrm{nn}}`$ the decrease in the glue energy by uniformly shrinking the cluster is exactly balanced by the increase in the pair energy. This energy balance occurs at larger $`r_{\mathrm{nn}}`$ for the global minimum for a number of reasons. Firstly, the decrease in the glue energy as the cluster is shrunk (this increases the $`\overline{\rho }`$ values of the surface atoms) is smallest for the global minimum because it has the fewest number of atoms in the surface of the cluster and has some atoms with high coordination number in the cluster interior (see below). Secondly, the non-linear increase in the pair energy as the cluster is shrunk is largest for the global minimum. This mainly stems from the nearest-neighbour contribution and reflects both the greater number of nearest neighbours for the global minimum and its greater $`\sigma _r`$.

The average nearest-neighbour distance also shows a systematic dependence on the coordination number. For energetic (see Section II.1) and geometric reasons the atoms with higher coordination numbers have larger average nearest-neighbour distances. For example, the nearest neighbours of the sixteen-coordinate atom in the centre of the global minimum are at an average separation of 2.933 Å, whereas the corresponding distance for the 12-coordinate atoms is 2.725 Å.

Although the analysis in this section has been done here for a single case study, similar analyses for other cluster sizes have revealed the same principles at work.

## V Larger clusters

Although polytetrahedral clusters are dominant in the size range considered in the last section, one would expect that clusters at sufficiently large size would adopt the bulk structure, which for aluminium is fcc. This assumes that the current glue potential can correctly reproduce the bulk structure. However, to the best of my knowledge, the relative energetics for different crystal structures have not been compared for this potential; for a somewhat similar aluminium empirical potential the fcc crystal is only marginally lower in energy than an A15 Frank-Kasper phase.Mishin et al. (1999)

That for Al61, and all the other clusters we examined in detail, $`E_{\mathrm{bulk}}^{\mathrm{fcc}}`$ is lower than that for the polytetrahedral and icosahedral clusters (Table 2) suggests that the system will converge to the correct bulk structure. To further understand the size evolution of the cluster structure I have reoptimized a series of Mackay icosahedra, fcc truncated octahedra and Marks decahedra, which are usually the optimal shapes for these morphologies.Raoult et al. (1989) The energies of these sequences are compared in Figure 7(a) along with the global minima obtained in the Section III. I have not attempted to construct a sequence of larger polytetrahedral clusters, because one cannot easily deduce from the global minima how the polytetrahedral growth sequence should continue. Instead, we also plot $`E_{\mathrm{ave}}`$, the fit to the energies of the global minima.

From Figure 7(a) it is clear that decahedral structures are never competitive for this potential. This is slightly unusual as it is often found that for materials which are close-packed in bulk, there is a progression from icosahedral to decahedral to close-packed structures as the size increases. Instead, the icosahedral and fcc lines cross at $`N520`$ suggesting that beyound this size fcc clusters will be lower in energy for the majority of sizes. The $`E_{\mathrm{ave}}`$ curve also seems to suggest that these fcc clusters will be global minima, but one has to be somewhat careful. Firstly, one should not necessarily trust this function outside the size range to which it has been fitted. Secondly, it represents the average energy of the polytetrahedral global minima, whereas the other curves in Figure 7(a) are lower bounds to the energies of the clusters of that type. For example, $`E_{\mathrm{ave}}`$ only dips below the interpolated energies of the complete Mackay icosahedra for a small range of $`N`$, even though polytetrahedral structures are the global minima for the vast majority of the size range considered in Section III.

To get a better idea of the causes of this size evolution I have plotted in Figures 7(b) and (c) the glue and pair contribution to the energies of these sequences of structures. It is immediately clear, in agreement with the analysis of Al61, that the polytetrahedral global minima are stabilized by their particularly low pair energies. Furthermore, the global minima that are based on Mackay icosahedra at $`N=54`$–55 and 147 clearly stand out from this trend.

It is also clear that the crossover from icosahedral to fcc structures is driven by the increasing advantage provided by its lower pair energy. By contrast, the glue energy of the icosahedron only becomes greater than that for the truncated octahedra at $`N6500`$. $`E_{\mathrm{glue}}`$ more strongly reflects the greater number of nearest neighbours associated with the Mackay icosahedra and is less sensitive to the strain inherent to the icosahedra. By contrast, as noted in section III and as is usual for pair potentials,Doye and Wales (1997) the nearest neighbour contribution to the pair energy is reduced due to the strain in the pair distances for the icosahedra.

## VI Experimental and theoretical comparisons

Although I have focussed on the general principles that can be drawn from the current results, given that the potential has been successfully used in a wide variety of applications, with in some cases excellent agreement with experimental data,Sandberg et al. (2002) one would hope that this system would also be a reasonable model for real aluminium clusters. However, one should remember that to correctly predict the energetics of cluster structure requires that the potential can successfully describe a wide variety of surface and bulk properties for a variety of structural types, and so is an extremely stringent test of a potential.

Although aluminium clusters have been much studied experimentally, there is little reliable information about the geometric structure in the size range of the global minima considered in Section III. The difficulty is that features in the experimental data can originate from the electronic shell structure rather than the geometric structure. Indeed, following the successful interpretation of many of the properties of alkali metal clusters in terms of the electronic shell structure, aluminium, a nearly free-electron metal, was a natural system to see how widely applicable these ideas might be. Indications of this electronic shell structure were found from the ionization potentials,Schriver et al. (1990) the magic numbers in mass spectra,Schriver et al. (1990); Pellarin et al. (1993) cluster ion mobilitiesJarrold and Bower (1993) and photoelectron spectroscopyLi et al. (1998); Akola et al. (1999, 2000).

However, a significant number of the features in these properties could not be explained thus, and so these were usually attributed to the effects of the geometric structure. Noteworthy of these is the magic number at $`N=55`$, whose most obvious interpretation is in terms of a Mackay icosahedron. No convincing explanations, however, have been given for the other ‘geometric’ features. For example, except for $`N=55`$, the magic numbers do not correspond to those of any of the usual structural forms for clusters, but neither do they correspond to the new polytetrahedral magic numbers found here.

Theoretically, aluminium clusters in this size range have been recently studied using two empirical potentials, namely the Murrell-MottramLloyd and Johnston (1998); Lloyd et al. (2002) and GuptaBailey et al. (2003) forms. For both these system icosahedral structures predominate for the global minima, although with fcc and decahedral structures competing at sizes intermediate between the complete Mackay icosahedra. I have also reoptimized a series of icosahedral, decahedral, fcc and polytetrahedral clusters for the potential of Mishin et al., whose $`\varphi (r)`$, $`\rho (r)`$ and $`U(\overline{\rho })`$ functions resemble those for the current glue potential and were fit to a wide range of bulk data.Mishin et al. (1999) In this case polytetrahedral structures again predominate.

There have also been many electronic structure calculations of small aluminium clusters using a variety of methods.Jones (1991); Yi et al. (1991); Cheng et al. (1991); Debiaggi and Caro (1992); Yang et al. (1993); Akola et al. (1998, 1999, 2000); Ahlrichs and Elliott (1999) For small clusters there has been some success in identifying the structure for $`N<20`$ by comparing calculated photoelectron spectra with experiment.Akola et al. (1999, 2000) However, for somewhat larger clusters there is little consensus among the results, partially because the computational expense prevents one rigorously searching for the global minimum—usually only a selection of clusters are compared—but also because different methodologies simply give different results. For example, there is disagreement over whether for Al55 the Mackay icosahedronYang et al. (1993) or the fcc cuboctahedronCheng et al. (1991); Ahlrichs and Elliott (1999); Akola et al. (2000) is lower in energy, and still others have found disordered structures to be lower in energy than both.Yi et al. (1991) A similar range of results have been found for Al147.Debiaggi and Caro (1992); Ahlrichs and Elliott (1999)

For large aluminium clusters there are clear experimental indicators of the geometric structure. For $`N>250`$ (up to at least $`N\mathrm{10\hspace{0.17em}000}`$) there are set of regular spaced magic numbers in the mass spectra.Lermé et al. (1992); Martin et al. (1992); Pellarin et al. (1993); Baguenard et al. (1994) Initially, this result was interpreted in terms of electronic-shell effects,Lermé et al. (1992); Pellarin et al. (1993) but Martin et al. soon showed that the magic numbers were in fact due to octahedra, with each additional magic number corresponding to the addition of a complete overlayer to a single face of the octahedra.Martin et al. (1992); Näher et al. (1993) Growth simulations further added to the plausibility of this interpretation.Valkealahti and Manninen (1994); Valkealahti et al. (1995); Valkealahti and Manninen (1998)

For large clusters one should always ask whether the observed structures occur because they are most stable or because of the growth kinetics. This is especially the case here, because the most stable fcc clusters are usually more spherical truncated octahedra,van de Waal (1989) rather than octahedra. Indeed as the life times of the clusters are increased, the magic numbers disappear,Baguenard et al. (1994) indicating their kinetic origin. However, it is probably reasonable to conclude that fcc structures are most stable in this size range, even if the appearance of octahedra is a kinetic effect that results from the nucleation of a new face being slower than the subsequent growth.

My estimate that fcc structures become lowest in energy beyond $`N500`$ is therefore reasonably close to experiment. This contrasts with the Gupta potential where fcc structures only become more stable beyond $`N1250`$ and the Murrell-Mottram potential where the crossover occurs at even larger sizes.Turner et al. (2000) However, one should remember that these energetic crossover sizes are only strictly valid at $`T=0`$K. If one compares the relative free energies, the crossover sizes can be sensitively dependent on the temperature, mainly due to the differences in the vibrational entropy of the competing structures.Doye and Calvo (2001); Doye and Wales (2002b)

## VII Conclusion

I have found that the global minimum of aluminium clusters described by a glue potential to be predominantly polytetrahedral for $`N<160`$. The structure of the larger clusters can be described in terms of their disclination networks, which resemble those of the square-triangle Frank-Kasper phases. This result is somewhat surprising given that ordered polytetrahedral structures are only found for alloys in bulk, but it reflects the greater tendency for polytetrahedral order in clusters. The bulk behaviour should of course emerge for sufficiently large clusters, and our results suggest that for $`N>500`$ clusters with the bulk fcc structure predominate.

An important aspect of this work has been to relate the observed structures back to the interatomic interactions. Polytetrahedral structures are favoured for a combination of reasons. Firstly, the clusters present surfaces with a low surface energy and high average cooordination number, so will be potentially favoured at small sizes. Secondly, the glue function, which gives rise to most of the cohesive energy of the clusters, is relatively insensitive to the strain in the pair distances and the high coordination numbers that are inevitable in polytetrahedral structures. Thirdly, as with the modified Dzugutov potential,Doye and Wales (2001) the two minima in the pair potential most closely match the distances in polytetrahedral structures.

As the first two of these reasons are true for most metal potentials, it is not so surprising that there is experimental evidence of polytetrahedral structures for monatomic cobalt,Dassenoy et al. (2000) strontium,Wang et al. (2001) and perhaps some rare-earthBréchignac et al. (1993) clusters. Nor would it be surprising if the number of theoretical and experimental examples were to increase. My results furtherDoye and Wales (2001) illustrate the types of polytetrahedral structures that one might expect for such clusters, along with a series of magic numbers that will be useful in comparisons with experiment.

Whether the preference for polytetrahedral structures in the current model carries over to real aluminium clusters is not yet clear. Because of the lack of easily interpretable experimental data on the geometric structure of aluminium clusters with less than 190 atoms, it is hard to draw any firm conclusions about the realism of the structures predicted by the glue potential in this size range. Moreover, ab intio calculations add little further insight because of the lack of consensus. However, comparisons to the experimental data for larger clusters suggest that the glue potential is more realistic than the other empirical potentials for which calculations have been performed.

The properties of these clusters can also provide useful insights into the behaviour observed for this potential in other contexts. For example, studies on aluminium nanowires have found polytetrahedral structures for those wires with sufficiently small cross-section.Gülseren et al. (1998) Some of these involve only icosahedral local coordination, but others clearly involve disclinations. For example, the wire that Gülseren et al. denote as A7 has a single disclination running along the centre of the wire. Similarly, from the current results it should be unsurprising that there is strong polytetrahedral ordering in the aluminium melt described by this potential.Liu et al. (2001)


format: text

text:
# Defect-mediated turbulence in systems with local deterministic chaos

## Abstract

We show that defect-mediated turbulence can exist in media where the underlying local dynamics is deterministically chaotic. While many of the characteristics of defect-mediated turbulence, such as the exponential decay of correlations and a squared Poissonian distribution for the number of defects, are identical to those seen in oscillatory media, the fluctuations in the number of defects differ significantly. The power spectra suggest the existence of underlying correlations that lead to a different and non-universal scaling structure in chaotic media.

Weakly driven, dissipative, pattern-forming systems often exhibit the spatiotemporal chaos in the form of defect-mediated turbulence, where the dynamics of a pattern is dominated by the perpetual motion, nucleation and annihilation of point defects (or dislocations) Coullet et al. (1989). Examples can be found in striped patterns in wind driven sand, electroconvection in liquid crystals Rehberg et al. (1989), nonlinear optics Ramazza et al. (1992), fluid convection Morris et al. (1993); Daniels and Bodenschatz (2002), autocatalytic chemical reactions Ouyang and Flesselles (1996) and Langmuir circulation in the oceans Haeusser and Leibovich (1997). These results suggest that the dynamics of these very different systems can be characterized by a universal description which is based only on the defect dynamics.

A statistical description of defect mediated turbulence based on a simple model for the defect dynamics was given by Gil, *et al.* Gil et al. (1990). Treating the defect pairs as statistically independent entities, the nucleation rate for pairs of defects was taken to be independent of the number of pairs $`n`$ and, based on the topological nature of the defects, the annihilation rate was taken to be proportional to $`n^2`$. This directly led to a squared Poissonian distribution for the probability distribution function (PDF) of $`n`$. Gil, *et al.* found that simulations of a spatiotemporal chaotic state of the complex Ginzburg-Landau equation (CGLE), which is the prototype of oscillatory media, agreed with their prediction. Rehberg, *et al.* Rehberg et al. (1989) measured the PDF of $`n`$ for defect-mediated turbulence in electroconvection of nematic liquid crystals and found it to be consistent with the predicted squared Poissonian distribution. Later, Ramazza, *et al.* Ramazza et al. (1992) investigated a defect turbulent state in optical patterns and found that their data were not conclusive. Very recently, Daniels and Bodenschatz Daniels and Bodenschatz (2002) studied the defect-mediated turbulent state of undulation chaos in inclined layer convection of a fluid and found that the observed pair nucleation and annihilation rates agree with the theoretical predictions in Gil et al. (1990). They also derived the PDF of $`n`$ for the case that single defects can leave and enter through the boundaries and found agreement with their measurements. When boundary effects are negligible, this PDF reduces to the squared Poissonian distribution.

The experimental and theoretical studies of defect-mediated turbulence have focused exclusively on first order statistics like the PDF of $`n`$ and creation and annihilation rates. Moreover, theoretical investigations have been carried out only for media with underlying oscillatory dynamics, generically described by the CGLE. Yet, in many cases the local dynamics may differ from simple oscillatory behavior; instead, complex-periodic or even chaotic attractors may exist (see, e.g., Ref. Goryachev et al. (2000)). Chemically reacting systems, notably the Belousov-Zhabotinsky (BZ) reaction, are known to exhibit deterministic chaos Scott (1991). This leads to a variety of new spatiotemporal states Park and Lee (1999). In this Letter, we show that defect-mediated turbulence can exist in media where the underlying local dynamics is chaotic and that second-order statistics can be used to distinguish between different media. This implies that a universal description of defect-mediated turbulence cannot encompass second-order or higher correlations.

Our focus is on systems where the dynamics of the spatially homogeneous system, described by ordinary differential equations, has a deterministic chaotic attractor; hence, at least three phase space variables are required in contrast to the two-variable descriptions of simple oscillatory media. A specific example of such a chaotic system is the Willamowski-Rössler (WR) reaction-diffusion model Willamowski and Rössler (1980) given by

$`_t𝐜(𝐫,t)`$ $`=`$ $`𝐑[𝐜(𝐫,t)]+D^2𝐜(𝐫,t),`$ (1)
where $`R_1=\kappa _1c_1\kappa _1c_1^2\kappa _2c_1c_2+\kappa _2c_2^2\kappa _4c_1c_3+\kappa _4`$, $`R_2=\kappa _2c_1c_2\kappa _2c_2^2\kappa _3c_2+\kappa _3`$ and $`R_3=\kappa _4c_1c_3+\kappa _4+\kappa _5c_3\kappa _5c_3^2`$. These rate equations were derived from the mass action kinetics of a reaction scheme with quadratic kinetics where certain pool species are taken to be fixed. Here $`c_i(𝐫,t)`$ is the local concentration of species $`i`$ at site $`𝐫`$ in a two-dimensional space of size $`L^2`$ with periodic boundary conditions. The parameters $`\kappa _{\pm j}`$ are rate coefficients that contain the concentrations of the pool species that are fixed to maintain the system out of equilibrium. The diffusion coefficients $`D`$ of all three species are taken to be equal. Even though the WR model is very simple, it exhibits a phenomenology Goryachev and Kapral (1996a) with many features in common with those observed in chemical experiments on the BZ reaction Park and Lee (1999). Consequently, this model may be expected to capture the qualitative features of chemical systems whose chaotic attractors arise from a period-doubling cascade. The chaotic attractor for a certain set of rate coefficients is shown in the left panel of Fig. 1.

A defect is characterized by its integer topological charge $`m_{top}`$ which is defined by Mernin (1979)

$$\frac{1}{2\pi }\varphi (𝐫,t)𝑑𝐥=\pm m_{top},$$
(2)
where $`\varphi (𝐫,t)`$ is the local phase and the integral is taken along a closed curve surrounding the defect. A topological defect corresponds to a point in the medium where the local amplitude is zero and the phase is not defined. Typically only topological defects with $`m_{top}=\pm 1`$ are observed. One-armed spiral waves with such defects at their centers are the only stable spiral waves for the CGLE Hagan (1982). The phase has to be defined in order to apply the notion of defects to chaotic media. This is not a trivial issue and for many chaotic systems it is impossible to introduce a phase field. The rather simple shape of the WR chaotic attractor (see left panel of Fig. 1) which arises from a period-doubling cascade admits a simple definition of the phase and the WR model belongs to the class of chaotic-oscillatory media Pikovsky et al. (2001). We chose $`\varphi (𝐫,t)=\mathrm{arctan}[(c_2(𝐫,t)c_2^0)/(c_1(𝐫,t)c_1^0)]`$ with $`(c_1^0,c_2^0)=(8.0,9.0)`$ as the center of rotation. For such simple chaotic attractors, the particular choice of a phase variable for chaotic systems does not influence the results Josić and Mar (2001). In the center panel of Fig. 1, the phase field of the WR medium is shown for a certain set of parameters numerics . The topological defects can be identified as the termini of the white equiphase contour lines. As for the CGLE in the defect-mediated turbulent state Chaté and Manneville (1996); Aranson and Kramer (2002), the defects in the WR system for these parameter values rarely emit waves. They behave as passive objects and are merely advected by the surrounding chaotic fluctuations. The right panel in Fig. 1 shows that the correlation function $`C(|𝚫𝐫|)=Re[A(𝐫,t)\overline{A}(𝐫+𝚫𝐫,t)]_{𝐫,t}`$, where $`A(𝐫,t)`$ is the complex amplitude in the $`(c_1,c_2)`$-plane with respect to the center of rotation, decays exponentially with a very short characteristic length scale, verifying the existence of a turbulent state. This is the most typical characteristic for defect-mediated turbulence and has been found in the CGLE as well Coullet et al. (1989).

The fluctuations in the number of pairs of topological defects shown in the left panel of Fig. 2 provide further evidence that a defect-mediated turbulent state can exist in chaotic-oscillatory media. The total number of defects in the medium is exactly twice the number of pairs because the net topological charge is conserved and equal to zero due to the periodic boundary conditions. Hence, topological defects can only be created and annihilated in pairs of opposite topological charge. One can easily derive a PDF $`p(n)`$ for the number of defects provided that the defects are statistically independent entities Gil et al. (1990); Daniels and Bodenschatz (2002). In the stationary state and for periodic boundary conditions, the master equation reduces to $`p(n)=p(n1)c(n1)/a(n)`$ where $`c(n)`$ and $`a(n)`$ are the creation and annihilation rates, respectively. Provided that $`c(n)=c=\mathrm{const}`$ and $`a(n)=an^2`$, $`p(n)(c/a)^n/(n!)^2`$. The right panel of Fig. 2 shows that the PDF for the WR model agrees with the predicted form of a squared Poissonian distribution reasonably well. Moreover, the assumptions leading to this distribution seem to be justified: the annihilation rate scales approximately with $`n^2`$ and the creation rate is approximately independent of $`n`$ as can be deduced from Fig. 3.

The fluctuations in the number of defects $`n`$ do not have the properties one would expect if they arose from independent random events. If this were the case the power spectrum $`S_L(f)=lim_T\mathrm{}1/2T\underset{T}{\overset{T}{}}𝑑tn(t)\mathrm{exp}^{i2\pi ft}^2`$ for a system of size $`L`$ would have a Lorentzian shape. The power spectrum does not have this form for the defect-mediated turbulence in the WR reaction-diffusion system. Figure 4 shows that $`S_L(f)1/f^\gamma `$ for intermediate frequencies with an exponent $`\gamma `$ that is far from the value $`\gamma =2`$ expected for a Lorentzian shape. For $`k_2=0.072`$, we find $`\gamma =1.43`$ and, for $`k_2=0.075`$, $`\gamma =1.60`$. Although in both cases the system exhibits defect-mediated turbulence, the exponents are significantly different from each other. This implies that different chaotic-oscillatory media can have different second-order statistics. To confirm that these results are not specific to the WR system, we have also studied the autocatalator reaction-diffusion system Peng, Scott and Showalter (1990); Kapral and Wu (1997) which also has a chaotic attractor arising from a period-doubling cascade but with very different Lyapunov spectra. We find a power-law decay with values of $`\gamma `$ that are similar to those reported here for the WR model footnote . These results suggest that the power spectrum of $`n(t)`$ may exhibit power-law decay with non-trivial exponents for defect-mediated turbulent states when the local dynamics exhibits deterministic chaos.

For large enough system sizes, $`S_L(f)L^2`$ for all $`f`$; thus, the power spectrum for large systems can be considered to be the superposition of the power spectra of subsystems. Such behavior is expected in view of the short correlation length. This implies that the standard deviation $`\mathrm{\Sigma }=(n^2n^2)^{1/2}`$ scales as $`\sqrt{n}`$, which is proportional to $`L`$. Such a scaling is a consequence of the law of large numbers and has been observed for the CGLE as well. It follows that $`\gamma `$ and the low-frequency cutoff are independent of $`L`$. The low-frequency cutoff is due to the fact that $`n(t)`$ is bounded. Its location depends on the density of defects in the medium. For lower densities the cutoff moves to lower frequencies.

To understand how the non-trivial correlations in $`n(t)`$ at intermediate time scales arise, we have analyzed the series of waiting times between consecutive creation events and consecutive annihilation events separately. In both cases, the waiting times are exponentially distributed and statistically uncorrelated as for a random walk. This implies that the correlations in $`n(t)`$ are due to the interaction of creation and annihilation events.

To compare our results for the power spectrum of $`n(t)`$ in chaotic media with those in oscillatory media, we have simulated the CGLE Coullet et al. (1989); Aranson and Kramer (2002)

$$_tA=A+(1+i\alpha )^2A(1+i\beta )|A|^2A,$$
(3)
in a domain of size $`L^2`$ with periodic boundary conditions numerics\_cgl . Here, $`A`$ is the complex amplitude field and $`\alpha `$ and $`\beta `$ are control parameters. The power spectrum of the CGLE in the defect-mediated turbulent state has a power-law decay for intermediate frequencies and the scaling of $`S_L(f)`$ with $`L`$, as well as the dependence of the low-frequency cutoff on the density of defects, agree with those in the WR medium. Yet, as Fig. 5 shows, the exponent $`\gamma =1.9`$ which is very close to that expected for a Lorentzian form if the dynamics of $`n(t)`$ is a simple random walk in a bounded domain. It is very different from the values observed for the chaotic media we analyzed. Thus, different media can have different second-order statistics for $`n(t)`$ depending on the underlying local dynamics. This argues against a universal description of defect-mediated turbulence.

We have shown that defect-mediated turbulence can arise both in oscillatory and chaotic-oscillatory media and, thus, applies to a broad range of systems. While most common diagnostic measures do not allow one to distinguish between oscillatory and various chaotic media, the fluctuations in the number of defects can be different for different media: for the CGLE they resemble the form of a bounded random walk; for the WR and autocatalator models pronounced non-trivial correlations exist. Our results may be tested experimentally on systems like the Belousov-Zhabotinsky reaction whose local temporal dynamics can be complex-periodic or even chaotic Park and Lee (1999), and where defect-mediated turbulence has been observed Ouyang and Flesselles (1996).

This work was supported in part by a grant from the Natural Sciences and Engineering Council of Canada.


text:
# Hole–defect chaos in the one–dimensional complex Ginzburg–Landau equation

## I Introduction

The formation of local structures and the occurrence of spatiotemporal chaos are among the most striking features of pattern forming systems. The complex Ginzburg–Landau equation (CGLE)

$$_tA=A+(1+ic_1)_x^2A(1ic_3)|A|^2A$$
(1)
provides a particularly rich example of these phenomena. The CGLE is the amplitude equation describing pattern formation near a Hopf bifurcation CH ; AK , and exhibits an extremely wide range of behaviors as a function of $`c_1`$ and $`c_3`$ CH ; AK ; chao1 ; phasediagram ; mvh ; maw ; mvhmh .

Defects and holes are local structures that play a crucial role in the intermediate regime between laminar states (small $`c_1`$, $`c_3`$) and hard chaos (large $`c_1`$, $`c_3`$). Isolated defects occur when $`A`$ goes through zero, where the complex phase $`\psi :=\mathrm{arg}(A)`$ is no longer defined. Homoclinic holes are localized propagating “phase twists”, which are linearly unstable. As illustrated in Fig. 1, holes and defects are intimately connected: defects can give rise to “holes”, which may then evolve to generate defects, from which further holes can be born, sometimes generating self–sustaining patterns. For more details see Refs. mvh ; mvhmh .

The aim of our paper is to understand and model the spatiotemporally chaotic hole–defect behavior of the 1D CGLE, built on the local interactions and dynamics of the holes and defects. Given the strength of the initial phase twist that generates a hole, and the wavenumber of the state into which it propagates, the hole lifetime $`\tau `$ turns out to be the key feature. Surprisingly, the initial phase twist and invaded state play very different roles. For hole–defect chaos, we will show that the defect profiles, which constitute the phase twist initial condition for the resulting daughter holes, show rather little scatter for fixed $`c_1`$ and $`c_3`$. Changes in $`c_1`$ and $`c_3`$, however, are encoded in changes in the defect profiles, and thus lead to changes in the typical lifetimes of the daughter holes. We then demonstrate that the chaos does not result from variations in defect profiles. It rather follows from the sensitivity of the holes to the states they invade, since the passage of each hole disorders the background wavenumber profile leading to disordered background states. It is the self–disordering action of the holes that is primarily responsible for the spatiotemporal chaos.

With these insights, we can then construct a simplified lattice model for hole–defect chaos, which both reproduces the correct qualitative behavior as $`c_1`$ and $`c_3`$ are varied and which captures the correct mechanism (propagating, self–disordering holes). Our initial findings on this subject can be found in Ref. mvhmh , where we introduced the concept of self–disordering, and outlined a simplified lattice model. However in this paper, we investigate the subject in considerably greater depth, and, in particular, provide much more conclusive evidence for the correctnesss of the self–disordering hypothesis.

We now briefly summarize the structure of the paper. Topics discussed already in earlier work mvh ; mvhmh  are dealt with rather briefly. We start in Section II by describing hole–defect dynamics on a local scale. In Section III, we then use this knowledge to investigate disordered hole–defect dynamics, and conclusively show that it is the coupling of the holes to a self–disordered background which is the dominant mechanism for spatiotemporal chaos. This concept is then illustrated by a minimal lattice model for hole–defect dynamics in Section IV, before we draw our conclusions in Section V.

## II Hole–defect dynamics

We begin by studying the hole lifetime $`\tau `$ as a function of the initial conditions (Fig. 2). This study motivates the central question of this paper: how does $`\tau `$ depend on the initial conditions and on the external wavenumber, and which of these dependencies is most important for spatiotemporal chaos? We then study general properties of defect profiles, and demonstrate that in hole–defect chaos the profiles of defects show rather little scatter. We also show how the lifetimes of “daughter” holes born from a typical defect vary with $`c_1`$ and $`c_3`$. Taken together, the data presented here forms direct evidence for the heuristic picture of hole–defect dynamics developed in Refs. mvh ; mvhmh .

### II.1 Incoherent homoclons

In full dynamic states of the CGLE, one does not observe the unstable coherent homoclinic holes unless one fine–tunes the initial conditions (see below). Instead evolving incoherent holes which either decay or start to grow towards defects occur mvh ; mvhmh . Let us consider the short–time evolution of an isolated incoherent hole propagating into a regular plane wave state. Holes can be seeded from initial conditions like mvhmh :

$$A=\mathrm{exp}(i[q_{ex}x+(\pi /2)\mathrm{tanh}(\gamma x)]).$$
(2)
The two essential parameters $`\gamma `$ and $`q_{ex}`$ represent, respectively, the initial conditions from which the incoherent hole is born and the background wavenumber of the state into which the hole propagates. In this context, a single parameter $`\gamma `$ is sufficient to scan through different initial conditions, since the coherent holes have just one unstable mode mvh .

A detailed contour plot of the lifetime of an initial incoherent hole as a function of $`\gamma `$ and $`q_{ex}`$ is shown in Fig. 2. These results were obtained using a semi–implicit numerical integration of the CGLE, with space and time increments $`dx=0.25`$ and $`dt=0.01`$. As expected, three possibilities can arise for the time evolution of the initial peak: evolution towards a defect (upper right part of Fig. 2), decay (lower left part of Fig. 2), or evolution arbitrary close to a coherent homoclinic hole (the boundary between these two regions). These possibilities, together with an illustrative sketch of the phase space, are shown in the inset of Fig. 2. The rather simple and monotonic behavior of $`\tau `$ with $`q_{ex}`$ and $`\gamma `$ is somewhat of a surprise, and this reinforces our simple phase space picture; no other solutions seem to be relevant in this region of phase space.

Since homoclons are neither sinks nor sources, Fig. 2 can be interpreted as follows: for a right–moving homoclon an incoming wave with positive wavenumber tends to push the homoclon more quickly towards a defect; previously we have referred to this as “winding up” of the homoclon. Similarly, an incoming wave with negative wavenumber “winds down” a right–moving homoclon, possibly even preventing the formation of a defect mvh .

### II.2 Defects

We now study the defect profiles themselves in more detail. In Fig. 3a we show complex plane plots of $`Re(a)`$ vs. $`Im(A)`$ just before, close to, and just after a defect. As can be seen, there is no singular behavior whatsoever: the real and imaginary parts are smooth functions of $`x`$ and $`t`$, even at the time of defect formation. However, when transforming to polar coordinates, a singularity manifests itself at the defect, where $`|A|0`$. This can also be seen from the $`q`$-profiles shown in Fig. 3b-d. In fact, it is straightforward to show that the maximum value of the local phase gradient $`q_m`$ diverges as $`(\mathrm{\Delta }t)^1`$ at a defect mvhmh , where $`\mathrm{\Delta }t`$ is the time before defect nucleation (see Fig. 3b of Ref. mvhmh ) error .

In Fig. 3(e-g), we overlay complex plane plots of $`A`$ around $`10^3`$ defects obtained from numerical simulations of the CGLE in the chaotic regime. Surprisingly (see Fig. 3e), defect profiles of the interior chaotic states are dominated by a single profile in the hole–defect regime, similar to fixing $`\gamma `$ in Fig. 2. This provides a strong indication that hole–defect chaos does not come from scatter in the defect profiles. For large enough $`c_1`$ and $`c_3`$, where holes no longer play a role, and where hard defect chaos sets in maw , the profiles show a much larger scatter (Fig. 3(f,g)).

### II.3 Defect $``$ holes

Suppose a hole has evolved to a defect; what dynamics occurs after this defect has formed? As Fig. 3(b-d) shows, defects generate a negative and positive phase–gradient peak in close proximity. The negative (positive) phase gradient peak generates a left (right) moving hole, and analogous to what we described in Fig. 2, the lifetimes of these holes depend on the initial peak and on $`q_{ex}`$. Hence the defect profile acts as an initial condition for its daughter holes, as can also clearly be seen in Fig. 1.

We now examine the fate of these daughter holes in the well–defined case where the initial defect is generated from the divergence of a right–moving, near–coherent homoclon in a $`q_{ex}=0`$ background state (see Fig. 4a). We then define $`t_1`$ and $`t_2`$ as the lifetimes of the resulting daughter holes. When a daughter hole does not grow out to form a defect, its lifetime diverges. In Fig. 4b we plot $`t_1`$ and $`t_2`$ for $`c_1=0.6`$ as a function of $`c_3`$. The initial hole that formed the first defect has a lifetime of at least $`60`$, and we have checked that a further increase of this time does not change $`t_1`$ and $`t_2`$ appreciably. When both $`t_1`$ and $`t_2`$ are infinite, no defect sustaining states can be formed, and the final state of the CGLE is in general a simple plane wave. When only $`t_1`$ is finite, isolated zigzag states are formed; such states have been discussed in Ref. mads , and we will see some examples below. When both $`t_1`$ and $`t_2`$ are finite, and of comparable value, more disordered states occur. We will later use this data on $`t_1`$ and $`t_2`$ to calibrate our minimal lattice model for spatiotemporal chaos.

Hence, we see that changes in $`c_1`$ and $`c_3`$ not only lead to changes in the defect profiles, but also modify the lifetimes of the resulting daughter holes. However, for fixed $`c_1`$ and $`c_3`$, we have seen that the defect profiles, which act as initial conditions for the daughter holes, show rather little scatter, at least in the hole–defect regime. In the next section, we will build on this knowledge to unravel the causes of hole–defect spatiotemporal chaos.

## III Mechanism of hole–defect chaos

In this paper and in earlier work mvhmh , we have argued that the principal cause of the spatiotemporally chaotic behavior in the 1D CGLE is the movement of holes through a self–disordered background. Clearly, as we can see from Fig. 2, a disordered background wavenumber $`q_{ex}`$ will give rise to varying hole lifetimes and thus to disordered hole–defect dynamics. In this section, we explicitly demonstrate the correctness of this mechanism by modifying the CGLE dynamics in two ways:

Model (i): Fixed defect profile. Whenever a defect occurs, this defect is replaced by a standardized defect profile (obtained from an edge defect). Here the dynamics will be chaotic, showing the irrelevance of the scatter in defect profiles.

Model (ii): Background between holes $``$ plane wave with $`q=0`$. At each timestep, the background between any two holes is replaced by a plane wave with wavenumber zero. Here no chaos will occur, illustrating the crucial importance of the self–disordered background.

In case (i), the size of the replaced defect profile was five centered around the defect; in case (ii) the background was defined to be all regions where $`|A|>0.95`$. Our results are substantially independent of the exact defect size or cutoff value. In both cases, it is crucial to ensure that no jumps in the phase occur at the edges of the replaced regions. This can be achieved by phase matching the replaced region (either defect profile or plane wave) at the left boundary, while the state to the right of the replaced region is multiplied by a phase factor to enforce phase–continuity at the right edge. We take open boundary conditions (i.e. $`_xA=0`$) and only study the behavior far away from these boundaries.

In Fig. 5 we show an example of the dynamics and spreading of a localized perturbation for the full CGLE (Fig. 5a–b) and for the “fixed defect” model (i) (Fig. 5c–d), both for $`c_1=0.6,c_3=1.4`$. For both models, we took as an initial condition a defect rich state, which after a few timesteps shows the typical hole–defect dynamics. At $`t167`$ we applied a local perturbation of strength $`10^6`$ to the middle gridpoint (corresponding to $`x=500`$) and followed the evolution of both the perturbed and the unperturbed systems in order to follow the spreading of perturbations. For the full CGLE (Fig. 5a–b), the perturbation spreads along with the propagation of the holes. We note that the initial growth of the perturbation manifests itself in slight “shifts” of the spatial and temporal positions of the defects. In particular when two holes collide, a strong amplification of the perturbations is observed.

We can now compare this with the above fixed defect profile model (i) (Fig. 5c,d). Clearly, the replacement of the defects does not destroy the chaotic behavior of the system, as confirmed by the spreading of a localized perturbation (Fig. 5d), which propagates in a similar fashion to the full CGLE (Fig. 5b). This strongly indicates that variation in the defect profiles is not contributing in a major way to the spatiotemporally chaotic behavior of the full CGLE. We should also point out one subtlety here: due to the discretization of time, the times at which defect profiles are replaced are also discretized, and one may worry whether this destroys the chaotic properties of the model. However, we have performed simulations for a smaller timestep ($`dt=0.001`$) and found no qualitative difference. As we will see, this issue of discretization will play a more important role in the lattice model discussed in section IV.

Turning now to model (ii), where laminar regions of the CGLE are replaced by $`q=0`$ plane waves, we see that the disorder is destroyed. This is illustrated in Fig. 6, where we show examples of model (ii) dynamics for $`c_3=1.4`$ and $`c_1=0.6,0.7,0.8`$. Clearly chaos is suppressed for $`c_1=0.6`$ and 0.7 (Fig. 6a–b), with zigzag type patterns being especially dominant. For $`c_1=0.8`$ (Fig. 6c), the initial dynamics do appear disordered, but after a transient the system evolves to a regular zigzag state (Fig. 6d).

From the behavior of models (i) and (ii) we conclude that the self–disordered background is an essential ingredient for hole–defect chaos, while scatter in the defect profiles is not.

## IV Lattice model

To further justify and test our picture of self–disordered dynamics, we will now combine the various hole–defect properties with the left–right symmetry and local phase conservation to form a minimal model of hole–defect dynamics. The model reproduces regular edge states, spatiotemporal chaos and can be calibrated to give the correct behavior as a function of $`c_1`$ and $`c_3`$. An earlier version of the model was presented in Ref. mvhmh . However, as will become clear, we have now modified and improved the model, and also for the first time made direct comparisons with the full CGLE.

From our earlier analysis (see also Ref. mvhmh ), we see that the following hole–defect properties must be incorporated in the model: (i) Incoherent holes propagate either left or right with essentially constant velocity. (ii) Their lifetime depends on $`c_1`$, $`c_3`$, and the wavenumber of the state into which they propagate. When the local phase gradient extremum diverges, a defect occurs. (iii) Each defect, in turn, acts as an initial condition for a pair of incoherent holes.

In our lattice model we discretize both space and time, and take a “staggered” type of update rule which is completely specified by the dynamics of a $`2\times 2`$ cell (see Fig. 7). We put a single variable $`\varphi `$ on each site, where $`\varphi `$ corresponds to the phase difference (the integral over the phase–gradient $`q`$) across a cell, divided by $`2\pi `$. Local phase conservation is implemented by $`\varphi _l^{}+\varphi _r^{}=\varphi _l+\varphi _r`$, where the primed (unprimed) variables refer to values after (before) an update. Holes are represented by active sites where $`|\varphi |>\varphi _t`$; here $`\varphi `$ plays the role of the internal degree of freedom. Inactive sites are those with $`|\varphi |<\varphi _t`$, and they represent the background. The value of the cutoff $`\varphi _t`$ is not very important as long as it is much smaller than the value of $`\varphi `$ for coherent holes, and is here fixed at $`0.15`$. Without loss of generality we force holes with positive (negative) $`\varphi `$ to propagate only from $`\varphi _l`$ ($`\varphi _r`$) to $`\varphi _r^{}`$ ($`\varphi _l^{}`$).

The details of the translation of these rules into the model can be found in the caption of Fig. 7 and in Ref. mvhmh , with one exception. A “defect” is formed when $`\varphi _l>\varphi _d`$. Here we have adopted two alternative schemes. In the simplest case (A) (studied before in Ref. mvhmh ), we take $`\varphi _r^{}=\varphi _{ad}`$, and $`\varphi _l^{}=\varphi _d1\varphi _{ad}`$. Here we completely fix the new holes. The factor $`1`$ reflects the change in the total winding number associated with the defects. Notice that the overall winding number does not change by exactly $`1`$. This is because $`\varphi _l+\varphi _r`$ is usually slightly different from $`\varphi _d`$. As we will discuss below, to avoid breaking this “phase conservation” we have also studied case (B), where we take $`\varphi _r^{}=\varphi _{ad}`$, and $`\varphi _l^{}=\varphi _l+\varphi _r1\varphi _{ad}`$. Here some (small) scatter in the defect profiles is allowed, but the change in the overall winding number is now strictly $`1`$.

The model does contain a large number of parameters, $`g`$, $`\varphi _n`$, $`\varphi _d`$ and $`\varphi _{ad}`$. We will first discuss the role of $`g`$ and the difference between the two defect rules (A) and (B).

In order for the model to reproduce the correct lifetime dependence of edge–holes in hole–defect states, the coupling of the holes to their background, $`g`$, should be taken negative (although its precise value is unimportant). For $`g=0`$ the hole lifetime $`\tau `$ becomes a constant, independent of the $`\varphi `$ of the state into which the holes propagate, and moreover, the dynamical states are regular Sierpinsky gaskets (Fig. 8a). For $`g<0`$ both the appropriate $`\tau `$-divergence and disorder occur (Fig. 8b), illustrating the crucial importance of the coupling between the holes and the self–disordered background.

It turns out that the model with defect rule (A) is not strictly chaotic; sufficiently small perturbations do not always chaotically spread. However, with defect rule (B) implemented, infinitesimal perturbations do spread (see Fig. 8c,d). To understand this difference consider the fate of a tiny, localized perturbation. Holes will sweep past and be influenced by this perturbation, but since time is discrete, a sufficiently small perturbation will not lead to a change in the time at which a hole evolves to a defect. We have found that after a number of holes have passed over such a perturbation, it can actually be absorbed, so that no chaotic amplification occurs. It is therefore the combination of the discreteness of time and the fixed defect profiles which do not, strictly speaking, lead to chaos. By lowering $`\lambda `$, this problem is diminished, but this makes the model far less effective computationally. Alternatively, we have found that defect rule (B) also circumvents this problem; perturbations can now never be absorbed, due to the nature of the defect rule (B). In this case the defect profile is not entirely fixed, but its scatter is still rather small: for $`\lambda =0.1`$ a typical scatter is of the order of 5%, and this diminishes as $`\lambda `$ is decreased. Therefore we can conclude that, in the continuous time limit of the lattice model, the scatter of the defect profiles is not necessary to obtain chaos. In the remaining part of the paper, we will use model (B) only.

The self–disordering can be very clearly observed in the minimal model, since its update rules unambiguously specify which sites are “background” and which are ‘active”. Two snapshots of the evolution shown in Fig. 8c, are plotted in Fig. 9. These snapshots clearly demonstrate how, after sufficient time has passed, the “inactive” background between the holes has become completely disordered.

The essential parameters determining the qualitative nature of the overall state are $`\varphi _n`$, $`\varphi _d`$ and $`\varphi _{ad}`$. These parameters determine the amount of phase winding in the core of the coherent holes with $`q_{ex}=0`$ ($`\varphi _n`$), and in the new holes generated by the defects ($`\varphi _d`$, $`\varphi _{ad}`$). When varying the CGLE coefficients, these parameters change too, leading to qualitatively different states. In particular, they determine the times $`t_1`$ and $`t_2`$ that we already studied for the full CGLE in section II.3. We found that when $`\varphi _n`$ and $`\varphi _{ad}`$ are both decreased, $`t_1`$ and $`t_2`$ roughly remain the same. We have therefore kept $`\varphi _{ad}=0.65`$, and varied $`\varphi _n`$ and $`\varphi _d`$ to tune the values of $`t_1`$ and $`t_2`$. Notice that the symmetry (or asymmetry) of the defect profile depends on $`\varphi _d1`$; a value of $`\varphi _d<1`$ typically promotes zigzag patterns. We have determined the appropriate values of $`\varphi _n`$ and $`\varphi _d`$ for three concrete cases, tabulated in Table I. Notice that the parameters for $`c_3=1.5`$ precisely correspond to those used in Fig. 8.

As can be seen in Fig. 10, the agreement between the simple model and the CGLE is satisfactory, although clearly the CGLE displays richer behavior. Note that in the full CGLE, small perturbations of the background wavenumber evolve in a non–trivial manner. For example, a nonzero average background wavenumber introduces a drift of the phase perturbations in the background between the holes torcinidrift . Since this phase dynamics is much slower than the hole–defect dynamics, we have chosen to ignore it in the grid model, and this accounts for the difference between Fig. 10a and Fig. 10b.

Finally, we emphasize that the grid model allows us to disentangle the mechanism of hole–defect chaos, by enabling us to completely control the behavior of defects and the coupling between holes and the laminar background. The grid model also has the advantage of being possibly the simplest model that captures the essence of the self–disordered hole-defect spatiotemporal chaos. We also emphasize that we have, for the first time, carried out a detailed comparison between the full CGLE and the grid model, both in our analysis of the spreading of perturbations and in the calibration of the grid model as a function of $`c_1`$ and $`c_3`$. Given the simplicity of the model, the agreement with the full CGLE is striking.

## V Conclusion

In conclusion, we have studied in depth the dynamics of local structures in the 1D CGLE. We have presented strong evidence that the origin of the chaotic behavior in the 1D CGLE lies in the self–disordering action of the holes, rather than in the scatter of the defect profiles. Using this insight, we have then developed a minimal lattice model for spatiotemporal chaos, which, despite its simplicity, reproduces the essential spatiotemporally chaotic phenomenology of the full CGLE.

How general are these results? We conjecture that there are two crucial properties needed for hole–defect type chaos: propagating saddle–like coherent structures (the holes) and a “conserved” field (the phase field). Of course, the phase is not strictly conserved here due to the occurrence of defects, but $`(𝑑x\psi )\text{mod }2\pi `$ is a conserved quantity. It is this conservation that is weakly broken in our grid model for defect rule (A), but is preserved for defect rule (B). Only the latter is truly chaotic. The conservation is also the underlying reason why an evolving hole leaves an inhomogeneous and self–disordered trail behind. Without such a conserved field, there is no reason for “self–disordering” to occur, and the holes then typically will exhibit a fixed lifetime, leading to Sierpinsky gasket type patterns as is often the case in reaction–diffusion models rd . A related scenario appears to occur in the periodically forced CGLE chatepikovsky . Conserved fields of the type described here may be expected more generally for systems undergoing a Hopf bifurcation, and can therefore be expected to also occur in Ginzburg–Landau type equations including higher order terms, and also in experiments. We have argued in Ref. grr , that saddle–type structures like the homoclons here, may be much more general. This leads us to believe that the type of dynamics described here is not an artefact of the pure CGLE, but could be far more widespread.

Our work opens up the possibility for quantitative studies of hole–defect and homoclon dynamics, states which have recently been observed in various convection experiments garnier ; willem . We hope that our simple picture will advance these experimental studies of space–time chaos into the quantitative realm. Local dynamics of the type studied here, such as the dependence of lifetime on initial conditions grr , or measurements of quantities like the daughter hole lifetimes $`t_1`$ and $`t_2`$ should be accessible in experiment, thereby circumventing the difficulties normally associated with characterizing fully–developed chaotic states.

Finally, we mention another commonly observed type of spatiotemporal chaos occurring in systems when a periodic state undergoes a certain symmetry breaking bifurcation faivre . Mathematically, such systems may be described by a complex Ginzburg–Landau equation, coupled to a phase field. Such models are sometimes referred to as $`A\varphi `$ models aphi . Theoretically, the role of holes and defects has not yet been studied in great detail for these systems, but the main ingredients for hole–defect chaos of the type described here appear to be present. We hope that our work will encourage further studies in this area.

M.H. acknowledges support from Stichting FOM and from The Royal Society.