Biaxial nematics of hard cuboids in an external field

Abstract

By computer simulation, we model the phase behaviour of colloidal suspensions of board-like particles under the effect of an external field and assess the still disputed occurrence of the biaxial nematic (N_B) liquid crystal phase. The external field promotes the rearrangement of the initial isotropic (I) or uniaxial nematic (N_U) phase and the formation of the N_B phase. In particular, very weak field strengths are sufficient to spark a direct I–N_B or N_U–N_B phase transition at the self-dual shape, where prolate and oblate particle geometries fuse into one. By contrast, forming the N_B phase at any other geometry requires stronger fields and thus reduces the energy efficiency of the phase transformation. Our simulation results show that self-dual shaped board-like particles with moderate anisotropy are able to form N_B liquid crystals under the effect of a surprisingly weak external stimulus and suggest a path to exploit low-energy uniaxial-to-biaxial order switching.

---

It is well established that anisotropic colloidal particles can self-assemble into a plethora of fascinating liquid crystal (LC) phases in a solvent. Onsager showed that mere excluded volume effects can force hard-core particles to align along a common director at sufficiently large concentrations.¹ The resulting LC phases found at the thermodynamic equilibrium and their structural properties strongly depend on the particle geometry. In particular, prolate particles tend to orient along their major axis, while oblate particles along their minor axis. Although this tendency is regularly observed in systems of uniaxial particles, such as disks, whose orientation is determined by a single unit vector, it is less predictable for biaxial particles, such as cuboids, whose orientation is fully determined by two unit vectors. For instance, slightly oblate hard board-like particles (HBPs) have been shown to orient along their major axis and thus form prolate (rather than oblate) smectic LC phases.² This unusual arrangement was observed in suspensions of HBPs with length-to-thickness ratio L* ≡ L/T = 12 and width-to-thickness ratio , a geometry where oblate and prolate shapes fuse into one.

Such an exclusive particle architecture, referred to as self-dual shape, was predicted to favour the formation of the biaxial nematic (N_B) phase in systems of HBPs with rounded³ or square⁴ corners. Nevertheless, these theoretical predictions were made within the context of the Zwanzig model, which does not allow free rotations of particles and restricts their orientations to only six. Computer simulations that explored the phase behaviour of freely rotating HBPs highlighted the challenge of observing stable N_B phases, even when an element of size-dispersity is incorporated.⁵ The very recent and insightful simulation study by the Dijkstra's group showed that monodisperse HBPs, with a geometry close or equal to the self-dual shape, are able to stabilise the N_B phase only if their anisotropy is significantly large, with L* ≥ 23.⁶ Nevertheless, the highly uniform colloidal board-like particles synthesised by Nie and coworkers were not observed to assemble into the N_B phase in a wide spectrum of aspect ratios, between L* = 20 and 180, and very close to the self-dual shape.⁷ A direct observation of a stable N_B phase was reported almost ten years ago in dispersions of purely repulsive and quasi dual-shaped goethite particles,⁸ whose stability was shown to depend on the extremely large particle size dispersity.⁴ These computational and experimental findings unambiguously indicate that self-dual shaped HBPs cannot self-assemble into an N_B phase, unless their anisotropy and/or size dispersity are especially relevant. Vroege and coworkers investigated the effect of a magnetic field on the phase behaviour of quasi-dual-shaped polydisperse goethite particles in suspension and observed a biaxial-to-uniaxial nematic transition above a certain magnetic field strength.⁹ To explain the origin of these experimental results, the same authors formulated a mean-field theory to calculate the phase diagram of HBPs in the presence of a magnetic field. Despite the good qualitative agreement between theory and experiments, the former neglects the effect of size polydispersity, being crucial to stabilise the N_B phase,⁸ and is limited by the strong approximations imposed by the Zwanzig model, which was shown to be not especially accurate to describe the phase behaviour of HBPs.^2,5

In this work, we perform Monte Carlo (MC) simulations of freely-rotating monodisperse HBPs in the presence of an external field that promotes a phase transition from the isotropic (I) or uniaxial nematic (N_U) phase to the N_B phase. Because our main interest is gaining an insight into the energetics associated to this process, the field strength is a simulation parameter that assumes values between ε_f* ≡ ε_fβ = 0.1 and 3, with β the inverse temperature. At the same time, we propose an alternative route for the formation of the N_B phase that does not necessarily imply extreme anisotropies⁶ or particle size dispersities,⁸ nor the addition of depletants.¹⁰ In the following, we discuss the main aspects of the model and simulation methodology applied. The interested reader is referred to ESI† and our previous works for additional details.^2,5

We have simulated fluids of HBPs with L* = 12 and 1 ≤ W* ≤ 12, including , which gives the above mentioned self-dual shape. The particle thickness, T, is our system unit length and is kept constant. The unit vectors, x̂, ŷ, and ẑ, aligned along W, T, and L, respectively, define the particle orientation in space. Interactions are described by a hard-core potential, which only depends on the distance between pairs of particles. Consequently, an MC move is rejected if an overlap occurs and accepted otherwise. The separating axes method, introduced by Gottschalk et al.¹¹ and later modified to study suspensions of tetragonal parallelepipeds,¹² has been employed to assess the occurrence of overlaps between HBPs. The interaction between HBPs and external field is described by a potential energy term that favours the alignment of the particle intermediate axis, x̂, with the direction of field, indicated by ê, and reads

where N = 2000 indicates the number of HBPs in the system. The field direction, ê, is always aligned with the same box axis, and the field itself coupled to the particle axis x̂ in all the systems studied. The latter could have been alternatively coupled to ŷ or ẑ, but this choice would have promoted the formation of, respectively, oblate (N_U⁻) or prolate (N_U⁺) uniaxial nematics, which are normally observed without the application of an external field.² We stress that there is no intention here to mimic a real electric field, whose effect on a colloidal suspension depends on particle polarizability, charge and dielectric constant difference with that of the solvent. The effective electric field perceived by a particle would be the sum of the external field and the field due to the induced and spontaneous dipoles. Dipole–dipole interactions are often large and cannot trivially be ignored as they would determine the stability of the phases observed. In the experimental work on the formation of the N_B phase in suspensions of goethite particles, the driving external field was magnetic rather than electric.⁸ In this situation, one can assume that the magnetic coupling between neighbouring particles is small^9,13 and thus the interactions between a particle and the external field will be the predominant effect. Rather than on the physical nature of the external stimulus, which we keep as simple as possible, we focus on its resulting effect, that is the ability of reorienting and aligning HBPs as well as the formation of uniaxial and biaxial LC phases. How electric and magnetic fields can be employed to control the internal organisation and self-assembly of colloidal suspensions has been explored in recent theoretical and simulation studies on spherical^14–17 and anisotropic^18–21 particles, including HBPs within the Zwanzig model approximation.⁹ The field-particle potential in eqn (1) can be regarded as an effective potential that promotes the alignment of the particle intermediate axis with ê. Therefore, we assume the inter-particle interactions to be predominantly governed by excluded-volume effects, ignoring any enthalpic contribution that could arise from the presence of the field.

The external field has been applied to suspensions of HBPs that were either in dense I phases, very close to the border with the nematics, or deeply in the N_U⁻ or N_U⁺ phases, at the packing fraction η ≡ Nv/V = 0.34, with v and V the particle and box volume, respectively. At this packing fraction, the nematic phase is stable for the whole range of particle geometries studied here. In particular, N_U⁺ and N_U⁻ phases are observed, respectively, at 1 ≤ W* < 3.46 and 3.46 < W* ≤ 12, as shown in the phase diagram available in ref. 2. By contrast, the packing fraction of the systems in the I phase changes with W*, between η = 0.20 and 0.31. All simulations have been run in the canonical ensemble, where number of particles, volume and temperature are kept constant. After equilibrating a number of I and N_U phases, the field was switched on and new equilibrium states were achieved.

To identify the long-range order of the phases at the initial and final equilibrium states, we calculated the nematic order parameter and nematic director coupled to each of the three particle axes. In particular, we applied the standard procedure of diagonalising the following traceless symmetric second-rank tensor²²

where i indicates a generic HBP, [small klambda, Greek, circumflex] = x̂, ŷ, ẑ is the particle unit orientation vector and I the second-rank unit tensor. The resulting eigenvalues (S_2,W, S_2,T, S_2,L) and associated eigenvectors (m̂, p̂, n̂) identify, respectively, the order parameters and nematic directors coupled to each particle axis. If the order is very weak, as in the I phase, the eigenvalues vanish and the eigenvectors are meaningless. However, if a preferential direction of alignment exists, this is identified by the largest positive eigenvalue of Q^λλ and its corresponding eigenvector. For instance, the N_U⁺ phase is characterised by a large value of S_2,L and a preferential alignment of the particle ẑ axes along the director n̂. The three tensors in eqn (2) can be applied to estimate three additional order parameters, B_2,W, B_2,T and B_2,L, coupled, respectively, to m̂, p̂ and n̂, that identify the occurrence of biaxiality.²³ However, it is sufficient to monitor the fluctuations of the unit vectors that are perpendicular to the nematic director associated to the largest eigenvalue of Q^λλ.^24–26 If S_2,L is the largest eigenvalue, then these fluctuations are quantified by the biaxial order parameter B_2,L = (m̂·Q^xx·m̂ + p̂·Q^yy·p̂ − m̂·Q^yy·m̂ − p̂·Q^xx·p̂)/3. In this case, S_2,L = 1 and B_2,L = 0 would indicate perfect alignment along n̂ only and the occurrence of the N_U⁺ phase, while S_2,L = 1 and B_2,L = 1 would indicate perfect alignment along three directors and the formation of the N_B phase. We stress that alternative definitions of order parameters, which explicitly incorporate the particle axes connectivity, are available.²⁷ Finally, to unambiguously discard the occurrence of smectic LC phases, we have calculated the pair distribution functions along the relevant nematic directors, which do not show any significant periodicity as expected from the relatively low packing.

In the light of these preliminary considerations, we now discuss our simulation results, which depend on three main factors: (i) the initial state of the system, being either isotropic or uniaxial nematic, (ii) the field strength, ε_f, and (iii) the particle geometry, which only depends on W*. In Fig. 1 and 2, we show the phases observed, in the new equilibrium state, when a field of strength 0.1 ≤ ε_f* ≤ 3 is applied to, respectively, I or N_U phases of HBPs with 1 ≤ W* ≤ 12. The values of the uniaxial and biaxial nematic order parameters that validate these results are provided in the ESI.† The state points at ε_f = 0 are included as a reference and identify I phases in Fig. 1, and N_U⁺ (W* < 3.46) or N_U⁻ (W* > 3.46) phases in Fig. 2. The criterion adopted in both figures to classify the new equilibrium phases essentially depends on the value of the uniaxial and biaxial order parameters and is summarised in Table 1.

Fig. 1 Equilibrium LCs resulting from the application of the external field given in eqn (1), with strength ε_f, to an isotropic phase of HBPs of length-to-thickness ratio L* = 12 and width-to-thickness ratio 1 ≤ W* ≤ 12. Packing fractions are in between η = 0.20 and 0.31 (see ESI† for details). Circles, empty and solid squares, empty and solid triangles, and diamonds indicate, respectively, I, N_U⁺, N_U⁻, N_B⁺, N_B⁻, and N_B phases. Shaded areas are guides for the eye. The vertical dashed line refers to the self-dual shape, where .

Fig. 2 Equilibrium LCs resulting from the application of the external field given in eqn (1), with strength ε_f, to uniaxial nematic phases (η = 0.34) of HBPs of length-to-thickness ratio L* = 12 and width-to-thickness ratio 1 ≤ W* ≤ 12. Circles, empty and solid squares, empty and solid triangles, and diamonds indicate, respectively, I, N_U⁺, N_U⁻, N_B⁺, N_B⁻, and N_B phases. Shaded areas are guides for the eye. The vertical dashed line refers to the self-dual shape, where .

Table 1 Range of uniaxial and relevant biaxial order parameters for isotropic (I), uniaxial (N_U) and biaxial (N_B) phases observed in suspensions of HBPs. The superscripts ⁺ and ⁻ indicate, respectively, prolate and oblate phase symmetry

S 2,L	S 2,T	S 2,W	B 2,L or B2,T	Phase
a S 2,T and/or S2,W should be in the range specified. b S 2,L and/or S2,W should be in the range specified.
0.00–0.20	0.00–0.20	0.00–0.20	—	I
0.40–1.00	0.00–0.35	0.00–0.35	0.00–0.30	NU^+
0.00–0.35	0.40–1.00	0.00–0.35	0.00–0.30	NU^−
0.40–1.00	0.35–1.00	0.35–1.00	0.20–0.35	NB^+
0.35–1.00	0.40–1.00	0.35–1.00	0.20–0.35	NB^−
0.40–1.00	0.40–1.00	0.40–1.00	0.35–1.00	NB

---

According to this criterion, two weak biaxial phases, similar to those previously found by the Patras group in suspensions of spheroplatelets,^28–30 have been identified. To distinguish weak and strong biaxial phases, we look at the magnitude of the biaxial order parameter that is associated to the largest uniaxial order parameter. If the latter is, for example, S_2,L, then the former is B_2,L. In a strong biaxial phase, simply labelled as N_B, the HBPs are highly orientationally ordered, with their three unit vectors x̂, ŷ and ẑ oriented along three mutually perpendicular directions. This significant degree of orientational order is reflected in the large values of all the uniaxial and biaxial order parameters. By contrast, the weak biaxial phases, labelled as N_B⁺ and N_B⁻, show a modest, but not negligible biaxial order parameter, between 0.20 and 0.35, and at least two pronounced uniaxial order parameters, above 0.40, which give them a prolate or oblate character, respectively.

An external field with strength ε_f* < 1 is able to reorient prolate and oblate HBPs along a common direction to form N_U⁺ and N_U⁻ phases, respectively (see Fig. 1). Interestingly enough, at the self-dual shape, even a very gentle field is sufficient to spark a direct I–N_B⁻ transition and then, at ε_f* = 0.3, the formation of the N_B phase, which becomes dominant at larger field strengths across the complete set of particle geometries. Similar tendencies are noticed when the same field is applied to a uniaxial nematic phase. In this case, however, the I–N_B phase transition is detected at an almost insignificant field strength, ε_f* = 0.1, with no intermediate formation of weak biaxial phases (see Fig. 2). Snapshots showing the phase order before and after the application of a relatively weak field to I and N_U phases are shown in Fig. 3.

Fig. 3 Direct I-to-N_B (top) and N_U⁺-to-N_B (bottom) phase transitions induced by an external field with strength ε_f* = 0.3. Particles have a relative length L* = 12 and relative width W* ≈ 3.46. The packing fraction of the isotropic and uniaxial nematic phase is, respectively, η = 0.307 and 0.340. The vertical arrow indicates the direction ê of the applied field. The particles are colored by RGB code where the red intensity is designed by the scalar product |z·n|, the green intensity by |y·p| and the blue intensity by |x·m|. This means that if the three orientations of the particles are parallel to the nematic director vector associated to each orientation, the particles will have white color. In a perfect biaxial phase, all particles will be white.

Therefore, regardless of the initial system configuration, the dual-shape is found to be the optimal particle geometry to form the N_B phase in terms of energy costs as the strength of the field applied at is at its minimum. In spite of being relatively low, such energy costs are crucial to achieve the desired orientational order with no increase of the translational order, which remains negligible and would otherwise imply the onset of biaxial smectic, rather than nematic, LCs. As a matter of fact, the calculation of the pair distribution functions along the nematic director, g_‖(r), excludes the formation of positionally ordered phases (see ESI†). At particle geometries that are different, but still relatively similar, to the self-dual shape, stabilising the N_B phase becomes more energetically demanding, whether the initial phase is I or N_U. For instance, an I phase incorporating HBPs with W* = 4 requires a field strength of approximately ε_f* = 1. At larger anisotropies, stronger fields are required, up to ε_f* = 1.5 at W* = 12. As expected, the energy to reorient HBPs originally in uniaxial nematics into biaxial nematics is sensibly lower. At W* = 4, this is approximately 50% of that needed to align them in the isotropic phase.

In summary, our results suggest that it is indeed possible to observe biaxial nematics in colloidal suspensions of monodisperse HBPs. To this end, one needs to apply an external field whose strength is comparable to the particle thermal energy. In particular, the self-dual shape is the most appropriate particle geometry for energy-efficient I–N_B and N_U–N_B phase transitions. Although these conclusions agree well with past theoretical works,^3,4,27,31de facto they indicate that this agreement is the consequence of restricting particle orientations and/or neglecting the occurrence of positionally ordered phases, which artificially enhance the stability of the N_B phase. Mere excluded volume effects were shown to be insufficient to observe biaxial nematics in monodisperse systems of HBPs,² unless extremely long particles are employed.⁶ Applying an external field to monodisperse freely-rotating HBPs produces a scenario that the above mentioned theories, despite their strong assumptions, had predicted: the N_B phase can be stabilised at moderate particle anisotropies and a direct I–N_B transition, so far only detected experimentally in systems of significantly size dispersed goethite particles,⁸ observed. The external field induces the desired alignment of particles with no effect on their positional order an at a relatively low energy cost. While promoting the particle alignment along one direction, this field does not prevent them from rotating freely, especially around the remaining two directions. By gradually increasing the field strength, the particles, regardless their geometry, assume a more and more narrow distribution of orientations, which eventually end up fluctuating around three main orthogonal directions, being the signature of the formation of the N_B phase.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

AC acknowledges project P12-FQM-2310 funded by the Junta de Andalucía-FEDER and C3UPO for HPC facilities provided. AC also acknowledges grant PPI1719 awarded by the Pablo de Olavide University for funding his research visit to the School of Chemical Engineering and Analytical Science, the University of Manchester. EMR acknowledges the Malaysian Government Agency Majlis Amanah Rakyat for funding his PhD at the University of Manchester. DC and AP acknowledge financial support from EPSRC under grant agreement EP/N02690X/1.

Notes and references

1. L. Onsager, Ann. N. Y. Acad. Sci., 1949, 51, 627.
2. A. Cuetos, M. Dennison, A. Masters and A. Patti, Soft Matter, 2017, 13, 4720.
3. M. P. Taylor and J. Herzfeld, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 3742–3751.
4. S. Belli, A. Patti, M. Dijkstra and R. van Roij, Phys. Rev. Lett., 2011, 107, 148303.
5. A. Patti and A. Cuetos, Mol. Simul., 2018, 44, 516.
6. S. Dussi, S. Tasios, T. Drwenski, R. van Roij and M. Dijkstra, Phys. Rev. Lett., 2018, 120, 177801.
7. Y. Yang, G. Chen, S. Thanneeru, J. He, K. Liu and Z. Nie, Nat. Commun., 2018, 9, 4513.
8. E. van den Pol, A. V. Petukhov, D. M. E. Thies-Weesie, D. V. Byelov and G. J. Vroege, Phys. Rev. Lett., 2009, 103, 258301.
9. A. B. G. M. Leferink op Reinink, S. Belli, R. van Roij, M. Dijkstra, A. V. Petukhov and G. J. Vroege, Soft Matter, 2014, 10, 446–456.
10. S. Belli, M. Dijkstra and R. van Roij, J. Phys.: Condens. Matter, 2012, 24, 284128.
11. S. Gottschalk, M. C. Lin. and C. Manocha, Comput. Graph., 1996, 30, 171–180.
12. B. S. John and F. A. Escobedo, J. Phys. Chem. B, 2005, 109, 23008–23015.
13. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, 1993.
14. A.-P. Hynninen and M. Dijkstra, Phys. Rev. Lett., 2005, 94, 138303.
15. A.-P. Hynninen and M. Dijkstra, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 051402.
16. F. Smallenburg and M. Dijkstra, J. Chem. Phys., 2010, 132, 204508.
17. F. Smallenburg, H. R. Vutukuri, A. Imhof, A. van Blaaderen and M. Dijkstra, J. Phys.: Condens. Matter, 2012, 24, 464113.
18. B. W. Kwaadgras, M. Verdult, M. Dijkstra and R. van Roij, J. Chem. Phys., 2011, 135, 134105.
19. B. W. Kwaadgras, R. van Roij and M. Dijkstra, J. Chem. Phys., 2014, 140, 154901.
20. T. Troppenz, A. Kuijk, A. Imhof, A. van Blaaderen, M. Dijkstra and R. van Roij, Phys. Chem. Chem. Phys., 2015, 17, 22423–22430.
21. A. Azari, J. J. Crassous, A. M. Mihut, E. Bialik, P. Schurtenberger, J. Stenhammar and P. Linse, Langmuir, 2017, 33, 13834–13840.
22. R. Eppenga and D. Frenkel, Mol. Phys., 1984, 52, 1303–1334.
23. P. J. Camp and M. P. Allen, J. Chem. Phys., 1997, 106, 6681–6688.
24. P. I. C. Teixeira, M. A. Osipov and G. R. Luckhurst, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 061708.
25. M. P. Allen, Liq. Cryst., 1990, 8, 499–511.
26. P. J. Camp, M. P. Allen and A. J. Masters, J. Chem. Phys., 1999, 111, 9871–9881.
27. B. Mulder, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 360–370.
28. D. D. Peroukidis and A. G. Vanakaras, Soft Matter, 2013, 9, 7419–7423.
29. S. D. Peroukidis, A. G. Vanakaras and D. J. Photinos, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062508.
30. S. D. Peroukidis, Soft Matter, 2014, 10, 4199–4207.
31. J. P. Straley, Phys. Rev. A: At., Mol., Opt. Phys., 1974, 10, 1881–1887.

---

Footnote

† Electronic supplementary information (ESI) available: Simulation details; pair distribution functions of a biaxial nematic phase; uniaxial and biaxial order parameters. See DOI: 10.1039/c8sm02283f

---