## Abstract

The purpose of the present paper is to report on the first computational study of the dynamical and rheological response of a self-assembled diamond mesophase under Couette flow in a ternary mixture composed of oil, water and an amphiphilic species. The amphiphilic diamond mesophase arises in a wide range of chemical and biological systems, and a knowledge of its rheological response has important implications in materials science and biotechnological applications. The simulations reported here are performed using a kinetic lattice–Boltzmann method. Lyotropic liquid crystals exhibit characteristic rheological responses in experiments that include shear-banding and a non-Newtonian flow curve as well as viscoelasticity under oscillatory shear. Their behaviour under steady and oscillatory shear is correctly reproduced in our simulations. On cessation of shear, as the morphology returns to the diamond phase, the relaxation of the stress response follows a stretched-exponential form for low initial strain rates.

## 1. Introduction

Complex fluids composed of oil, water and amphiphilic species can undergo self-assembly to form triply periodic, bicontinuous structures at specific concentrations of the components and/or under particular temperature and pressure conditions (González-Segredo & Coveney 2004; Conn *et al*. 2006). The most commonly studied cubic phases are the bicontinuous gyroid (G, spacegroup *Ia3d*), diamond (D, spacegroup *Pn3m*) and primitive (P, spacegroup *Im3m*) phases. Here, we will consider ternary amphiphilic mixtures symmetric in the oil and water composition. The periodic bicontinuous structures consist of a three-dimensional network of interpenetrating yet non-intersecting channels filled with oil and water separated by an amphiphilic monolayer (Chin & Coveney 2006) in which the hydrophobic and hydrophilic ends of the amphiphile particles are oriented towards oil and water, respectively. The stability of these structures arises due to low bending energy as a result of the characteristic saddle-shaped surfaces of zero spontaneous curvature and the reduction of surface tension owing to the presence of amphiphiles at the interfaces (Lynch & Spicer 2005). The first bicontinuous structures were observed experimentally in soap–water systems and later in other systems including biological lipids. These phases are ubiquitous in nature; for example, they are observed in biomineralization processes and in lung tissue, being formed by both low molecular weight surfactants and amphiphilic polymers (Schwarz & Gompper 2002; Lynch & Spicer 2005). More recently, triply periodic minimal surface-containing cubic phases have been finding applications in consumer products, food technology (Mezzenga *et al*. 2005*b*), in materials science for the synthesis of mesoporous nanomaterials with interesting structural and electronic properties (Schwarz & Gompper 2002), in biosensing and ultrafiltration, protein crystallization (Grabe *et al*. 2003) as well as drug delivery applications in biotechnology (Shearman *et al*. 2006). Furthermore, elucidation of the mechanism and kinetics of self-assembly and phase transitions in cubic phases should provide insight into cellular processes, an active area of experimental investigation recently (Squires *et al*. 2005; Angius *et al*. 2006; Conn *et al*. 2006; Seddon *et al*. 2006).

The non-equilibrium behaviour of the amphiphilic cubic phases and their response to external fields such as an applied shear flow are commonly observed in experiments but not well understood theoretically (Chen *et al*. 1997). The rheological response of these fluids arises as a result of coupling between the microscopic dynamical processes governed by local interactions and the macroscopic hydrodynamic flow field. Applied shear may induce macroscopic transitions to morphologically and/or orientationally distinct self-assembled states; it can induce order–disorder transitions as well as order–order transitions from one phase with long-range order to a distinct phase also with long-range order (Chen *et al*. 1997; Mezzenga *et al*. 2005*a*,*b*; González-Segredo *et al*. 2006; You *et al*. 2007). In addition to oscillatory and steady shear, other flow fields, such as roll casting and extensional flow, are used in the processing of various self-assembled systems (Chen *et al*. 1997). Applied fields can be used to process materials, for example leading to defect-free assembly of mono-domains and mono-crystalline samples (Williams & Mackintosh 1994) and in the selection of particular orientations of the amphiphilic liquid crystal relative to the applied shear flow. A knowledge of the variational response of liquid crystals to a range of processing conditions opens up possibilities for the rational design and fabrication of functional materials for switchable optical, transport and electronic properties and of ordered nanostructures such as nanopores, nanosheets and nanofibres with desired directional properties (Fahmi *et al*. 2006). However, in order to exploit the characteristic response of amphiphilic liquid crystals to various types of external fields, one has to understand the physical basis of the effects of processing conditions such as the shear frequency, strain, temperature and pressure. A universal description applicable to a wide range of amphiphilic liquid crystals may be independent of the microscopic characteristics of the amphiphilic system and instead depend on the mesoscopic response of the system to the applied shear field.

The self-assembly of ternary amphiphilic mixtures into the gyroid, diamond, primitive and sponge phases has been studied previously using our lattice–Boltzmann (LB) approach (González-Segredo & Coveney 2004; Chin & Coveney 2006; Giupponi *et al*. 2006; González-Segredo *et al*. 2006; Saksena & Coveney 2008). Rheological simulations of the gyroid phase under Couette flow for 64^{3} and 128^{3} lattice-site models exhibit non-Newtonian viscoelastic behaviour (Giupponi & Coveney 2006; Giupponi *et al*. 2006; González-Segredo *et al*. 2006), as is widely observed in rheological experiments of self-assembled amphiphilic systems (Bonn *et al*. 1998; Pitzalis *et al*. 2000; Porcar *et al*. 2003; Mezzenga *et al*. 2005*a*; Pouzot *et al*. 2007); in these LB simulations, the gyroid mesophase is found to possess enhanced viscosity as compared with the sponge and lamellar mesophases. The cubic diamond phase is observed to form quite commonly in amphiphilic fluid mixtures (Squires *et al*. 2005; Conn *et al*. 2006; Seddon *et al*. 2006; Shearman *et al*. 2006) and some experimental studies of its rheological properties have also been reported in the literature (Pitzalis *et al*. 2000; Mezzenga *et al*. 2005*a*). However, the self-assembly and rheology of the diamond phase have not received as much attention in simulation studies (Marrink & Tieleman 2001, 2002; Saksena & Coveney 2008) as the gyroid phase. We note here that, although treatment of spontaneous self-assembly of amphiphilic species into mesoscale structures is currently beyond all-atom molecular dynamics simulations because of the large length and time scales involved, considerable progress has been made in the application of coarse-grained molecular dynamics (CGMD) simulations to such problems (Michel & Cleaver 2007; Klein & Shinoda 2008). Ellison *et al*. (2006) have simulated the spontaneous self-assembly of the gyroid phase for hard- and soft-repulsive pear-shaped particles using coarse-grained Monte Carlo and molecular dynamics methods, respectively. CGMD has also been used to treat interactions of self-assembled lipid membranes with embedded proteins (e.g. Ayton *et al*. 2007; Treptow *et al*. 2008).

In this work, we focus on the rheology of the self-assembled diamond mesophase simulated using the kinetic LB method under applied Couette flow. Briefly, the kinetic LB method takes into account the hydrodynamic flow under shear as well as quasi-molecular interactions that describe the coupling between the applied shear field and microscopic dynamical processes within the system. Thus, the method provides access to chemical detail as well as longer (continuum) time scales compared with microscopic simulation methods such as molecular dynamics. The kinetic LB method differs from other LB methods based on the relaxation of a suitably chosen free energy functional (Theissen *et al*. 1998; Lamura *et al*. 1999; Denniston *et al*. 2001, 2004; Cates *et al*. 2005; Yeomans 2006). In the kinetic LB method, the system dynamics and its long-time behaviour (including the approach to equilibrium) are emergent from the quasi-molecular interactions specified in the particular model: we do not *a priori* impose a form for the free energy functional that would govern the self-assembly dynamics or the rheological response.

The paper is organized as follows. In §2 we describe our LB method for simulation of ternary amphiphilic fluids, as well as the visualization and analysis tools and computational infrastructure used in this work. In §3, we describe results from our simulations of diamond phase self-assembly, while, in §4, we present rheological simulations under both steady and oscillatory shear. In §5, we present our conclusions from this work.

## 2. Modelling and simulation

In this section we describe the ternary amphiphilic LB fluid model and simulation methodology. We also discuss visualization and analysis tools and the computing infrastructure used in this work.

### (a) LB simulation method

We use a kinetic, mesoscale LB method that, in the limit of low Mach number, corresponds to the Navier–Stokes equations for nearly incompressible single-phase fluids. The amphiphilic species is modelled as a point dipole with one end having affinity for the oil species and the other end having an affinity for the water species. The strength of the amphiphile–amphiphile coupling, the repulsion strength between the oil and water species and the coupling of amphiphile to the oil and water species are governed by three independent interaction parameters. The kinetic LB model has been implemented in a Fortran 90 code called LB3D, which has been parallelized with a message passing interface using a domain decomposition strategy. The code has been described in detail previously (Love *et al*. 2003; Giupponi *et al*. 2006) and here we briefly mention salient points of the model. The evolution of the LB single-particle velocity distribution functions, and belonging to component *σ* (oil or water) and amphiphile (*s*), respectively, and with discrete velocity *c*_{i}, is given by the following set of equations:(2.1)(2.2)where *τ*^{α} is the relaxation time for the species *α*=(*σ*, *s*). The terms , , and are matrix elements of the collision operators that result from mean-field interactions among different fluid components (Nekovee *et al*. 2000). The local equilibrium distribution functions, , are calculated by shifting the average velocity by an increment proportional to the net interaction force on the corresponding fluid species,(2.3)The time evolution of amphiphilic dipoles ** d**(

**,**

*x***) is given by(2.4)(2.5)where**

*t***is the average dipole vector per site, whose relaxation is governed by the Bhatnagar–Gross–Krook (BGK) equation (2.5). is the average dipole vector at site**

*d***prior to collision and is calculated according to equation (2.4). This simple model of a ternary amphiphilic mixture gives rise to self-assembly into a range of periodic structures and complex transition behaviour under steady shear, which is consistent with experiments and simulations (Chin & Coveney 2006; Giupponi**

*x**et al*. 2006; González-Segredo

*et al*. 2006; Saksena & Coveney 2008). These simulations allow us to directly observe the evolution of the initial mesophase as it undergoes morphological transitions under shear and to determine how the interactions between the component species affect the dynamics. In Couette flow simulations, shear is applied in the

*xz*-plane using the Lees–Edwards boundary condition (Giupponi

*et al*. 2006) with the shear velocity in the

*z*-direction and the shear gradient in the

*x*-direction. For a particle subject to Couette flow under Lees–Edwards boundary conditions, the new position and velocity along the flow direction are given, respectively, by the following:(2.6)(2.7)where

*Δ*

_{z}≡

*U*Δ

*t*;

*U*is the shear velocity;

*u*

_{z}is the

*z*-component of the average fluid velocity

**; and**

*u**N*

_{x}is the system length in the

*x*-direction. Additionally, an interpolation scheme suggested by Wagner & Pagonabarraga (2002) is used as

*Δ*

_{z}is not generally a multiple of the lattice site. For the cubic

*L*

^{3}lattice site systems we use in these simulations, the relationship between shear rate and applied shear velocity is given by . The minimum shear velocity we have used in our simulations is

*U*=0.05. The maximum applied shear rate is determined by the requirement that the shear velocity should be much smaller than the speed of sound in the D3Q25 lattice used here (González-Segredo

*et al*. 2006) for the stability of the LB algorithm.

The components of the pressure tensor are computed as the sum of a kinetic term and a virial mean-field term (González-Segredo *et al*. 2006)(2.8)The summation in the second term on the right-hand side runs over all combinations of the three fluid components, oil (*r*), water (*b*) and amphiphile (*s*). Since the interaction matrix is symmetric with only the nearest neighbour interactions considered, the virial term reduces to(2.9)Under applied Couette flow, the negative of the off-diagonal component of the pressure tensor in the shear plane, −*P*_{xz} above, is equal to the stress response, *σ*_{xz}, of the sheared mesophase to the applied strain.

### (b) Mapping of lattice units to physical units

We determine the mapping of simulation units to physical units by comparing the simulated diamond phase self-assembly with the lyotropic diamond phase self-assembly observed in experiments (Squires *et al*. 2005; Rappolt *et al*. 2006). Based on the experimental lattice parameter size of the order of 10 nm, the simulation resolution is approximately equal to 2 nm per lattice unit; this mapping is also in agreement with the previous simulations of the gyroid-phase assembly (González-Segredo & Coveney 2004). In experiments on diamond-phase self-assembly, for example, as induced by changing the concentration of water or on application of a thermodynamic jump, the self-assembly times are found to be of the order of minutes; comparing this with simulations, we estimate that each simulation time step corresponds to approximately 4–40 ms of physical time. Similarly, the mass in lattice units corresponds to approximately six surfactant molecules (of the order of 10^{−24} kg), based on observations of particle distribution functions around a single lattice site, which approximate the physical liquid crystal lattice site. According to the above mapping between lattice units and physical units, the applied strain rates in the Couette flow simulations are in the range of 0.01–0.1 s^{−1} in physical units with the frequency of applied oscillatory shear between 0.001 and 0.2 rad s^{−1}. All simulation results in §§3 and 4 are reported in lattice units.

The mapping of lattice to physical units using the observed behaviour of the system, as described here, is likely to lead to optimistic estimates of the length and time scale covered by the simulation. One might attempt a more precise mapping by comparing dimensionless groups such as the Reynolds and capillary numbers. However, for complex fluids such as those in the present work, full resolution of all the length, time and energy scales is not generally possible with the LB method and, at best, one expects to be able to correctly treat the hierarchy of scales in order to reproduce the physics of interest (Cates *et al*. 2005).

### (c) Isosurface visualization and topological characteristic

The interface between the oil and water phases is defined by the implicit surface of zero colour order parameter, *ϕ*(** r**)=

*ρ*

^{r}(

**)−**

*r**ρ*

^{b}(

**)=0, for**

*r***∈**

*r*^{3}(here, we use superscripts

*σ*=

*b*for water species and

*σ*=

*r*for the oil species). As described by Chin & Coveney (2006), the continuum-order parameter field on the discrete simulation lattice is estimated by trilinear interpolation. In order to graphically visualize the interface, a polygon mesh approximation is generated using the marching cubes algorithm implemented using our own isosurfacing software. The marching cubes algorithm detects changes in the sign of

*ϕ*(

**) corresponding to a transition from an oil-rich site to a water-rich one, indicating the point where the link between the sites intersects the isosurface. The marching cubes algorithm proceeds by considering the sign of**

*r**ϕ*(

**) at the vertices of each cube on the lattice to identify a set of polygons approximating the oil–water interface within the cube. The resulting polygon mesh representation of the isosurface of zero order parameter provides direct visualization of the interface between the oil and water phases. Given the polygonal mesh representation of a surface, we can calculate its topological Euler characteristic,**

*r**Χ*

^{surf}=

*V*−

*E*−

*F*, where the surface has

*V*vertices,

*E*edges and

*F*polygon faces. The Euler characteristic depends only on the topology of the surface and not on its shape or the particular polygon mesh representation (Kinsey 1993). Hence, the topological Euler characteristic of the surface or, equivalently, the genus

*g*, related to

*Χ*by

*Χ*=2(1−

*g*), is a useful metric to identify and compare cubic surfaces and has been used in this work. The Euler characteristic per unit cell is –8 for the gyroid surface, –16 for the diamond surface and –4 for the primitive surface (Schwarz & Gompper 1999). Our isosurfacing software runs on a single processor and has been used to analyse lattice sizes of up to 256

^{3}lattice sites.

The morphology of the self-assembled cubic phases is verified with a numerically exact, ideal structure generated using the periodic nodal approximation (PNA; Schwarz & Gompper 1999; Squires *et al*. 2005). The PNAs for the three simulated cubic surfaces are

gyroid surface is ,

diamond surface is , and

primitive surface is .

The topological characteristics of the simulated self-assembled mesophase can be compared with the numerically generated surface, using the corresponding PNA on a lattice of the same size, in order to compare and contrast the simulated morphologies.

### (d) Grid infrastructure

The simulations reported in this work were performed on high-performance transatlantic computing grid infrastructure. The computations were done on the petascale machine Ranger, at the Texas Advanced Computing Centre (http://www.tacc.utexas.edu/resources/hpcsystems/), with a theoretical peak of 0.58 petaflops and the Cray XT4 machine, HECToR, in the UK (http://www.hector.ac.uk), with a theoretical peak of 59 teraflops. The parallel kinetic LB code, LB3D, used in this work demonstrates excellent scaling on both resources, notably performing at a parallel efficiency of close to unity on up to 32 768 cores on Ranger (Saksena *et al*. in press). The simulation code was hosted in the application-hosting environment that enables distribution of simulations across grid resources (Coveney *et al*. 2007). Data transfer from the computational resources to local machines at UCL for post-processing and visualization was performed over the JANET-switched lightpath network (Coveney *et al*. in press; http://www.ja.net/services/lightpath/index.html), which is connected to the Starlight-switched router facility in the USA (http://www.startap.net/starlight/ABOUT/). The disk performance of computational resources is important in determining the overall performance of the LB3D code in our simulations; this is due to the large volume of data produced in this work (approximately equal to 3 TB). We found good scaling of I/O performance on up to 256 cores, the maximum used in this work. Weak-scaling benchmarks on Ranger, however, indicate that parallel scaling of our typical I/O requirements degrades when tens of thousands of cores are used; for example, performance decreases by 10 per cent, instead of remaining constant, on going from 1024 cores to 32 768 cores while keeping a constant load of 64^{3} lattice sites per core.

## 3. Simulation details

Here we describe LB simulations of the diamond phase in a ternary amphiphilic mixture on two different lattice sizes—64^{3} and 128^{3}—that are chosen to determine the scaling of the rheological response with system size as well as any finite-size effects. The overall computational cost of simulating the entire rheological flow curve increases with lattice size: the computational cost of each simulation scales as *L*^{3} where *L* is the lattice dimension; additionally, the number of simulation time steps required to reach steady state under shear also increases with lattice size (∼*L*^{2}).

In our LB simulations, the particle distribution functions in the initial random quench, *f*^{α}(0) (where *α* can be oil (*r*), water (*b*) or amphiphile (*s*)), are distributed uniformly as 0≤*f*^{α}(0)≤*F*_{α}, with the maximum density, *F*_{α}≤1.0, specified as input to the LB3D simulation code. The initial configuration is then allowed to evolve under the collision–advection scheme of the LB method in the presence of interaction forces (see equations (2.1)–(2.5)) for asymptotically long times. In simulations performed at specific values of the amphiphilic fluid composition and interaction parameters of the LB model, the ternary amphiphilic mixture undergoes self-assembly into the corresponding periodic mesophase, irrespective of the lattice sizes used in the simulation. Briefly, the interaction parameter *g*_{br} governs the interaction between the non-amphiphilic water (blue, *σ*=*b*) and oil (red, *σ*=*r*) species. A positive *g*_{br} is used to model repulsion between the two species (Chen *et al*. 2000; Love *et al*. 2003). A negative value of the water–amphiphile coupling parameter, *g*_{bs}<0, favours assembly of periodic structures in which the headgroups of the amphiphilic particles are oriented towards water-rich regions and the tail groups are oriented towards oil-rich regions. The amphiphilic coupling parameter *g*_{ss} controls the interaction between amphiphile particles. For *g*_{ss}<0, it describes attraction between two tail groups and repulsion between a head and a tail group, while, for *g*_{ss}>0, the coupling parameter models repulsion between two headgroups, as might arise in ionic surfactants (Bhattacharya & Mahanti 2001). For *g*_{ss}>0, we find that, despite the repulsion between headgroups, the amphiphile particles still align themselves at the oil–water domain interfaces with headgroups pointing towards the water domains; this arises due to their directional affinity to oil and water species and the repulsive interaction (*g*_{br}>0) between the oil and water species.

Previously, we have performed systematic parameter sweeps to extend the mesophase diagram explored by the kinetic LB simulation method (Saksena & Coveney 2008). We reported a wide range of mesomorphic behaviours obtained by systematically varying the interaction parameters, *g*_{br} and *g*_{ss}, at a fixed species concentration. In these kinetic LB simulations, we were able to observe mesomorphic variations in agreement with experimental studies of lyotropic solutions as (i) the amphiphilic strength *g*_{ss} was varied and (ii) the oil–water repulsive interaction *g*_{br} was varied.

### (a) Cubic-phase self-assembly

For the same concentration of oil, water and amphiphilic species and different values of *g*_{br}, *g*_{bs} and *g*_{ss}, we observe the self-assembly of the primitive, gyroid and diamond mesophases (Saksena & Coveney 2008). The parameters corresponding to the three phases are listed in table 1. The morphology of the self-assembled cubic domains is verified with an ideal structure generated on a lattice of the same size and with the same number of unit cells using the PNA (Schwarz & Gompper 1999; Squires *et al*. 2005) by means of (i) visual comparison and (ii) comparison of the topological Euler characteristic (Gandy *et al*. 1999; Chin & Coveney 2006) with the PNA domain (Squires *et al*. 2005; Saksena & Coveney 2008). In figure 1, we visualize the three-dimensional isosurfaces of the simulated self-assembled gyroid, diamond and primitive domains and the corresponding PNAs. We have selected sub-domains of the self-assembled cubic mesophases which are closest to the ideal PNA sub-domain. During the simulated self-assembly process, there are small changes in local alignment of the cubic phase with respect to the simulation axes and spontaneous formation of defects, which result in differences between the simulated system morphology and the corresponding ideal PNA structures that are idealized static representations, perfectly aligned with the *x*, *y* and *z* simulation axes.

In the rest of this paper, we will focus on the diamond phase that undergoes self-assembly in our LB simulations at parameters specified in table 1 (system II). The salient morphologies that are encountered as the system self-assembles into the diamond phase on a 128^{3} lattice-site system, from a random quench, are shown in figure 2. Systems of both sizes took between 13 000 and 16 000 simulation time steps to self-assemble into the diamond phase from the initial random quench. For the 64^{3} lattice, we observe a self-assembled ideal diamond phase with an Euler characteristic value equal to that of the corresponding PNA surface, while, on the 128^{3} lattice, the self-assembled diamond phase contains defects and its Euler characteristic differs from that of the corresponding PNA surface.

After the self-assembly of the diamond phase, we follow the diamond phase dynamics on the 64^{3} and 128^{3} lattices over the next 50 000 time steps and find that the system undergoes transitions to a short-lived, intermediate state with primitive domains; this reconverts to a purely diamond morphology within 1500 time steps. Finite-size effects play a major role in the defect dynamics and the subsequent phase transition to primitive domains. In particular, in the simulations on smaller 64^{3} lattices, the primitive sub-domains grow to occupy the entire simulation cell; that is, we see a spontaneous ideal diamond→ideal primitive transition prior to reconversion to the diamond phase. By contrast, in simulations on the larger 128^{3} lattice, transient primitive-phase sub-domains nucleate around defects and coexist with the diamond and sponge phases before annealing to a lattice with only diamond domains. In figure 3, we plot the dynamical evolution of the Euler characteristic *Χ* for the 64^{3} and 128^{3} systems; the ideal line shows the Euler characteristic value of a perfect diamond mesophase with 10^{3} unit cells on a 64^{3} lattice-site system and 20^{3} unit cells on a 128^{3} lattice-site system. The troughs in the plot of the Euler characteristic correspond to the short-lived intermediate state comprising primitive sub-domains that form and disappear as the simulation progresses. These transitions are hydrodynamically driven (as is confirmed by direct correlations between the fluid velocity field and these non-equilibrium phase transitions; data not shown). Asymptotically, the system is eventually expected to reach an equilibrium state, where such transitions are absent (although it may be noted that interacting LB models neither satisfy detailed balance nor, in turn, support an H-function; Boghosian *et al*. 2001). However, experimentally, it is known that D and P cubic phases often coexist and thus one may expect the barrier to their interconversion to be small (Czeslik *et al*. 1995; Winter *et al*. 1998; Conn *et al*. 2006).

In figure 4, we show isosurface plots of the ternary system in and around the diamond→(primitive+diamond+sponge)→diamond transition for the 128^{3} lattice-site system. This sequence of isosurfaces illustrates the morphological changes that occur as the diamond phase converts to the intermediate primitive phase. This particular morphological sequence corresponds to the sixth peak in figure 3*b*. We note that although the genus and Euler characteristic are important metrics for classification and verification of topologies, they do not furnish a sufficient condition to guarantee the existence of a particular cubic phase. However, for systems that are already known to be in a particular cubic morphology, for example the diamond phase in the present case, the topological parameters provide a succinct metric by means of which to follow the morphological dynamics of the periodic mesophase, highlighting points in the morphological landscape which are worthy of more detailed investigation.

The relatively high temporal and spatial resolution of the dynamics of the D→P transition that we obtain from our simulations provide an opportunity to analyse the phase transition mechanism. It is notable that these kinetic LB simulations allow us, for the first time, to study cubic→cubic transitions at a high resolution. Preliminary investigation into the origin of the formation of primitive sub-domains indicates that growth of primitive sub-domains is nucleated in the vicinity of defects in the diamond phase.

We turn next to consider the non-equilibrium properties of the diamond phase under an applied Couette flow.

## 4. Rheology

In this section we report the rheological response of the diamond phase to applied Couette flow. In the LB simulations of Couette flow, a linear shear velocity profile is applied using Lees–Edwards boundary conditions as described in §2. Our simulations are able to correctly simulate the shear-banding effect and non-Newtonian responses of liquid crystals under steady shear as observed in experiments (Bonn *et al*. 1998; Pitzalis *et al*. 2000; Rey & Denn 2002; Smith *et al*. 2002; Porcar *et al*. 2003; Mezzenga *et al*. 2005*a*). The simulated viscoelastic response under an oscillatory shear is also in very good agreement with analogous experiments on the diamond phase (Pitzalis *et al*. 2000; Mezzenga *et al*. 2005*a*).

### (a) Approach to steady-state Couette flow and shear banding

The number of simulation time steps required for steady-state Couette flow to be established depends on the lattice size, with a shear velocity profile occurring in approximately 1000 time steps for the 64^{3} lattice and within 4000 time steps for the 128^{3} lattice. The evolution of the velocity profile from the unsheared system to steady state for the two system sizes is shown in figure 5. The nonlinear velocity profile that is established here at steady state is evidence of a backflow effect (Jewell & Sambles 2002; Smith *et al*. 2002; Marenduzzo *et al*. 2003): the movement of the fluid components themselves under the influence of the applied shear velocity as well as local surface reordering causes a modification of the linear velocity profile under Couette flow. This is a mesoscopic phenomenon involving an interplay between the externally applied macroscopic velocity flow field and microscopic dynamical processes governed by the interspecies interactions. The nonlinear velocity profile is accompanied by the appearance of distinct liquid crystal morphology (texture) in regions along the shear velocity gradient within the simulation cell. Such shear banding has been observed in our rheological simulations of the diamond phase for both system sizes. This non-Newtonian effect has also been observed in experiments and other simulations of liquid crystals under shear (Rey & Denn 2002; Marenduzzo *et al*. 2003).

In figure 6, we show the shear-banded texture in the *xz* (shear)-plane and the *xy* (transverse)-plane for the 64^{3} and 128^{3} lattice-site system, at a shear velocity *U*=0.11 in lattice units. The isosurface visualization of the interfacial structure in the 64^{3} and 128^{3} lattices is shown here.

In contrast to the morphology at the centre of the simulation cell, the shear-banded texture in the high shear velocity regions comprises a smectic lamellar phase, oriented in the *transverse* direction relative to the applied shear velocity (figure 6). Such an orientation relative to the applied velocity field is particularly striking as one would expect it to be unstable and rather a parallel or perpendicular lamellar phase to exist at a steady state. However, intermediate smectic domains aligned transverse to the shear flow have indeed been observed in experimental studies of di-block copolymers, during the transition from an isotropic mixture to the parallel/perpendicular lamellar phase under shear (Chen *et al*. 1997). Similarly, in our simulations, transverse lamellar domains are formed when the diamond channels are subjected to shear flow in the *z*-direction as the oil and water channels get reoriented under the flow field. We conjecture that this transverse orientation is reinforced due to backflow effects and hence it persists, an effect also observed in twisted nematic devices where a short-lived rotation of the director field occurs in the ‘wrong’ direction upon switching (Smith *et al*. 2002; Marenduzzo *et al*. 2003). There are also examples, in numerical studies, of the development of transverse orientational patterns in liquid-crystalline polymer systems under steady shear where the interactions are modelled by the Marrucci–Greco potential or by the simpler Leslie–Ericksen theory (Rey & Denn 2002). Backflow effects will, of course, be significant only when the applied rate of external shear is low enough to be able to couple to intrinsic fluid dynamical processes such as surface reordering; we would expect that, at high shear rates, shear flow will become the dominant factor and break the ordering of lamellae transverse to it. This is indeed what is observed in our simulations beyond a critical value of the strain rate: at such strain rates greater than the critical value, the transverse smectic phase is ruptured at steady state as we will describe in §4*b*. Although many experimental reports exist in the literature on the non-Newtonian rheology of amphiphilic cubic liquid crystals (Pitzalis *et al*. 2000; Porcar *et al*. 2003; Mezzenga *et al*. 2005*a*; Pouzot *et al*. 2007), there are very few papers reporting morphological transitions *induced* by shear (for example, Porcar *et al*. (2003) discuss the shear-induced transition from the sponge to the lamellar phase in an amphiphilic solution), although shear-induced morphological transitions in self-assembled block copolymer systems have been studied in more detail (Fredrickson & Bates 1996; Park & Char 2005; Darling 2007). Experimental verification of the formation of such a transverse smectic phase at low strain rates when a self-assembled diamond morphology is subjected to steady shear is thus not currently available; we hope our work will inspire further laboratory studies.

### (b) Stress response

Once steady shear had been established, we averaged the stress response over 5000 and 3000 time steps for the 64^{3} and 128^{3} lattice-site systems, respectively. The average stress-response versus strain-rate profiles and viscosity versus strain rate for the two system sizes are plotted in figures 7 and 8.

For both lattice sizes, we performed multiple runs to generate the flow curve starting from different initial diamond phase configurations, which have been aged for different lengths of time. At a shear velocity close to *U*≈0.13, we observe dramatic shear thinning behaviour, i.e. the viscosity of the system reduces sharply as the strain rate is increased. However, at this strain rate, the value of the stress response itself exhibits variations between systems sheared from different initial configurations in the diamond morphology, indicative of hysteresis.

Shear thinning associated with a stress plateau in the flow curve has been widely observed experimentally in a range of self-assembled systems and morphologies such as lyotropic lamellar, multi-lamellar and hexagonal phases, thermotropic liquid crystals and block copolymers. We are unaware of analogous experiments on bicontinuous cubic liquid crystals. Here, we draw comparison with phenomenologically related experiments on other lyotropic self-assembled systems. In experiments involving Couette flow, flow curves have been reported, which are in qualitative agreement with our simulations over the entire range of strain rates simulated. As examples, we mention here experiments comprising fixed strain-rate rheological measurements of (a) the micellar bcc cubic crystal in aqueous triblock copolymers (Eiser *et al*. 2000) and the lyotropic liquid-crystalline lamellar ‘onion’ phase (Bonn *et al*. 1998), in which flow curves with a stress plateau together with a sharp drop in viscosity with increasing strain rate have been observed. The existence of the stress plateaux in these experiments has been attributed to the coexistence of high- and low-viscosity phases, which are stable on the lower and higher shear rate branches, respectively. The occurrence of hysteresis in the experiments in the vicinity of the stress-response drop has been attributed to the fact that the system was bistable, with a small energy barrier between the two states, thus resulting in two values of the viscosity for a given strain rate (Bonn *et al*. 1998). The hysteresis in our simulations would also arise if the system were indeed bistable in the vicinity of this critical value of the strain rate corresponding to a morphological transition in the sheared system. In our diamond phase rheology simulations, the bistable states are (i) the periodic (smectic+diamond) morphology found at lower strain rates and (ii) the disordered morphology, with oil and water channels running parallel to the shear flow, which is found at higher strain rates as described below.

At the transition strain rate (figure 7), the diamond phase undergoes a morphological transition, in which the oil and water channels aligned transverse to the shear flow are disrupted and structure at steady state under shear is no longer composed of the periodic (transverse smectic+diamond) phase. Figure 9 shows the steady-state morphology at strain rates beyond the morphological transition; this can be compared with the steady-state morphology at a lower strain rate corresponding to *U*=0.11 in figure 6. For the 64^{3} lattice systems and to a lesser extent for the 128^{3} lattice systems, the disruption in periodicity extends into the centre of the simulation cell.

In the steady-shear experiments of the lamellar ‘onion’ phase (Bonn *et al*. 1998), the disappearance of the shear-banded morphology was detected at high strain rates by phase contrast microscopy and was attributed to either (i) the actual disappearance of the shear bands or (ii) inability to detect shear bands as they became too small. Our results suggest the former to be the more likely explanation for the disappearance of the shear bands at high shear rates in these experiments (Bonn *et al*. 1998). We do not observe any strain-rate-dependent variations in the size of the shear bands in simulations of diamond phase rheology.

To the best of our knowledge, there are no experimental results in the literature on the amphiphilic diamond mesophase under steady shear which exhibit the stress plateau and with which we could make a quantitative comparison. However, in the phenomenologically similar steady-strain rheological experiments of the biphasic onion phase (Bonn *et al*. 1998), analogous extreme shear-thinning behaviour was observed at , while in the steady-strain rheological experiments investigating the conversion of the *L*_{3} sponge to the lamellar phase, the transition occurred at with a decrease in viscosity (Porcar *et al*. 2003); it was also noted that the critical shear rate could be reduced by a factor of 8 by halving the lipid volume fraction.

To summarize, for the two system sizes simulated, we find agreement in the profile of the flow curve comprising a non-Newtonian shear-thinning regime that includes a striking drop in viscosity. For the 64^{3} lattice system, we observe subtle finite-size effects on closely following the shear-induced morphological transition dynamics from the initial unsheared configuration. In contrast to the gradual morphological evolution under shear that is observed in the larger systems, for the 64^{3} lattice-site system, the initial diamond morphology completely breaks up into a random quench as soon as shear is applied; the shear-induced assembly of the transverse smectic+diamond phase is observed subsequently. This indicates that the long-range correlations in the initial diamond liquid crystal present in the smaller 64^{3} lattice-site system are coupled to the periodic boundary conditions: a slight perturbation in the boundary condition causes the liquid-crystalline morphology to break up.

Finally, we find that the steady-state flow curves for the two system sizes have the same shape but differ quantitatively. For the 128^{3} lattice-site system, the sharp drop in viscosity occurs at a strain rate that is close to half of that for the 64^{3} lattice-site system. There is no theoretical basis for such a scaling of the critical strain rate with system size; however, it is notable that non-overlap of Couette flow curves measured for different Couette cell gap widths and non-trivial correlation of Couette cell gap width with wavelength of vorticity structures have been observed in experiments on amphiphilic fluids (Ganapathy *et al*. 2008).

### (c) Relaxation of the stress response

We also analysed the time dependence of the relaxation of the stress response once the shear flow field is extinguished. In previous LB and lattice-gas simulations, domain growth kinetics during self-assembly of amphiphilic ternary mixtures followed a stretched-exponential functional form (Emerton *et al*. 1997; Nekovee *et al*. 2000). The viscoelastic response of the gyroid phase was also found to be consistent with stretched-exponential kinetics (Giupponi *et al*. 2006). Following earlier work, we fit the relaxation kinetics to the following form:(4.1)We used Matlab (http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml) to perform a least-squares fit of the relaxation curves of the 128^{3} diamond phase lattice sheared at various initial strain rates to the stretched exponential. The particular sheared system we study here corresponds to the system whose flow curve is plotted in figure 8 with legend *t*_{0}=43 000. The relaxation curve and corresponding least-squares fit for two initial strain rates are shown in figure 10, while the inset shows the relaxation curves for a range of initial strain rates. The fitted parameters *A*, *B*, *C* and *D* for the strain rates corresponding to different shear velocities, *U*, are listed in table 2. The value of the stretched-exponential exponent *D*, and hence relaxation times, increases with initial applied strain rate and remains close to 1 for the highest strain rates considered here, corresponding to pure exponential relaxation dynamics. The time dependence of the stress relaxation after the cessation of shear flow gives a good indication of the time scales involved in the intrinsic relaxation processes of the self-assembled mesophase. The internal stress response of the system approaches zero and, accompanying this, re-adopts the diamond morphology within 5000 time steps. This is similar to experimental observations of a sheared gyroid phase where the cubic phase recrystallizes after the cessation of shear albeit with a different orientation (Olsson & Mortensen 1995).

### (d) Oscillatory shear and viscoelasticity

Oscillatory shear is simulated by modulating the applied shear velocity according to the equation(4.2)where *ω*/2*π* is the frequency of oscillation and shear velocity, *U*, is set to 0.1. Here we follow the approach taken by Giupponi *et al*. (2006) in which such an oscillating shear was applied to the gyroid mesophase to measure its viscoelastic response. The simulations were run for at least three complete oscillations over a range of shear frequencies, *ω*. We calculated the storage and loss moduli by fitting the stress-response profile to the equation(4.3)where and *G*′(*ω*), *G*″(*ω*) are the storage and loss moduli, respectively. The fitted values for *G*′ and *G*″ versus *ω* are plotted in figure 11.

At low frequencies, the loss modulus dominates, indicating a viscous response, while at high frequencies the system exhibits a solid-like elastic response. The frequency of oscillation at which the storage and loss moduli curves intersect signals the transition from liquid-like behaviour to solid-like behaviour; it gives an estimate of the characteristic stress relaxation times in the system and provides clear evidence of the viscoelastic behaviour of the simulated liquid-crystalline diamond phase. The form of the dynamic moduli versus strain oscillation frequency is in qualitative agreement with experimental studies of the shear rheology of the diamond phase formed in lyotropic liquid crystals (Pitzalis *et al*. 2000; Mezzenga *et al*. 2005*a*). In figure 11, the crossover between liquid-like and solid-like behaviours occurs at *ω*≈0.015, which corresponds to a relaxation time *τ*≈419 in lattice units and 17 s in physical units according to the mapping in §2*b*. The observed relaxation time is also in good agreement with that obtained for the diamond phase in the experiments of Mezzenga *et al*. (2005*a*). It should be noted that, in the experiments of Pitzalis *et al*. (2000), the diamond phase exhibits a dominant storage modulus and is essentially elastic except for the lowest shear frequencies. We plot the stress-response curve versus simulation time in figure 12 at *ω*=0.015 where the diamond phase exhibits a viscoelastic response with a phase lag *ϕ* between 0 (purely elastic) and *π*/2 (purely viscous).

## 5. Conclusions

In this paper, we have reported the first simulations of a ternary amphiphilic diamond mesophase characterizing the rheological response under an applied Couette flow. Under shear, the diamond phase exhibits many of the characteristic phenomena of sheared liquid-crystalline systems observed in experiments (Bonn *et al*. 1998; Pitzalis *et al*. 2000; Rey & Denn 2002; Smith *et al*. 2002; Mezzenga *et al*. 2005*a*), including a non-Newtonian flow curve and shear banding. Non-Newtonian behaviour arises in complex amphiphilic fluids because of the coupling between the length and time scales of molecular ordering, determined by microscopic interactions and the applied macroscopic flow field. It is notable that such rich phenomenology can be accessed using a minimal amphiphilic model simulated on a three-dimensional lattice with a discrete set of velocity vectors using the LB method. At a critical value of the strain rate, we observe extreme shear thinning corresponding to a steep drop in viscosity accompanied by a morphological transition in which the ordered (smectic+diamond) phase is broken up under shear and channels are formed oriented in the direction of shear. We observe hysteresis in the rheological response of the amphiphilic fluid at the transition strain rate and beyond, particularly in the 64^{3} lattice simulations. Finally, we note that the flow curves for the 128^{3} and 64^{3} systems exhibit the same shape but do not overlap. The rheological responses of the two system sizes differ in part due to the presence of multi-domain defects in the larger 128^{3} lattice-site systems.

Finite-size effects are found to be pronounced on following shear-induced morphological transitions in the smaller (64^{3}) systems where, under shear, the initial diamond morphology exhibits a striking break-up into a random quench before assembling into the steady-state shear-induced morphology. For larger lattices, the transition to steady-state morphology under shear occurs gradually. This highlights the importance of using a sufficiently large model in simulations of dynamical processes involving local interactions and in which long-range correlations arise due to hydrodynamics, inherent periodicity of the system or are induced by an externally applied field. However, despite the difference in the dynamical approach to steady state, we obtain the same shear-induced steady-state morphology for both smaller and larger system sizes.

On comparing various functional forms (Emerton *et al*. 1997; Nekovee *et al*. 2000) for the relaxation of the stress response on the cessation of shear, as the diamond phase reforms, the best fit is obtained by a stretched-exponential function for systems initially sheared at low strain rates. The stretched-exponential form has previously been found to describe domain growth kinetics in simulations of ternary amphiphilic fluids and in experimentally observed micellar self-assembly (Wilcoxon *et al*. 1995). Stretched-exponential kinetics are characteristic of glassy systems with multiple relaxation times. In our simulations, at high strain rates, as channels develop along the flow direction there is a transition from a high-viscosity glassy steady state to a low-viscosity steady state, for which the stretched exponential gives a poorer fit. The stress-response relaxation provides insight into the probable functional form of the kinetics of intrinsic dynamical processes such as the growth of defects and rearrangement of the domains in the unsheared system.

## Acknowledgments

We are grateful to EPSRC for funding this research through the grant EP/E045111/1, which also provided access to HECToR. This work was also supported by US National Science Foundation MRAC grants TG-DMR070013N and TG-DMR070014N. We are grateful for DECI allocations on the DEISA grid (http://www.deisa.org) provided by the EU (FP6 project 508830 and FP7 project 222919). We would like to thank Davide Marenduzzo for useful discussions.

## Footnotes

- Received November 21, 2008.
- Accepted March 4, 2009.

- © 2009 The Royal Society