## Abstract

We investigate a model of hard pear-shaped particles which forms the bicontinuous Iad structure by entropic self-assembly, extending the previous observations of Barmes *et al.* (2003 *Phys. Rev. E* **68**, 021708. (doi:10.1103/PhysRevE.68.021708)) and Ellison *et al.* (2006 *Phys. Rev. Lett.* **97**, 237801. (doi:10.1103/PhysRevLett.97.237801)). We specifically provide the complete phase diagram of this system, with global density and particle shape as the two variable parameters, incorporating the gyroid phase as well as disordered isotropic, smectic and nematic phases. The phase diagram is obtained by two methods, one being a compression–decompression study and the other being a continuous change of the particle shape parameter at constant density. Additionally, we probe the mechanism by which interdigitating sheets of pears in these systems create surfaces with negative Gauss curvature, which is needed to form the gyroid minimal surface. This is achieved by the use of Voronoi tessellation, whereby both the shape and volume of Voronoi cells can be assessed in regard to the local Gauss curvature of the gyroid minimal surface. Through this, we show that the mechanisms prevalent in this entropy-driven system differ from those found in systems which form gyroid structures in nature (lipid bilayers) and from synthesized materials (di-block copolymers) and where the formation of the gyroid is enthalpically driven. We further argue that the gyroid phase formed in these systems is a realization of a modulated splay-bend phase in which the conventional nematic has been predicted to be destabilized at the mesoscale due to molecular-scale coupling of polar and orientational degrees of freedom.

## 1. Introduction

This article addresses the spontaneous formation of one of nature's most symmetric, most complex and most ordered structures—the self-assembled bicontinuous gyroid structure. Further, it does so in the context of a system which involves no attractive forces. Therefore, the formation of this gyroid structure does not happen despite entropy, but because of it: the long-range order emerges precisely because this highly ordered structure maximizes entropy. While this runs counter to the commonly held, but flawed, notion that maximal entropy is always associated with increasing levels of disorder [1–3], in fact, several examples of entropy-driven ordering have been previously identified. Foremost among these are the celebrated ordering transition in the hard sphere equilibrium fluid [4,5], beautifully realized in colloidal crystals by Pusey and van Megen [6] and the longstanding debate about the corresponding hard disk system [7,8], as well as the emergence of nematic or smectic order in simple liquid crystal models like spherocylinders [9–11], ellipsoids [12], helices [13] and other multi-sphere objects [14,15]. See also the recent review by Manoharan [16].

This paper analyses the formation of the double gyroid in equilibrium ensembles of hard pear-shaped particles, thereby extending earlier work by Barmes *et al.* [17] and Ellison *et al.* [18]. Pear-shaped particles are tapered versions of spheroids, best thought of as prolate ellipsoids with a wider ‘blunt’ end and a narrower ‘sharp’ end (see figure 1). For appropriate parameter values, equilibrium ensembles of such pear-shaped particles adopt the so-called double gyroid structure [19–21]. The double gyroid is a bicontinuous geometry that is best visualized as two intertwined, identical, periodic and highly symmetric network-like or labyrinth-like domains (see [22] for animation). The common interface between these domains is a triply periodic, saddle-shaped surface, with symmetry group including operations that exchange the two labyrinth-like domains. The pears adopt this diffusive phase in an arrangement that fills space fairly uniformly, at fluid-like densities. The blunt ends of the pears point into the labyrinthine domains, with pear positions adopting a distribution of distances from the (hypothetical) minimal surface and pear orientations adopting characteristic angles around the minimal surface normal directions. The pear blunt ends, hence, can straightforwardly be subdivided into two subsets, each occupying one or the other of the labyrinthine domains, whereas the locations of the sharp ends are staggered near the minimal surface (see figure 2*b*).

In various forms, the double gyroid is a commonly observed phase in a variety of self-assembled soft matter systems [22–24]. In particular, systems of amphiphilic molecules—with energetic terms that favour local segregation—adopt the double gyroid: in the inverse lipid/water phase, a warped bilayer membrane separates two aqueous domains [25–27] (see figure 2*c*). Also, dry monovalent soaps form *Q*^{230}_{I} double gyroid phases when their hydrophobic moieties form the network-like domains in the absence of water [28]. The second important classification of molecules that adopt double gyroids are copolymeric melts [29], in particular, diblock copolymers. These form the double gyroid phase at volume fractions of approximately 33% of the minority component [30,31] (see figure 2*d*). The double gyroid has also recently been observed in experimental investigation of dendrimer systems [32] and simulations using sphero-symmetric potentials with short-range attractions and longer-range repulsions [33–35].

The importance of entropy for the self-assembly of amphiphilic systems has long been recognized and is implicit in both the molecular shape concept [23,36] and the Helfrich formalism [37,38]. Indeed, soft matter physics is generally concerned with systems in which entropy plays a significant role, that is, where interaction terms are typically of the same order as the thermal energy *k*_{B}*T*. This is certainly the regime in which the bicontinuous double gyroid phase is formed by amphiphilic molecules. However, all of the systems mentioned in the previous paragraph also have a clear enthalpic component, evident in the amphiphilic (segregating) nature of their constituent molecules. This is a significant difference from the hard pear-shaped particles studied here. One goal of this article is, therefore, to contribute to a deeper understanding of the mechanisms that lead to the formation of double gyroid phases. Are the mechanisms in amphiphilic systems fundamentally different from those exhibited by pear-shaped particles? Is the mechanism by which the gyroid phase increases the entropy in the pear-shaped system the same as that by which this happens in the amphiphilic systems? This article intends to provide data that will contribute to the resolution of these questions.

We also note that there is one particular way in which the behaviour of these pear-shaped systems is distinguished from those of other entropy-driven ordering processes such as those of hard spheres and liquid crystals. The onset of the double gyroid phase creates crystallographic order on a larger length-scale (‘mesoscale’) than the particle or molecular scale, forming unit cells whose dimension is one or more orders of magnitude larger than the particle or molecule size. This contrasts with the ordered phases of hard spheres and liquid crystals where the crystallographic parameter corresponds to the translation from one particle to its nearest neighbour (and hence to the molecular scale itself). While the formation of emergent mesoscale order by self-assembly is by now a well-established concept, a clear discrimination between the energetic and entropic contributions is an only partially understood question and one to which we aim to contribute here. Ultimately, understanding how entropy can be used to create ordered structures on the technologically important mesoscale (typically nanometers) is an essential step towards low-energy sustainable nano-engineering. ‘Only entropy comes easy’ [Anton Chekhov, *ca* 1900].

This paper is structured as follows. In §2, we present the main simulation methodology and results that underpin this work and, thus, present a phase diagram showing the density–tapering parameter phase diagram of hard pears, including the extent of its gyroid phase region. Following this, in §3, detailed structural analysis is undertaken, through which insights are gained into the nature of the hard pear gyroid and how it differs from those exhibited by different classes of experimental system. Finally, in §4, the entropy-driven gyroid phase formed by the hard pear system is appraised in the round. Through this, a compelling argument is made that it is a realization of a modulated splay-bend phase in which the conventional nematic has been predicted to be destabilized at the mesoscale due to molecular-scale coupling of polar and orientational degrees of freedom.

## 2. Phase diagram

In this section, we extend previous reports on the simulation of systems of hard pear-shaped particles. Barmes *et al.* showed that for long, tapered particles with aspect ratio *k* = 5 and *k*_{θ} = 5 the nematic, bilayer smectic and crystalline phases are formed [17], see figure 1 for the definitions of *k* and *k*_{θ}. For shorter particles, however, this initial study only found a ‘domain ordered’ arrangement on compression of systems of 1000 tapered particles. Subsequently, by simulating 10 000 particles with *k* = 3 and *k*_{θ} = 3.8, Ellison *et al.* showed that, on compression from the isotropic fluid, the system was actually entering a gyroid arrangement with long-range, three-dimensional periodicity [18]. The same study also showed spontaneous development of a gyroid structure on expansion of a high density, interdigitated bilayer smectic, but the reverse transition was not observed.

We here use molecular dynamics (MD) and Monte Carlo (MC) simulation techniques to undertake a systematic simulation study of the ordering behaviour of *k* = 3 hard pears for a range of tapering angles. The interaction potential used to represent the hard core interaction for these pears is the modified version of the purely repulsive Weeks–Chandler–Andersen potential (WCA) [39], which Barmes *et al.* introduced. This is an implementation of the parametric hard Gaussian overlap (PHGO) approximation which is, in turn, based on the generalized Gay–Berne interaction [40,41]. The PHGO approach is founded on the observation that contacts between convex particles can locally be approximated by those between appropriately chosen ellipsoids. To maintain adherence to this, particle convexity is achieved here by using a core description defined by two continuously differentiable Bézier curves [17]. Each pear–pear interaction can then be described via an expansion based on multiple hard core ellipsoid potentials, weighted to mimic the Bézier descriptions at arbitrary relative orientation. Full details of the PHGO approach and its application to hard pear systems are given in [17,42].

For this part of the study, simulations are performed on systems of *N* = 3040 particles with *k*_{θ} between 2 and 6 (tapering angle between 28.1° and 9.5°) within a cubic box with three-dimensional periodic boundary conditions. The lower boundary of *k*_{θ} = 2 is set to ensure that there is no particle concavity. Full simulation sets are performed using MD in the canonical (NVT) ensemble, with a time step Δ*t* = 0.0015 and the dimensionless temperature *T* set to 1, since phase behaviour at fixed density is independent of the temperature for hard-particle systems [43]. Additionally, MC simulation sets using the same parameters are run for *k*_{θ} = 2.2, *k*_{θ} = 3.8 and *k*_{θ} = 5.4 and for other systems as outlined in subsequent sections. The MC translation step and the rotation step are set as Δ*s* = 0.05 and Δ*ϕ* = 2.0°, respectively. The results of these two simulation sets show no significant differences. The global density of the system
2.1which corresponds to the packing fraction, is defined by the pear volume *V*_{pear}, calculated numerically using a mesh of the particle's surface, and the volume of the simulation box *V*_{box}. Here, and subsequently, lengths are expressed in the units of the pear width at its waist *w*_{p}.

For each value of *k*_{θ} an initial, crystalline ordered configuration generated at low density (*ρ*_{g} = 0.28) is run to erase all configurational memory before being compressed to a starting density, *ρ*_{g} = 0.44. Subsequently, a sequence of small compression steps are imposed (see symbols in figure 3); each of which entails an equilibration run of 1.5 × 10^{6}Δ*t* and a production run of 10^{6}Δ*t*. Compressions are made up to *ρ*_{g} = 0.67, which is found to be a crystalline state for all *k*_{θ}. Expansion sequences are performed in an equivalent, but reverse, manner from each *ρ*_{g} = 0.67 state. The resultant phase diagram is shown in figure 3. To examine the sensitivity of the phase diagram simulations with different system sizes from *N* ≈ 1750 to *N* ≈ 10 000 are performed as well. These show only modest changes in the phase diagram in regard to system size and box shape. We, therefore, conclude that the very small crystallographic moduli at play here mean that commensurability effects, while obviously present, are surprisingly weak.

Phase identification is based on four main observables. The first of these was the excess pressure
2.2which is based on the virial and dependent on the distance *r*_{ij} of the centre positions and the forces *f*_{ij} between particles *i* and *j*. All forces *f*_{ij} are repulsive, since the PHGO model is a soft-repulsive particle model [17,42]. The excess pressure, however, has to be treated with caution, as we deal with a potential close to the hard core limit (maximally 2% overlap), in which *p*_{ex} is constant. The second key observable is the nematic order parameter
2.3where *u*_{i} denotes the orientation vector of the *i*th particle and the system director, *n*, is the eigenvector associated with the largest eigenvalue of the tensor
2.4The standard deviation of local orientations
2.5based on the scalar product *α*_{i} = *u*_{i} · *u*_{j} of the orientation vectors of nearest-neighbour particles *i* and *j*, and the diffusion coefficient *D* are the other main observables. *D* is determined from the diffusive-regime slope of the mean squared displacement characteristics such as those shown in figure 4. These observables also serve to confirm system equilibration through their constancy under consideration of fluctuations within production runs (see figure 5). Additionally, systems are analysed using cluster identification algorithms applied to the pear blunt ends, performed to enable structural characterization of phases. From this, six distinct phases are identified—isotropic, nematic, smectic, gyroid, solid smectic and solid gyroid—as well as narrow biphasic or hysteretic regions (marked in grey in figure 3) between isotropic and ordered fluid phases. Since the resultant phase diagram, figure 3, can readily be divided into three sections with regard to the particle tapering parameter *k*_{θ}, details of observable characterization are now given in three separate paragraphs §§2.1–2.3 corresponding to these three sections in the phase diagram.

### 2.1. Strong particle tapering: 2.0 < *k*_{θ} < 2.4

For small tapering parameter values, between *k*_{θ} = 2.0 and *k*_{θ} = 2.4, systems have low orientational order at low densities, *ρ*_{g} = 0.49 (*k*_{θ} = 2.0) and *ρ*_{g} = 0.52 (*k*_{θ} = 2.4), and cluster algorithms cannot identify any full bilayer structures. However, randomly oriented bilayer-like clusters of interdigitating particles, such as those apparent in figure 6 (centre), are observed as the ordered phase is approached from low density.

On compression from the isotropic, the system exhibits smectic lamellar and solid-smectic phase behaviours. While the excess pressure is an effective indicator of the transition from isotropic to smectic lamellar (figure 5), the main signal of this transition is the adoption of high-orientational order parameter (see figure 5*b*). As depicted in figure 6*a*, in the smectic phase, layers of interdigitating bilayer sheets or leaflets are formed in which all particles are orientationally aligned either parallel or antiparallel to one another. At intermediate densities, mobility remains diffusive, particle motion being dominated by in-leaflet diffusion but also involving occasional ‘flip-flop’ of pear-shaped particles between neighbouring sheets. A second transition, between smectic and solid-smectic, is characterized by a steep drop in mobility (see figure 5*d*) as well as features in the excess pressure and order parameter characteristics.

### 2.2. Intermediate particle tapering: 2.4 < *k*_{θ} < 4.5

Pears with intermediate tapering, 2.4 < *k*_{θ} < 4.5, form an isotropic fluid for small values of *ρ*_{g}. However, as the density is increased for these systems, the gyroid arrangement is adopted rather than the bilayer smectic phase. The isotropic–gyroid transition occurs at densities between *ρ*_{g} = 0.52 (*k*_{θ} = 2.4) and *ρ*_{g} = 0.55 (*k*_{θ} = 4.5), the systems adopting a super-structure characterized by two interpenetrating channel networks (see figure 7). The interdigitating curved bilayers indicate a smectic order on the gyroid minimal surface (note related work on possible nematic ordering in block copolymers [45]). In this phase, the pear-shaped particles are able to traverse the simulation box isotropically, through both in-leaflet diffusion and leaflet to leaflet flip-flop. In addition to being bicontinuous and triply periodic, clusters identified from the positions of the blunt ends of the pears display nodes with three branch junctions, as is characteristic for the double gyroid. These clusters are, indeed, readily identifiable with the medial surface of the gyroid. This medial surface (or axis) is a geometric construction that produces a centred skeleton of the original shape (see [46] for a recent review). For the case of bicontinuous structures, it represents a generalized line graph that also provides a robust definition of local domain (or channel) size and hence relates to questions of chain stretching frustration and geometric homogeneity [44,47,48,49]. For an object defined by its bounding surface, for every surface point *p* with surface normal vector *n*(*p*), the corresponding medial surface point is defined as *p* + *d*(*p*) · *n*(*p*), where *d*(*p*) is the medial surface distance function which describes the distance from *p* to the corresponding centre of the channel (for the analysis here *p* is a point on the gyroid minimal surface and *p* + *d*(*p*) · *n*(*p*) the corresponding point on the medial surface skeleton). The structure of this simulated gyroid phase is analysed in greater detail in §3.

At the phase transition from the isotropic to the gyroid phase, the orientational order parameters of these intermediate particle tapering systems remain low (figure 5*b*) in sharp contrast to what is seen for strong and weak particle tapering. However, the associated drop in the standard deviation of local orientations *σ*(*α*) (figure 5*c*) suggests that the systems have adopted short-range orientational order. The transition from the isotropic to gyroid phase is also indicated by a feature in the excess pressure (figure 5*a*), which coincides with the point of inflection of *σ*(*α*). Entering the gyroid phase is also associated with a decrease in the gradient of the diffusion coefficient with respect to density (figure 5*d*). A second change in the diffusion characteristic is seen at *ρ*_{g} = 0.62, above which the pear particles no longer traverse the simulation box in the course of a production run and, so, the system is characterized as solid. Owing to the kinetic nature of this method of ascertaining solidification, the transition between the diffusive and solid phase is not defined distinctly and, consequently, is indicated as a dotted line in the phase diagram. Compared to the crystallization of the smectic bilayers, the gyroid systems remain diffusive for higher densities. This may be an indication that a kinetically inaccessible solid lamellar arrangement is the true stable phase at these high densities.

Systems of pears simulated at the phase boundary between the smectic and gyroid phases, i.e. those with *k*_{θ} = 2.4, cannot be unambiguously assigned to either phase since both run cycles (isotropic–smectic–solid and isotropic–gyroid) are observed on re-running simulation sequences. At this apparent transition region between the smectic and gyroid phases, some configurations show long-lived coexistence between regions of parallel and curved bilayers and commensurability with the periodic boundary conditions is an issue. To investigate this behaviour further a second-phase diagram is generated by performing simulation sequences with changing *k*_{θ} and constant *ρ*_{g} (see figure 8). Even though this procedure of changing the shape of particles is hardly possible for experiments on hard particles, we follow this idea which is inspired by the ‘molecular shape concept’ in lipid self-assembly [50,51]. In binary lipid–water mixtures, many facets of the phase diagram can be understood simply by the recognition that changes in pressure, salt concentration or temperature translate to changes in the shape of the individual molecules. Starting from smectic configurations with *k*_{θ} = 2.0 the tapering parameter is increased by Δ*k*_{θ} = 0.1 steps until *k*_{θ} = 6.0. Afterwards the tapering parameter is decreased again until the particles take their original shape. The disposition of the phases is similar compared to the phase diagram in figure 3. However, the grey area indicates major hysteresis effects between the smectic/nematic and gyroid phase, which mainly reduce the parameter space where the gyroid phase forms. The hysteresis effect is more dominant for higher densities which suggests again that for these high densities lamellar structures are similarly stable as the gyroid.

### 2.3. Weak particle tapering: *k*_{θ} > 4.6

For hard pear systems with *k*_{θ} > 4.6, the simulations exhibit four different phases over the density range 0.46 < *ρ*_{g} < 0.65. In addition to the solid-crystalline, smectic-lamella and isotropic phases formed by particles with small tapering parameter values, these weakly tapered pears also adopt nematic order between the isotropic and smectic regions (see figure 6). This window of nematic phase stability integrates straightforwardly into the phase diagram, as can be seen from the excess pressure and diffusion characteristics in figure 5*a*,*d*, neither of which distinguish between the nematic and the gyroid. We note from the phase diagram that the density range of smectic phase stability narrows as *k*_{θ} increases. This behaviour is expected since, for particles with which represent ellipsoids, there is no smectic phase and systems exhibit transitions directly from the solid to the nematic phase [12].

Given that geometrical theories exist for phases, which are observed in lipid/water (shape parameter *s* [36]) and copolymer systems (strong segregation limit [52]), we compare the phase diagrams of lyotropic [23,50], AB diblock copolymer [30] and pear-shaped particle systems (see figures 3 and 8). It is apparent that the phase diagrams exhibit both similarities and differences. Lipids with a shape parameter close to *s* = 1 and copolymers with an A-monomer fraction *f* = 0.5 generate lamellar macrostructures, similar to weakly tapered hard pear-shaped particles. By decreasing *s*, *f* and *k*_{θ}, respectively, and making the particles more and more cone shaped, all three systems first transition into the bicontinuous gyroid phase. However, whereas close to *s* = 0 and *f* = 0 lipid/water and copolymer systems form cylindrical hexagonal phase structures and eventually spherical micelles, hard pear-shaped particles with small *k*_{θ} generate smectic (aka lamellar) structures again, which are, as already mentioned, characteristics for cylindrically shaped lipids and copolymers with *s* = 1 and *f* = 0.5.

## 3. Geometric analysis of the gyroid structure

The simulations presented in §2 demonstrate that the gyroid structure is formed by hard pears with tapering angles of between 12.4° and 23.5°. To draw further comparison between lipid/water, copolymer systems and the behaviour captured in our simulations; in the following sections we characterize the unit cell of the hard-pear gyroid structure, interrogate the microstructure of its minimal surface and perform Voronoi analysis to examine its underlying correlations.

### 3.1. Crystallographic lattice parameter and number of particles per unit cell

The triply periodic feature of the gyroid dictates that its structure is determined by its periodicity. As a consequence, the simulations of the gyroid phase presented in the previous sections are subject to the commensurability issues encountered by any self-assembled system with crystallographic periodicity. In the thermodynamic limit, the lattice parameter *a* and the number *N*_{TUC} of particles within a translational unit cell result from thermodynamic equilibration. Fluctuations in the number of molecules in a unit cell or the size of that cell can then be accommodated since the total number of unit cells and the total number of molecules are infinite. In a simulation, by contrast, only a finite system volume can be modelled. In the hypothetical case where *a* and *N*_{TUC} are known *a priori*, the obvious choices would be to impose a simulation box with dimensions that are integer multiples of the translational unit cell, or an alternative (such as the one based on (110) and (001) directions described below) which represents a larger, differently oriented translational periodicity of the crystal structure. Where a simulation box is not both correctly oriented and commensurate with *a* and *N*_{TUC}, however, its behaviour would be distorted from that of the infinite system. The system might respond by forming defects, adopting a geometric structure that is not thermodynamically stable (but stabilized by the imposed simulation box), or attaining a modified version of the true equilibrium structure.

In practice, *a* and *N*_{TUC} are not known *a priori*, and the challenge lies in detecting rigorous estimates for these parameters from simulations of finite systems—and indeed determining beyond reasonable doubt that an observed geometric phase corresponds to the equilibrium (i.e. infinite) structure, rather than being stabilized by the finite simulation box. Even when the geometry is known, determining *a* and *N*_{TUC} from simulations is not straightforward. This problem has been considered in the context of cluster crystals [53] and a low-density low-temperature gyroid phase [35]; however, the methodology applied therein is not transferrable to our hard pear systems.^{1} In the context of the results summarized in figure 7, we note that, for *N* = 3040 particles with *k*_{θ} = 3.8, a cubic simulation box with edge length 20.84 *w*_{p} was required to form eight unit cells of the gyroid in a 2 × 2 × 2 arrangement. However, from repeated compression sequences performed on a series of further simulated systems at the same number density (*ρ*_{g} = 0.55), we find that equivalent behaviour is obtained for a range of between 3000 and 3200 particles within the simulation box. To investigate these boundary-condition effects further, we therefore perform additional simulations with systems of 10 000 particles in an ensemble where the simulation box is able to adapt its three edge lengths independently while maintaining fixed total volume.

Using this approach to compress hard pear systems with 3.2 < *k*_{θ} < 4.4 to density *ρ*_{g} = 0.555 does not achieve mono-domain gyroids in all cases. In particular, these large systems with 2.4 < *k*_{θ} < 3.2 prove prone to adopting poly-domain curvy bilayer structures, the particles of which cannot be readily assigned to two independent channels using cluster analysis, even after 2 × 10^{7} time steps. Further, the crystallographic (100) direction of the gyroids formed with particles 3.2 < *k*_{θ} < 4.6 are generally not aligned to the (1,0,0)-direction of the simulation box. However, by performing fast Fourier transformation of the density profile of the largest identifiable gyroid domain and calculating the mean peaks of the resultant three-dimensional scattering pattern, it proves possible to reliably determine the appropriate reciprocal lattice vectors of the FCC lattice in Fourier space. From these, the lattice vectors of the gyroid BCC lattice, the volume and consequently the number of particles within each unit cell *N*_{TUC} can be determined for each *k*_{θ}. A representative structure of such a system (*k*_{θ} = 3.8) is shown in figure 1.

This analysis of scattering functions is used to determine the *k*_{θ}-dependence of the mean and variance of the unit cell parameters (see figure 9). These show that the tapering parameter has relatively little influence on the gyroid unit cell, all systems yielding a particle number of 379 ± 11. While the corresponding scattering patterns show that some of the analysed cells are slightly elongated towards the (111) direction, the results of these Fourier analysed systems are in very good agreement with those from the *N* = 3040 simulations described in §2 (see figure 7). Additionally, single unit cell systems with *N* = 380 are simulated as depicted in figure 11 where both channels contain a similar number of particles. Finally, we generate systems of pears within a cuboidal simulation box with edge length ratio and *N* = 2*N*_{TUC} = 760, which display both channel systems. In these, the *x*-axis points into the crystallographic (110), *y* into the (−110) and *z* into the (001) directions of the gyroid.

Estimating the number of particles within the unit cell as *N*_{TUC} = 379 within the gyroid phase at a density *ρ*_{g} = 0.555, the cubic lattice parameter is *a* = 10.4 in units of the width of the pears. It is instructive to compare the number *N*_{TUC} to that found in gyroid-forming lipid systems. Using the fact that the surface area of the gyroid minimal surface is *S* = 3.0915 × *a*^{2} and the average area of a single chain lipid such as monoolein is 37 Å^{2} (at 25°C [56]), we can estimate the number of lipid molecules in an cubic gyroid phase with lattice parameter *a* = 140 Å to be *N* = 2(3.0915 × 140^{2}/37^{2}) ≈ 89 [57].

### 3.2. The gyroid minimal surface

Having determined the unit cell size and shown, in figures 1 and 7, that pear-shaped particles can be assigned to a channel system by cluster analysis of their blunt ends, we now consider another important aspect of the observed gyroid phase—the characterization and analysis of its minimal surface. In binary lipid/water systems, lipids form sheets of bilayers such that the surfactant–water interfaces envelope the gyroid minimal surface. The pears form bilayers as well, which have the same topology as the gyroid. However, unlike lipid/water systems, the pear bilayers have to fully occupy space, such that the distance between the interpenetrating bilayer sheets has to be able to accommodate variable pore radii. As a result, the bilayer thickness cannot be assumed to be constant, and the distance between each pear centre and the hypothetical interface which optimally bisects that bilayer has to be determined.

As a corrective we define a local particle–particle distance measure, called the local bilayer staggering length Δ, as the local distance between the centres of the two interpenetrating sheets. For this calculation, we consider the longitudinal distribution functions *g*(*z*) of the double unit cell systems at a density of *ρ*_{g} = 0.555 to avoid possible errors caused by the deformed gyroid in the 10 000 particles system. Here, *z* is the relative distance distribution between two pears measured along a particle rotational symmetry axis, and the calculation is restricted such that only pears within a cylinder of radius 0.9*w*_{p} around that axis are taken into account. This limiting radius is applied to ensure that pears from the same leaflet are excluded from the calculation. The resultant profiles are given in figure 10.

We extract the local bilayer staggering length from the first peak attributed to the mean relative distance of two next neighbouring pears of interpenetrating sheets, measured along their rotational symmetry axes. The location of this peak shifts to larger *z* with an increase in *k*_{θ} (the error bars in figure 10 indicate the full width at half maximum of the first peak in *g*(*z*) rather than the measurement error). This shift implies that particles with smaller tapering parameter, and consequently a higher tapering angle, interpenetrate more deeply, hence implying a smaller distance between these sheets. The width of the first peak shows that, for all systems, Δ is not constant within the gyroid. This is analysed in further detail in §3.3. The *g*(*z*) curves are terminated before their second peak, which would correspond to the distance between the two bilayers. This is because curvature of the sheets introduces unacceptable levels of uncertainty in the data at this range of *z*. The other noteworthy trend in figure 10 is that the peak heights drop, and the tail at intermediate *z* grows with an increase in *k*_{θ}. This indicates that reducing the pear tapering angle broadens the distribution of observed stagger distances. Finally, we recall that these observations are made in the context that, as noted in §3.1, the overall unit cell size does not change with the tapering parameter.

### 3.3. Voronoi tessellation

In this section, we investigate the relationship between particle packing and minimal surface curvature in our simulated gyroid phase. The gyroid surface is characterized by its mean curvature with the principal curvatures *κ*_{1} and *κ*_{2}. Consequently, the Gauss curvature *K* = *κ*_{1}·*κ*_{2} ≤ 0 and, more precisely, its absolute value |*K*| can be used to quantify curvature. Two approaches are used to relate this measure to simulation configurations. First, multiple configurations are taken from simulations of freely self-assembled single unit cell structures generated by *N* = 380 particles. Second, unit cells are constructed using an MC algorithm which induces the gyroid structure by artificially restricting the pears to the gyroid surface. In figure 11*a*, the channel systems of the restricted and translated self-assembled unit cell are separated by the gyroid minimal surface parametrized by the Enneper–Weierstrass representation [58,59]. We note that the two systems are qualitatively indistinguishable.

The local density
3.1around the *i*th pear is proportional to the inverse volume of the Set Voronoi cell *V*_{i} of particle *i*, see figure 11*c*. The Set Voronoi cell is defined as the space containing all points which are closer to the surface of a given particle than to any other particle in the system [60]. To calculate correlations between these cells and the gyroid, the minimal surfaces shown in figure 11 are triangulated and tessellated according to their intersections with the configuration-derived Voronoi diagrams. As a result, characteristics of the Voronoi diagram and the gyroid are assigned to every point/triangle on the minimal surface. In figure 12, the local Gauss curvature and the mean (from 1000 configurations) volume of intersecting Voronoi cells are, respectively, depicted on two gyroid minimal surfaces. From this, it is apparent that highly curved regions at the necks of the gyroid tend to be intersected by cells with higher volume, whereas more tightly packed particles reside in lower curvature zones (particularly the nodes). This can be interpreted as meaning that higher curvature required lower particle density since this avoids the restrictions otherwise associated with leaflet interdigitation. To achieve a more quantitative measure of this effect, the area *A*_{0} occupied by each Voronoi cell is summarized in a plot of 〈*V*/*A*_{0}〉 against |*K*|—an anti-correlation is observed (see figure 13*a*). Similarly, another anti-correlation is found when the Gauss curvature is, equivalently, plotted against the distance between the gyroid surface and the pear positions (see figure 13*b*). This means that the pears remain further away from the minimal surface at the nodes than at the necks. These two anti-correlations are mutually consistent since low curvature is compatible with high interpenetration.

Comparing these findings for our pear systems with those for lipids and polymers, differences are apparent. The molecular geometry is described by the ‘surfactant parameter’ *v*/*A*_{0} · *l*, where *v* is the effective surfactant chain volume and *l* is the chain length [23]. By invoking the relationship between molecular shape and resulting interface curvatures, one can then express the so-called Steiner's formula as
3.2Lipids and copolymers which form the gyroid surface are conventionally sketched as cones (*v*/*A*_{0} · *l* < 1), whereas in the lamella phase the molecules are considered to be cylinders (*v*/*A*_{0} · *l* = 1) [23,36]. This is illustrated in figure 2. Plotting Steiner's formula in figure 13 and analysing the shape of the tapered Voronoi cell of a single pear (see figure 11*c*) show that the pear particles have a surfactant parameter greater than 1. However, while molecular flexibility means that it is often feasible for lipids and copolymers to have differing surfactant parameters in opposing leaflets, this is patently not the case in the systems studied here. Owing to the interdigitation and fixed particle shape in our hard pear systems (see figure 2), it is also necessary for pear blunt ends to point into the opposing channels. This contrasts with Steiner's theorem, leading to the poor agreement between simulation data and equation (3.2) and leading us to conclude that the mechanisms behind the formation of the gyroid by lipids/copolymers and hard pears may be different.

## 4. Conclusion

We have simulated and numerically analysed systems of hard pear-shaped particles capable of forming the Iad gyroid phase. Through this, we have confirmed that attractive interactions are not necessary for the formation of this structure and that the gyroid can be stabilized by purely entropic effects. The phase diagram obtained here indicates that particles with a range of tapering parameters, corresponding to tapering angles of between 12.4° and 23.5°, are able to form the gyroid phase. Also structural analysis has been used to determine key characteristics of the gyroid, such as the unit cell size and occupancy and the local bilayer staggering length. Correlations between the curvature of the gyroid surface and the Voronoi cells of the pears have shown that a more open structure is adopted in areas of higher curvature and that the gyroid's range of channel widths is accommodated by variation in leaflet interpenetration.

It has been shown that gyroids assembled by our pear systems differ from those generated by lipids in terms of the surfactant parameters appropriate to particles in the bilayer. Additionally, we note that the phase boundaries found here for the gyroid are fundamentally different from those seen in many experimental systems where gyroid phases are present. Conventionally, the gyroid is sandwiched between planar lamellar and hexagonal phases, and its stability is argued in terms of curvature elasticities. This, though, is difficult to reconcile with the phase diagram of figure 3, in which the gyroid borders isotropic and nematic fluids for which there is no curvature elasticity. A possible explanation for this is offered by recent arguments from the Selinger group [61] that sufficiently strong bulk splay-bend coupling between polar and orientational degrees of freedom can destabilize the nematic with respect to supra-molecular modulations (i.e. periodic structures). These arguments, in turn, hark back to the classic paper of Dozov [62] in which the central ideas of the twist-bend nematic were set out. Given that our hard pear systems clearly possess steric coupling between molecular-scale splay and bend, there appears to be a strong argument that the gyroid region observed here is indeed a realization of a modulated splay-bend phase predicted by Dozov and Selinger.

In general, it is interesting to further investigate if it is possible to self-assemble other cubic phases purely entropically—like the polar blue phase, which is predicted for liquid crystals [61], or other minimal surface phases (note recent work on quenched polymeric phases [63–65], including simulations of bicontinuous structures [64,65]). Here it is advisable to be guided by biological systems again and to introduce a second component such as an oligomer or a solvent of hard spheres, which might promote a change in curvature. Alternatively, approaches based on polydisperse systems and active particles may prove effective in influencing mesoscopic structure. The aim here would be to partition space into cells according to a given target desired structure with, preferably, small variation in total volume. This idea is the entropic equivalent to designing potentials which induces particle clustering [66]. In pursuing such a goal, the approach of §3.3 of restricting particles to a given minimal surface and tracking their behaviours, might prove useful.

A recurrent question in the context of gyroid-like phases is chirality, discussed in particular in terms of optical properties [67–69] and in terms of an observed, but unexplained enantiomeric imbalance in butterfly gyroid nanostructures [70,71]. While the hard-pear system described here adopts only the achiral gyroid phase with symmetry , it may also inform on the issue of chiral generalizations. While synthetic self-assembly protocols exist to generate chiral single gyroids with symmetry *I*4_{1}32 [72] or solid replicas thereof [68], these do not break the chiral symmetry, that is, left-handed and right-handed enantiomers occur with equal probability. In these systems, the two network-like domains are chemically distinct, say A and C, yet the probability for A and C to be the right-hand gyroid network is equal at 50 : 50. (Note that more complicated chiral gyroid-like arrangements have been observed in simulations [73] or analysed in terms of geometric free energy concepts [74].) The open question is what is required to make the A moiety adopt the, say, right-handed network with a higher probability than the left-handed one. It has been demonstrated that molecular twist in copolymeric components can affect the mesoscopic structural chirality (and enantiomeric type) of the self-assembled nanostructure [75]. The gyroid phase described here is interesting in the sense that the pears naturally subdivide into two groups, each occupying one of the two labyrinthine domains. It is, though, conceivable that a small adaptation of the particles to embed chiral character may lead to an adapted mesoscale geometry (where the labyrinthine domain that ‘matches’ one particle enantiomeric type is different from the other, e.g. higher density); only one of the two enantiomers of the *I*4_{1}32 gyroid would then form. Further, if a mixture of both enantiomeric particle types was considered, it is conceivable that a chiral microphase separation may result, with RH particles occupying the RH domain and LH particles the LH domain, again producing only one of the two enantiomers of the *I*4_{1}32 single gyroid. While such demixing is unlikely, due to entropy, it is not impossible, and so would be a fascinating topic for future research.

## Competing interests

We declare we have no competing interests.

## Funding

We thank the German Academic Exchange Service and Universities Australia for travel funding through a collaborative grant scheme. We also thank the Cluster of Excellence ‘Engineering of Advanced Materials’ (EAM) and the DFG through the ME1361/11-2 grant for funding.

## Acknowledgements

We gratefully acknowledge John Seddon's help in calculating the number of lipids in bilayer lipid gyroid phases and Mark Lukas for a careful reading of this manuscript. P.W.A.S. acknowledges a Murdoch University Postgraduate Research Scholarship.

## Footnotes

One contribution of 17 to a theme issue ‘Growth and function of complex forms in biological tissue and synthetic self-assembly’.

↵1 The method of Mladek

*et al.*[53] relies on Widom's test particle insertion method [54,55] to determine the chemical potential, but sampling efficiency for this approach is very poor for short-ranged repulsive potentials, particularly at the packing fractions of interest in this work.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.