## Abstract

We describe the first dynamical simulations of domain growth during the self-assembly of the gyroid mesophase from a ternary amphiphilic mixture, using the lattice Boltzmann method. The gyroid is a chiral structure; we demonstrate that, for a symmetric amphiphile with no innate preference for left- or right-handed morphologies, the self-assembly process may give rise to a racemic mixture of domains. We use measurements of the averaged mean curvature to analyse the behaviour of domain walls, and suggest that diffusive domain growth may be present in this system.

## 1. Introduction

Emergent phenomena occur when many elements, each with very simple properties, display unexpectedly complex behaviour when gathered together in large numbers (Skår & Coveney 2003; Coveney & Fowler 2005). Such behaviour gives rise to a vast number of patterns and structures in nature. A good example of this is provided by amphiphile molecules, which, despite having a very simple structure in general, can spontaneously assemble into a wide variety of ordered structures; in particular, the membranes of living cells are constructed from these molecules.

One of these structures is called the gyroid mesophase, whose properties, once assembled, have been studied experimentally and theoretically. However, little is understood about the process of self-assembly. Using a lattice Boltzmann model which can capture both surfactant dynamics and fluid flow effects, we are able to simulate the self-assembly of such a phase. In particular, we demonstrate below that the self-assembly process of the gyroid mesophase can give rise to imperfections or defects in its structure.

## 2. Amphiphile mesophases

The word *amphiphile* comes from the Greek words amfi, for ‘both’, and filein, ‘to like’, and means a kind of molecule that is composed of two parts, each attracted to a different species. Soap is a commonly encountered amphiphile: soap molecules are composed of a polar head group which is hydrophilic, or attracted to water, and a long hydrophobic or lipophilic tail group which is repelled by water and attracted to oil. In a mixture of oil, water and soap, the soap molecules are strongly attracted to the interface between the oil and the water, allowing them to sit in the minimal-energy configuration with the head group facing towards water and the tail group facing into the oil. This tendency to migrate towards interfaces is why such molecules are commonly termed *surfactants*, or surface-active agents. Surfactants have immense industrial importance, mainly due to their tendency to sit at interfaces and thereby reduce surface tension; this is the principle by which all detergents work.

While a solution of surfactant in water is easily prepared in the kitchen sink, such solutions show very rich and complicated behaviour, summarized by, for example, Gompper & Schick (1994) and Langevin (1999), or the reviews of Seddon & Templer (1993, 1995). In a very dilute solution, amphiphile molecules exist as monomers, but beyond a certain critical concentration, they may assemble into a variety of morphologies, such as spherical, wormlike, or branching micellar clusters, or bicontinuous structures, comprised of interpenetrating networks of oil and water channels separated by a monolayer of amphiphile. The networks may be isotropic sponges with no long-range order, but can also form liquid crystalline structures with translational symmetry. This behaviour depends on many variables, such as temperature, surfactant volume fraction, non-polar chain length, polar head-group strength (which, in turn, depends on factors such as pH and ionic concentration) and the form of the interaction between tail groups, between head and tail groups and between head groups.

### (a) The Canham–Helfrich Hamiltonian

A surfactant molecule in a monolayer will have a variety of lateral forces acting on it. Repulsion between tail groups and between head groups acts to push molecules apart; surface tension effects act to pull them together. The effects of these forces can often be difficult to separate, and the nature of the individual forces is still poorly understood (Seddon & Templer 1995).

A useful idealization of this system is to think of it as being composed of an oil-filled region, a separate water-filled region and a surfactant membrane separating the two. The oil and water regions may have a highly tortuous morphology, giving rise to a correspondingly complicated shape for the surfactant membrane. However, the complicated interactions between molecules can now be described in terms of their net effect on the surface, and, in particular, their effect on its curvature.

Consider a hypothetical symmetric amphiphile, whose head and tail groups are of the same size and exert the same magnitude of forces. At an interface, the ideal configuration of amphiphile molecules is planar; by symmetry, there is no innate preference for the interface to curve in one direction or another. If one bends the surface, then head groups will be pushed together and tail groups pulled apart, or vice versa; this incurs an energy cost. At any point on a surface, maximum and minimum curvatures can be found, corresponding to the minimum and maximum radii of curvature, respectively; the sign of the curvature indicates its direction. These are called the principal curvatures and ; the mean and Gaussian curvatures are defined as and , respectively. The free-energy cost of curvature at a point on an interface was characterized by Canham and Helfrich (Canham 1970),(2.1)where *κ* is called the *splay curvature modulus* and is the *saddle splay modulus*. This is integrated over the interfacial surface to find the Hamiltonian for the complete system. The Helfrich Hamiltonian, which describes the energetics of a fictional surface of zero thickness (Schröder *et al*. 2004), is still far from the full story for surfactants; packing frustration and hydrocarbon chain stretching are important, for example, in the inverse hexagonal wormlike micellar phases (Seddon & Templer 1993). However, it is sufficient to predict many of the possible mesophase morphologies.

### (b) Bicontinuous triply periodic minimal surface phases

The existence of bicontinuous phases was first suggested by Scriven (1976). In these phases, the amphiphile surface must both extend throughout space and minimize its curvature energy. These constraints are satisfied by a class of structures known as triply periodic minimal surfaces (TPMS). A ‘minimal surface’ has zero mean curvature and negative Gaussian curvature at all points; a TPMS is a minimal surface with cubic symmetry, repeating in the *X*, *Y* and *Z* directions. A brief history of the TPMS is given by Schwarz & Gompper (2002); until 1970, only five TPMS were known, discovered by H. A. Schwarz and his students in the late nineteenth century. Schoen (1970) discovered 12 more, describing them in a NASA technical report. Little attention was paid to the new surfaces, until their existence was rigorously proven by Karcher, who discovered the existence of others in addition (Karcher & Polthier 1996). Further theoretical details of various aspects of periodic minimal surfaces are described in the literature (Mackay 1985; Góźdź & Hołyst 1996; Karcher & Polthier 1996; Klinowski *et al*. 1996; Hyde & Ramsden 2000; Hyde & Schroeder 2003; Lord & Mackay 2003; Enlow *et al*. 2004).

Three surfaces of particular interest in surfactant morphologies are called the *P*, *D* and *G* surfaces (figure 1). *P* and *D* were discovered by H. A. Schwarz; *G*, the ‘gyroid’, was discovered by Schoen. Fogden & Hyde (1999) give a detailed analysis of the *P*, *D* and *G* surfaces. An analytical form for these three surfaces (among others) is known in terms of the Enneper–Weierstrass representation (Fogden & Hyde 1999; Gandy & Klinowski 2000), which maps the fundamental patch of each surface into the complex plane. It can be shown (Gandy & Klinowski 2000) that the *P*, *D* and *G* surfaces all have the same Weierstrass representation, except for a single parameter *θ*, called the Bonnet angle. Variation of *θ* maps the surfaces on to one another, and is called a Bonnet transformation. The Enneper–Weierstrass representation of these surfaces is unfortunately a little cumbersome to work with, since it involves several elliptic integrals. Nonetheless, it has been used to derive many properties of the surfaces, such as the area, Euler characteristic and Gaussian curvature, per unit cell. These have also been derived numerically for the gyroid and related surfaces of constant mean curvature (Große-Brauckmann 1997), by discretizing a surface and permitting it to evolve, under volume constraints, to *H*=0.

## 3. The gyroid

Of these three surfaces, the ‘gyroid’ *G* is perhaps the most interesting. It is the only known TPMS, which is balanced (i.e. the two labyrinths can be mapped onto each other through a Euclidean transformation) while containing no straight lines; it is also the only known TPMS composed entirely from triple junctions. The gyroid has space group symmetry ; the unit cell consists of 96 copies of a fundamental surface patch, related through the symmetry operations (Gandy & Klinowski 2000) of this space group. Channels run through the gyroid labyrinths in the (100) and (111) directions; passages emerge perpendicular to any given channel as it is traversed, the direction at which they do so gyrating down the channel, giving rise to the ‘gyroid’ name (Große-Brauckmann 1997). The labyrinths are chiral, so that the channels of one labyrinth gyrate in the opposite sense to the channels of the other. Looking down the (111) direction of a gyroid shows a distinctive ‘wagon wheel’ pattern (Avgeropoulos *et al*. 1997), which has been observed experimentally in transmission electron micrographs (TEM) of gyroid phases (Shefelbine *et al*. 1999). Gyroids have been observed in many experimental systems, and are usually regarded as the most commonly occurring of the cubic TPMS geometries. They have been seen in triblock copolymers (Shefelbine *et al*. 1999) and lipid–water mixtures (Seddon & Templer 1993, 1995; Mariani *et al*. 1996; Czeslik & Winter 2002). Attempts have been made (Chan *et al*. 1999) to construct a ceramic nanostructured film with the gyroid morphology through the use of a self-assembling polymer template, raising the interesting possibility of fabricating chirally selective porous media. The possibility of using TPMS morphologies for photonic crystals is under active investigation (Urbas *et al*. 2002).

There have been hints that lyotropic liquid crystals may play a part in biological processes such as lipid digestion (Patton & Carey 1979; Rigler *et al*. 1986), and offer insights into cell membrane properties and dynamics (Stier *et al*. 1978; Rand 1981; Prestegard & O'Brien 1987). Donnay & Pawson (1969) suggested that periodic minimal surfaces could be found in nature, and pointed at the microscopic structure of sea urchin skeletons (Nissen 1969) as a possible example; micrographic evidence has since been emerged to suggest that cubic TPMS phases may be present in certain plant cell organelles, and Landh (1995) suggests, with micrographic evidence, that gyroid structures may exist in the endoplasmic reticulum. Schwarz & Gompper (1999, 2000) examined several TPMS morphologies, and suggested that, for oil–water symmetric systems giving rise to zero spontaneous mean curvature, the gyroid was the most stable structure.

## 4. Dynamical simulations

Despite the significant progress (Schwarz & Gompper 1999) made in free-energy models, which through analysis and Monte Carlo simulation give a good understanding of equilibrium properties and phase stability, comparatively little is known about the dynamics of the cubic phases, or indeed mesophases in general; free-energy functionals will not give information about the dynamics of a system far from equilibrium. The limitation to equilibrium states is a severe one, for several reasons. First, real-world systems may contain boundaries or imperfections, which prevent ‘perfect’ equilibrium phases from forming. Second, rheological or otherwise dynamical properties of a mesophase are non-equilibrium by definition, and so cannot be captured by such treatments. Many of the mechanical properties of metals and electrical properties of semiconductors are determined primarily by the nature of defects and impurities, so an understanding of the non-equilibrium behaviour is essential for the investigations of material properties. Mesoscale techniques, such as dissipative particle dynamics or lattice Boltzmann, since they can handle kinetic descriptions, offer a non-equilibrium alternative to the free-energy descriptions. Groot & Madden (1998) performed dissipative particle dynamics simulations of diblock copolymer melts, which allowed reproduction of several well-known morphologies, and also showed the dynamical pathways through which they were reached. In particular, they showed the existence of a ‘gyroid-like’ phase with many triple junctions. However, the phase was not exactly of cubic symmetry, possibly because of finite size effects and frustration. Various course-grained molecular dynamics based mesoscale simulation methods have been developed in the past few years, with applications *inter alia* to surfactant self-assembly (Marrink *et al*. 2004).

Using an early version of the LB3D code, Nekovee & Coveney (2001) showed that it was possible to use the lattice Boltzmann model of Chen *et al*. (2000) to simulate the assembly dynamics of lamellar and bicontinuous phases of binary water–surfactant systems; indeed, they showed the existence of a *P*-surface cubic phase of the model. Then, while investigating the effect of the presence of surfactant on spinodal decomposition, González-Segredo & Coveney (2004*a*,*b*) discovered the existence of a gyroid cubic phase. Importantly, use of the lattice Boltzmann method allowed modelling of the dynamics of the self-assembly of the gyroid. These lattice Boltzmann simulations of the gyroid phase demonstrated several important features. First, it had been suggested (Prinsen *et al*. 2002) that simulation of cubic phases might require a model with more complicated amphiphiles than symmetric dimers, and also that long-range interactions might be required. The lattice Boltzmann simulations, which used symmetric amphiphiles with short-range interactions, refuted this. Second, analysis of the X-ray structure factor of the simulation data showed the existence of oscillating modes at long times; visualization of the simulation output strongly suggested that these were due to Marangoni effects (Scriven & Sternling 1960; Grotberg & Gaver 1996).

Theissen *et al*. (1998) describe a lattice Boltzmann model of amphiphilic systems which is based around a Ginzburg–Landau free-energy functional (Gompper & Schick 1990; Schwarz & Gompper 2000); however, it should be noted that this model lacks an explicit surfactant density, making the assumption that the surfactant is all adsorbed onto the oil–water interface. It is not expected, therefore, that this model would reproduce Marangoni effects at the interface. Lamura *et al*. (1999) describe a Ginzburg–Landau model with explicit surfactant density, but their model, in turn, lacks explicit surfactant orientation. The model used in LB3D has both explicit density and orientation for surfactant particles.

The gyroid self-assembly simulations used initial conditions of a randomized mixture of oil, water and surfactant, with the oil and water in 1 : 1 ratio. Phase separation would occur rapidly, within the first thousand time-steps or so. After phase separation, the system would form a ‘molten gyroid’ phase, consisting of many oil or water rods surrounded by surfactant; these rods would join up to form triple junctions, giving rise to a gyroid morphology. Towards the end of these simulations, after around 30 000 time-steps, the Marangoni oscillations were observed. These simulations were originally performed to investigate the dynamics of surfactant-limited phase separation; discovery of the gyroid phase was an interesting side effect. Owing to this, the simulations of González and Coveney were not optimally suited to examining the gyroid phase, in several ways. First, they were limited in size: most of the simulation work was done with 64^{3} systems. These are perfectly adequate for phase-separation studies; however, it was observed that while at early times there might be several different gyroid regions with different orientations, at long times these would join up to form a single gyroid grain spanning the entire simulation grid, a so-called ‘perfect gyroid’. In the real world, cubic mesophases are far from perfect: despite careful experimental procedures (Hajduk *et al*. 1994), artificially created gyroid materials may not consist of a perfect gyroid repeating throughout space, but rather of many gyroid grains with different orientations. In addition, there may exist dislocations and defects within these grains, analogous to the defects found in ordinary solid crystals (Feynman *et al*. 1964; Kittel 1996). Any investigation of the gyroid domains or the nature of their interactions would require a simulation of sufficient physical size to accommodate several domains for the majority of the simulation.

Second, they were limited in time: the time-scales probed (maximum of 35 000 time-steps) offered few hints as to how transient the effects observed were. Phase separation was observed to happen extremely quickly, over a 1000 time-steps or so, and formation of gyroid morphology on a slower but comparable time-scale. The behaviour of gyroid defects, on the other hand, takes place on a much slower scale, at least tens of thousands of time-steps, and these time-scales were not probed. With these limitations in mind, a new set of simulations was performed specifically to look at the nature of defect dynamics in the gyroid liquid crystal phase.

## 5. The TeraGyroid simulations

It was clear that larger simulations were required; how much larger was not so clear, since defects had not been observed in the unambiguous absence of finite size effects. Defect dynamics simulations were performed as part of a larger project (Blake *et al*. 2005) in which a transatlantic supercomputing grid was constructed. As part of this project, a very substantial amount of computing power became available for use. The simulations of most interest used the parameter sets (González-Segredo & Coveney 2004*a*) called numbers 8 and 9, which are shown in table 1.

Simulations were run at sizes of, among others, 64^{3}, 128^{3} and 256^{3}, for durations listed in table 2. If a stable gyroid state was reached, then a simulation was terminated; the systems which did not reach this state ran for as long as possible under the resource constraints at the time. During each simulation, the order parameter field was stored to disk, as well as the surfactant density field . This was performed every time-steps, for values of also shown in table 2. The output rate was chosen according to several constraints; it had to be sufficiently rapid to allow observation of the gyroid dynamics, but limits were imposed by local disk quotas, which would often shift due to data from other concurrent jobs.

## 6. Analysis techniques

### (a) The structure function

A commonly examined property of mesoscale systems is the *structure function* , defined as the spherically averaged Fourier transform of the order parameter ,(6.1)This is very closely related to the ‘structure factor’ measured in small angle X-ray scattering (SAXS) experiments, which are commonly used in experimental examination of mesophases (Avgeropoulos *et al*. 1997; Laurer *et al*. 1997; Sakurai *et al*. 1998; Shefelbine *et al*. 1999), and has been an important tool in the analysis of many mesoscale simulations (Velasco & Toxvaerd 1993; Rybka *et al*. 1995; Wagner & Yeomans 1998; González-Segredo & Coveney 2004*b*). To measure the dominant length-scale of a system, the reciprocal space first moment of is usually calculated,(6.2)The corresponding length-scale is then .

### (b) Polygonal approximation of the interface

Since the interfacial properties are so crucial to the behaviour of the gyroid phase, it is important to be able to extract the location of the two-dimensional interfacial surface from the three-dimensional simulation data. In a continuum system, the interface between the oil and water phases can be defined as the implicit surface , for . In contrast, the simulation output data consists of the values of *ϕ* at discrete sites on the lattice; it is not quite so straightforward to define an interfacial surface because, in general, while can be positive or negative, corresponding to the lattice site containing a majority of red or blue particles, respectively, it is almost never exactly zero. However, a continuum order parameter field can still be estimated from the discrete simulation output, by linear interpolation; this interpolated function is single-valued and continuous across the simulated region, and can be regarded as the simulated approximation to the order parameter field; the isosurface is then well defined.

Generation of a polygon mesh approximation to this surface can be achieved through several means. The most popular method is known as marching cubes, this technique was developed by Lorensen & Cline (1987). In marching cubes, each point on the lattice is assigned a value of 1 if is ‘negative’, or unambiguously less than zero, and 0 if is ‘positive’, i.e. greater than zero or within floating-point error of zero. If a positive site is adjacent to a negative site, then changes sign somewhere along the link between the two sites, and therefore the link between the sites intersects the contour surface; the position of these intersections can be found easily through interpolation. The marching cubes algorithm proceeds by considering each cube of eight sites in turn, calculating a number between 0 and 255, representing the state of these sites, and then using this number to find a set of polygons approximating the interface shape inside the cube. The polygon sets can be stored in a lookup table, making the algorithm extremely fast.

Unfortunately, the algorithm is not guaranteed to produce a surface with the same topology as the implicit interpolated surface . In fact, it is possible for marching cubes to produce ‘non-manifold’ surfaces containing holes; this arises because some of the 256 possible states of a cube have more than one way of being tiled with polygons, and the tilings are topologically different.

Several approaches have been suggested to correct this problem. One such technique is to tile the space with tetrahedra rather than cubes (Carneiro *et al*. 1996; Treece *et al*. 1999), which produces a topologically consistent surface without holes, but at the expense of speed, memory and geometrical accuracy (Lewiner *et al*. 2003). Another approach is the recursive dividing cubes algorithm (Cline *et al*. 1988), which has similar performance penalties, particularly for large datasets.

Chernyaev (1995) observed that it was possible to extend the marching cubes approach to produce polygon surface meshes which are not only topologically consistent (i.e. do not have any holes), but also topologically correct (i.e. homeomorphic to the surface ). An implementation of this technique was described by Lewiner *et al*. (2003), who used lookup tables to determine which tests to perform, resulting in a very efficient algorithm to produce a triangle mesh which is both topologically consistent and topologically correct.

### (c) Fluid properties at the interface

By applying this algorithm to the simulation output datasets, a polygonal approximation to the interface between the two fluid components can be found. This is useful not only for direct visualization of the system, but also for calculating many interesting properties.

First, the interfacial area of the system can be calculated, simply by adding up the areas of all of the triangles in the surface. Without knowing the location of the interface, any space-varying property , such as surfactant density, can only be averaged over the bulk of the three-dimensional system. However, once the interface has been polygonized, the variation of *α* over the interface can now be found, by interpolating *α* to the surface vertices and taking an area-weighted average. Consider a single triangle , area , vertices . The value of *α* averaged over the triangle is(6.3)This relationship is exact if *α* varies linearly over the surface of the triangle. Integrating *α* over the entire closed polygon mesh *M* is equivalent to summing the values integrated over each individual triangle,(6.4)Hence, the average of *α* over the surface is(6.5)

### (d) Interfacial curvature

Any point on the surface has a unit normal , defined as . The shape operator is a rank-2 tensor defined as .

The Gaussian and mean curvatures can then be found from the determinant and half the trace of the shape operator, respectively. More usefully, the surface-averaged value of the squared mean curvature, , can also be calculated; this quantity is physically interesting because it appears as a term in the Helfrich Hamiltonian. The contribution to the surface-averaged curvature from each lattice site can also be calculated. This curvature field has low values where there is no interface, or where the interface has a minimal surface structure, and it has a high value where the interface is highly curved; consequently, it can be used to distinguish non-minimal-surface defect regions from the bulk gyroid phase (figure 2).

### (e) Interfacial Euler characteristic

Consider a two-dimensional closed surface tiled with a polygon mesh; each polygon is joined to its neighbours by a straight edge, and the edges meet at point vertices. Suppose that the surface has *V* vertices, *E* edges and *F* polygon faces; the Euler characteristic or Euler number of the surface is defined as . A surprising and very powerful result of elementary topology (Kinsey 1993; Nakahara 1993; Armstrong 1997) is that the Euler characteristic depends only on the topology of the surface, and not on its shape or the manner in which it is tiled with polygons. More precisely (Kinsey 1993), if and are compact, connected surfaces without boundary, then is topologically equivalent to if and only if , and either both are orientable or both are non-orientable. A surface is non-orientable if it contains Möbius bands; all of the surfaces under consideration here are orientable, so determination of topological equivalence reduces to comparison of the Euler characteristic.

The Euler number is additive; a system with *n* surfaces each with Euler number has total Euler number . The surface meshes produced by the marching cubes algorithm are composed entirely of triangles; this property simplifies calculation of the Euler characteristic. Each triangle on the mesh has three edges; each edge on the mesh joins exactly two triangles. Hence, the number of edges is 3/2 times the number of triangles, .

The Euler number of a surface is directly related to another topological quantity called the genus *g*, which is often regarded as simply the number of holes in the surface. For an orientable surface (Kinsey 1993), . The Euler number is an extremely useful concept when trying to understand the output of mesoscale simulations. Imagine a lattice Boltzmann simulation of a spherical-droplet phase, where the droplets are defined as regions for which the order parameter . Since the Euler number of a sphere is 2, and since it is additive, the number of droplets in the system is simply half the Euler number of the *ϕ*=0 surface. Moreover, the Euler number, being a topological invariant, is insensitive to fluctuations in the droplet shape or size, provided that each droplet remains simply connected.

Recently, the Euler number was used by Hołyst & Oswald (1998) to characterize fluctuating wormhole-like passages (Goos & Gompper 1993) appearing in Monte Carlo simulations of lamellar phases; similar wormholes have been observed in a lattice gas dynamical model (Boghosian *et al*. 1996). In a system with periodic boundary conditions, a membrane spanning both the *x* and *y* directions is topologically equivalent to the surface of a torus, and therefore has Euler number zero; thus, a lamellar stack of membranes will also have Euler number zero. However, if a wormhole forms between two lamellae, the Euler number will suddenly jump to −2, hence the Euler number is for *n* such wormholes. The Euler number has also been used in examinations of TPMS surfactant morphologies (Schwarz & Gompper 2000) and experimental MRI data (Gerig *et al*. 1999). A gyroid surface has Euler characteristic −8 per unit cell.

## 7. Analysis of the simulation data

So far, both the techniques used to simulate assembly of the gyroid mesophase and the techniques used to analyse the resulting data have been presented. The following sections present the actual analysis of the generated data, and a discussion of the physics involved.

### (a) Initial phase separation

The simulated system is initialized to a very highly mixed and disordered state, often called a ‘quench’ in the literature. The system undergoes spinodal decomposition (Bray 1994; Kendon *et al*. 1999; Chin & Coveney 2002; González-Segredo *et al*. 2003), and rapidly divides into interpenetrating regions of bulk water and bulk oil, separated by a layer of surfactant. If the system contained only water and oil, each component would eventually form a single large domain; however, this situation would leave insufficient surface area to allow all of the surfactant to sit at the oil–water interface. Therefore, the surfactant has the effect of limiting the domain growth (Love *et al*. 2001).

Figure 3 shows the dominant length-scale for the first 5000 time-steps of three simulations using parameter set 8. This shows that the length-scale rapidly increases as the oil and water components separate, but stays relatively constant after around 1000 time-steps. Figure 4 shows renderings of the oil–water interface at *ϕ*=0 over the first 500 time-steps, indicating that the rapid increase in structure factor corresponds to the formation of a mesh of interconnecting channels, all of roughly the same width. Visual inspection of the interface shows that no gyroid structure is present during the channel formation.

### (b) Gyroid formation: parameter set 9

Figure 5 shows the mean curvature *H* of the *ϕ*=0 interface, averaged over the whole surface. The most obvious feature of this graph is the magnitude of the curvature fluctuations in each system. The 256^{3} system, being the largest, had the most interfacial area over which to average; therefore had the lowest amount of fluctuation in this system, with a standard deviation of 8.41×10^{−5} inverse lattice units. The 128^{3} system had fluctuations roughly twice as large, with a standard deviation of 1.70×10^{−4}, and the 64^{3} a little over twice as large again, at 3.86×10^{−4} inverse lattice units.

The next obvious feature of figure 5 is that, while the 128^{3} and 256^{3} systems fluctuated about zero for the entire simulation duration, the 64^{3} simulation dropped to a negative-mean-curvature geometry after around 10^{5} time-steps. Visual inspection showed that this corresponded to the system collapsing into a stable state, forming a regular structure which filled the entire simulation lattice, apart from two point defects.

Looking along the *y*-axis, the structure reconnected with itself through the periodic boundary conditions, moving half a unit cell in each of the *x* and *z* directions after traversing the lattice. There were seven unit cells in the *x* and *z* directions, and six and a half in the *y* direction; a topologically perfect gyroid with this many unit cells would have an Euler characteristic of −2548. The simulated structure had due to the presence of the point defects. The skewing of the ‘gyroid’ in the *x* and *z* directions and consequent rhombohedral (rather than cubic) liquid crystal structure meant that the oil–water interface did not form a zero-mean-curvature gyroid structure, but a skewed gyroid with, on average, negative mean curvature. It should be noted that the production of a gyroid with this orientation is a reassuring sign that the production of a gyroid mesophase is not a simulation artefact due to some anisotropy induced by the lattice, in which case such misalignment would not be expected.

Figure 6 shows the *ϕ*=0 isosurface, coloured red on the oil-majority side and blue on the water-majority side, at time-step 2.5×10^{5}. The region shown lies at the interface between the two domains; the domain wall is shown running vertically through the centre of the diagram. Careful examination of the channels running through the diagram shows that, while the two domains contained gyroids with roughly the same orientation, the domains had opposite chirality: red channels spiral clockwise away from the viewer in the right-hand side domain, but anticlockwise away from the viewer in the left-hand side domain. Since the model used is oil–water symmetric, there should be no preference for, say, oil channels to be left-handed, so the existence of chiral domains is to be expected.

Figure 2 shows volume renderings of the curvature contribution, highlighting regions with high interfacial curvatures. At time-step 10 000, the gyroid structure had not emerged, but by time-step 40 000, multiple gyroid domains were visible. These domains were observed to merge, until only two domains were left in the system, around time-step 200 000. One of the domains then collapsed over the next hundred thousand time-steps, leaving behind two columnar defect regions which persisted for as long as the simulation continued.

Figure 7 shows a volume rendering of the order parameter field at time-step 5×10^{5}, with the transfer function chosen to highlight regions of high oil density and regions of high water density in red and blue, respectively. By this time, the domain walls had collapsed into a pair of dislocation lines running right through the system. The green lines show a Burgers loop (Kittel 1996) around each dislocation; it can be seen that the Burgers vectors for each dislocation are equal in magnitude and opposite in direction.

### (c) Gyroid formation: parameter set 8

Figures 8 and 9 show the surface-averaged mean curvature and surface-averaged squared mean curvature, respectively, for simulations using parameter set 8. The fluctuations in the mean curvature again dropped in magnitude with increased system size, just as observed with parameter set 9. On closer examination, the fluctuations for the 128^{3} simulation dropped in magnitude around time-step 2.5×10^{5}, with a similar halt in the evolution of the other bulk and interface properties. The domain wall visualizations (figure 10) show that the system produced two domains at late times, and one of these disappeared around time-step 2.5×10^{5}, leaving no dislocations or other major defects in the system. It was verified by direct inspection that the last two domains were again of opposite chirality. At the end of the simulation, time-step 378 500, the system was found to contain a gyroid with 15 unit cells in each direction, with an Euler characteristic of −26 828, corresponding to 3353.5 gyroid unit cells. A 15^{3} periodic gyroid would be expected to contain 3375 unit cells; the discrepancy is attributed to persistent point defects, which appeared as the domains formed, up to time-step 10^{5}, and remained unchanged for the rest of the simulation.

The 256^{3} simulation, limited only by available computer time, ran until time-step 57 500, at which point it was still composed of multiple chiral gyroid domains. The domains were observed to interpenetrate in a manner not unlike the morphology of single-component fluid domains during spinodal decomposition.

Figure 11 shows the averaged squared mean curvature plotted against time for 128^{3} and 256^{3} simulations of parameter sets 8 and 9. Both of the systems which used parameter set 8 appeared to take on fairly close power-law behaviour (seen as straight lines on the log–log graphs) after time-step 10^{4}, when the system recognizably contained gyroid regions (see figure 12). The parameter set 9 systems also did this, but took longer to do so; the value of for the set 9 systems was somewhat larger than that of the set 8 systems, until around time-step 5×10^{4}, at which point they appeared to show the same behaviour. The 128^{3} systems eventually deviated from the power-law behaviour as finite-size effects set in, around time-step 2×10^{5} for set 8 and 3×10^{5} for set 9.

Figure 13 shows renderings of the domain walls and defects of the parameter set 9 system at points before and just after it reached power-law behaviour, with the defects for the parameter set 8 system plotted for comparison. During the period when the set 8 system showed power-law scaling but the set 9 system did not, the latter appeared to have a rather higher density of defective regions, shown as dark regions in the volume renderings. After time-step 50 000, both systems appeared clearly separated into domains.

The datasets for for each simulation were trimmed to only include the time-steps after 10 000, when gyroid structures form, for the parameter set 8 simulations, and the time-steps after 50 000, after clear domain formation, for the parameter set 9 simulations; they were also restricted to the time-steps before the sudden drop in and stabilization of the squared mean curvature for the 128^{3} simulations, to exclude the corresponding finite size effects.

A relation of the form was then fitted to each trimmed dataset; the fitted values of the exponent *n* are given, along with the relevant time periods, in table 3. The fitted exponents together had a mean of −0.481 and standard deviation of 0.037. This is quite close to an exponent of ; the data points to which these fits were made are plotted together in figure 11, along with a power law for comparison. A possible reason for power-law scaling and for why the exponent should be −1/2 is given below.

Consider a large system of volume *V*, containing many gyroid domains. Let the length-scale of a typical single domain be , and suppose that the gyroid length-scale grows as a power law of time, , just as domains often do, e.g. in spinodal decomposition. The surface area of a single domain will then scale as , and the volume of a single domain as .

The total number *N* of domains in the system is , so . The total surface area *A* of domains is proportional to both the number of domains and the typical surface area of a single domain, so therefore . If the domain walls all have roughly the same width *λ*, then the total volume of the domain walls scales as . Therefore, the contribution of the domain walls to any extensive property of the system would also be expected to scale as .

A gyroid is a minimal surface, so one would expect the contribution to the average squared curvature of the system from the well-formed gyroid regions inside domains to be small. On the other hand, the domain walls are non-gyroidal, and have non-zero mean curvature; if there is a roughly constant contribution to per unit volume of domain wall, then could also be expected to scale as . If this is the case, then that would further suggest that the length-scale of the gyroid domains indeed scales as , indicative of random-walk or diffusive behaviour.

## 8. Conclusions

Simulations of the self-assembly of a surfactant gyroid mesophase were performed; for the first time, the simulations were sufficiently large to clearly show the formation and evolution of multiple gyroid domains.

All systems underwent a rapid phase separation, during which oil and water phases formed a network of channels, separated by surfactant. After the phase separation, the systems began to form a gyroid mesophase. At late times, when clear domain walls had formed, power-law scaling of the average squared mean curvature of the system was observed, with exponent close to −1/2. This may indicate diffusive scaling of the length-scale of the domains themselves.

Analysis of simulated mesophases is non-trivial, due to the difficulty of defining a gyroid order parameter. In a ferromagnetic system, one can use the magnetization as an order parameter; all points within a ferromagnetic domain will have the same value of . In an ordinary fluid phase transition such as spinodal decomposition, the density of one fluid at a point can be subtracted from the density of the other, and normalized to give a scalar order parameter which varies from +1 in one fluid to −1 in the other, and which is zero in disordered, mixed regions. This kind of localized order parameter is sufficient to characterize the rapid separation of components before a mesophase is formed; however, after the initial demixing, further ordering of the system is morphological in nature, as gyroid regions are formed. Such a localized order parameter cannot describe this sort of ordering, since it is defined on too small a length-scale, and is not morphological in nature.

A set of tools was developed for calculating a polygonal approximation to the *ϕ*=0 surface of a simulated mesophase, in a topologically consistent manner; calculations of useful quantities such as the topologically invariant Euler characteristic, total interfacial area and interfacially averaged squared mean curvature were then possible. The average squared mean curvature was observed to vary as a power law with exponent approximately −1/2 over the régime where domains are clearly formed but finite size effects have not yet set in; this may indicate random-walk motion of the domain walls. Direct visualization of the interface shows the existence of chiral domains, as should perhaps be expected for a symmetric amphiphile which has no intrinsic chirality. Once domains effectively disappeared due to reaching the simulation system size, the formation of line dislocations analogous to those in solid crystals was observed.

These simulations have opened up several avenues for further work. Even disregarding the physics involved, the techniques for characterization of lyotropic liquid crystal phases could still be improved, by looking at, for example, localized contributions to the interfacial curvature or Euler characteristic. Gyroid domain walls were observed to have a high density of points where both the local order parameter *ϕ* and the local order parameter gradient were small; while this may be due to the structure of the *ϕ*=0 surface inside domain walls, an improvement of this technique beyond a simple visual tool could lead to an algorithm for spatial localization of domain walls, and, consequently, the ability to make more quantitative analyses of the domain behaviour.

The simulations raise a number of physical questions. Those performed so far have all had ‘oil’/‘water’ symmetry and a symmetric amphiphile, a situation which may be hard to realize exactly in experiment. The simulations produced chiral domains; it would be interesting to see if these could ever be produced experimentally, since the existence of such domains would affect the feasibility of producing gyroidal chirally selective filters. Observations of gyroid chirality would be difficult to make using TEM or SAXS techniques, but might be possible using, for example, transmission electron microtomography (Laurer *et al*. 1997).

Finally, it is interesting to note that the lattice Boltzmann method allowed the bridging of two different length-scales in these simulations: the scale (approx. five lattice sites) of the single-fluid-component domains, and the scale (greater than or approx. 50 lattice sites) of gyroid domains.

## Acknowledgements

We thank EPSRC for funding this research under RealityGrid grant GR/R67699 and GR/R94916/01. Most of the scientific data on which this paper is based were obtained during the TeraGyroid project (http://www.realitygrid.org/TeraGyroid.html), jointly funded by NSF and EPSRC, in which we were able to use over 6000 processors from a federated computational grid formed from the US TeraGrid and the UK national supercomputing centres at HPCx and CSAR, to generate and visualize the results over a very short period of time. We are grateful to Thomas Lewiner for helpful discussions and triangulation code.

## Footnotes

- Received December 24, 2005.
- Accepted May 11, 2006.

- © 2006 The Royal Society