Simulation of spherulite growth using a comprehensive approach to modeling the firstorder isotropic/smecticA mesophase transition
Abstract
A comprehensive modeling and simulation study of the firstorder isotropic/smecticA transition is presented and applied to phase diagram computation and twodimensional spherulite growth. An approach based on nonlinear optimization, that incorporates experimental data (from 12CB, dodecylcyanobiphenyl), is used to determine physically realistic model parameters. These parameters are then used in conjunction with an optimized phase diagram computation method. Additionally, a timedependent formulation is presented and applied to the study of twodimensional smecticA spherulite growth. These results show the growth kinematics and defect dynamics of nanoscale smecticA spherulite growth in an isotropic phase with an initially radial layer configuration.
1 Introduction
Liquid crystallinity and other forms of selforganization are key phenomena both technologically and in Nature. Liquid crystalline order ranges from liquids that show some degree of orientational order to those that, in addition, show various types of translational order. The myriad of types of material that exhibit this behavior range from traditional low molecular mass molecules, currently used in display technology, to biological membranes composed of phospholipids [1]. To date, the main focus of liquid crystal research has been on the simplest of the class of mesophases where some degree of orientational order is present: the nematics. Research on mesophases that also show a degree of translational order, including smectics and columnar liquid crystals, has been less abundant. Recognizing the increasing importance of these mesophases, particularly in biological systems, the need for practical methods to access the time and length scales at which these phenomena occur becomes more important.
Experimental work in this general field has made great progress in the basic understanding of translationally ordered liquid crystals [2, 3, 4, 5, 6]. Nonetheless, it is currently infeasible to access much of the dynamic phenomena of translational phaseordering processes. Recent experimental work has begun to address these issues [4, 5], but theoretical approaches are currently the only way to access the length scales (nanometers) and time scales (nanoseconds) at which liquid crystal dynamics occur. The use of highorder models in conjunction with advanced numerical simulation techniques has shown a great deal of promise for theoretical study [7, 8, 9, 10]. Recent computational advances have allowed for the possibility of simulation in greater detail than ever before.
Utilizing a highorder Landaude Gennes type model of the firstorder isotropic/smecticA mesophase transition [11, 7], the objectives of this work are:

to present a comprehensive approach to modeling and simulation of the firstorder isotropic/smecticA transition.

to determine phenomenological model parameters through incorporation of experimental data (from 12CB, dodecylcyanobiphenyl).

to efficiently compute a phase diagram predicted by the model and parameter set.

to study the twodimensional growth kinematics and defect dynamics of an initially radial smecticA spherulite in an isotropic matrix.
This approach builds upon previous work [12], which incorporates experimental data into the phenomenological model and derives equations for phase diagram computation. A timedependent formulation [9, 10] and the nanoscale growth of an initially radial spherulite are presented. This work is organized as follows: a brief background on relevant types of liquid crystalline order is given (Section 2.1), the model and simulation approach are explained (Sections 2.32.6), and simulation results are presented and discussed (Section 3).
2 Background and theory
2.1 Liquid crystalline order
Liquid crystalline phases or mesophases are materials which exhibit partial orientational and/or translational order. They are composed of anisotropic molecules which can be disclike (discotic) or rodlike (calamitic) in shape. Thermotropic liquid crystals are typically purecomponent compounds that exhibit mesophase ordering most greatly in response to temperature changes. Lyotropic liquid crystals are mixtures of mesogens (molecules which exhibit some form of liquid crystallinity), possibly with a solvent, that most greatly exhibit mesophase behavior in response to concentration changes. Effects of pressure and external fields also influence mesophase behavior. This work focuses the study of calamitic thermotropic liquid crystals which exhibit a firstorder mesophase transition.
An unordered liquid, where there is neither orientational nor translational order (apart from an average intermolecular separation distance) of the molecules, is referred to as isotropic. Liquid crystalline order involves partial orientational order (nematics) and, additionally, partial translational order (smectics and columnar mesophases). The simplest of the smectics is the smecticA mesophase, which exhibits onedimensional translational order in the direction of the preferred molecular orientational axis. It can be thought of as layers of twodimensional fluids stacked upon each other. Schematic representation of these different types of ordering are shown in Figure 1.
2.2 Order parameters and the phenomenological model
Theoretical characterization of mesophase order is accomplished using order parameters that adequately capture the physics involved. These order parameters typically have an amplitude and phase associated with them. In order to characterize the partial orientational order of the nematic phase, a second order symmetric traceless tensor can be used [11]:
(1) 
where are the eigenvectors of Qtensor, which characterize the average molecular orientational axes, and are scalars which represent the extent to which the molecules conform to the average orientational axes [13, 14, 15]. Uniaxial order is characterized by and , which correspond to the maximum eigenvalue(and its corresponding eigenvector) of the Qtensor, . Biaxial order is characterized by and , which correspond to the lesser eigenvalues and eigenvectors, .
The smecticA mesophase has onedimensional translational order in addition to the orientational order found in nematics. Characterizing this mesophase can be accomplished through the use of primary (orientational) and secondary (translational) order parameters together [16]. This is accomplished using the tensor order parameter (1) and the complex order parameter [11]:
(2) 
where is the phase, is the scalar amplitude of the density modulation. The density wave vector, which describes the average orientation of the smecticA density modulation, is defined as . The smectic scalar order parameter characterizes the magnitude of the density modulation, and is used in a dimensionless form in this work. In the smecticA mesophase the preferred orientation of the wave vector is parallel to the average molecular orientational axis, .
A Landaude Gennes type model for the first order isotropic/smecticA phase transition is used that was initially presented by Mukherjee, Pleiner, and Brand [11, 7] and later extended by adding nematic elastic terms [17, 18, 8]:
(3)  
where is the free energy density, is the free energy density of the isotropic phase, terms 15 are the bulk contributions to the free energy, terms 67 are couplings of nematic and smectic order; both the bulk order and coupling of the nematic director and smectic densitywave vector, respectively. Terms 810 are the nematic and smectic elastic contributions to the free energy, respectively. The order parameters are defined in (12), is temperature, / are the hypothetical second order transition temperatures for isotropic/nematic and isotropic/smecticA mesophase transitions (refer to [19] for more detail), and the remaining constants are phenomenological parameters. Further explanation and justification for the use of this highorder model can be found in [10].
2.3 Homogeneous free energy and parameter determination
Following past work [12], in order to compute the phase diagram from the free energy equation (3) a homogeneous uniaxiallyordered volume assumption can be used, resulting in a simplified pointvolume free energy density:
(4) 
where is the magnitude of the wave vector and is the equilibrium layer spacing. Note that a different definition of the Qtensor (1) is used in this work than in ref. [12].
A major complication of applying this model is the determination of a suitable set of model parameters. In order to overcome this challenge, a nonlinear programming formulation is derived which then allows to the application of nonlinear solution methods. Utilizing experimental data [19, 20, 21] and minimization criteria of the homogeneous free energy (4) [7, 22], a nonlinear optimization problem is formulated:
(5) 
where the objective function minimizes the change in the layer spacing between the unknown bulk transition value, , and the known value at some minimum valid temperature for the model, . Each constraint is evaluated at the corresponding temperatures , the bulk transition temperature, and , the minimum valid temperature for the model parameter set.
2.4 Phase diagram determination
Utilizing the homogeneous free energy (4) and model parameters determined from the nonlinear optimization solution in Section 2.3, a phase diagram for the system can be computed. Figure 2 shows a schematic of the phase diagram of a firstorder isotropic/smecticA phase transition. Coexistence regions are present, due to the firstorder nature of the transition, where both the isotropic and smecticA phase are either stable or metastable.
Minimization of the homogeneous free energy (4) is computationally challenging due to the presence of three degrees of freedom. This can be simplified by parameterization of the free energy equation [7] using minimization invariants and the assumption that the smecticA phase exists. The resulting free energy equation, parameterized as a function of the nematic scalar order parameter (see (1)), is:
(6) 
where is the nematic scalar order parameter with the assumption that the smecticA phase is stable/metastable. The minima of (2.4) are easily found in that they are the roots of a polynomial, . Once the minima are determined, the validity of the assumption that the smecticA phase exists (at the specified temperature) can be tested by i) verifying that the computed and the corresponding values of / are real, and ii) verification of the minimization criteria of the full homogeneous free energy (4). The equations to determine and from are derived from the application of the minimization criteria to the homogeneous free energy equation (4):
(7) 
2.5 Timedependent formulation
The LandauGinzburg timedependent formulation[23] is used to capture the kinetics of the phase transition. Due to the higher order derivative term in the free energy functional, a higher order functional derivative must be used. Additionally, in order to utilize standard numerical solution techniques, the complex order parameter (2) is separated into its real and imaginary contributions:
(8) 
The general form of the timedependent formulation is as follows [23]:
(9)  
(10) 
where / is the rotational/smectic viscosity, is the heterogeneous free energy density (3), and the volume. As previously mentioned, a higher order functional derivative must be used due to the secondderivative term in the free energy (3):
(11) 
where corresponds to the order parameter.
The use of the full Qtensor in this timedependent model does not neglect biaxiality as was assumed in Sections 2.32.4. Substituting (8), the free energy (3), and high order functional derivative (11) into the timedependent formulation (9) yields the closed set of simulated equations:
where the asterisk denotes an nondimensionalized value, the superscript denotes the symmetric/traceless portion of a tensor, and is the ratio of the smectic and rotational viscosities. The nondimensionalized model parameters are as follows:
(13) 
where is the simulationspecific imposed length scale.
2.6 Simulation conditions and computational approach
A square computational domain with imposed length scales of (approximately 25 layers) and (approximately 50 layers) were used in two separate simulations. Referring to Figure 3a, both symmetry and Neumann boundary conditions were used to simulate bulk conditions. Symmetry conditions for the Qtensor (1) must take into account vector symmetry, which results in the following boundary conditions [24]:
(14) 
The initial condition for both simulations was a smecticA spherulite in an initially radial layer configuration (see Figure 3b). The radius of the spherulite was initially set to , a value on the order of the layer spacing at (Section 2.4), [20]. The initial value used for , , and the layer spacing correspond to the homogeneous values at , determined from the computed phase diagram. The Heaviside step function was used to generate the initial spherulite. The constraint that the spherulite does not impinge on the domain boundaries was verified postsimulation.
A commercial package, Comsol Multiphysics, was used to solve the timedependent model (2.5). Quadratic Lagrange basis functions were used for the Qtensor variables and quartic Hermite basis functions used for the complex order parameter components. Standard numerical techniques were utilized to ensure convergence and stability of the solution. This platform does not support adaptive mesh refinement, thus a uniform mesh was used with a density of approximately nodes/ for the 25 layer simulation and nodes/ the 50 layer simulation. Previous simulations using this model and numerical method have shown good agreement with both past experimental and theoretical findings [9, 10]. Additionally, exhaustive past work using this numerical method and the Landaude Gennes model for the firstorder isotropic/nematic phase transition [25, 26, 27, 28] has served to further validate this simulation approach.
3 Results and discussion
3.1 Phase diagram computation
The algorithm for phase diagram computation presented in Section 2.4 was implemented using a highlevel numerical programming language [29]. The resulting phase diagram is shown in Figure 4; refer to the phase diagram schematic, Figure 2, for a detailed explanation of the features.
The experimentally observed bulk transition temperature and nematic scalar order parameter [19] are exactly matched using the parameter determination method presented in Section 2.3. Additionally, the smectic layer spacing on the order of [20] is also well reproduced below the lower stability limit of the isotropic phase, . Above this temperature the model predicts a layer spacing trend that increases approaching the superheated stability limit of the smecticA phase.
3.2 Twodimensional spherulite growth
Simulation results from the growth of an initially radial spherulite are shown in Figure 5 from the 50 layer simulation (see Section 2.6). Past work on the isotropic/nematic mesophase transition [14, 30, 25, 26, 27, 28, 24] and smecticA filament growth [31] provide a great deal of insight into the two general growth processes observed: shape dynamics [15] and selfsimilar growth [32]. As will be shown, these simulation results definitely show that the growth of an initially radial textured smecticA spherulite in an isotropic phase follows a similar topological growth process as observed in its nematic counterpart [25, 26, 27, 28].
Shape dynamic growth
The shape dynamic growth regime is dominated by transient texturing and interfacial forces. When the spherulite radius is on the order of the smectic coherence length:
(15) 
longrange energy effects (gradients in molecular and layer orientation) are dominant over shortrange energy effects (gradients in the bulk order parameters and ). As the spherulite radius increases, interfacial anchoring affects the core of the spherulite less and shortrange energy gradients become dominant in this region. This results in increased texturing and defect dynamics [13, 14]. Referring to Figure 6, this process is observed for an initially radial smecticA spherulite.
The shape dynamic growth regime of this spherulite morphology has been found to have topologically similar dynamics compared to the growth of initially radial nematic spherulites [25, 26, 27, 28]. In agreement similar nematic growth simulations [25, 26, 27, 28], the initially imposed isotropic core is observed to undergo a selfselected transition to a disclination, shown in Figure 7a. This initial morphology minimizes elastic longrange energy through the formation of a highenergy defect. As the spherulite radius exceeds the smectic coherence length (15), longrange elastic gradients become energetically favorable compared to gradients in shortrange order. Subsequently, the disclination splits into a pair of disclinations [33, 14]. This is a topologically equivalent texture that minimizes shortrange energy, where the energy of a defect is proportional to its the square of its strength [13]. The smecticA layer configurations corresponding with this splitting process is shown in Figure 6.
The 25 layer simulation, with a more refined mesh, was used to both verify mesh independence of the 50 layer simulation (see Section 3.2.2) and obtain a more refined view of the smectic disclinations observed in this work. In order to identify disclination defects, the degree of biaxiality (see Section 2.2) can be computed as follows [34, 35]:
(16) 
The core of the disclination is shown in Figure 7a as was determined from the 25 layer simulation. The simulated structure of this disclination is both topologically correct and its biaxial halo/uniaxial center structure wellreproduced, compared to its nematic counterpart [36, 37, 38, 35, 39]. The spherulite length scale of the 25 layer simulation was not large enough to observe the splitting of the disclination into a pair of , thus Figure 7b shows the partially biaxial core structure of the computed in the 50 layer simulation. The core of this defect is also in agreement with its nematic counterpart, with a solely biaxial core [33].
Selfsimilar growth
The radius of the spherulite versus time was determined for both simulations. Figure 8a shows this data and a comparison of the results from the 50 layer and more refined 25 layer simulations show good agreement. Deviation in the early stages of growth result from the more refined 25 layer simulation capturing the formation and core structure of the disclination more accurately. The convergence of the evolution of the two spherulite radii establishes the accuracy of the 50 layer simulation results beyond the radius of convergence.
Following the shape dynamic growth regime, the spherulite forms a scaleindependent, or selfsimilar, shape [32]. Spherulite growth obeys a power law in this growth regime, where is both experimentally observed [40, 41] and theoretically predicted [26, 42] to approach for phaseordering processes quenched below the lower stability limit of the isotropic phase. A power law fit of the final two points of the spherulite radius evolution (Figure 8, 50 layer simulation) shows that the simulation results in this work correctly predict a selfsimilar growth regime, where is found. Additionally, the results from Figure 8 show that the transition from the initial shape dynamics regime to the selfsimilar growth regime occurs at a relatively small length scale compared to the growth of initially radial textured nematic spherulites [26].
Figure 9 shows the order parameter profiles during the beginning of the selfsimilar growth regime (Figure 5iii). Good agreement is observed between the bulk values of the order parameters (determined in Section 3.1) and those found in the dynamic simulation. The imposition of layer compression/expansion and curvature from the spherulite morphology and interfacial anchoring result in a decrease of the smecticA order from that of the ideal homogeneous layer configuration. Additionally, a local increase in smecticA order is observed in the region immediately outside of the spherulite core where layer compression/expansion is minimized by the presence of the defect structure within the core (See Figure 6iii). This phenomenon is an example of the competition between gradients in long and shortrange order. The existence of such high smecticA ordering in the presence of substantial layer curvature is in agreement with both theoretical and experimental observations that layer bending is a low energy distortion in the smecticA mesophase [11, 43].
4 Conclusions
A comprehensive approach to modeling the firstorder isotropic/smecticA phase transition was presented and applied to phase diagram computation and growth of an initially radial smecticA spherulite in an isotropic phase. Summarizing the results determined from this work are as follows:

Shape dynamics in the early spherulite growth period were found in agreement with the past theoretical work on a similar system, the firstorder isotropic/nematic transition [14, 30, 25, 26, 27, 28, 24], where the initial imposed isotropic core forms a selfselected disclination. As the spherulite radius increases, this disclination splits into two disclinations, which is governed by a competition of long and shortrange ordering in the presence of interfacial anchoring [13, 14] (Figures 57).

A selfsimilar spherulite growth regime was observed following the initial shape dynamic growth regime where a power law fit of the spherulite radius was shown to approach (Figure 8), in agreement with past experimental [40, 41] and theoretical [26, 42] studies of mesophase spherulite growth under a “deep” quench.

The general structure of both and smectic disclinations were found in agreement with studies of their nematic counterparts [33, 36, 37, 38, 35, 39]. Note that due to computational constraints, the mesh density used in this simulation was adequate only to the extent of resolving the general structure of the defect core and not resolution with full detail.
The simulation results presented show that the use of a highorder phenomenological model and experimentally based model parameters results in a substantially more complete reproduction of the physical smecticA system. Though current computational resources restrict this initial work to twodimensions, these results show a strong correlation to past experimental [33, 36, 37, 38, 35, 39, 40, 41] and theoretical observations [14, 30, 25, 26, 27, 28, 42, 24]. These promising results show that threedimensional simulation of this model could be used to study the formation and dynamics of smecticA spherulites at length and time scales currently inaccessible via experimental study.
Acknowledgments
This work was supported, in part, by a grant from the Natural Science and Engineering Council of Canada.
References
 Fisch, M. R. Liquid Crystals, Laptops and Life, volume 23 of World Scientific Series in Contemporary Chemical Physics. World Scientific, (2004).
 Harrison, C., Adamson, D. H., Cheng, Z., Sebastian, J. M., Sethuraman, S., Huse, D. A., Register, R. A., and Chaikin, P. M. Science 290(5496), 1558–1560 (2000).
 Harrison, C., Cheng, Z., Sethuraman, S., Huse, D. A., Chaikin, P. M., Vega, D. A., Sebastian, J. M., Register, R. A., and Adamson, D. H. Phys. Rev. E 66(1), 011706 Jul (2002).
 Ruiz, R., Sandstrom, R., and Black, C. Advanced Materials 19(4), 587–591 (2007).
 Michel, J.P., Lacaze, E., Goldmann, M., Gailhanou, M., de Boissieu, M., and Alba, M. Physical Review Letters 96(2), 027803 (2006).
 Kleman, M. and Friedel, J. Reviews of Modern Physics 80(1), 61 (2008).
 Mukherjee, P. K., Pleiner, H., and Brand, H. R. European Physical Journal E: Soft Matter 4, 293–297 (2001).
 Biscari, P., Calderer, M., and Terentjev, E. Phys Rev E Stat Nonlin Soft Matter Phys 75(5), 051707 May (2007).
 Abukhdeir, N. and Rey, A. Solid State Phenomena 139, 135–140 (2008).
 Abukhdeir, N. and Rey, A. New Journal of Physics 10(6), 063025 (17pp) (2008).
 de Gennes, P. and Prost, J. The Physics of Liquid Crystals. Oxford University Press, New York, second edition, (1995).
 Abukhdeir, N. and Rey, A. In Modelling and Simulation. International Association of Science and Technology for Development (IASTED), ACTA Press, (2007).
 Rey, A. and Denn, M. Annual Review of Fluid Mechanics 34(1), p233 – (2002).
 Yan, J. and Rey, A. D. Phys. Rev. E 65(3), 031713 Feb (2002).
 Rey, A. D. Soft Matter 3, 1349 – 1368 (2007).
 Toledano, J.C. and Toledano, P. The Landau Theory of Phase Transitions: Application to Structural, Incommensurate, Magnetic, and Liquid Crystal Systems (World Scientific Lecture Notes in Physics). World Scientific Pub Co Inc, (1987).
 Brand, H. R., Mukherjee, P. K., and Pleiner, H. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics 63, 061708/1–061708/6 (2001).
 Mukherjee, P. K., Pleiner, H., and Brand, H. R. The Journal of Chemical Physics 117(16), 7788–7792 (2002).
 Coles, H. J. and Strazielle, C. Molecular Crystals and Liquid Crystals 55, 237–50 (1979).
 Urban, S., Przedmojski, J., and Czub, J. Liquid Crystals 32, 619–624 (2005).
 DasGupta, S. and Roy, S. Physics Letters A 306, 235–242(8) (2003).
 Prostakov, A. E. Crystallography Reports (Translation of Kristallografiya) 47, 856–858 (2002).
 Barbero, G. and Evangelista, L. R. An Elementary Course on the Continuum Theory for Nematic Liquid Crystals (Series on Liquid Crystals , Vol 3). World Scientific Publishing Company, (2000).
 Abukhdeir, N., Soule, E., and Rey, A. Arxiv preprint arXiv:0805.3539 (2008).
 Wincure, B. and Rey, A. D. The Journal of Chemical Physics 124(24), 244902 (2006).
 Wincure, B. and Rey, A. Continuum Mechanics and Thermodynamics 19(1), 37–58 June (2007).
 Wincure, B. and Rey, A. Nano Letters 7(6), 1474–1479 (2007).
 Wincure, B. and Rey, A. D. Liquid Crystals 34(12), 1397–1413 (2007).
 Eaton, J. Gnu Octave Manual. Network Theory Ltd., (2002).
 Sharma, D. and Rey, A. D. Liquid Crystals 30(3), p377 – (2003).
 Rey, A. D. and Abukhdeir, N. M. Journal of NonNewtonian Fluid Mechanics 153(23), 177–182 August (2008).
 Sethna, J. P. Statistical Mechanics: Entropy, Order Parameters and Complexity. Clarendon Press, (2007).
 Schopohl, N. and Sluckin, T. J. Phys. Rev. Lett. 59(22), 2582–2584 Nov (1987).
 Kaiser, P., Wiese, W., and Hess, S. J. NonEquilib. Thermodyn 17, 153–169 (1992).
 Kralj, S., Virga, E. G., and Zumer, S. Phys. Rev. E 60(2), 1858–1866 Aug (1999).
 Hudson, S. D. and Larson, R. G. Phys. Rev. Lett. 70(19), 2916–2919 May (1993).
 Sonnet, A., Kilian, A., and Hess, S. Phys. Rev. E 52(1), 718–722 Jul (1995).
 Sigillo, I., Greco, F., and Marrucci, G. Liquid Crystals 24(3), p419 – 425 (1998).
 Andrienko, D. and Allen, M. P. Phys. Rev. E 61(1), 504–510 Jan (2000).
 Dierking, I. Appl. Phys. A: Mater. Sci. Process. 72, 307–310 (2001).
 Dierking, I. and Russell, C. Physica B (Amsterdam, Neth.) 325, 281–286 (2003).
 Huisman, B. A. H. and Fasolino, A. Physical Review E (Statistical, Nonlinear, and Soft Matter Physics) 76(2), 021706 (2007).
 Kleman, M. and Lavrentovich, O. Soft Matter Physics: An Introduction. Springer, New York, first edition, (2003).