{'fraction_non_alphanumeric': 0.03266250101436339, 'fraction_numerical': 0.014160512862127729, 'mean_word_length': 4.527472527472527, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo simulation of binary mixtures of hard colloidal cuboids

We perform extensive Monte Carlo simulations to investigate the phase behaviour of colloidal suspensions of hard board-like particles (HBPs). While theories restricting particle orientation or ignoring higher ordered phases suggest the existence of a stable biaxial nematic phase, our recent simulation results on monodisperse systems indicate that this is not necessarily the case, even for particle shapes exactly in between prolate and oblate geometries, usually referred to as self-dual shape. Motivated by the potentially striking impact of incorporating biaxial ordering into display applications, we extend our investigation to bidisperse mixtures of short and long HBPs and analyse whether size dispersity can further enrich the phase behaviour of HBPs, eventually destabilise positionally ordered phases and thus favour the formation of the biaxial nematic phase. Not only do our results indicate that bidisperse mixtures of self-dual shaped HBPs cannot self-assemble into biaxial nematic phases, but they also show that these particles are not able to form uniaxial nematic phases either. This surprising behaviour is also observed in monodisperse systems. Additionally, bidisperse HBPs tend to phase separate in coexisting isotropic and smectic phases or, at relatively large pressures, in a smectic phase of mostly short HBPs and a smectic phase of mostly long HBPs. We conclude that limiting the particle orientational degrees of freedom or neglecting the presence of positionally ordered (smectic, columnar and crystal) phases can dramatically alter the phase behaviour of HBPs and unrealistically enlarge the region of stability of the biaxial nematic phase.


Introduction
Colloids are dispersions of solid particles or liquid droplets, between 1 nm and 1 µm in size, evenly suspended in a fluid. Their motion, basically controlled by the stochastic collisions with the fluid molecules, is driven by a thermal energy of the order of few kT per particle. Very interestingly, the interparticle forces controlling their phase behaviour and dynamics, including electrostatic and excluded volume interactions, are the same forces that determine the behaviour of atoms and molecules. As such, the possibility to synthesise particles in a variety of sizes and shapes and control their interactions, makes colloids ideal model systems to understand the physical laws underpinning structure and dynamics of atomic and molecular systems. This remarkable similarity provides a unique opportunity to unravel a number of processes at the molecular scale that are too fast to be detected by conventional microscopy. In particular, investigating how anisotropic colloidal particles organise themselves can contribute to unveil the phase and aggregation behaviour of molecular liquid crystals (LCs), enhance their fundamental understanding and thus optimise the design of functional electro-optical devices [1]. Nevertheless, the scientific relevance of colloids goes well beyond their use as mere model systems. The impressive advances in the synthesis of anisotropic 2. Methodology

Model
In Figure 1, we provide a visual representation of the particles studied in this work. They are freelyrotating hard rectangular parallelepipeds (cuboids) of thickness T , width W , and length L. The thickness T is the unit length, while the reduced width and length are, respectively, W * = W/T and L * = L/T . Following the theoretical work by Belli et al, we make use of a geometrical parameter, s, to describe the degree of particle size dispersity [13]. The dimensions of the short and long HBPs are then defined as a function of this parameter. For long HBPs, thickness, width and length are, respectively, T 1 = T (1 + s), W 1 = W (1 + s) and L 1 = L(1 + s). Similarly, for short HBPs, the dimensions are T 2 = T (1 − s), W 2 = W (1 − s) and L 2 = L(1 − s). By definition, the reduced length and width of the two species are the same: L * = L 1 /T 1 = L 2 /T 2 = L/T and W * = W 1 /T 1 = W 2 /T 2 = W/T . In order to compare our simulation results with the above mentioned theoretical predictions, we assign to these two aspect ratios the same values as those in Ref. [13], namely L * = 9.07 and W * = 2.96. This particular geometry was chosen to reproduce the experimental observations by Vroege and coworkers, who found a remarkably stable N B phase in systems of polydisperse boardlike particles [8]. Finally, the size-dispersity index is set to s = 0.2, a value that provides the richest phase behaviour among the binary mixtures investigated in Ref. [13]. While our particles are identical to those designed by Belli, they are completely free to assume any possible orientation, rather than only the six allowed by the Zwanzig model. This difference is particularly relevant when assessing the formation of uniaxial and biaxial phases, as we have recently noticed in systems of monodisperse boardlike particles [10]. Long and short board-like particles with thickness, width and length given by (T 1 , W 1 , L 1 ) and (T 2 , W 2 , L 2 ), respectively. The reduced length L * = L 1 /T 1 = L 2 /T 2 and width W * = W 1 /T 1 = W 2 /T 2 are the same for both species. The size-dispersity index of the particles shown in this figure is s = 0.2 (see text).

Monte Carlo simulations
We perform Monte Carlo simulations of systems containing between N = 2000 and N = 4000 HBPs, where N = N 1 + N 2 is the sum of long (N 1 ) and short (N 2 ) HBPs. Our cuboids interact via a hard-core potential and, as such, their path towards the equilibrium state, which only depends on their geometry and system packing, can be investigated in terms of free volume available, or, equivalently, by checking and discarding the occurrence of particle overlaps. To this end, we applied the separating axes method described by Gottschalk et al [20] and adapted by John and Escobedo to study the phase behaviour of tetragonal parallelepiped particles [21]. We refer the interested reader to the Appendix of Ref. [21] for additional details.
All the simulations were run in the isobaric-isothermal ensemble, which constrains a constant number of HBPs into a periodic simulation box of volume V = L x ×L y ×L z , with L x , L y , and L z the box dimensions. To study phase coexistence, we decided to apply a direct coexistence simulation method, where the two species are initially separated by an interface and thus occupy different regions of the simulation box. In both regions, the particles are perfectly aligned, so to create a dense biaxial nematic phase consisting of short or long HBPs. Undesired finite size effects have been taken into account by employing elongated boxes, where L x is at least four or five times larger than the other two dimensions. This choice ensures a bulk behaviour in the regions far from the interface and allows us to impose an isotropic pressure gradient to a system that, due to the presence of an interface, is actually affected by an anisotropic pressure tensor. However, as long as the box is sufficiently elongated, this pressure tensor, scaling with L −1 x , can be safely neglected [32]. The systems in such an initial configuration, containing the two species at different relative concentrations, in the range 0 ≤ x 1 = N 1 /N ≤ 1, have been gradually expanded or compressed to the target reduced pressure P * = βP T 3 , where β is the inverse temperature. Moreover, to better identify the boundary of the isotropic monophasic regions, we also compressed systems of randomly mixed HBPs. This also allowed us to confirm that the interface imposed in our direct coexistence simulations was actually the one with the lowest free energy. Each MC cycle consisted of N attempts of displacing and/or rotating randomly selected particles, plus an attempt to modify the three box lengths independently. Translational and rotational moves as well as volume changes were accepted if no overlap was detected.
To address equilibration, we calculated the uniaxial (S 2 ) and biaxial (B 2 ) order parameters as well as the packing fraction η = (N j v j )/V , with v j = T j × W j × L j the volume of a particle of species j = 1 or 2. The systems were considered to be at equilibrium when η, S 2 and B 2 achieved a steady value within reasonable statistical fluctuations. To weight the contribution of each species to the global order of the system, the order parameters have been calculated for short and long HBPs separately. In particular, to determine the nematic order parameter and nematic director associated to each particle axis, the following traceless symmetric second-rank tensor has been diagonalized: whereλ λ λ=x,ŷ,ẑ refers to the unit orientation vectors of the generic particle i along W , T ,and L, respectively, I is the second-rank unit tensor, and the angular brackets indicate ensemble average. Three eigenvalues and three eigenvectors are obtained from diagonalization of Q λλ , with the largest positive eigenvalue and the associated eigenvector determining, respectively, the nematic order parameter S 2 and the uniaxial nematic director [22]. In particular, the nematic order parameter S 2,L is the largest eigenvalue of the tensor Q zz and reads where the eigenvectorn is the nematic director associated to the preferential orientation of the particle axisẑ. Similar expressions can be obtained to determine the remaining two uniaxial order parameters, S 2,W =m·Q xx ·m and S 2,T =l·Q yy ·l, withm andl their respective nematic directors. By contrast, the biaxial order parameter associated to the nematic directorn is calculated as follows [23] In principle, B 2,W and B 2,T can be obtained by similar expressions. However, to assess phase biaxiality is enough to monitor the biaxial order parameter associated to the principal particles axis, defined as the axis displaying the largest uniaxial order parameter [24][25][26]. Therefore, here we only consider the biaxial order parameter associated to the main nematic director and refer to it as B 2 . Additional details are available in our recent work on monodisperse HBPs [10]. Finally, to calculate the particle composition in each coexisting phase, we run N V T simulations within the miscibility gap at x 1 = 0.5 and, in some cases to improve statistics, along the same tie-line at x 1 = 0.2, 0.4, 0.6 and 0.8. We then measured the particle composition by dividing the box in layers of cross section A c = L y × L z and height L 1 parallel to L x .

Results
Details of the binary mixtures investigated in this work are reported in Table 1, where we provide a summary of the ordering observed in the monophasic and biphasic regions. The P * − x 1 (pressurecomposition) phase diagram of these systems is given in Figure 2. Table 1. Details of the systems studied in this paper, consisting of binary mixtures of short and long HBPs. For comparison, we report the reduced pressure P * , phases at equilibrium a , composition of long HBPs x 1 , packing fraction η, and uniaxial and biaxial order parameters. Superscripts (1) and (2) refer to the phases at equilibrium in the biphasic region of the phase diagram. Four main regions can be identified: (i ) the region of stability of the I phase, where both species are well-mixed together across the complete range of compositions; (ii ) the region of stability of a Sm phase that incorporates mostly long HBPs, referred to as Sm 1 ; (iii ) a dark gray shaded area being the region of I/Sm 1 coexistence; and (iv ) a light gray shaded area indicating equilibrium highly ordered phases. The empty squares indicate overall system concentrations, set to x 1 = 0.4, 0.5 or 0.6. Regardless of these values, the resulting coexistence concentrations in the I and Sm phases are generally in good agreement, which is of the order of the error bars shown in Figure 2. The reason why we are showing an overall concentration rather than an other is only due to the fact that some Sm phases display a number of defects whose annihilation time is especially long. Although the presence of these defects does not seem to particularly influence the concentrations at coexistence, in Figure 2 we decided to include those overall concentrations referring to Sm phases whose defects have been (almost) completely annihilated within our simulation time. At increasing pressure, the systems undergo a transition from the I phase to a two-phase region where I phases rich in short HBPs coexist with Sm phases rich in long HBPs. At larger pressures, we observe a region where two Sm phases coexist: Sm 1 , rich in long HBPs, and Sm 2 , mostly composed by short HBPs. The pure-component systems, indicated by the vertical lines x 1 = 0 and x 1 = 1, experience transitions from the I phase to the Sm phase with no evidence of intermediate N U phases. This result is rather remarkable as theoretical works on spheroplatelets [11,16] and HBPs [10,13] did identify the presence of stable N U phases at particle geometries equal or close to the self-dual shape, where the relation √ LT = W between particle dimensions holds. Nevertheless, as we could establish in our work on monodisperse HBPs with L * = 9 and L * = 12, the region of stability of the N U phase for this particle shape is so small that pinning down its boundaries becomes very challenging, especially due to the weak first-order character of the I-to-N phase transition of this particular particle geometry. As a general tendency, the transition pressure from the monophasic to the biphasic region strictly depends on the relative content of short and long HBPs and spans from P * ≈ 0.328 for x 1 = 0.0 to P * ≈ 0.094 for x 1 = 1.0. In terms of packing fraction, the monophasic-to-biphasic phase transitions occur at approximately 0.35 < η < 0.40, the same interval at which monodisperse self-dualshaped HBPs [10] and monodisperse hard spheroplatelets [19] experience an I-to-N phase transition. Although the theory of Ref. [13] does not predict the existence of an I/Sm biphasic region, it locates the I-to-N transition of binary mixtures of HBPs at η ≈ 0.28, a packing fraction at which our binary mixtures are still isotropic.
At pressures generally below P * < 0.1, depending on particle composition, a single I phase is found, where short and long particles are evenly mixed throughout the simulation box. Figure 3 displays a typical configuration of an I phase at P * = 0.093 and composition x 1 = 0.8. While the configuration in the top frame shows both particle species, that in the middle frame only consists of short HBPs to provide an unambiguous evidence of their homogeneous distribution thoughout the simulation box. We stress that the initial configuration of this mixture consisted of two separated phases, each incorporating either short or large HBPs. The average composition profiles along the longest box dimension for long and short HBPs are given in the bottom frame and confirm the presence of a well-mixed I phase. Although some small clusters of long HBPs are observed in the top frame, the long-ranged order of this system is very weak, with the pair correlation functions (not shown here) decaying fast to 1 at length scales comparable to L 1 . The phase diagram in Figure 2 unveils the coexistence between a Sm 1 phase, rich in long HBPs, and an I phase, rich in short HBPs. A typical configuration of this intriguing phase coexistence is provided in Figure 4 for a binary system of 4000 HBPs at P * = 0.13 and x 1 = 0.6. The corresponding average composition profiles, provided in the bottom frame, unveil the formation of a Sm phase rich in long HBPs with x 1 ≈ 0.85 (black circles) and an I phase rich in short HBPs with x 1 ≈ 0.28 (red triangles). Interestingly enough, the former incorporates short HBPs in the interlayer spaces as well as within the smectic layers. Transverse interlayer particles had already been observed in Sm phases of monodisperse spherocylinders [27][28][29][30][31], but their occurrence had always been considered particularly rare, especially for the prohibitive free-energy barriers to reach the transverse state [29]. Here, it seems that the probability of observing transverse inter-layer HBPs is not so low, although we have not assessed this point quantitatively. In particular, short HBPs accommodate themselves in the spaces available between long HBPs and give rise to a very peculiar Sm phase with prolate-like layers, formed by short and long HBPs oriented along their main axis, and oblate-like layers, incorporating mainly short HBPs oriented along their minor axis. Transversely oriented long HBPs can in principle also be observed, but, as also found in monodisperse systems of spherocylidners, their occurrence is expected to be rare. To quantify these observations, we calculated the longitudinal pair correlation function, g (r ·n), between particles at distance r projected along the nematic directorn. It is convenient to calculate these functions in cubic rather than elongated boxes, where a single phase is observed. To this end, once that the boundaries of the miscibility gap had been determined, we equilibrated a number of Sm phases in the region of the P * −x 1 phase diagram where monophasic equilibrium is expected. In Figure 5, we show g (r·n) for long and short HBPs in a pure Sm 1 phase at P * = 0.15 and x 1 = 0.95. The black solid line, g 11 , suggests that the relative distance between long HBPs is d/T ≈ 12 in the direction ofn, corresponding to a particle length (L 1 /T = 10.884) plus the interlayer spacing, approximately given by the short-particle thickness (T 2 /T = 0.8). The non-zero values of g 11 at intermediate distances between contiguous smectic layers confirm that long HBPs can also be observed in the interlayer region, where they can assume a transverse orientation. The red dashed line (g 22 ) provides the spatial correlations between short HBPs. Although in this case the peaks are significantly smaller than those of g 11 , a well-defined pattern can still be recognised. While the main peaks of this pattern cannot unambiguously resolve whether short HBPs are located in the prolate-like or in the oblate-like layers, the minor peaks in between clearly indicate that they are most probably located in both of them. This conclusion is confirmed by the analysis of g 12 , given by the green dot-dashed line in Figure 5, which describes the distribution of short HBPs with respect to long HBPs. The broadness of its main peaks denotes a relatively wide in-layer positional distribution of short HBPs, which, due to their reduced length as compared to the layer thickness, are free to fluctuate in the direction ofn and around the center of mass of the layers. The secondary peak at exactly d/2T indicates that short HBPs are indeed laying in the interlayer spacing or, equivalently, forming oblate-like smectic layers. At larger pressures, above P * = 0.32, two coexisting Sm phases, one rich in long HBPs (Sm 1 ) and the other in short HBPs (Sm 2 ), are observed. An example of this two-phase equilibrium is shown in Figure 6, where 2000 cuboids at P * = 0.35 and x 1 = 0.2 give rise to two separate layered structures. Both phases present a number of unresolved defects, which are expected to be annihilated at time scales larger than our simulation time. While short HBPs are again observed in the bulk of Sm 1 , it is very difficult to see long HBPs in the bulk of Sm 2 , as the average composition profile in the bottom frame of Figure 6 points out. A simple visual inspection of the top frame of this figure suggests that the minor axes of long HBPs appear to be significantly aligned and thus that the Sm 1 phase might possess a degree of biaxiality. These observations are actually not limited to the region of the phase diagram where two Sm phases coexist, but also at lower pressures, where an I/Sm coexistence has been detected. To address this point, we have calculated the order parameter defined in Eq. (2) and (3).
Typical uniaxial and biaxial order parameters are shown in Figure 7 for a binary mixture at x 1 = 0.4 and 0.1 ≤ P * ≤ 0.32. For the sake of clarity, it is important to note that both sets of parameters have been calculated separately for each species and thus provide a measure of the order resulting from the orientation of the HBPs regardless of the phase they belong to. Nevertheless, since March 16, 2018 Molecular Simulation paper˙V2 the phases observed in the biphasic region predominantly consists of either short or long HBPs, to a very good approximation S 2 and B 2 can also be employed to measure their orientational order. In particular, the left frame of Figure 7 shows the uniaxial order parameter of the long HBPs (black circles) abruptly increasing at P * > 0.125, where the transition from the I to the two-phase region is observed. By contrast, short HBPs (red squares) are very weakly ordered and persist in the I phase up to relatively large pressures. Above P * = 0.25, these particles start to form small oriented clusters in the I phase, which can slightly enhance the value of S 2,L . Additionally, similarly to the organisation observed in the bulk of the Sm 1 phase in Figure 4, some of them succeed to diffuse through the layers of the Sm phase and contribute to further increase their global longrange orientational order. As far as the biaxial order parameter B 2 is concerned, we notice that long HBPs do show evidence of an appreciable degree of biaxiality at P * > 0.25 (right frame of Figure 7). Biaxial Sm phases had also been observed in monodisperse systems of slightly oblate HBPs [10]. By contrast, no significant biaxiality arises from short HBPs.

Conclusions
In summary, we have performed Monte Carlo simulations to investigate the phase behaviour of binary mixtures of short and long HBPs. Upon increasing density, these systems experience a phase transition from a pure I phase, where the two species are homogeneously mixed, to a two-phase region where I phases of mostly short HBPs coexist with Sm phases of mostly long HBPs. At relatively larger densities, a separation between a long-HBP-rich Sm phase (Sm 1 ) and a short-HBP-rich Sm phase (Sm 2 ) is observed. The positionally ordered Sm 1 phase incorporates short HBPs within its layers, aligned along the nematic director, and in the inter-layer spacing, aligned perpendicularly to the nematic director. Transverse interlayer particles had been previously observed in Sm phases of monodisperse spherocylinders, although their occurrence was rather rare and usually restricted to one or two transverse particles per interlayer spacing [27][28][29][30][31]. Although we have not assessed this specific phenomenon quantitatively, here we have detected a significant number of interlayer HBPs oriented perpendicularly to the nematic director, as our pair correlation functions highlight. This intriguing organisation unveils oblate-like smectic layers embedded in a prolate-like smectic phase of mainly long HBPs. By contrast, the Sm 2 phase is mostly made of short HBPs which are preferentially ordered along the nematic director. Interestingly, and perhaps surprisingly, our HBPs do not form N U phases, even at x 1 = 0 and x 1 = 1. This is in agreement with our recent results on monodisperse HBPs, whose phase diagram at W ≈ √ LT , the geometry explored here, eventually shows a very small region of stability for the N U phase. Finally, we do not find any significant evidence of the formation of the N B phase, as predicted by Onsager-type theory within the Zwanzig model [13] and by variational cluster expansion theory applied to HBPs with their long axes fully oriented [14]. We conclude that restricting the orientational degrees of freedom by imposing a full or partial alignment of particles can have a significant impact on the extension of the stability region of the N B phase. We are aware that the complex free-energy landscape of these systems, which is not fully captured by Onsager's or variational cluster theories, might also play an important role. The experimental observation of especially stable N B phases were performed on colloidal dispersions of board-like goethite particles with L/W W/T , exactly as in the present work, but a size polydispersity of 20%-25% in all directions was applied [8]. This significant size-dispersity, as compared to the bi-dispersity studied in this work, is most probably a strategic ingredient to circumvent phase separation and destabilise the Sm phase in favour of the N B phase. Preliminary simulation results in our group indicate that this is the case, but further investigation is currently ongoing.