Linear aggregation and liquid-crystalline order: comparison of Monte Carlo simulation and analytic theory

Linear aggregation and liquid-crystalline order: comparison of Monte Carlo simulation and analytic theory

Tatiana Kuriabova Department of Physics and Liquid Crystal Materials Research Center, University of Colorado at Boulder Institute for Complex Adaptive Matter, University of California, Davis    M. D. Betterton Department of Physics and Liquid Crystal Materials Research Center, University of Colorado at Boulder    Matthew A. Glaser Department of Physics and Liquid Crystal Materials Research Center, University of Colorado at Boulder
July 5, 2019
Abstract

Many soft-matter and biophysical systems are composed of monomers which reversibly assemble into rod-like aggregates. The aggregates can then order into liquid-crystal phases if the density is high enough, and liquid-crystal ordering promotes increased growth of aggregates. Systems that display coupled aggregation and liquid-crystal ordering include wormlike micelles, chromonic liquid crystals, DNA and RNA, and protein polymers and fibrils. Coarse-grained molecular models that capture key features of coupled aggregation and liquid-crystal ordering common to many different systems are lacking; in particular, the role of monomer aspect ratio and aggregate flexibility in controlling the phase behavior are not well understood. Here we study a minimal system of sticky cylinders using Monte Carlo simulations and analytic theory. Cylindrical monomers interact primarily by hard-core interactions but can stack and bind end to end. We present results for several different cylinder aspect ratios and a range of end-to-end binding energies. The phase diagrams are qualitatively similar to those of chromonic liquid crystals, with an isotropic-nematic-columnar triple point. The location of the triple point is sensitive to the monomer aspect ratio.We find that the aggregate persistence length varies with temperature in a way that is controlled by the interaction potential; this suggests that the form of the interaction potential affects the phase behavior of the system. Our analytic theory shows improvement compared to previous theory in quantitatively predicting the I-N transition for relatively stiff aggregates, but requires a better treatment of aggregate flexibility.

I Introduction

Self-assembly of aggregates is ubiquitous in soft matter and biophysics Hamley (2007). Aggregation requires only a pool of monomers with some type of attractive interaction. Reversible aggregation occurs when monomers interact via relatively weak, non-covalent attractions that are typical of soft matter; in this case the aggregates are in equilibrium with the pool of monomers. Reversible aggregation is also called equilibrium polymerization or supramolecular polymerization Ciferri (2002). While different geometries of aggregates are possible (depending on the nature of the attraction), in many cases anisotropic interactions favor linear or filamentous aggregates. These rod-like aggregates can form liquid-crystal (LC) phases; the liquid-crystal order then couples to the aggregation, often promoting the formation of longer aggregates in the LC phase Taylor and Herzfeld (1993); Ciferri (1999). The liquid-crystal ordering of the aggregates can occur even when the monomers alone do not form liquid-crystal phases.

Since the ingredients required for coupled aggregation and LC ordering are so simple, this basic physics occurs for a variety of systems, including worm-like micelles and microemulsions, chromonic liquid crystals, nucleic acids, proteins, and protein assemblies. Some of the first work on coupled aggregation and LC order was inspired by experiments on sickle-cell hemoglobin protein, a fiber-forming protein important in sickle-cell anemia Herzfeld and Briehl (1981); Hentschke and Herzfeld (1989); later work has considered other peptides and proteins that form fibers or fibrils Aggeli et al. (2003); Ciferri (2007); Lee (2009, 2009); Park et al. (2009); Jung and Mezzenga (2010). Other early work considered worm-like micelles, formed either of amphiphilic molecules in water or microemulsions of water and oil which are stabilized by amphiphilic molecules. If the micelles are relatively stiff they can display liquid-crystal phases Khan (1996); these systems have been the subject of extensive experimental study Charvolin et al. (1979); Hendrikx and Charvolin (1981); Hendrikx et al. (1983); Boden et al. (1985, 1990); Quist et al. (1992); Leaver and Holmes (1993); Furo and Halle (1995); Johannesson et al. (1996); von Berlepsch et al. (1996); Angelico et al. (2000); Pereira et al. (2000); Fischer and Callaghan (2001); Goodchild et al. (2007); Kuntz and Walker (2008). Liquid-crystal phases have also been observed in folic acid salts Ciuchi et al. (1994), guanosine derivatives Franz et al. (1994); Spindler et al. (2002), and nucleosome core particles Mangenot et al. (2003). Recent work on chromonic liquid crystals and nucleic acids has renewed interest in coupled aggregation and LC order. Chromonic LCs are formed from relatively flat, disk-like dye molecules. When the dye molecules are suspended in water, hydrophobic interactions drive the stacking and formation of rod-like aggregates that can form LC phases Tiddy et al. (1995); Harrison et al. (1996); Lydon (1998); von Berlepsch et al. (2000); Lydon (2004); Horowitz et al. (2005); Nastishin et al. (2004, 2005); Prasad et al. (2007); Edwards et al. (2008); Park et al. (2008); Joshi et al. (2009); Wu et al. (2009). Hydrophobically-driven end stacking also causes short pieces of DNA and RNA to assemble into aggregates and form LC phases Nakata et al. (2007); Zanchetta et al. (2008, 2008).

Extensive work on the analytic theory of aggregation and LC ordering has been done since the early work of Herzfeld and Briehl Herzfeld and Briehl (1981). Much of the initial work focused on perfectly rigid aggregates Herzfeld and Briehl (1981); Gelbart et al. (1985); McMullen et al. (1985); Herzfeld (1988); Herzfeld and Taylor (1988); Hentschke and Herzfeld (1989); Taylor and Herzfeld (1990, 1991). However, when aggregates are assumed perfectly rigid and excluded volume interactions are treated in the second-virial approximation, a “nematic catastrophe” occurs in which the aggregates in the nematic phase grow infinitely long Odijk (1987); van der Schoot and Cates (1994). Later work therefore emphasized adding aggregate semiflexibility to the models Hentschke (1991); Hentschke and Herzfeld (1991); van der Schoot and Cates (1994); van der Schoot (1994, 1995, 1996); Odijk (1996); Kindt and Gelbart (2001); Lu and Kindt (2004, 2006). Analytic theory qualitatively reproduces results of experiments and simulations, including a first-order I-N transition and longer aggregates in the nematic phase. Some work has also predicted an isotropic-nematic-columnar triple point Taylor and Herzfeld (1991). However, quantitative agreement between simulation and theory remains lacking, particularly for aggregates with higher flexibility Lu and Kindt (2004, 2006).

Surprisingly little simulation work has addressed coupled aggregation and LC order. Some papers have focused on the molecular details that promote aggregate formation in specific systems and considered single aggregatesBast and Hentschke (1996, 1996); Tieleman et al. (2000); Mohanty et al. (2006); Yakovlev and Boek (2007); Turner et al. (2010). Simulating larger systems of interacting aggregates is too computationally costly for atomistic simulation and coarse-grained models are required. Edwards, Henderson, and Pinning did pioneering Monte Carlo simulations of disks formed from hard spheres, with interactions that promoted disk stacking. They observed isotropic, nematic, columnar, and crystalline phases Edwards et al. (1995). A similar approach was used by Maiti et al. to simulate formation of columnar aggregates in chromonic liquid crystals. This work focused on conditions required to find columnar aggregates and didn’t study the phase behavior Maiti et al. (2002). Rouault used a coarse-grained 2D lattice model of wormlike micelles to study how varying the stiffness affects the ordering transition; stiffer molecules display higher orientational ordering Rouault (1998).

More recently, a number of simulation papers have considered hard spheres with an added anisotropic interaction that promotes linear aggregation. Hentschke and coworkers performed molecular dynamics and Monte Carlo simulations of spheres with anisotropic Lennard-Jones potential Fodi and Hentschke (2000); Ouyang and Hentschke (2007). They observed isotropic, nematic, and columnar phases, and observed the narrowing and eventual disappearance of the nematic phase in a I-N-C triple point as the temperature is increased. Similarly, Chatterji and Pandit added an anisotropic interaction to hard spheres via a 3-particle potential that favors linear aggregation Chatterji and Pandit (2003). They observe a first-order I-N transition in Monte Carlo simulations. Some of the most detailed simulation work of a simplified model was done by Lü and Kindt, who performed Monte Carlo simulation of spheres with sticky “patches” that allow assembly into linear aggregates Lu and Kindt (2004, 2006). Their simulations did not show the columnar LC phase but did see the I-N transition.

The key physical ingredients required for aggregation coupled to LC order are (1) monomers with an anisotropic attractive interaction that promotes the formation of long, thin filaments and (2) excluded volume interactions between aggregates that promote liquid-crystalline order at high aggregate density. While other physical effects such as electrostatic interactions are clearly important in some systems Aggeli et al. (2003); Mangenot et al. (2003); Kroeger et al. (2007); Kuntz and Walker (2008); Jung and Mezzenga (2010), a minimal physical picture which incorporates filament formation and excluded-volume interactions is a valuable starting point to understand common features of the diverse experimental systems.

In this work we develop a coarse-grained model of aggregation-induced LC order, and study the model using Monte Carlo simulations and analytic theory. Our goal is to ultimately compare the simulation and theoretical results to experimental data on chromonic and DNA liquid crystals, so we have chosen the form of the model and the approximations used to facilitate this comparison. Motivated by recent work on nucleic acids and chromonic liquid crystals, we consider molecules with anisotropy both in shape and interactions. Therefore, our monomer is a cylinder of length and diameter , so the aspect ratio . The cylinders are sticky: if two cylinders stack so that their circular ends are sufficiently close together, they experience an attractive interaction. Other than this binding energy, the cylinders experience only a hard-core repulsion. This coarse-grained model captures the main physical effects governing linear aggregation and mesophase formation in these systems, namely excluded volume interactions and end-to-end association. For the appropriate range of binding energy, density of cylindrical monomers, and temperature, the Monte Carlo simulations of sticky cylinders find linear aggregates and a range of phases, including isotropic, nematic, columnar liquid crystal, and columnar crystal phases.

The hard cylinder system (corresponding to zero binding energy between cylinders) has been studied previously in simulations by Blaak et al. Blaak et al. (1999) and by Ibarra-Avalos et al. Ibarra-Avalos et al. (2007), but to our knowledge this is the first simulation study of hard cylinders with attractive interactions.

In simple (non-aggregating) LC systems, the rod aspect ratio is a primary determinant of possible phases and the location of the phase transitions. The role of monomer aspect ratio in controlling the phase behavior of aggregating LC systems is not well understood. Our sticky-cylinder model enables systematic investigation of the dependence of aggregation and phase behavior on monomer aspect ratio for systems ranging from disk-like chromonic liquid crystals () to duplex DNA oligomers ().

In section II we discuss the simulation model and the Monte Carlo simulation methods, including cluster moves we used to improve configurational sampling. Section III outlines the analytic theory, the free energy minimization, and the I-N phase coexistence equations. In section IV we present the simulation results, including the phase behavior, order parameter and pair correlation functions, aggregate length distributions, and aggregate persistence length. The comparison of analytic theory and simulation for the I-N transition is presented in section V, with a comparison both to our simulations and to the work of Lü and Kindt Lu and Kindt (2004). Section VI is the conclusion. Appendix A describes some calculation details of the analytic theory.

Ii Simulation methods

The monomers are hard cylinders of length and diameter with short-range attractive interactions between cylinder ends. The short-range attraction between cylinder ends is a generalized ramp potential,

(1)

where is the distance between the centers of the circular end faces of two neighboring cylinders. An exponent of corresponds to a linear ramp potential, while the limit corresponds to a square well potential. We chose and , parameter values that we empirically found promote the formation of linear (as opposed to branched) aggregates; larger values of and give rise to branched networks.

This model depends on three dimensionless parameters: cylinder aspect ratio , packing fraction (or, alternatively, dimensionless pressure ), and dimensionless binding energy . Here is the pressure, is the inverse temperature, is the cylinder volume, and is the volume per particle (inverse number density). Unless otherwise noted, the cylinder diameter is the unit of length.

We investigated systems of sticky cylinders with three aspect ratios, , 1, and 2. We carried out Monte Carlo simulations over a range of pressures for binding energies ranging from to .

The persistence length of sticky cylinder aggregates is implicitly determined by the cylinder-cylinder pair interaction potential; no explicit bond angle bending potential is included. An implicit dependence of on the nature and strength of effective intermolecular interactions is a key feature of nucleic acid and chromonic aggregates, which may influence the temperature dependence of the I-N coexistence curve. In this respect our model is more specific than that of Lü and Kindt Lu and Kindt (2004), which included a bond angle bending potential, enabling the binding energy and persistence length to be varied independently.

ii.1 Monte Carlo Simulations

We carried out Monte Carlo (MC) simulations of sticky cylinder systems with periodic boundary conditions, with fixed monomer number , pressure , and temperature Frenkel and Smit (2001). Our simulations employed several types of MC move, including small displacements and rotations of single particles and changes in simulation cell volume and shape under an applied pressure. To improve configurational sampling, we also used cluster moves (see below), and flip moves ( reorientation of individual cylinders). An MC cycle includes trial single-particle displacement/rotation moves, one trial volume-changing move, and (in cases where these are used) trial cluster moves and/or trial flip moves.

We simulated systems with between and cylinders. For each aspect ratio and binding energy we performed both an expansion run (a series of simulations for decreasing pressure starting from a high-density columnar crystal) and a compression run (a series of simulations for increasing pressure starting from a low-density isotropic fluid). For each value of the pressure, we simulated MC cycles, and measured thermodynamic and structural properties during the final MC cycles. To accommodate crystalline and LC phases of arbitrary symmetry without defects or strain, we utilized a fully flexible simulation cell, except in cases where this led to highly elongated or deformed simulation cells (e.g., at low densities).

ii.1.1 Cluster moves.  

Long aggregates undergo slow effective translational and rotational diffusion in simulations that utilize only single-particle MC moves. To improve configurational sampling of aggregates, we used cluster moves: simultaneous displacements and rotations of groups of neighboring cylinders. Similar cluster moves have been used in simulations of ionic fluids Orkoulas and Panagiotopoulos (1994). Clusters are defined by proximity: any two cylinders with belong to the same cluster. The cluster cutoff can be varied to control the average cluster size . For the simulations described here, we adjusted to maintain for all pressures, binding energies, and aspect ratios.

A trial cluster MC move consists of the following five steps. 1. Select a root cylinder at random, and find all cylinders in the cluster containing the root cylinder. 2. Generate a trial displacement and rotation of the cylinders in the cluster. The cluster is rotated as a rigid body about an axis passing through the center of the root cylinder. We performed two types of rotation move with equal probability: spin moves, for which the axis of rotation is the axis of the root cylinder, and tilt moves, for which the axis of rotation is a randomly selected axis perpendicular to the root cylinder axis. 3. Reject the trial move if cylinder-cylinder overlaps are found in the new configuration. 4. Reject the trial move if the number of cylinders in the cluster has changed (i.e., if cylinders not in the original cluster are members of the cluster in the new configuration), as such moves violate detailed balance. 5. If no overlaps are found, and if the cluster size is unchanged, accept the move with probability , where is the change in potential energy.

Cluster moves were used in all simulations of cylinders having a nonzero binding energy (). Flip moves were previously introduced by Blaak et al. Blaak et al. (1999) to improve configurational sampling for hard cylinders at high density, in the vicinity of a partially ordered (cubatic-like) intermediate phase. We used flip moves for small binding energy () for the same reason.

Cluster moves moderately improve sampling efficiency, as measured by the cylinder orientational decorrelation rate. The orientational autocorrelation function is

(2)

where is a unit vector along the axis of cylinder , and the angle brackets denote an average over all time origins (with “time” measured in units of MC cycles). In Figure 1 we show the orientational autocorrelation function for an example simulation in the isotropic phase. In this simulation, cylinders assemble into aggregates containing about 10 cylinders on average, which leads to slow effective rotational diffusion. Note that decays exponentially for large . Cluster moves reduce the orientational correlation time by approximately a factor of two, from MC cycles in a simulation with only single-particle moves (dashed curve) to MC cycles in a simulation incorporating cluster moves (solid curve). This decrease occurs despite the fact that the acceptance rate for cluster moves is rather low (around in this case). Performing cluster moves per MC cycle adds approximately to the computational cost of the simulation, so the gain in efficiency is around . Similar performance was found for other parameter values. More substantial improvements in sampling efficiency will likely require more sophisticated schemes such as configurational bias Monte Carlo Frenkel and Smit (2001) or the cluster cleaving method of Whitelam and Geissler Whitelam and Geissler (2007, 2008).

Figure 1: Orientational autocorrelation function for an isotropic fluid of sticky cylinders with aspect ratio and binding energy at a pressure of . Results are shown for simulations with (solid curve) and without (dashed curve) cluster MC moves. Cluster moves decrease the orientational decorrelation time by approximately a factor of two.

ii.1.2 Overlap test.  

The test for cylinder-cylinder overlaps is more time-consuming than the corresponding test for hard spheres, spherocylinders or ellipsoids. We use a procedure similar to that used in previous simulation studies of hard cylinders Blaak et al. (1999); Ibarra-Avalos et al. (2007), and that consists of three tests: 1. Spherocylinder overlap: test for overlaps between the enclosing spherocylinders (obtained by adding hemispherical endcaps to hard cylinders). If the spherocylinders don’t overlap, then neither do the hard cylinders. 2. Disk-disk overlap: if an overlap is found in the first step, test for overlaps between the flat end faces of the cylinders. 3. Disk-cylinder overlap: if no overlaps are found in the second step, test for overlaps between the flat end faces of one cylinder and the cylindrical surface of the other.

If the enclosing spherocylinders overlap, but no cylinder-cylinder overlaps are found, then the two cylinders may interact via the short-range attractive potential (Eq. (1)). Computing the attractive interaction between sticky cylinders requires minimal additional effort, because the distance between the centers of the circular end faces of the cylinders is available from the spherocylinder overlap test.

Iii Theory

To compare the I-N transition in simulation and theory, we extended the work of van der Schoot and Cates van der Schoot and Cates (1994). This analytic theory considers cylindrical monomers (of length and diameter ) that reversibly assemble into linear aggregates. The form of the free energy is determined from three key assumptions. First, the aggregates are treated as nearly rigid cylinders. The ratio of aggregate length to the persistence length is assumed small, . Second, the aggregates interact via steric repulsion. Only pairwise interactions are taken into account within the second virial approximation, but the Parsons-Lee approximation is used to extend the theory to higher density. Third, the assembly of monomers into aggregates is driven by an energetic penalty for monomers to have free ends. We assume that is independent of the aggregate length .

The assumption of nearly rigid aggregates does not hold in our simulations: the results shown in this paper have a ratio of aggregate lengths to persistence length . We discuss below the limitations of the theory in this flexibility regime.

The work of van der Schoot and Cates considered two-particle interactions in the second virial approximation. Since this approximation is accurate only in the dilute limit, we used the Parsons-Lee approximation to extend the theory to the regime of higher packing fractions Parsons (1979). The Parsons-Lee method effectively inlcudes contributions from higher virial coefficients in the interaction free energy.

The total free energy of the system depends on the density of monomers that belong to a subset of aggregates of contour lengths between and and have orientations within the solid angle and with respect to some specified direction . The monomers have length , diameter , aspect ratio , and volume . The free energy per unit volume (in units of ) is

(3)

All length integrals (over ) are from zero to infinity and all angular integrals (over ) are over of solid angle. The first term is the free energy of a mixture of ideal gases of aggregates. Aggregates of different lengths and orientations are treated as different chemical species. The second term describes semiflexibility of aggregates; it is a perturbative correction to the orientational degrees of freedom (implicitly included in the first term) due to the flexibility of aggregates Khokhlov and Semenov (1982). The third term describes the excluded-volume interaction between pairs of monomers with relative orientations described by . The prefactor reflects the Parsons-Lee correction for higher virial coefficients Parsons (1979); Cinacchi et al. (2004); Cuetos et al. (2007) based on the Carnahan-Starling expression for the correlation function of hard spheres Carnahan and Starling (1969),

(4)

where is the packing fraction of monomers. In the usual second-virial approximation . The last term in Eq. (3) is the free energy of polymerization, which is proportional to the total number of aggregates in the system. In writing Eq. (3) we assumed that in equilibrium the distribution of monomers is spatially homogenous and that the conformational degrees of freedom of aggregates decouple from the interaction degrees of freedom.

For a closed system the total number of monomers is fixed, leading to the normalization condition

(5)

To determine the I-N phase diagram we take the standard approach of requiring mechanical and diffusive equilibrium of the two phases at coexistence. We calculate the osmotic pressure and chemical potential per monomer in each phase. Coexistence requires

(6)

The solution of Eqs. (6) provides information about the packing fractions and of monomers at the coexistence boundaries, the mean aggregation number (and distribution of aggregation numbers), as well as the order parameter , where is the angle between the monomer’s axis and the director , and averaging is performed over all monomers in the system.

In the sections below we outline the calculations and resulting equations for the isotropic and nematic phases as well as the numerical methods.

iii.1 Isotropic phase

In the isotropic phase the monomers are distributed uniformly both in space and in orientation, so independent of angle. The free energy functional Eq. (3) for the total free energy of the system in the isotropic phase reduces to

(7)

where the monomer aspect ratio and the monomer packing fraction .

Functional minimization of the free energy is subject to the normalization condition Eq. (5), which can be handled with a Lagrange multiplier :

(8)

This results in a length distribution of aggregates that is, as assumed, independent of orientation

(9)

with

(10)

This solution corresponds to an exponential distribution for the density of aggregates as a function of aggregate length :

(11)

The chemical potential per monomer is given by

(12)

where is the derivative of with respect to . The osmotic pressure is

(13)
Figure 2: Instantaneous configurations from MC simulations of sticky cylinders with binding energy in various phases: (a) columnar LC phase, ; (b) nematic phase, ; (c) isotropic fluid phase, . Cylinder color encodes orientation: the RGB color index of cylinder is ), where is a unit vector along the cylinder axis.
Figure 3: Equation of state for measured on expansion from a high-density columnar crystal state (symbols). The horizontal lines indicate the approximate location of first-order phase transitions.

iii.2 Nematic phase

In the nematic phase the system is spatially homogenous but orientationally ordered due to alignment of aggregates along a spontaneously chosen direction . The system preserves azimuthal symmetry with respect to rotations about . Therefore, in spherical coordinates the density is independent of the polar angle and

(14)

with , where is the angle to the director . As in the previous section, we minimize the system’s free energy with respect to the orientation-dependent density function ,

(15)

This results in the integro-differential equation

(16)

The interaction kernel is

(17)

Solving Eq.(16) for the function that minimizes the free energy is a difficult task that requires either numerical solution or introduction of a trial function. We chose a mixed approach that starts with a trial function and uses numerical solution of the resulting algebraic equations. We adopted the trial function proposed in reference van der Schoot and Cates (1994):

(18)

where the parameter is given by Eq. (10) and is an orientation-dependent aggregation number.

The free energy of the nematic phase as a functional of becomes

(19)

The bars denote dimensionless variables: and . The next step is functional minimization of the free energy to determine an integro-differential equation for . We find a series solution for by expanding in Legendre polynomials, which turns the equation for into a system of algebraic equations that we solve numerically. For details of the calculations, see Appendix A. Here we summarize our results for the chemical potential per monomer and osmotic pressure in the nematic phase:

(20)
(21)

where are integrals involving that are defined in Eqs. (33)-(35).

Figure 4: Phase diagrams of the sticky cylinder system as a function of binding energy and packing fraction (left) and binding energy and pressure (right), for aspect ratios (a), (b), and (c). A variety of phases are observed, including isotropic (), nematic (), columnar liquid crystal (), columnar crystal (), and cubatic-like () phases. Shaded areas correspond to regions of two-phase coexistence.

The mean aggregation number in the nematic phase is

(22)

where

(23)

is the distribution function of aggregates of length in the nematic phase. The order parameter of aggregates of length is

(24)

where is the second-order Legendre polynomial. The order parameter (averaged over all monomers in the system) is

(25)
Figure 5: Hysteresis in the equation of state for with (left) and (right). The horizontal lines indicate the approximate location of first-order phase transitions observed on expansion from a high-density columnar crystal state (solid lines) and on compression from a low-density isotropic fluid state (dashed lines).
Figure 6: Nematic order parameter as a function of pressure for and .

Iv Simulation results

In Figure 2 we illustrate the phase behavior with sample snapshots of the simulations. Instantaneous configurations from simulations of the sticky-cylinder system with binding energy are shown. This figure displays the typical phase sequence observed upon expansion from a high-density columnar crystal for large binding energy, which includes the columnar LC phase, the nematic phase, and the isotropic fluid phase.

To determine phase boundaries, we determined the locations of plateaus in the equation of state obtained on expansion from a high density columnar crystal state. Such plateaus indicate first-order phase transitions. Phase identification was based on structural properties such as the pair-correlation function and nematic order parameter, as discussed below. For example, the equation of state for and various values of is shown in Figure 3, with the approximate location of first-order phase transitions indicated by horizontal solid lines.

Phase diagrams of the sticky-cylinder system for aspect ratios , 1 and 2 are shown in Figure 4. All three phase diagrams share the same key features. First, both nematic () and columnar LC () phases are observed for large . The range of nematic stability (in terms of pressure or packing fraction) increases with increasing binding energy. Second, an isotropic-nematic-columnar triple point is observed at intermediate binding energies. For binding energies below the triple point, there is a direct transition from the isotropic phase to the columnar LC phase.

We also find a columnar crystal () phase at high density, and the system exhibits a cubatic-like phase for small . The cubatic-like phase was previously described by Blaak et al. Blaak et al. (1999) in the hard cylinder system. As the focus of this paper is on LC phases, we have not investigated the crystalline or cubatic phases in detail, and the phase boundaries involving these phases should be regarded as schematic.

Comparison of the equation of state obtained on compression from an isotropic fluid state with that obtained on expansion from a high-density columnar crystal state reveals considerable hysteresis (Figure 5). This hysteresis becomes increasingly pronounced with decreasing , and compression from the low-density isotropic state for small typically leads to a high-density disordered (jammed) state, likely due to the dominance of hard-core interactions in this limit. For large , the I-N transition occurs at relatively low packing fractions, which facilitates annealing into a well-ordered nematic state starting from an isotropic state. Because we generally observed slow kinetics of formation of an ordered state starting from disordered (isotropic) initial configuration, we used isotherms obtained on expansion (e.g., Figure 3) to construct the phase diagrams in Figure 4. This implies that the phase boundaries shown are lower limits (both in density and pressure), and should be considered provisional. In future work we will refine the phase diagram using free energy calculations.

Figure 7: Instantaneous configurations (left) and pair distribution functions in planes containing the nematic director (, center) and perpendicular to the nematic director (, right) for and . (a) Columnar LC phase at , obtained by expansion from a high-density columnar crystal. There are no long-range correlations in the positions of cylinders, as shown by , and a well-defined hexagonal lattice of columns is apparent in . (b) Nematic LC phase at , obtained by compression from a low-density isotropic fluid. Again, there are no long-range correlations in the positions of cylinders, as shown by . The plot of reveals isotropic and rapidly decaying positional correlations in the plane perpendicular to , indicative of nematic order.

iv.1 Order parameter and correlation functions

We measured nematic order by calculating the tensor

(26)

where and refer to cartesian components, indexes the sum over monomers, and is the Kronecker delta function. The nematic order parameter is the largest eigenvalue of the average order tensor , and the nematic director is the corresponding eigenvector. Typical results for as a function of pressure are shown in Figure 6, for and . The order parameter jumps abruptly from a small value () in the isotropic phase to in the nematic phase, and exhibits a further abrupt jump to in the columnar LC phase. We list values of the nematic order parameter in the nematic or columnar liquid crystal phases at coexistence ( or ) with the isotropic phase in Table 1. The values of are generally large (), indicating a high degree of orientational order in all LC phases investigated.

To distinguish between the various orientationally ordered phases we computed the three-dimensional pair distribution function and examined positional correlations as a function of pair separations parallel and perpendicular to the nematic director. The pair distribution function is defined as the normalized probability of finding a pair of particles with separation ,

(27)

where is the average number density of particles and is the Dirac delta function.

Figure 7 shows example plots of pair distribution functions. The nematic director defines the direction (), so is in a plane containing the nematic director and is in a plane perpendicular to the nematic director. For this set of simulations, and . The columnar LC phase shown in Figure 7(a) is characterized by a well-defined hexagonal lattice of columns and by the absence of long-range correlations in the positions of cylinders. The nematic LC phase shown in Figure 7(b) exhibits isotropic and rapidly decaying positional correlations in the plane perpendicular to .

Figure 8: One-dimensional pair distribution function along a line parallel to the nematic director in the columnar LC phase (blue squares) and the nematic phase (red circles) for and .

We computed the positional correlation length along the nematic director in both the nematic and columnar LC phases from the one-dimensional pair distribution function along a line parallel to the nematic director, . Figure 8 shows typical results for the columnar LC and nematic phases for and . We find exponential decay of correlations to a constant in both phases (in the nematic phase approaches 1 for large , while in the columnar LC phase we estimated that approaches 3.5 for large ). The correlation length determined by fitting the exponential decay is in the nematic phase at and in the columnar LC phase at . In the same simulations, the mean aggregate lengths are 44 in the nematic phase and 352 in the columnar LC phase. This gives a ratio of aggregate length to correlation length that is large: in the nematic phase and in the columnar LC phase. (For more discussion of aggregate lengths, see Section IV.2 below.) In other words, we find that the correlation length can differ from the mean aggregate length by 1-2 orders of magnitude. This in turn suggests that it may not be feasible to determine the mean aggregation number from X-ray scattering experiments that measure correlation lengths, such as the experiments of references Park et al. (2008); Joshi et al. (2009).

Figure 9: Mean aggregation number as a function of pressure (left) and aggregation number distributions (right) for and , measured on expansion from a high-density columnar crystal state. The nematic phase displays fast exponential decay for small crossing over to slow exponential decay for large (bi-exponential distribution). The small peak near in the nematic aggregation number distribution is a finite-size artifact, associated with aggregates that connect to themselves across the periodic boundaries.

iv.2 Aggregate Statistics

We measured the aggregation number distribution function and its first moment, the mean aggregation number , in the isotropic, nematic, and columnar LC phases. Here is the aggregation number (the number of monomers in an aggregate). Aggregates are defined by proximity: any pair of cylinders with endpoint separation smaller than belongs to the same aggregate.

The use of a finite simulation volume imposes a cutoff on the aggregate-length distribution. This finite-size effect is particularly evident in the limit of large binding energies, where the mean aggregate length can exceed the maximum linear dimension of the periodic simulation cell. To ameliorate this effect, the initial columnar crystal configurations used for expansion runs were constructed with an oblique orientation of columns with respect to the edge vectors of the periodic simulation cell, so that a given column traverses the periodic box multiple times before intersecting itself. With this expedient, we are able to measure aggregation number distributions over a wide range of conditions, although small finite-size artifacts associated with aggregates that connect to themselves across the periodic boundaries are observed in many of the expansion runs (see Figure 9).

We illustrate the typical aggregation behavior in Figure 9, which shows the mean aggregation number as a function of pressure (left) and aggregation number distributions in the isotropic and nematic states (right) for and , measured upon expansion from a high-density columnar crystal state. The mean aggregation number jumps discontinuously at the I-N transition, and increases strongly with increasing pressure in the nematic phase. This behavior is qualitatively consistent with that observed in previous simulation and theoretical studies of semiflexible linear aggregates. The distribution has an exponential dependence on in the isotropic phase, but displays biexponential behavior in the nematic phase, i.e., fast exponential decay for small crossing over to slow exponential decay for large . Such biexponential behavior has been noted previously by Lü and Kindt Lu and Kindt (2004), who argued that the short length scale behavior of comes from a distinct population of short aggregates with relatively low orientational order. Biexponential behavior in the nematic phase is also predicted by the analytic theory; see Figure 11. In the simulations we also observed a biexponential aggregation number distribution in the columnar LC phase.

Mean aggregation numbers in the isotropic phase at coexistence () and in the nematic or columnar LC phase at coexistence with the isotropic phase ( or ) are listed in Table 1. The mean aggregation number at coexistence tends to increase with increasing , although the trend in the nematic or columnar LC phase at coexistence isn’t strictly monotonic. We note that the mean aggregation numbers at coexistence decrease with increasing cylinder aspect ratio for a given , possibly because the corresponding packing fractions at coexistence decrease with increasing aspect ratio.

0.5 12 1920 14.9 0.66 9.96 24.08 - 0.704 -
10 1680 13.5 0.67 6.91 18.75 - 0.729 -
8 1680 11.4 0.68 4.06 - 45.90 - 0.908
6 1440 9.2 0.69 2.74 - 27.44 - 0.895
4 1440 7.2 0.70 2.10 - 15.85 - 0.872
2 1200 5.4 0.70 1.78 - 14.96 - 0.869
0 1200 - - - - - - -
1 12 1440 31.4 1.15 5.65 14.98 - 0.761 -
10 1200 27.3 1.16 4.15 16.38 - 0.841 -
8 1200 22.3 1.17 2.83 - 58.42 - 0.947
6 1200 17.7 1.17 1.89 - 12.95 - 0.903
4 960 13.3 1.18 1.55 - 11.04 - 0.902
2 960 - - - - - - -
0 720 - - - - - - -
2 12 1440 62.6 2.14 3.80 11.88 - 0.840 -
10 1200 51.0 2.15 2.54 11.15 - 0.873 -
8 1200 44.4 2.16 1.82 7.99 - 0.870 -
6 960 32.2 2.16 1.43 - 19.58 - 0.965
4 960 24.3 2.17 1.22 - 9.31 - 0.956
2 720 - - - - - - -
0 720 - - - - - - -
Table 1: Properties of the sticky-cylinder system measured in MC simulations. Columns contain the cylinder aspect ratio , the dimensionless binding energy , the number of particles , the measured persistence length , the mean intra-aggregate bond length , the mean aggregation number in the isotropic phase at coexistence , the mean aggregation number in the nematic or columnar LC phase at coexistence with the isotropic phase ( or ), and the nematic order parameter in the nematic or columnar phase at coexistence with the isotropic phase ( or ).

iv.3 Persistence Length

As discussed above, the persistence length of sticky-cylinder aggregates is determined by cylinder-cylinder interactions. We measured the persistence length from the bond orientational correlation function,

(28)

where is a unit vector along the th bond (the vector joining the centers of two neighboring cylinders) in an aggregate, the contour length , and is the average bond length. In the isotropic phase, decays exponentially with increasing , with a decay constant that can be identified with the persistence length , . By fitting this exponential decay to the simulation data, we determined the persistence length. For a given and , we find that the calculated persistence length is nearly independent of pressure within the isotropic phase, which gives us confidence that this procedure measures the intrinsic bending rigidity of aggregates. Table 1 lists the persistence length and average bond length as a function of and , where both quantities are averaged over isotropic state points.

The persistence length increases approximately linearly with increasing binding energy, and is also approximately proportional to cylinder aspect ratio for a given binding energy. Because of this, the reduced persistence length is reasonably well approximated as a universal linear function of , independent of cylinder aspect ratio, as shown in Figure 10. The solid line in Figure 10 is the best linear fit to the data, . The ratio of persistence length to monomer length approaches a nonzero constant in the limit of zero binding energy. This indicates that there is a significant entropic contribution to the persistence length due to hard-core interactions between neighboring cylinders in an aggregate.

Figure 10: Ratio of persistence length to monomer length as a function of for aspect ratios (circles), (squares), and (triangles). The data for all three aspect ratios is well approximated by a universal linear function, , which represents the best linear fit to the data (solid line).

We expect that the precise way that varies with to depend strongly on the form of the intermolecular potential. This dependence can be estimated by simple dimensional analysis arguments. When two bound cylinders interact, thermal fluctuations will typically cause their separation to vary to about 1 above their minimum interaction energy . Therefore a typical separation . This suggests a characteristic angle defined by the distance and the cylinder radius : . The typical angle is therefore . This gives an estimate of the persistence length, based on the relation and assuming , of . In other words, the persistence length should vary linearly with , with a slope of 2, as observed in simulations.

For other forms of the potential, this dependence would change. Using the same dimensional analysis argument, we would predict that for an exponent of in the interaction potential (Eqn. (1)), the scaling is . In the limit of a square well () the persistence length would become independent of . Because the phase behavior of aggregates is strongly dependent on the persistence length, this implies a strong dependence of the phase behavior on the form of the monomer-monomer interaction potential, in addition to the obvious dependence on the binding energy. In the future it would be interesting to investigate the interplay of these two effects, focusing how sticky-cylinder phase behavior varies with the form of the interaction potential.

The aggregates in this study are relatively flexible: the ratio of persistence length to mean aggregate length for the nematic phase at coexistence ranges from for and to for and .

V Theory-simulation comparison

We compared our analytic theory results both to the simulations of Lü and Kindt Lu and Kindt (2004) and to our simulations. These comparison illustrate the key role of aggregate flexibility in determining the theory-simulation agreement.

Figure 11: Theoretical predictions of aggregate length and order parameter distributions. Top, aggregate length distributions at I-N coexistence of the isotropic (I, dashed curves) and nematic (N, solid curves) phases with and (subscripts). Note that in this figure the predictions of the simple second-virial approximation and the Parsons-Lee approximation overlap. Bottom, length-dependent aggregate order parameter at I-N coexistence in the nematic phase with and (subscripts). Solid curves: simple second-virial approximation. Dashed curves: Parsons-Lee approximation. In both cases, the association constant , for comparison with the simulations of Lü and Kindt Lu and Kindt (2004).
Figure 12: Theoretical predictions of mean aggregation number in the isotropic and nematic phases at I-N coexistence. (A), . (B), . Comparison of theory (solid and dashed curves) to simulation results of Lü and Kindt Lu and Kindt (2004) (points). Curves denote borders of regions of I-N coexistence. Solid curves: simple second-virial approximation. Dashed curves: Parsons-Lee approximation. Points and bars: simulation results of Lü and Kindt Lu and Kindt (2004). Note that for (B) Lü and Kindt gave ranges of measured values, which we plot using horizontal bars, while for (A) Lü and Kindt gave only mean values, which we plot using points.

First, we calculated the aggregate length distributions and order parameter as a function of length to compare with the simulation results of Lü and Kindt (Figure 11). Lü and Kindt considered two persistence length values, and Lu and Kindt (2004). Note that is the persistence length divided by the monomer length. (In Lü and Kindt’s work, the monomers are spheres.) For this figure, we used the the association constant .

A significant and unexpected result is the calculated biexponential aggregate length distribution in the nematic phase. This distribution arises naturally from the free-energy minimization (Figure 11, top). To our knowledge, our work is the first to find this biexponential distribution in analytic theory. Lü and Kindt’s initial analytic work found only a single exponential distribution Lu and Kindt (2004); in later work they found that assuming a biexponential distribution improved the theory-simulation agreement Lu and Kindt (2006).

Overall we find good quantitative agreement between the predictions of the analytic theory for aggregate length distributions and order parameter as a function of aggregate length (compare Figure 11 to Figure 2 of reference Lu and Kindt (2004)). The aggregate length distribution is exponential in the isotropic phase and biexponential in the nematic phase. The nematic phase also shows significantly longer aggregates. Stiffer aggregates also tend to be longer in the nematic phase. The average alignment depends strongly on aggregate length in the nematic phase, with short aggregates only weakly aligned and a rapid increase in the length-dependent order parameter with aggregate length.

Figure 13: Comparison of nematic order parameter at I-N coexistence predicted by theory (dashed curves) to simulation results of Lü and Kindt Lu and Kindt (2004) (points). Persistence lengths and indicated by subscripts.

Next we compared our predictions for the aggregate length and nematic order parameter at I-N coexistence to the simulations of Lü and Kindt for a range of monomer binding energies . In Figure 12 we plot the average aggregation number in the isotropic and nematic phases at I-N coexistence for a range of values of . The curves are predictions of our analytic theory and the bars and points are values from the simulations of Lü and Kindt Lu and Kindt (2004). Note that for the shorter persistence length Lü and Kindt gave ranges of measured values, which we plot using horizontal bars, while for the larger persistence length Lü and Kindt gave only mean values, which we plot using points. (The same holds for the two figures below in which we compare to Lü and Kindt’s work). This result illustrates the strong enhancement of aggregation in the nematic phase. The theory clearly captures the correct qualitative trends of the simulations: as the binding energy increases, the aggregates become longer in both the isotropic and the nematic phases. Quantitatively agreement is reasonable for the larger persistence length, but for the shorter persistence length the predicted values of are off by approximately 50%. This illustrates the decreased quantitative accuracy of the theory as the aggregates become less flexible.

The order parameter in the nematic phase at I-N coexistence is shown in Figure 13. Here our calculations for the order parameter agree well with the simulation results. As expected, the stiffer, longer aggregates show a higher average order parameter that is less sensitive to the binding energy.

Figure 14: Comparison of theory (solid and dashed curves) to simulation results of Lü and Kindt Lu and Kindt (2004) (points). (A), . (B), . Curves denote borders of regions of I-N coexistence. Solid curves: simple second-virial approximation. Dashed curves: Parsons-Lee approximation.
Figure 15: Comparison of I-N coexistence in theory (curves) and simulation (points). Top, aggregation number. Botton, phase diagarm Curves denote borders of regions of I-N coexistence; gray: onset of I-N coexistence; black: upper limit of I-N coexistence. Solid curve: . Dashed curves: . Dotted curves: . In doing the calculations we used the empiricial persistence length found from simulations . Triangles: . Squares: . Circles: .

Our final comparison with the simulation results of Lü and Kindt is the phase diagrams shown in Figure 14. The analytic theory agrees well with the simulation results and less well with the simulation results. For the analytic theory gives errors of up to about 50% in the predicted volume fraction of the coexistence region. Perhaps surprisingly, we find that the simple second-virial approximation (solid curves in Figure 14) are closer to the simulation results than the results that include the Parsons-Lee approximation (dashed curves).

In our simulations, the ratio of persistence length to monomer length is significantly less than in the work of Lü and Kindt; we find values in the I-N coexistence range (Table 1). The ratio of mean aggregate length to persistence length in this range (Table 1). Because our analytic theory assumes nearly rigid aggregates () we expect that the theory will not show good quantitative agreement with the simulations in this regime. We note that in our simulations the persistence length is not a control parameter but is determined by the balance of binding energy and entropy of the aggregates; we found empirically that (figure 10). We used this empirical relationship in the analytic theory.

In Figure 15 we show the aggregation number and volume fraction at the borders of I-N coexistence. While the qualitative trends in the simulations are correctly captured by theory, there is significant quantitative disagreement. The analytic theory predicts significantly larger aggregation numbers and significantly lower volume fractions than found in the simulations. The biggest discrepancy is with the aggregation numbers, where the analytic theory predictions are larger than found in simulation by factors of 2–5. For the volume fraction, the discrepancy is a factor of 1/3–4. However, we do note that although the position of the phase boundaries is not correctly predicted by our theory, their slope is. We note that the slope is affected by the variation of the persistence length with : including this variation in the analytic theory is necessary to reproduce the slope seen in the simulations.

Figure 16: Comparison of order parameter in theory (curves) and simulation (points) in the nematic phase at I-N coexistence. Solid curve: . Dashed curves: . Dotted curves: . In doing the calculations we used the empiricial persistence length found from simulations . Triangles: . Squares: . Circles: .

We also find significant quantitative discrepancy in the order parameter in the nematic phase at coexistence (Figure 16). The theory predicts lower values of than found in the simulations by almost a factor of 2.

Vi Conclusion

Recent experiments on the coupled aggregation and liquid-crystal ordering of chromonic liquid crystals and short double-stranded nucleic acids motivated us to consider coarse-grained simulation and theory of this problem. In nucleic acid systems in particular, the monomer aspect ratio can be varied by varying the number of base pairs per monomer; however, the effects of varying monomer aspect ratio on the aggregation/LC ordering problem have not been extensively studied in theoretical work.

In this paper we used both Monte Carlo simulation and analytic theory to study aggregation and liquid-crystal ordering of a simple model of sticky cylinders. The model assumes hard cylinders (of length and diameter , so the aspect ratio of the monomers can easily be varied) that experience attractive interactions when they stack end to end.

We determined approximate phase diagrams by using calculations of the equation of state in expansion and compression runs; we then studied the aggregate length distributions, order parameter, and correlation functions in the different phases. We find isotropic, nematic, columnar LC, columnar crystal, and cubatic-like phases are possible for this system. This family of phases is similar to the results of Hentschke and colleagues Fodi and Hentschke (2000); Ouyang and Hentschke (2007) but does show differences with the work of Lü and Kindt, who only observed isotropic and nematic phases Lu and Kindt (2004). Similarly to previous work we observe the disappearance of the nematic phase (as the monomer binding energy drops or temperature increases) at an isotropic-nematic-columnar triple point Taylor and Herzfeld (1991); Fodi and Hentschke (2000); Ouyang and Hentschke (2007); Park et al. (2008); Edwards et al. (2008). These phase diagrams bear a strong qualitative resemblance to those of chromonic liquid crystals, which show a strong decrease in nematic phase width with increasing temperature Park et al. (2008); Edwards et al. (2008). The location of the I-N-C triple point appears to be sensitive to the monomer aspect ratio, occurring at lower for larger .

In the simulations we studied three aspect ratios: , 1, and 2. The same qualitative phase behavior is present for all aspect ratios. The key differences are that, first, we only observe the cubatic-like phase for . This is consistent with previous work on hard cylinders Blaak et al. (1999). Second, we do not observe the columnar crystal phase for the lowest aspect ratio of 0.5. Naïvely this might be understood because having columns made of a larger number of shorter monomers makes the increase in entropy allowed by disorder along the columns more favorable. Third, as noted above, the nematic phase persists to lower when the aspect ratio is larger. Qualitatively, this makes sense. We find that the mean aggregate length at I-N coexistence increases with increasing monomer aspect ratio, which suggests (as in Onsager theory) that the density at coexistence should also decrease with increasing monomer aspect ratio. We do observe this decrease, although the trend is not extremely strong. We also find that the persistence length increases linearly with increasing monomer aspect ratio; therefore aggregate stiffness increases with increasing aspect ratio (assuming fixed binding energy). Both increasing aggregate length and increasing aggregate stiffness tend to favor nematic order, and therefore tend to broaden the nematic region.

Our results give evidence that the phase diagram depends on the the form of the interaction potential, not just the binding energy between monomers. The phase behavior is sensitive to the flexibility of the aggregates, and the aggregate persistence length varies with temperature in a way that is controlled by the interaction potential. This suggests that the temperature dependence of the I-N phase boundary as well as the location of the I-N-C triple point are influenced by the implicit dependence of on . This is an interesting direction for future work.

Our simulation results demonstrate a remarkable difference between the correlation length along aggregates and the length of aggregates in the nematic and columnar LC phases. Indeed, the aggregate length is typically 10-100 times larger than the density correlation length along the columns. This has important implications for attempts to measure aggregate length distributions: X-ray scattering experiments, which measure correlation lengths, may not be able to probe aggregate lengths Park et al. (2008); Joshi et al. (2009).

In our computations the aggregate persistence length is determined from the interaction potential of the monomers; we found a universal (independent of ) linear dependence of on .

Our extensions of previous theoretical work suggest that quantitative agreement between analytic theory and simulation for the I-N transition of aggregates is within reach. Starting with the relatively simple second-virial model of aggregation and nematic order of van der Schoot and Cates van der Schoot and Cates (1994), we added (1) the Parsons-Lee approximation for higher-density systems, and (2) expansion of the angular trial function in Legendre polynomials and numerical solution of the resulting algebraic equations. This combination significantly improves the quantitative agreement between the analytic theory and the simulations of Lü and Kindt Lu and Kindt (2004).

Our analytic theory predicts the biexponential distribution of aggregate lengths in the nematic phase. Previous analytic work of Lü and Kindt assumed a biexponential distribution in the nematic phase, and based on this assumption found good agreement with their simulation results Lu and Kindt (2006). Here the biexponential distribution arises directly from the free-energy minimization and does not have to be added to the theory as an assumption.

The comparison between simulation and theory highlights the role of aggregate flexibility in controling the phase diagram. In our analytic theory, the aggregates are treated as nearly rigid. Therefore, our theoretical results match the simulation results best where the aggregate persistence length is long (). In this regime, the quantitative agreement between theory and simulation is remarkably good, particularly in the prediction of the phase boundaries, the order parameter, and the biexponential aggregate length distribution in the nematic phase. When comparing to the simulations of Lü and Kindt with shorter aggregate persistence length () or our simulations (), the quantitative agreement between analytic theory and simulation decreases. Future improvement to the theory will require improved treatment of aggregate flexibility.

A major direction for future research is the comparison of our results to experimental results on chromonic and nucleic-acid liquid crystals, and refinement of both the simulation model and the analytic theory based on the comparison. In the future, improvements to the simulation model and the analytic theory may allow true quantitative agreement between experiment, simulation, and theory.

Appendix A Free energy minimization in the nematic phase

Here we outline the free energy minimization and numerical solution for the function introduced in Eq. (18). Plugging the trial function Eq. (18) into the normalization condition Eq. (5) we arrive at the normalization condition for

(29)

Minimizing the free energy (Eq. (19)) with respect to subject to the normalization constraint gives

(30)

leading to the equation for ,

(31)

Integrating Eq. (A) over and using the normalization condition we obtain the equation for the Lagrange multiplier :

(32)

where the three integrals are

(33)
(34)
(35)

To solve equation (A) and determine the function that minimizes the free energy we seek a series solution: expansion of in Legendre polynomials turns equation (A) into a system of coupled algebraic equations that can straightforwardly be truncated and solved numerically. We expand in even Legendre polynomials (note that the nematic phase is symmetric under , so only even Legendre polynomials are allowed):

(36)

Plugging this series expansion into Eq (A) gives the equations for the coefficients , after using the orthogonality condition and completeness relation for Legendre polynomials:

(37)

The dimensionless parameters and are

(38)
(39)

The coefficients come from the expansion of the interaction kernel

(40)

in Legendre polynomials,

(41)

with and

(42)

The coefficients are

(43)

and the integral of 3 Legendre polynomials is

(44)

This integral is known exactly van der Schoot and Cates (1994): it is equal to 0 when or or . In all other cases

(45)

The three integrals in Eqs. (33)-(35) become, after plugging in the series expansion Eq. (36),

(46)
(47)
(48)

We truncated and solved the system of algebraic equations (32), (37) and (46)-(48) numerically in Matlab. We typically kept 25-30 expansion coefficients for the function , depending on how well the numerical solution converged.

In terms of the expansion coeficients the mean aggregation number in the nematic phase is

(49)

and the order parameter is