Self-assembly of Freely-rotating Polydisperse Cuboids: Unveiling the Boundaries of the Biaxial Nematic Phase

Colloidal cuboids have the potential to self-assemble into biaxial liquid crystal phases, which exhibit two independent optical axes. Over the last few decades, several theoretical works predicted the existence of a wide region of the phase diagram where the biaxial nematic phase would be stable, but imposed rather strong constraints on the particle rotational degrees of freedom. In this work, we employ molecular simulation to investigate the impact of size dispersity on the phase behaviour of freely-rotating hard cuboids, here modelled as self-dual-shaped nanoboards. This peculiar anisotropy, exactly in between oblate and prolate geometry, has been proposed as the most appropriate to promote phase biaxiality. We observe that size dispersity radically changes the phase behaviour of monodisperse systems and leads to the formation of the elusive biaxial nematic phase, being found in an large region of the packing fraction vs polydispersity phase diagram. Although our results confirm the tendencies reported in past experimental observations on colloidal dispersions of slightly prolate goethite particles, they cannot reproduce the direct isotropic-to-biaxial nematic phase transition observed in these experiments.

The first known theory on entropy-driven phase transitions was proposed by Onsager in his seminal work dating back to 1940s [1]. Onsager demonstrated that systems of infinitely long hard rods exhibit isotropic-to-nematic phase transition as a result of mere volume effects. Later on, experiments and computer simulations showed that entropy-driven phase transitions can also lead to the formation of positionally ordered liquid crystals (LCs), such as smectic and columnar phases [2,3]. The equilibrium structures stemming from the self-assembly of colloidal particles are especially determined by the architecture of their building blocks. In particular, biaxial particles, such as bent-core and cuboidal particles, have been reported to form biaxial nematic (N B ) LCs [4][5][6]. First theorised by Freiser in 1970 [7], the N B phase has attracted widespread attention for being a promising candidate to be engineered into next generation liquid crystal displays. In contrast to uniaxial nematics (N U ), where long-range orientational order exists only along one direction, the N B phase possesses three orthogonal directors and hence two distinct optical axes that can pave the path to highperformance displays [8][9][10]. Despite having been extensively studied over the past 50 years, the stability of the N B phase still remains an open question, predominantly, but not only, at the molecular scale. Answering this question is hampered by the fact that the N B phase tends to be metastable with respect to other morphologies, such as the N U and Sm phases [11].
About a decade ago, Vroege and coworkers reported on the first experimental evidence of the N B phase in a system of colloidal goethite (roughly board-like) particles [6]. The stability of this N B phase was ascribed to the particles' quasi self-dual shape, a geometry in between oblate and prolate, and to their significant size disper-sity, which hinders the formation of the Sm phase. This key work has reignited recent interest, sparking numerous theoretical, experimental and computer simulation studies on hard board-like particles (HBPs) [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] and other biaxial geometries [24,[28][29][30][31][32]. Theoretical and computational studies on monodisperse systems have suggested that self-dual-shaped particles exhibit a higher tendency to form biaxial nematics. However, these studies applied rather strong approximations, limiting the particles orientation to six orthogonal directions [11][12][13], freezing the rotation of the particle long axes [33], or neglecting the occurrence of positionally ordered LC phases [34][35][36][37]. Our recent theoretical work and computer simulations of freely-rotating HBPs suggested that these approximations might artificially magnify the stability of the N B phase [22,23]. We note that stable N B phases have been found in systems of cuboids with rounded corners (spheroplatelets) with length-to-thickness ratio L * ≡ L/T > 9 [16,17] and in systems of especially elongated HBPs (L * ≥ 23) [24]. Nevertheless, experiments on highly uniform and monodisperse colloidal cuboids with 15 ≤ L * ≤ 180 did not report the formation of the N B phase, which, for self-dual-shaped cuboids, was found to be pre-empted by the biaxial smectic (Sm B ) phase [25]. This lack of agreement between experiments, theory and simulations keeps the discussion on the ability of HBPs to form the N B phase still alive.
Size dispersity has been identified as a key ingredient to destabilise the Sm phase and thus promote the formation of biaxial nematics. In particular, the effect of size dispersity in systems of HBPs was investigated by Onsager's theory within the restricted-orientation (Zwanzig) approximation [13,18]. This theory provided an elegant and solid explanation on the origin of the N B stability ex-perimentally observed in colloidal dispersions of goethite particles [6]. Nevertheless, its conclusions were strongly determined by the use of the Zwanzig model, which only allows six orthogonal particle orientations and cannot describe the phase behaviour of cuboids accurately, as recent simulations and theory have indicated [22,23]. Consequently, fully unlocking the particle rotational degrees of freedom is of paramount importance to ascertain the impact of polydispersity on the phase behaviour of HBPs and accurately map the boundaries of the N B phase. While it is extremely challenging to formulate a theory that simultaneously incorporates particle size dispersity and unrestricted orientations, molecular simulation can provide an insightful contribution to shed light on this combined effect. To this end, we have performed Monte Carlo (MC) simulations of freely-rotating HBPs with Gaussian size distribution peaked at L * = 12. The particle thickness, T , is the system unit length and is the same for all HBPs, whereas L * changes with standard deviation σ L ⟨L * ⟩, where 0.05 ≤ σ L ≤ 0.30 measures the particle length dispersity. Finally, the width-to-thickness ratio, W * ≡ W /T = L * , sets the self-dual shape for all particles. Our systems consist of N p = 2000 to 3000 HBPs that are initially arranged in cubic or rectangular boxes with periodic boundaries and are equilibrated in the isothermal-isobaric ensemble. Phase transitions have been assessed by expansion or compression of a perfect biaxial nematic phase at an extensive range of pressures. To ensure that the equilibrium configurations were independent of the initial configurations, expansion of Sm phases and compression of isotropic (I) phases have also been carried out. We have also simulated significantly larger systems, with N p = 6000, to discard the occurrence of finite-size effects, where these could especially influence the symmetry of the phases observed, that is at the I − N transition. Generally, up to 10 7 MC cycles were needed to equilibrate the systems, with a cycle consisting of N p attempts of displacing and/or rotating randomly selected particles and one trial volume change. Because the force field employed here only consists of a hard-core potential, these moves were accepted if no overlaps were detected, according to the separating axes theorem [38,39]. Systems were considered at equilibrium if packing fraction (η) and uniaxial (S 2 ) and biaxial (B 2 ) order parameters achieved steady values within reasonable statistical uncertainty. In particular, η ≡ ∑ where V is the box volume and v i the volume of a generic particle i. The calculation of S 2 and B 2 was done by diagonalisation of the traceless second-rank symmetric tensor Q λλ = ⟨∑ where I is the second-rank unit tensor,λ i =x,ŷ andẑ is the particle unit orientation vector and angular brackets denotes ensemble average. Isotropic, perfect uniaxial and perfect biaxial phases are observed at (S 2 , B 2 ) = (0, 0), (1, 0) and (1, 1), respectively. Finally, to assess the long-range ordering of the phases at equilibrium, we have analysed the spatial correlations along the relevant phase directors and perpendicularly to them by computing the longitudinal, g ∥ (r ∥ ), and transverse, g ⊥ (r ⊥ ), pair distribution functions, where r ∥ and r ⊥ are the corresponding projections of the inter-particle distance. The interested reader is referred to the ESI † for further details on the calculation of order parameters and pair correlation functions. Bearing in mind these introductory considerations, we now report on the LC phases that polydisperse HBPs are able to form at equilibrium, with specific interest in the critical polydispersity that stabilises the N B phase and the extent of this stability. For consistency with our former work [26], we defined the axial symmetry of N and Sm phases according to the magnitude of the order parameters (see ESI † for additional details). The σ L −η phase diagram of polydisperse and self-dual-shaped HBPs is shown in Fig. 1. All points at σ L = 0 refer to monodisperse systems, obtained from previous simulation results [22] and included here as a reference. The diagram reveals the existence of uniaxial and biaxial LC phases at η ≳ 0.30. At lower density and regardless the size dispersity, we only detect I phases, in good agreement with the Zwanzig-based Onsager theory by Belli et al, who predicted stable I phases for η ≲ 0.28 and 0 ≤ σ L ≤ 0.4 [13]. We stress that these authors studied slightly prolate particles with W * ≈ 1.017 L * or ν ≡ L/W − W /T = 0.1. This slight deviation from the self-dual shape, for which ν = 0, might appear insignificant, but in fact it determines the oblate or prolate symmetry of the nematic phase, at least in monodisperse systems, as established by Mulder in the 1980s [35].
As far as the N U phase is concerned, its stability re-(a) (c) gion (yellow-shaded area in Fig. 1) is not especially influenced by the particle polydispersity, being the lower and upper boundaries, approximately constrained between η = 0.30 and 0.35, similar to those observed in monodisperse systems (σ L = 0). A slight difference is detected at 0.15 ≤ σ L ≤ 0.20, where four different phases seem to converge and the region of stability of uniaxial nematics expands up to η = 0.41 at σ L = 0.18. In former simulation studies on monodisperse systems of self-dualshaped HBPs with L * = 12, the N U stability was also observed to be relatively small [24], and even vanish in binary mixtures [23]. Within the N U domain, we came across nematic LCs with oblate and prolate symmetry, respectively labelled as N − U and N + U . The occurrence of the N − U phase is predominantly detected close to the I-N transition and also far from it for σ L ≥ 0.25. By contrast, the N + U phase is mostly observed at larger packing fractions for 0 ≤ σ L ≤ 0.25. Consequently, at a given polydispersity, increasing the system density can produce an inversion of the nematic phase symmetry. While prolate and oblate particles commonly tend to trigger, respectively, N + U and N − U phases [11,18,34,35,37,40], although a significant polydispersity can counteract this tendency [13], the self-dual shape is in principle not expected to exhibit a clear oblate or prolate nature, but rather an ambivalent one. By applying our Onsager-like theory [22], we observed that the free-energy difference between N + U and N − U phases is very small close to the I-N transition and not much larger at increasing η (see ESI †). Although this theory is strictly valid for monodisperse systems and should serve here as a mere qualitative guideline, it confirms the slightly larger stability of the N − U in the vicinity of the I-N transition, in very good agreement with the tendencies observed in our simulations. Very small free-energy differences between oblate and prolate symmetries had also been reported by Martínez-Ratón and co-workers, who applied density-functional theory to study the phase behaviour of nearly self-dual-shaped monodisperse HBPs [12].
Uniaxial and biaxial smectic phases are obtained at larger packing fractions. To distinguish them from the nematic phases, we calculated the density distribution function along the nematic director, g ∥ (r ∥ ), identified by the order parameters. Uniaxial smectic (Sm U ) phases (green area in Fig. 1) only exhibit a prolate symmetry (Sm + U ), in agreement with former simulations of monodisperse systems [22]. In particular, their layers, whose thickness is roughly 13T to 15T , are perpendicular to the average direction of the particle length. The g ∥ (r ∥ ), calculated along this direction, displays periodically peaked profiles of the type reported in Fig. 2 and no indication of structural order in the other directions. By contrast, biaxial smectics, found at η > 0.50 and 0.05 ≤ σ L ≤ 0.20 (blue area in Fig. 1), can present a prolate (Sm symmetry. The former is characterised by layers piling along the particle length, while the latter by layers piling along the particle thickness. Upon increasing polydispersity, the Sm B phases acquire a more and more defined structural identity, with a weak-to-strong biaxiality crossover at approximately 0.05 < σ L < 0.10 and a full biaxial character (B 2 > 0.6) at σ L ≥ 0.18. At this value of size dispersity, our simulations highlighted a particularly rich phase behaviour, unveiling an I → N U → N B → Sm B sequence of phases that is ex- emplary shown in Fig. 3. In qualitative agreement with theory and experiments, we found the N B phase to be stabilised by a substantial degree of particle size dispersity.
In particular, our simulation results indicate σ L = 0.18 as the critical polydispersity above which the N B phase can form. Close to the N U − N B phase boundary, we find nematics with a relatively weak (but non-negligible) biaxiality and a residual oblate or prolate character. These phases, here referred to as weak biaxial nematics, are characterised by a biaxial order parameter in the range 0.2 ≤ B 2 ≤ 0.30 and one predominant uniaxial order parameter granting them prolate (N In contrast with the experiments on goethite particles [6], we do not observe a direct I − N B phase transition here. This apparent lack of agreement deserves some comments. First of all, the cuboids studied in these experiments are not self-dual-shaped. Their shape parameter, ν = 0.1, indicates a prolate geometry, which in monodisperse systems is expected to promote an I − N + U , rather than I − N B , transition as predicted by theories spanning almost five decades [11,18,34,35,37,40]. However, the goethite particles employed by Vroege and co-workers are not monodisperse, but exhibit a polydispersity between 20% and 25% in the three directions. Because polydispersity leads to fractionation [41] and these authors studied the phase behaviour in capillaries, the longer particles tend to accumulate towards the bottom, where the N B phase was found, de facto increasing the shape parameter of this subset of particles to the effective value of ν = 0.6 [6]. This particle geometry, evidently prolate, is very different from the self-dual shape applied here and a quantitative analogy is therefore not directly possible. Onsager theory within the Zwanzig approximation does not predict a direct I−N B transition in systems of polydisperse HBPs with ν = 0.1, but suggests the existence of the N B phase in a wide region of the η − σ L phase diagram, including for σ L < 0.1 [13]. While it is known that restricting orientations can significantly enhance the stability of the N B phase, both Belli's theoretical work and our simulations do not report a direct I−N B phase transition, whose existence has never been unambiguously confirmed by off-lattice simulations spanning more than twenty years [16,24,30,37,[42][43][44][45]. Indeed, our recent MC simulations and generalised Onsager theory applied to freely-rotating monodisperse HBPs had even excluded the existence of the N B phase, also at the self-dual shape [22]. To the best of our knowledge, there are no theoretical works on freely-rotating polydisperse HBPs that might help resolve this conundrum. While the phase diagram in Fig. 1 presents relevant discrepancies with that proposed by Belli, both works agree very well on the key role of polydispersity in the stabilisation of the N B phase. This is especially evident at σ L ≥ 0.18, where the stability region of the N B phase widens, remarkably reducing that of Sm and N U phases. This is not surprising as a large size dispersity is expected to hinder the formation of layered structures due to the absence of a well-defined structural periodicity in the longitudinal direction. In summary, our MC simulations of freely-rotating HBPs have revealed a rich phase behaviour that is characterised by three key results: (i ) a significant degree of particle size dispersity is needed to stabilise the N B phase; (ii ) self-dual-shaped HBPs do not exhibit direct I − N B phase transition in the range of size dispersities studied here; (iii ) the ambivalent nature of the self-dual shape provides uniaxial nematics that, in a relatively wide region of the η − σ L phase diagram, might well be oblate or prolate. More specifically, a significant polydispersity (σ L ≥ 0.18) prevents the discretization of the space along the nematic director and thus enfeebles the stability of the Sm phase, practically enhancing that of the nematic phase. This result is in line with experiments and theory, but on a mere qualitative basis only. The lack of a direct I − N B phase transition might appear in evident disagreement with previous experimental observations, which anyway employed prolate rather than self-dual-shaped particles, but agrees with a relevant number of simulation studies that could not confirm its occurrence in a reasonable range of particle anisotropies. As far as the ambivalence of the self-dual shape is concerned, it is interesting to observe an oblate-to-prolate symmetry inversion in the N U domain upon increasing density. Our modified version of Onsager theory for monodisperse biaxial particles suggests that the free-energy difference between N The classification of our equilibrium phases was done by examining orientational and positional ordering. To measure the long-range orientational order, we diagonalised the following second-rank symmetric tensor: where i indicates a generic particle,λ =x,ŷ,ẑ denotes its unit orientation vector along itslength (L), width (W ) and thickness (T ), respectively, I is the second-rank unit tensor, and the angular brackets denote ensemble average. When diagonalised, the tensor Q λλ produces three eigenvalues S 2,L , S 2,W , and S 2,T and their corresponding eigevectorŝ n,m, andl. For example, the tensor Q zz is related to the largest eigenvalue S 2,L , and corresponding eigenvectorn, which provides alignment along the particle axisx. The calculation of the eigenvalues S 2,L , S 2,W and S 2,T , referred to as uniaxial order parameters, allows to distinguish between an isotropic phase, where all eigenvalues vanish, and an ordered phase, where at least one of the eigenvalues is significantly larger than zero. We arbitrarily set the formation of a uniaxial LC phase when one of the three uniaxial order parameters is at least 0.40 (see Table I). To assess the system biaxiality, the biaxial order parameter for each axes can be calculated using the same symmetric tensor. To this end, the following equation is applied: The other two biaxial order parameters, B 2,W and B 2,T , can be calculated from similar expressions. To determine B 2 , it is sufficient to monitor the fluctuations of axes perpendicular to the main nematic director. For instance, if S 2,L is the dominant uniaxial order parameter, it is sufficient to monitor B 2,L as it indicates the fluctuations along axesŷ andẑ in the planes ofm andl. Table I shows the criteria to determine the symmetry of the phases observed in this work and consistent with Ref [1]. In Tables II to VIII, we report the uniaxial and biaxial order parameters of HBPs with σ L = 0.00 to σ L = 0.30.        To distinguish nematic from smectic phases, we calculated the pair correlation function parallel to the nematic director(s), which formally reads: where ρ = N p /V , N p is the number of particles, V the box volume, δ is the Dirac delta. S is a surface resulting from the intersection of a sphere with radius half of the simulation box and a plane perpendicular to the nematic director at a distance r from the center of the sphere. To practically evaluate it in our simulations, we employed the following expression [3]: where N c is the number of sampled configurations. N (r ) is the number of times that the distance between two particles projected along the nematic director, r ,ij , is in the interval |r + ∆r |. Finally, V is the volume of this region, that can be calculated as with R the half of the shorter side of the simulation box. In Fig. I

III. SMECTIC LAYER THICKNESS
Tables IX to XIII show the average thickness of smectic layers, defined by τ * = τ /T of HBPs with σ L = 0.05 to σ L = 0.20.

IV. FREE ENERGY CALCULATION
To gain an insight into the competition between prolate and oblate symmetries close to I − N U transition, we have compared the Helmholtz free energy of both N + U and N − U phases. To this end, we applied Onsager-like theory developed in our work for monodisperse systems, where the Helmholtz free energy is obtained as a combination of an entropy of mixing-like term and a correction for many-body interactions via virial expansions [2]. We stress that we are using a theory developed for monodisperse systems to calculate the free energy of polydisperse systems and hence our results should only be taken as a qualitative estimation of the free energy difference between the N + U and N − U phases at a particular value of η. Interested readers are referred to Ref [2] for additional details. The free energy calculation of our polydisperse systems are performed in the context of free particle rotation and assumes that the dominant uniaxial order parameter in the prolate (S 2,L ) and oblate (S 2,T ) nematic phases have the same value. The four virial coefficients are obtained from the S 2 calculated by Monte Carlo simulations (reported in Table II to Table VIII). The phases we consider in this analysis are the N U phases above the I − N U phase boundary for each polydispersity. For instance, at σ L = 0.05, our simulation results indicate the formation of an N U phase at η = 0.310 with dominant uniaxial order parameter S 2 = 0.402. These values of η and S 2 are used to obtain the corresponding virial coefficients and the entropy of mixing-like term for both N + U and N − U phases. In Fig. 2, the reduced free energy for N + U and N − U are shown for each size dispersity. For monodisperse systems, where σ L = 0.00, the free energy of the N + U phase is lower than the N − U phase, suggesting that the system has a higher propensity to possess prolate symmetry, in agreement with Ref [2]. In polydisperse systems, the N U phases are generally more likely to form N + U phases as well. However, the difference in free energy is such that relatively small density fluctuations, of the order of the thermal energy, might direct the systems to follow the oblate or prolate route, so determining their final symmetry. For σ L = 0.25, we find that the system is more likely to form the N − U phase, a tendency that is consistent with our simulation results.
σ L βF N P FIG. 2: Free energy per particle, plotted as a function of length polydispersity, σL. The parameter β is the inverse temperature. Each point corresponds to the first NU phase that is observed above the I-N phase boundary. Solid and empty circles indicate, respectively, the free energy of prolate, F * N + ≡ βF N + /Np, and oblate, F * N − ≡ βF N − /Np systems as obtained from theory.
With the same theory, the free energies F N + and F N − of N U phases deep in the N U region are also calculated. Tables XIV and XV show the estimated values of these free energies in monodisperse and polydisperse systems, respectively. Similar to what was observed for N U phases just above the I-N phase boundary, the free energy of prolate and oblate phases deep in the N U region are also very close to one another, suggesting an easy tendency for the systems to undergo director inversions with density fluctuations.