**Experimental Conditions for the Stabilization of the Lyotropic Biaxial Nematic Mesophase**

## **Erol Akpinar <sup>1</sup> and Antônio Martins Figueiredo Neto 2,\***


Received: 1 February 2019; Accepted: 15 March 2019; Published: 19 March 2019

**Abstract:** Nematic phases are some of the most common phases among the lyotropic liquid crystalline structures. They have been widely investigated during last decades. In early studies, two uniaxial nematic phases (discotic, ND, and calamitic, NC) were identified. After the discovery of the third one, named biaxial nematic phase (NB) in 1980, however, some controversies in the stability of biaxial nematic phases began and still continue in the literature. From the theoretical point of view, the existence of a biaxial nematic phase is well established. This review aims to bring information about the historical development of those phases considering the early studies and then summarize the recent studies on how to stabilize different nematic phases from the experimental conditions, especially, choosing the suitable constituents of lyotropic mixtures.

**Keywords:** Lyotropic liquid crystals; uniaxial nematic phase; biaxial nematic phase; stabilization of nematic phases; micelle; surfactants

## **1. Introduction**

Lyotropic liquid crystalline phases may be encountered in mixtures of amphiphilic molecules and a solvent, with or without a salt, in particular regions of the multidimensional phase diagram (relative concentrations and temperature) of these mixtures [1,2]. From the symmetry point of view, two types of lyotropic nematic phases exist: uniaxial (two of them were identified, the discotic, ND, and the calamitic, NC, phases, *D*∞*<sup>h</sup>* symmetry), and one biaxial (NB, *D*2*<sup>h</sup>* symmetry). In the uniaxial phases we define a director represented by a vector <sup>→</sup> *n* that corresponds to the optical axis of the phase. The <sup>→</sup> *<sup>n</sup>* and <sup>−</sup><sup>→</sup> *n* states are indistinguishable. In the case of the NB phase, there are three orthogonal two-fold symmetry axes (<sup>→</sup> *l* , → *m* and <sup>→</sup> *n*, where <sup>→</sup> *<sup>n</sup>* <sup>=</sup> <sup>→</sup> *<sup>l</sup>* <sup>×</sup> <sup>→</sup> *m*) and two optical axes [3–5]. These phases are diamagnetic and high-intensity magnetic fields ( *H* ∼ 10*kG*) are used to orient them in actual experiments.

After the discovery of the first lyotropic mixture showing the biaxial nematic phase [1], intensive studies have begun to characterize its properties and structure [6]. These studies were mostly conducted to understand whether the biaxial phases are thermodynamically stable or whether a mixture of the two uniaxial phases. Many experimental and theoretical studies have confirmed the first hypothesis.

In the phase diagrams of lyotropic mixtures presenting the three nematic phases, the NB phase is mainly located between the other two uniaxial phases [1,2,7–15]. The transition from the uniaxial (ND or NC) to the biaxial phase is of second order, as predicted by the Landau-type mean-field theory [16,17]. Since the NB phase domain (at least in some phase diagrams) is located in between those from the ND and NC phases, some researchers claimed that the NB phase was thermodynamically unstable and even constituted by a mixture of both uniaxial phases [18–22]. Assuming that the ND and NC phases consist

of disc-like and cylindrical-like micelles, respectively, the NB phase should be composed by a mixture of these types of micelles. However, it was theoretically proved that disc-like and cylindrical-like objects could not coexist, i.e., a mixture of both objects (in the case of nematic phases, micelles) is not stable and phase separation occurs [23].

This review article deals with the recent studies related to the stabilization of different lyotropic nematic phases from the experimental point of view. Mainly, we will review some important studies obtained in our research laboratories, giving to experimentalists hints of mixture compositions that favor the presence of the NB phase in the phase diagram.

## **2. Background**

The first lyotropic nematic phase was reported towards the end of the 1960s [24]. In the following years, many studies were conducted to understand the properties of this phase. In these early studies, the existence of two types of nematic phases was determined and these phases were classified as Type I and Type II, according to the orientation of the phase director with respect to an external applied magnetic field [25–27]. Since the lyotropic mixture investigated was composed by amphiphilic molecules with carbonic chains, it was verified that the director of the Type I (Type II) phase aligns parallel (perpendicular) to the magnetic field direction. In the following years, Type I (Type II) phase was named NC (ND) [28]. At that time, it seemed natural to assume that Type I and Type II phases were composed of cylindrical-like and disc-like micelles [28–32], respectively, as shown in Figure 1.

**Figure 1.** Orientations of (**a**) disc-like and (**b**) cylindrical-like micelles perpendicular and parallel, respectively, to magnetic field direction, <sup>→</sup> *H*, along z. <sup>→</sup> *n* (phase director) represents the preferred alignment of the individual micelles with their local directors.

Freiser [16,33] theoretically predicted the existence of the NB phase for the first time [5,34–38] after Taylor et al. [39] reported a biaxial smectic C phase in thermotropic liquid crystals. However, at that time, there were no experimental results to support this theoretical prediction. In 1980, Yu and Saupe [1] experimentally showed the existence of the NB phase for the first time in the ternary mixture of potassium laurate (KL)/decanol (DeOH)/D2O via Nuclear Magnetic Resonance (NMR) and microscopic conoscopy studies, as shown in Figure 2. They constructed the partial phase diagram of this ternary mixture where the NB phase was an intermediate phase between the two uniaxial ND and NC phases. In 1985, Neto et al. [40] investigated the same ternary mixture in more details. They used the optical microscopy, laser conoscopy and X-ray diffraction techniques and found that the biaxial phase domain was relatively large (~15 ◦C) for appropriate relative concentrations of the mixture compounds.

**Figure 2.** Phase diagram of KL/DeOH/water mixture which was established for the first time in the Reference [1].

In 1982, Bartolino et al. [7] reported a partial phase diagram of the ternary mixture of sodium decylsulfate (SdS)/DeOH/H2O to contribute to the understanding of the biaxial nematic phases, as shown in Figure 3a. They showed that the uniaxial ND (NC) phase had positive (negative) birefringence, as shown in Figure 3b. Their results seem to show a discontinuity of the birefringence in the biaxial domain, however, to prove this claim it should be necessary measurements of the birefringences with more precision in temperature.

**Figure 3.** (**a**) Partial phase diagram of SDS/DeOH/water mixture given in the Reference [7]. (**b**) The birefringences of the uniaxial nematic phases for the sample of 35.8% SDS.

Finding the first biaxial nematic phase Yu and Saupe not only verified Freiser's prediction on the existing that phase, but also started interesting discussions about the nature and structure of the NB phase [1,36,41–47]. Theoretical and computer simulation studies of binary mixtures of rod-like and disk-like rigid particles were reported [19–21,48–57]. In 1984 Stroobants and Lekkerkerker proposed a model, based on the Onsager's theory, where rod-like and disc-like particles coexist (mixture of disc-like and cylindrical-like micelles, MCD model), giving rise to a biaxial nematic phase [21]. However, some theoretical and experimental studies showed that MCD-type models presented problems to describe the NB phase [58]. For instance, in 2000, Kooij and Lekkerkerker reported an interesting experimental study about the liquid crystal phase behavior of suspensions with mixtures of rod- and plate-like units. They verified that this type of mixture separates, i.e., the coexistence of rod- and plate-like particles is not stable [23]. In another study, Palffy-Muhoray et al. [22] investigated the phase behavior of a binary mixture of uniaxial nematic liquid crystals (rods and disks), in the absence of the external field, employing the Maier-Saupe mean-field theory. They concluded that the NB phase was not thermodynamically stable. In a recent study Martinez-Raton et al. showed that phase-separation of rod-plate mixtures and the stability of the biaxial phase depend on the aspect ratio and molar fraction of the rods [59]. Sharma et al. [60] showed that there could be a possibility to stabilize the biaxial phase in mixtures by increasing either the isotropic or the anisotropic parts of the interspecies interaction strength. The latter situation was also supported by other theoretical studies, i.e., biaxial phase may be favored in multi-component (polydisperse) mixtures [61] by the specific interactions (indeed, hydrogen bonding) between the unlike basic units [62,63]. However, for the symmetric [64–67] or asymmetric [23,68–73] binary rod-plate mixtures, considering the stability of the biaxial phase with respect to phase separation, the coexistence of rod-plate particles is not possible and the phase separation is favored. In summary, the separation of the rod and plate-like objects (in the case of the lyotropic systems, prolate or cylindrical-like and oblate or disc-like micelles) requires the rejection of the MCD-type model, as experimentally showed by Kooij and Lekkerkerker [23].

Although the discussions about the stability of the NB phase were going on, experimentalists tried to find new lyotropic mixtures presenting the NB phases. In 1985, Galerne et al. [8] introduced a new ternary mixture of rubidium laurate (RbL)/DeOH/H2O. They investigated the defect lines (disclinations) in both uniaxial and biaxial nematic phases and observed that the disclination lines in the biaxial region made "zig-zags" while the uniaxial counterparts exhibited a "flexible" behavior, i.e., the disclinations observed in the biaxial phase seem to be different from those in the uniaxial phases. However, some time later the same authors show that the same zig-zag disclinations were observed in uniaxial phases with different mechanism than in the biaxial [74–76]. Santos and Galerne investigated the uniaxial discotic nematic phase of KL/DeOH/D2O by Rayleigh-scattering technique [77] in the vicinity of the uniaxial to biaxial phase transition. They observed the dynamic of fluctuations of the biaxial order parameter close to the uniaxial to biaxial phase transition, at the uniaxial phase domain, and showed that the micelles in the uniaxial discotic nematic phase exhibit biaxial ordering.

In 1985, a reliable model was proposed by Neto and co-workers [41,78], assuming that the micelles in all the three nematic phases in lyotropics have orthorhombic symmetry. They investigated the uniaxial and biaxial nematic phases of the mixture KL/DeOH/D2O by X-ray diffraction and observed that the micelles in all the nematic phases are locally biaxial, arranged in a pseudo-lamellar structure. This means that in the three nematic phases, the micelles are similar from the point of view of the local symmetry. Their model, so-called "intrinsically biaxial micelles model, IBM", is mainly based on two pillars: (a) micelle symmetry and (b) orientational fluctuations of the micelles. According to the IBM model, the micelles have orthorhombic symmetry, as shown in Figure 4a, in all three different nematic phases and different orientational fluctuations are responsible for the stabilization of the three nematic phases. Figure 4a shows a sketch a micelle as an object of orthorhombic symmetry, with typical dimensions A', B' and C'. The axes of the local coordinate system fixed in the micelle are *α*, *β* and *γ*. These orthorhombic micelles exhibit different orientational fluctuations for the stabilization of the ND, NB or NC phases. The orientational fluctuations around the axis perpendicular to the largest micelle surface (along axis y) give rise to the ND phase, as shown in Figure 4b. In this situation, the micelle size A' approximately equal to B' (A'~B'). By changing temperature to reach the transition from

ND to NB, the micelle size along the symmetry axis *α* increases, and the small-amplitude orientational fluctuations along the three axes (x, y and z) lead to the NB phase. If the temperature is changed until obtaining the NC phase, the micelle dimension along *α* continues to increase, which favors the orientational fluctuations around the local α micellar axis (along the z axis).

**Figure 4.** (**a**) Sketch of the orthorhombic micelle in the framework of the IBM model (Reference [79]). The detergent amphiphilic bilayer is represented by C'. According to the IBM model, based on the different orientational fluctuations, the formations of (**b**) ND, (**c**) NB, and (**d**) NC phases. <sup>→</sup> *n* represents the director of the ND and NC phases."

A neutron scattering study reported by Hendrikx et al. [80] in 1986 supported the IBM model. In that study, authors experimentally showed that the micelles of the NC phase are statistically biaxial.

Until 1989, since the lyotropic mixtures exhibiting the biaxial nematic phases contained alcohol and surfactant, Oliveira et al. [9] achieved to find a new mixture of KL/decylammonium chloride (DACl)/water, alcohol-free. This new mixture was important because, as suggested by Saupe et al. [81], if a mixture has surfactant and alcohol, there is a possibility of a slow esterification occur in the mixture. In the study of the [9], researchers investigated the mixtures of KL/DACl/water and KL/DeOH/water, i.e., alcohol-free and with alcohol mixtures, respectively, to compare the physical-chemical stability of both types of mixtures. Their long-term studies were based on the birefringence and X-ray diffraction measurements. It was shown that the phase-transition temperatures, birefringences and microscopic structures changed with time in the case of the mixtures including alcohol. However, for the alcohol-free mixtures, especially the X-ray diffraction results of the nematic phases, revealed that the alcohol-free mixtures were more stable than the one with alcohol.

In 1991, Vasilevyskaya et al. [82] investigated a new lyotropic mixture exhibiting the biaxial nematic phase SdS/DeOH/H2O/Na2SO4. Their experimental study was based on the measurement of the birefringences of the nematic phases, as a function of the temperature. They observed that this new system had much smaller birefringence values with respect to the conventional lyotropic mixture KL/DeOH/H2O [1]. As a result of their study, they concluded that the phase transitions from biaxial to uniaxial arise from different orientational fluctuations of the micelles around their axes. Their conclusion is similar to those from the IBM model.

Experimental and theoretical studies stated that the biaxial nematic phase is an intermediate phase between two uniaxial nematic ones. Ho et al. [83] reported an interesting study related to a new lyotropic mixture of sodium lauroyl sulfate/hexadecanol (HDeOH)/water, exhibiting a biaxial nematic phase in the case of diluted solutions. They obtained this biaxial phase after applying the sonication process on a gel-like phase. They showed, surprisingly that a biaxial phase could be found in a different place in the lyotropic phase diagram, not as an intermediate phase between two uniaxial nematics. Another interesting study was published by Quist [84] to investigate novel lyotropic mixture of sodium dodecylsulfate SDS/DeOH/H2O via NMR spectroscopy. Based on the Landau theory, the phase transitions from uniaxial to biaxial nematic phase transitions are of second-order. However, Quist proposed that these phase transitions might be of first-order. The author attributed the first–order ND-NB and NB-NC phase transitions to a variation in the aggregate (micelle) shape. This result is really very interesting, because other experimental and theoretical studies showed that the uniaxial to biaxial transitions have to be of second-order. In addition, early [40] and recent [13,15] studies, especially obtained from the temperature dependence of the birefringence, showed that these phase transitions are of second-order, as predicted by mean-field theory [16,18].

From the diamagnetic-susceptibility anisotropy (Equation (1)) point of view, considering the diamagnetic susceptibilities along three two-fold symmetry axes being χ33, χ<sup>22</sup> and χ11, the biaxial nematics are identified with positive and negative Δχ, NB+ and NB−.

$$
\Delta \chi = \chi\_{33} - \frac{\chi\_{11} + \chi\_{22}}{2}. \tag{1}
$$

In the case of Δχ >0(Δχ < 0), the largest (smallest) diamagnetic susceptibility is larger (less) than the average of the diamagnetic susceptibilities of other axes and the biaxial phase aligns parallel (perpendicular) to the magnetic field direction [84,85]. Quist [84] reported a study with these two biaxial nematic phases, with positive and negative diamagnetic-susceptibility anisotropy in the phase diagram; however, they did not study an eventual phase transition between these NB+ and NB<sup>−</sup> phases. Moreover, from the symmetry point of view, there is no transition between them. This type of "phase transition" was investigated by de Melo Filho et al. [85] in the novel lyotropic mixture of tetradecyltrimethylammonium bromide (TDTMABr)/DeOH/H2O by measuring the deuterium quadrupolar splitting, via Nuclear Magnetic Resonance. In 2009, van den Pol et al. [36] reported a model-system to experimentally obtain the biaxial phase in colloidal dispersions of board-like particles. According to theory and simulations, the particles have the dimensions L/W ≈ W/T, where L, W and T are length, width and thickness, to obtain a biaxial nematic phase [48,50,86]. They used goethite particles as a model-system with L/W = 3.1 and W/T = 3.0, and investigated the properties of the model-system by small angle x-ray scattering, under the effect of external magnetic field. Their results indicated the evidence of the biaxial nematic phase. Indeed, their model system is similar to what is proposed by the IBM model in terms of the micelle shapes and, from this respect, the study in [36] supports the IBM model.

As reviewed above, despite the evidence of the stability of the biaxial nematic phase in lyotropics, some studies still reject the existence of this phase. Indeed, these controversies arise from the absence of enough number of lyotropic mixtures in the literature showing this phase and the missing thermotropic counterpart.

In the next sections, we will discuss how to stabilize the lyotropic biaxial nematic phase from the sample preparation point of view. In other words, we will show which experimental conditions have to be fulfilled to prepare a lyotropic mixture presenting this phase.

## **3. Recent Experimental Studies on the Stabilization of the Lyotropic Biaxial Nematic Phase**

Lyotropic mixtures are prepared by dissolution of a surfactant molecule in water (mainly). In some cases, a cosurfactant and/or an electrolyte are/is added to that binary solution. Indeed, while the former provides both the decrease in the repulsions between similarly charged headgroups of the surfactants at the micelle surfaces and the increase in the attractive forces (van der Waals) between the surfactant alkyl chains in the interior of the micelles, the latter one just reduces the repulsions between the headgroups. Because of this, modifications at the micelle surfaces and in the micelles shape anisotropy play a key role to obtain not only uniaxial nematic phases but also the biaxial one. Thus, we will review our recent experimental results focusing, mainly, in the modifications at micelle surfaces and in the micelles shape anisotropy.

## *3.1. Effects of Alkyl-Chain Lengths of Surfactant and Cosurfactant*

Surfactant molecules consist of two different parts, hydrophilic water-soluble head groups and hydrophobic water-insoluble alkyl tails or chains. By the same way, cosurfactants have water-insoluble parts and polar water-soluble parts attached. In general, long-chain alcohols, such as 1-decanol, are added to the lyotropic mixtures as cosurfactants.

In 2012, we studied the effect of the alkyl-chain length of the alcohols in the mixture of KL/K2SO4/alcohol/water [11]. We constructed a partial phase diagram by changing the number of carbon atoms in the alcohol-chain length from 8 (1-octanol) to 16 (hexadecanol), as shown in Figure 5. The laser conoscopy technique, which is a very accurate method to measure the birefringences of the nematic phases [87], was applied to determine the uniaxial to biaxial nematic phase transitions from the temperature dependences of the birefringences. The results indicated that there is a strong relation between the alkyl-chain length of the surfactant and that of the alcohol. As the alkyl-chain length of the alcohol gets longer (shorter), the NC (ND) phase is stable and at intermediate alkyl chain length there exists the possibility to obtain the NB one. In general, at constant mixture composition, the NB phase is stabilized when the *n* = *m* ± 2, where m (n) is the number of carbons in alkyl chain of the main surfactant (co-surfactant) molecule.

**Figure 5.** Partial phase diagram of KL/K2SO4/alcohol/water mixtures [11]. n is the number of carbon atoms in the alkyl chain of the alcohol molecules. I, 2P, MP and C represent an isotropic phase, two-phase region, multi-phase region and crystalline-like phase, respectively. Solid and dashed lines are only guides. The grey region shows the biaxial nematic phase domain.

In a more recent study [15], we have also investigated the relation between the alkyl-chain length of the surfactant and the alcohol in a different experimental design. In this case, we examined the effect of the surfactant alkyl-chain length on the stability of the different nematic phases, at constant alcohol chain-length and mixture composition. The mixtures chosen were potassium alkanoates/Rb2SO4/DeOH/water and sodium alkylsulfates/Na2SO4/DeOH/water. The experimental techniques employed were the laser conoscopy, polarizing optical microscopy and small-angle x-ray scattering. In both systems, it was observed that as the number of the carbon atoms in the surfactant alkyl-chain length increases (decreases), the stabilization of the ND (NC) phase

is favored. Here, again, for the intermediate level of the relation between the alkyl-chain length of the surfactant and the alcohol, the NB is stabilized. It was also observed that the *n* = *m* ± 2 rule is still applicable, except that the concentration of the surfactant, which is the main component of the lyotropic mixtures, is bigger than that of alcohol.

## *3.2. Effect of Interactions between Head Groups and Ions of Electrolytes*

Lyotropic systems are also known as micellar systems because their building blocks (or basic units) are micelles. The micellar systems consist of three main regions: (a) intermicellar region, (b) interfacial region and (c) micelle core. The intermicellar region includes water, some counterions of the surfactants and, if an electrolyte is added to the mixture, the ions of the electrolytes, and free surfactant molecules. The micelle cores consist of the hydrocarbon part, i.e., the alkyl chain of the surfactants. The interfacial region is a region at the micelle surfaces where the head groups of the surfactants, some counterions, ions of electrolytes and structured water are present. In the previous section, we have dealt with the modifications in the micelle core by examining the relative effect of the alkyl-chain lengths of the surfactant and of the alcohol. In this section, we discuss the effect of the modifications at the micelle surfaces, i.e., in the interfacial region.

The interaction of alkali ions with the head groups of the surfactant molecules on obtaining the different nematic phases, especially biaxial one, was investigated [88]. In this study, two different lyotropic mixtures were examined, KL/alkali sulfate/undecanol/water, as shown in Figure 6, and sodium dodecylsulfate (SDS)/alkali sulfate/dodecanol/water, as shown in Figure 7, where alkali sulfates are Li2SO4, Na2SO4, K2SO4, Rb2SO4 and Cs2SO4. In Figure 6 the optical birefringences are given. The refraction indices along the x, y and z axes of the laboratory coordinate system are *nx*, *ny* and *nz*, and the birefringences are defined as: *δn* = *nz* − *ny* ; <sup>Δ</sup>*<sup>n</sup>* <sup>=</sup> <sup>|</sup>*nx* <sup>−</sup> *nz*|. The magnetic field is applied along the z axis. In the ND phase, *nx* = *nz* = *ny*, in the NC phase *nx* = *ny* = *nz*, and in the NB phase *nx* = *nz* = *ny*.

**Figure 6.** Temperature dependences of the birefringences of the mixtures KL/alkali sulfate/UndeOH/water, including (**a**) Na2SO4, (**b**) K2SO4, (**c**)Rb2SO4 and (**d**) Cs2SO4, separately. (-) and (•) represent δn and Δn, respectively [88].

**Figure 7.** Temperature dependences of the birefringences of the mixtures SDS/alkali sulfate/DDeOH/water, including (**a**) Li2SO4, (**b**) Na2SO4, (**c**) K2SO4, (**d**)Rb2SO4 and (**e**) Cs2SO4, separately. (-) and (•) represent δn and Δn, respectively [88].

As seen in Figure 6, except for Na2SO4, three nematic phases were obtained for K2SO4, Rb2SO4 and Cs2SO4 and the biaxial temperature ranges were about the same, ~5–6 ◦C. The mixture with Na2SO4 has higher birefringence in the ND region (about twice higher) with respect to those from the other salts. Although we did not observe the biaxial regions for SDS mixtures, as shown in Figure 7, we extracted some useful information from those samples, which were in good agreement with the results of the KL mixtures. While the sample with Li2SO4 and Na2SO4 exhibited maximum birefringence of ~4.0 × <sup>10</sup>−<sup>3</sup> and ~5.0 × <sup>10</sup>−3, respectively, K2SO4, Rb2SO4 and Cs2SO4 presented similar birefringence values of ~6.0 × <sup>10</sup><sup>−</sup>3, in the ND phases of the lyotropic host mixture SDS/DDeOH/water. The ~50% increase in the birefringence value, going from Li2SO4 to K2SO4, was really significant because the birefringences are related to the micelle size and shape anisotropy. The higher the birefringences the bigger the micelle shape anisotropy. Moreover, the amount (or degree) of structured water around the ions of the salts provides a different electrostatic capability to them against the head groups of the surfactants in the micelle interface that plays an important role on the stabilization of the different nematic phases. Indeed, this effect arises from the chaotropic and kosmotropic properties of the ions/counterions. Small ions have relatively high surface-charge density [89] and they exhibit high tendency to be hydrated by a large amount of free water molecules (water-structuring or kosmotropic ions [90–92]) with the negatively greater enthalpy of hydration with respect to bigger ions [90–92]. On the contrary, big ions show opposite behaviour with respect to small ones and they are less hydrated (water-breaking or choatropic ions [90–92]). The kosmotropic and chaotropic properties of the ions are in close relation with the hydrodynamic radius of the ions, *R*. Comparing two ions with different physical (ionic) radius, *rs* and *rb* (*rs* < *rb*), and same electrostatic charge, as stated above, smaller ions are surrounded by a higher amount of structured water molecules with respect to the bigger ions, as shown in Figure 8. If two ions with *rs* = *rb* and *R*<sup>s</sup> ≈ *R*b, it is expected that they show similar kosmotropic or chaotropic character. This can be seen in our experimental results given in Figures 6 and 7. While the sequential ordering of the ionic radii of the alkali ions is Li<sup>+</sup> (0.59) < Na<sup>+</sup> (0.99) < K<sup>+</sup> (1.37) < Rb<sup>+</sup> (1.52) < Cs<sup>+</sup> (1.67), their hydrodynamic radii sequence is Li<sup>+</sup> (3.40) >Na+ (2.76) >K+ (2.32) ≈ Rb<sup>+</sup> (2.28) ≈ Cs+ (2.28),

where the numbers between parentheses are the ionic and hydrodynamic radius of the ions in Å [93]. In Figures 6 and 7, it was observed that K2SO4, Rb2SO4 and Cs2SO4 affected the phase topology of the nematic-host phases in the similar way because their alkali ions have almost same hydrodynamic radii. At first sight, it seems that the birefringences of the nematic phases could be related to the hydrodynamic radius of ions. Although, for instance, K+ ion has smaller hydrodynamic radius than Na+, the mixture with K+ has lower birefringences in the ND phase with respect to the one with Na+ in KL system (Figure 6a,b). However, in the case of SDS system, as shown in Figure 7b,c, there is an opposite situation for the same ions. Remembering that KL is kosmotropic and SDS is chaotropic type surfactants, the head groups of the surfactants KL and SDS interact with strongly kosmotropic Na<sup>+</sup> and chaotropic K+ ions, respectively, to produce strongly bound ion-pairs at the micelle surfaces, which leads to the formation of larger micelles. X-ray diffraction studies indicated that the larger the micelles the higher the birefringences [15]. Consequently, it would be better to conclude that the larger value of the birefringences is attributed to the formation of tightly bound ion-pairs at the micelle surfaces.

**Figure 8.** Comparative representation of ionic (2*r*) and hydrodynamic (2*R*) diameter of small kosmotropic (**a**) and the bigger chaotropic ions (**b**) with the same electrostatic charge. dw is the thickness of the water layer.

The results discussed in [88] may be considered as preliminary. For this reason, complementary investigations were needed to understand the effect of the specific interactions between the ionic species at the micelle surfaces on the stabilization of the lyotropic nematic phases, especially the biaxial. Sodium salts of some Hofmeister ions were added [12] to the host lyotropic mixture of dodecyltrimethylammonium bromide (DTMABr)/dodecanol (DDeOH)/water. The results showed that the interactions between the positively charged head group of DTMABr, trimethylammonium, and the Hofmeister series kosmotropic (SO4 <sup>2</sup>−, F−, Cl−), and chaotropic (Br−, NO3 −, I−, ClO4 −) anions play an important role on stabilizing the lyotropic nematic phases. The relatively strong interactions of the trimethylammonium head group, which has a chaotropic property, with the most chaotropic anions ClO4 − promoted the stabilization of just the ND phase. The higher the kosmotropic (chaotropic) character of the anions the bigger the NC (ND) temperature domain.

Although the results of [12] provided useful results about the stabilization the different nematic phases, it was necessary to generalize the role of specific interactions on the micelle's surfaces in terms of the kosmotropic and chaotropic properties of both the surfactant head groups and the counterions/ions. For this purpose, three surfactant molecules were selected, two of them with chaotropic head groups (one positively charged, the other negatively charged) and one with a negatively charged kosmotropic head group. The chaotropic surfactants were sodium dodecylsulfate (SDS) and tetradecyltrimethylammonium bromide (TTMABr), and the kosmotropic one potassium laurate (KL). Three mixtures were prepared, SDS/alkali salt/dodecanol/water, TTMABr/sodium salt/DeOH/water and KL/alkali salt/DeOH/water. The alkali salts (sodium salts) were Li2SO4, Na2SO4, K2SO4, Rb2SO4 and Cs2SO4 (Na2SO4, NaF, NaCl, NaBr, NaNO3, NaClO4 and NaSCN). The results [14] indicated that when the surfactant head groups and the ions bound to them at the micelle's surface exhibit the same property in terms of the chaotropic or kosmotropic character, the ionic species shows highly-strongly bound ion pairs (Figure 9a), which lead to the stabilization of the ND phase. If they form highly-loosely bound ion pairs (Figure 9b), the NC phase is stabilized. For the ion pairs formed by the intermediate level of the interactions between the head groups and the ions, the NB is favored, as shown in Figure 9c.

**Figure 9.** Formation of the ion pairs, considering the interactions between the positively charged kosmotropic surfactant head group and (**a**) kosmotropic, (**b**) chaotropic, (**c**) weakly chaotropic or kosmotropic counterion/ion. This phenomenon was described by Moriera and Froozabadi in Ref [94]. For a surfactant with chaotropic head group, similar situation is in question, except that it produces tightly bound ions pairs with a chaotropic ions and loosely bound ion pairs with a kosmotropic ions, as in the (**a**) and (**b**), respectively.

At this point, it would be better to discuss how the surfactant molecules aggregate in the micelles as a result of the formation of tightly/loosely bound ion-pairs at the micelle surfaces. When an ion is strongly bound to the head group of a surfactant, it screens the Coulombic repulsions between the head groups and the surfactant molecules are packed well in the micelles. In this case, the micelle surface curvature is flatter with respect to the one appeared with loosely interacting head group-ion pairs. This situation applies to strong kosmotrope-kosmotrope or chaotrope-chaotrope interactions between the ionic species. As the surfactant head groups and the ions presence the opposite kosmotrope/chaotrope characters, the repulsions between the head groups are screened less and the micelle surface curvature starts to increase. Thus, in the case of lyotropic nematic phases, from NC to ND phase passing through the NB phase by changing the temperature, the surfactant head groups are packed well at the micelle surfaces leading to the growth in the micelle size and the micelle surface curvature decreases.

## *3.3. Effect of Localization of Weak Electrolytes at Micelle's Surface*

Another effect of the modifications on the micelle's surface in the stabilization of uniaxial and biaxial nematic phases is the location of the dopant molecules (e.g., strong or weak electrolytes) on the micelle's surface. It is known that when a strong electrolyte is added to a lyotropic mixture, its oppositely charged ions with respect to the surfactant head groups produce highly bound ion pairs as a result of strong Coulombic attractions on the micelle's surface. This gives rise to the packing of the surfactant molecules within the micelles and the growth of the micelle size in the A'-B' plane of the orthorhombic micelles (Figure 4a). This process leads to the stabilization of ND and/or NB phase. In the case of weak electrolytes, as expected, the interactions of polar parts of these electrolytes with the surfactant head groups are not as strong as those of ions from the strong electrolytes. At first sight, it might be assumed that this situation is due to the difference in the solubilities of strong and weak electrolytes. We chose some weak electrolytes with –OH and –COOH polar parts [13] to investigate this hypothesis. The experimental results indicated that there exists no direct relation between the solubilities of those weak electrolytes in water and their effectiveness for stabilizing the uniaxial and biaxial nematic phases. In contrast to the solubility, the acidity constants, Ka, of the dopant molecules (i.e., weak electrolytes) emerged as an important parameter in obtaining different types of nematic phase. If a dopant molecule with pKa < 7 is chosen, it is highly possible to stabilize three nematic phases. However, in the case of pKa > 7 and exactly the same amount of dopant with pKa < 7, the mixture with the dopant molecule exhibit the same nematic phase type that the mixture without the dopant. In other words, the dopant with pKa > 7 does not change the type of the nematic phase of the host mixture. This situation was attributed to the different location of the dopant molecules in the micelles and their effectiveness of screening the repulsions between the ionic head groups of the surfactants.

## **4. Conclusions**

These recent studies have provided advances on the understanding of the biaxial nematic phase stabilization. We may conclude that the choosing the optimum mixture constituents is a key point to stabilize different nematic phases. Parameters as the relative alkyl chain lengths of surfactant and alcohols and occurrence of strong or weak interactions between the ionic species at the micelle surfaces have to be considered in the preparation of a lyotropic mixture aiming the obtaining of particular nematic phases. The following conclusions may be helpful for experimentalists to obtain the lyotropic mixture presenting uniaxial and/or biaxial nematic phases:


"KL was synthesized by neutralization of lauric acid with KOH". Here, the term "neutralization" does not include the control of the pH. However, Berejnov et al. showed that the neutralization of all lauric acid to give KL is a crucial point to have reproducible results for any mixture of the surfactant KL [95]. To do so, the reaction pH is kept at ~10.8. Otherwise, if the pH is less than 10.8, some amount of lauric acid may remain in the solution without neutralizing, and then excess of lauric acid and DeOH presence in the lyotropic mixture may cause the esterification of lauric acid. Melnik and Saupe proposed a slow esterification reaction between the alcohol and KL; however, most probably this reaction occurs between the alcohol and lauric acid, instead of the soap [81]. Thus, we can say that the experimentalists have to consider the purity of the KL to have chemically stable lyotropic mixtures.

Consequently, if we have a chance to give "a receipt" to experimentalists for preparation of a lyotropic mixture presenting the biaxial phase, the followings may be considered:


**Funding:** This study was funded by the Scientific and Technological Research Council of Turkey (TÜB˙ ITAK) [grant number: 113Z469 and 114Z031], CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) (465259/2014-6), FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo) (2014/50983-3 and 2016/24531-3), CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), INCT-FCx (Instituto Nacional de Ciência e Tecnologia de Fluidos Complexos), and NAP-FCx (Núcleo de Apoio à Pesquisa de Fluidos Complexos).

**Conflicts of Interest:** The authors declare no conflicts of interest.

## **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Magnetic Field and Dilution Effects on the Phase Diagrams of Simple Statistical Models for Nematic Biaxial Systems**

## **Daniel D. Rodrigues, André P. Vieira \* and Silvio R. Salinas \***

Instituto de Fisica, Universidade de Sao Paulo, Rua do Matao, 1371, Sao Paulo 05508-090, Brazil; daniel.dias.rodrigues@usp.br

**\*** Correspondence: apvieira@if.usp.br (A.P.V.); ssalinas@if.usp.br (S.R.S.)

Received: 19 June 2020; Accepted: 20 July 2020; Published: 22 July 2020

**Abstract:** We use a simple statistical model to investigate the effects of an applied magnetic field and of the dilution of site elements on the phase diagrams of biaxial nematic systems, with an emphasis on the stability of the Landau multicritical point. The statistical lattice model consists of intrinsically biaxial nematogenic units, which interact via a Maier–Saupe potential, and which are characterized by a discrete choice of orientations of the microscopic nematic directors. According to previous calculations at zero field and in the absence of dilution, we regain the well-known sequence of biaxial, uniaxial, and disordered structures as the temperature is increased, and locate the Landau point. We then focus on the topological changes induced in the phase diagram by the application of an external magnetic field, and show that the Landau point is destabilized by the presence of an applied field. On the other hand, in the absence of a field, we show that only a quite strong dilution of nematic sites is capable of destabilizing the Landau point.

**Keywords:** biaxial nematic transition; field behavior; diluted nematic systems

## **1. Introduction**

In a number of recent investigations [1–6], we have performed some calculations for simple statistical lattice models to characterize the biaxial structures in liquid crystalline systems. These statistical models are based on fully-connected Maier–Saupe pair interactions, with a restricted choice of orientational degrees of freedom. They are amenable to detailed calculations, and have been shown to account for most of the qualitative features of the nematic phase diagrams, including the well-known sequences of biaxial nematic, uniaxial nematic, and isotropic structures, as the temperature is raised. At a special value of temperature and of a parameter gauging the degree of biaxiality, model calculations point out the existence of a highly symmetric Landau multicritical point, with a direct transition from a biaxial nematic to an isotropic phase.

In this work, we revisit the same simplified lattice models with a view to perform detailed calculations to investigate the effects produced by an applied field and by the dilution of nematogenic elements. Besides assuming a generalized form of the Maier–Saupe interactions between nematogenic units on a lattice site, we still make a discrete choice of the orientations of the microscopic nematic directors. This special choice, which is reminiscent of an old work of Zwanzig to treat the Onsager model of rigid cylinders, has also been shown to lead to the well-known first-order transition between simple uniaxial nematic and isotropic phases.

The effects of a magnetic field on the phase diagram of nematogenic molecules exhibiting biaxial phases have been the subject of a number of investigations over the years. For micellar lyotropic systems such as potassium laurate/1-decanol/water mixtures, in which shape anisotropy of the micelles is presumably strongly dependent on the temperature and on the concentration of the various components, the existence of a biaxial nematic phase is known since the pioneering work of Yu and Saupe [7]. Besides the biaxial phase, these systems also exhibit rod-like and disk-like uniaxial nematic phases, and, in general, upon heating, these uniaxial phases exhibit a discontinuous transition to an isotropic phase. Nevertheless, the three nematic phases and the isotropic phase become equal at the Landau multicritical point. In the neighborhood of the Landau point, experiments performed for three different values of concentration [8] show that application of a magnetic field turns the transition between the rod-like uniaxial and the isotropic phases into a crossover, while preserving the discontinuous nature of the transition between the disk-like uniaxial and the isotropic phases.

For mineral lyotropic systems, the combination of an intrinsic magnetic moment and a negative anisotropy of the magnetic susceptibility leads to a very strong induced magnetic birefringence [9] and to phase separation into a uniaxial nematic and a biaxial nematic due to entropic effects associated with polydispersity [10].

In thermotropic systems, prompted by experimental results [11] for bent-core, intrinsically biaxial molecules with an essentially fixed shape, there have been studies on the combined effect of biaxiality and a magnetic field on the temperature of the first-order uniaxial to a paranematic phase transition, on the stability of the various phases, and on the presence of multicritical points [11–16]. In particular, within the phenomenological Landau–de Gennes theory, Mukherjee and Rahman [13] have identified the destabilization of the Landau point and its splitting into one critical and two tricritical points.

In the present work, we go beyond the analysis of a Landau expansion. As we start from a microscopic statistical model, we are able to perform calculations to analyze the global phase diagram, which includes considerations on the stability of the Landau multicritical point and also illustrations of the behavior of the uniaxial and biaxial order parameters as temperature and degree of biaxiality are changed.

On the other hand, our motivation for investigating the effects of dilution is purely theoretical, as we would like to check whether the phase diagram of our simple model survives in the presence of an ingredient that is inescapable in the real world. Our results show that this is indeed the case, making clear that representing a fluid by a lattice model does not introduce any artifacts in our obtained phase diagrams.

In Section 2, we define the statistical model in terms of a quite general and elegant two-tensor formalism, proposed [17] by Sonnet, Virga, and Durand (SVD), and which includes all physically reasonable pair interactions between nematogenic elements. With a discrete choice of the microscopic directors, this model Hamiltonian is restricted to a six-state model, which we dub SVD6. Furthermore, we choose a relation between two interaction parameters, the "geometric-mean condition", which preserves the main features of the phase diagrams, and which has been widely used in the area. We obtain a number of results for the phase diagram of this SVD6 model, in terms of temperature *T* and a parameter Δ, which comes from the geometric-mean condition and provides a form of gauging the degree of biaxiality. In particular, we locate the Landau multicritical point. In Section 3, we discuss the inclusion of an external field in the SVD6 Hamiltonian with the geometric-mean condition. The effects of a magnetic field come from an intricate competition between the magnetic free energy, the entropy and the interactions between nematogens, which can be fully accounted for by the elementary statistical model. In general, the presence of an external field breaks the symmetry necessary to stabilize the Landau point. Except for the very special circumstances of a fully isotropic diamagnetic susceptibility, the main result is the splitting of the transition lines and the disappearance of the Landau multicritical point. We then turn to the study of the effects of dilution. In Section 4, we use the SVD6 model, with the recourse to the geometric-mean condition, to show that the presence of dilution, as far as it is not a strong dilution, does not introduce any qualitative changes with respect to the initial phase diagrams. However, in the zero-field limit and for quite large dilutions, the Landau point turns into a first-order transition, which is similar to the behavior of the critical line in a strongly diluted Ising ferromagnet.

## **2. The SVD Model**

The Hamiltonian representing the interactions between the nematogens in the mean-field SVD model [17] can be written as

$$\mathcal{H}\_{\text{SVD}} = -\frac{A}{2V} \sum\_{\mu,\nu \in \{\text{x,y,z}\}} \left[ \left( \sum\_{i=1}^{V} q\_i^{\mu\nu} \right)^2 + 2\gamma \left( \sum\_{i=1}^{V} q\_i^{\mu\nu} \right) \left( \sum\_{i=1}^{V} b\_i^{\mu\nu} \right) + \lambda \left( \sum\_{i=1}^{V} b\_i^{\mu\nu} \right)^2 \right],\tag{1}$$

in which the parameter *A* > 0 gauges the interaction between nematogens, the dimensionless parameters *γ* and *λ* are responsible for a dependence of that interaction on the relative orientation of the nematogens, and *V* is the number of nematogens, while *q μν <sup>i</sup>* and *b μν <sup>i</sup>* , with *μ*, *ν* ∈ {1, 2, 3}, are the components of two traceless tensors associated with the *i*th nematogen. Explicitly, for a generic nematogen with principal axes along the unit vectors *n*ˆ 1, *n*ˆ <sup>2</sup> and *n*ˆ 3, we have

$$\mathbf{q} = \frac{3}{2}\hbar\_1 \otimes \hbar\_1 - \frac{1}{2}\mathbf{I} \qquad \text{and} \qquad \mathbf{b} = \frac{3}{2}\left(\hbar\_2 \otimes \hbar\_2 - \hbar\_3 \otimes \hbar\_3\right),\tag{2}$$

**I** being the 3 × 3 identity matrix.

For our simplified model, we restrict the directors to point along the Cartesian axes [18,19], so that each nematogen takes only six distinct states, defining a new Hamiltonian

$$\mathcal{H}\_{\text{SVD}} = -\frac{A}{2V} \sum\_{a \in \{x, y, z\}} \left[ \left( \sum\_{i=1}^{V} q\_i^{aa} \right)^2 + 2\gamma \left( \sum\_{i=1}^{V} q\_i^{aa} \right) \left( \sum\_{i=1}^{V} b\_i^{aa} \right) + \lambda \left( \sum\_{i=1}^{V} b\_i^{aa} \right)^2 \right]. \tag{3}$$

Labeling the states of a generic nematogen as *<sup>ζ</sup>*(*k*), *<sup>k</sup>* <sup>∈</sup> {1, 2, 3, 4, 5, 6}, the values of the diagonal components of the tensors **q** and **b**, as well as of the *z* component of the first principal axis are given in Table 1. Further restricting ourselves to the case *λ* = *γ*2, the widely used geometric-mean condition [20,21], Equation (3) can be written as

$$\mathcal{H}\text{svd}\mathfrak{a} = -\frac{A}{2V} \sum\_{a \in \{x, y, z\}} \left[ \sum\_{i=1}^{V} \left( q\_i^{aa} + \gamma b\_i^{aa} \right) \right]^2. \tag{4}$$

As further discussed in Section 3, the parameter *γ* appearing in the last expression can be interpreted as a simplified measure of the intrinsic biaxiality of the nematogens. The effect of *γ* on the interaction energies under different relative orientations of a pair of nematogens replaces the entropic effects due to the different excluded volumes associated with these orientations for a pair of hard biaxial polyhedra, which, in the case of strong shape anisotropy, can also lead to the spontaneous appearance of stable biaxial phases [22].

**Table 1.** Values of the diagonal components of the tensors *q* and *b* and of the *z* component of the principal axes of a generic nematogen in the six possible states *ζ*(*k*) of the restricted model.


*Crystals* **2020**, *10*, 632

## **3. Effects of a Magnetic Field**

In order to consider the effects of a magnetic field, we add to the Hamiltonian of our simplified model the interactions of the nematogens with a uniform and constant field *H*- = *Hz*ˆ applied along the *z* axis in the laboratory frame [23]. We then write

$$\mathcal{H} = \mathcal{H}\_{\text{SVD6}} - \frac{1}{2} \sum\_{i=1}^{V} \sum\_{j=1}^{3} \chi\_{j} \left( H \cdot n\_{j,i} \right)^{2} \,, \tag{5}$$

in which *χ<sup>j</sup>* gauges the coupling with the magnetic field when the *j*th principal axis *n*ˆ*<sup>j</sup>* of an object lies along the field direction, and the index *i* runs over all objects. In view of the fact that, even under the geometric-mean condition, there are in general different interaction energies associated with distinct relative orientations of the principal axes of a pair of nematogens, no symmetry should be expected between field effects associated with choices of the sets of parameters {*χ*1, *χ*2, *χ*3} such as {*χ*, 0, 0}, {0, *χ*, 0} and {0, 0, *χ*}.

The partition function *Z* of the model can be evaluated by introducing an auxiliary diagonal tensor **Q** = diag *Qxx*, *Qyy*, *Qzz* via three Gaussian identities which allow us to eliminate the quadratic terms arising from Equation (4), yielding

$$Z = \sum\_{\{\tilde{\zeta}\}} e^{-\beta \mathcal{H}} = \left(\frac{\beta AV}{2\pi}\right)^{\tilde{\beta}/2} \sum\_{\{\tilde{\zeta}\}} \int \left[dQ\_{aa}\right] \exp\left\{-\frac{\beta AV}{2} \sum\_{a} Q\_{aa}^2\right.\tag{6}$$

$$+ \sum\_{a,i} \beta A Q\_{aa} \left(q\_i^{aa} + \gamma b\_i^{aa}\right) + \frac{\beta H^2}{2} \sum\_{i,j} \chi\_j \left(n\_{j,i}^z\right)^2\right\},\tag{7}$$

in which [*dQαα*] = *dQxxdQyydQzz* and *β* = *T*−<sup>1</sup> is the inverse temperature (we use a system of units for which Boltzmann's constant *kB* = 1). Evaluating the sums over the nematogen index *i* and the states of each nematogen, we obtain

$$Z = \left(\frac{\beta AV}{2\pi}\right)^{3/2} \int \left[dQ\_{\rm rad}\right] \exp\left\{-\beta V\psi\left(\mathbf{Q}\right)\right\},\tag{7}$$

with

$$\psi\left(\mathbf{Q}\right) = \frac{A}{2} \sum\_{a} Q\_{aa}^{2} - T \ln 2 - T \ln \phi\left(\mathbf{Q}\right) \tag{8}$$

and

$$\begin{split} \boldsymbol{\Phi} \cdot \boldsymbol{\Phi} \left( \mathbf{Q} \right) &= \varepsilon^{\frac{3}{2} \beta A Q\_{xx} + \frac{1}{4} \beta (\chi\_2 + \chi\_3) H^2} \cosh \left[ \frac{3 \beta A \gamma}{2} \left( Q\_{yy} - Q\_{zz} \right) - \frac{\beta \left( \chi\_2 - \chi\_3 \right) H^2}{4} \right] \\ &+ \varepsilon^{\frac{3}{2} \beta A Q\_{yy} + \frac{1}{4} \beta (\chi\_2 + \chi\_3) H^2} \cosh \left[ \frac{3 \beta A \gamma}{2} \left( Q\_{zz} - Q\_{xx} \right) + \frac{\beta \left( \chi\_2 - \chi\_3 \right) H^2}{4} \right] \\ &+ \varepsilon^{\frac{3}{2} \beta A Q\_{xx} + \frac{1}{2} \beta \chi\_1 H^2} \cosh \left[ \frac{3 \beta A \gamma}{2} \left( Q\_{xx} - Q\_{yy} \right) \right]. \end{split} \tag{9}$$

The partition function can be evaluated by means of the steepest descent method. The equilibrium states correspond to those values **Q** of the tensor **Q** that minimize the functional *ψ* (**Q**). As these values are such that **Q** is traceless, we can write

$$
\langle \mathbf{Q} \rangle = \begin{pmatrix}
0 & -\frac{1}{2} \begin{pmatrix} S - \eta \end{pmatrix} & 0 \\
0 & 0 & S
\end{pmatrix},
\tag{10}
$$

in which *S* and *η* play the roles of order parameters associated with the uniaxial and biaxial nematic phases, respectively. Noticing that at the mean-field level **Q** can be identified with *q* + *γb*, we also have

$$
\langle \mathbf{Q} \rangle = \frac{3}{2} \left[ \langle \hbar\_1 \otimes \hbar\_1 \rangle + \gamma \left( \langle \hbar\_2 \otimes \hbar\_2 \rangle - \langle \hbar\_3 \otimes \hbar\_3 \rangle \right) \right] - \frac{1}{2} \mathbf{I}.\tag{11}
$$

Substituting Equation (10) in Equation (8) we obtain

$$\psi\left(\mathcal{S},\eta\right) = \frac{3A}{4}\mathcal{S}^2 + \frac{A}{4}\eta^2 - T\ln 2 - T\ln \phi\left(\mathcal{S},\eta\right),\tag{12}$$

with

$$\begin{split} \phi\left(S,\eta\right) &= e^{-\frac{3}{4}\beta A\left(S+\eta\right) + \frac{1}{4}\beta\left(\chi\_{2}+\chi\_{3}\right)H^{2}} \cosh\left[\frac{\beta A\Delta}{4}\left(\eta - 3S\right) - \frac{\beta\left(\chi\_{2}-\chi\_{3}\right)H^{2}}{4}\right] \\ &+ e^{-\frac{3}{4}\beta A\left(S-\eta\right) + \frac{1}{4}\beta\left(\chi\_{2}+\chi\_{3}\right)H^{2}} \cosh\left[\frac{\beta A\Delta}{4}\left(\eta + 3S\right) + \frac{\beta\left(\chi\_{2}-\chi\_{3}\right)H^{2}}{4}\right] \\ &+ e^{\frac{3}{4}\beta A\mathcal{S} + \frac{1}{2}\beta\chi\_{1}H^{2}} \cosh\left(\frac{\beta A\Delta}{2}\eta\right). \end{split} \tag{13}$$

Notice that in the above equations we introduced a biaxiality parameter Δ = 3*γ* so that we can compare our phase diagrams with the one obtained at zero field by Boccara, Medjani and De Seze [24]. If, for the sake of comparison, we assume that the traceless part of the inertia tensor and the *q* + *γb* tensor of an object are proportional to each other, we may interpret that, as Δ is increased, the shape of the nematogen changes from calamitic ("rod-like") uniaxial for Δ = 0, to strongly biaxial for Δ = 1, and finally to discotic ("disk-like") uniaxial for Δ = 3 (see the appendix in [4]). In this paper, we always work in the interval 0 ≤ Δ ≤ 3.

Minimizing Equation (12) we arrive at the state equations

$$\begin{aligned} \frac{\partial \psi}{\partial S} &= \frac{3A}{2} S - T \frac{\partial \ln \phi}{\partial S} = 0\\ \frac{\partial \psi}{\partial \eta} &= \frac{A}{2} \eta - T \frac{\partial \ln \phi}{\partial \eta} = 0 \end{aligned} \tag{14}$$

If there are multiple solutions to Equation (14), the equilibrium state corresponds to the solution with the smallest value of *ψ* (*S*, *η*). At high temperatures and at zero field, there is an isotropic, disordered solution *S* = *η* = 0. At lower temperatures, depending on the choice of Δ and *χH*2, there may appear uniaxial solutions, for which *η* = 0 = *S*, and biaxial solutions, for which *S* = 0 and *η* = 0. At zero field, due to the six-fold symmetry of the problem, there are degenerate solutions for the ordered phases, so that for every uniaxial solution (*S*, *η*) = (*S*0, 0) there are equivalent solutions (*S*, *η*) = −1 <sup>2</sup>*S*0, <sup>±</sup><sup>3</sup> <sup>2</sup>*S*<sup>0</sup> . Of course, these do not correspond to biaxial solutions, which are not three-fold but six-fold degenerate. In fact, for each pair of biaxial solutions (*S*0, ±*η*0), with *η* = 0, there are equivalent solutions −1 <sup>2</sup>*S*<sup>0</sup> <sup>∓</sup> <sup>1</sup> <sup>2</sup> *<sup>η</sup>*0, <sup>3</sup> <sup>2</sup>*S*<sup>0</sup> <sup>∓</sup> <sup>1</sup> <sup>2</sup> *η*<sup>0</sup> and −1 <sup>2</sup>*S*<sup>0</sup> <sup>±</sup> <sup>1</sup> <sup>2</sup> *<sup>η</sup>*0, <sup>−</sup><sup>3</sup> <sup>2</sup>*S*<sup>0</sup> <sup>∓</sup> <sup>1</sup> <sup>2</sup> *η*<sup>0</sup> . In the biaxial phase, comparison of Equations (10) and (11) in the *T* → 0 limit indicates that the zero-field values of *η* range from *η* = − max {Δ,(Δ + 3)/2} to *η* = max {Δ,(Δ + 3)/2}. When plotting the values of *S* and *η* in Figures 1–3, we select for each choice of Δ and *T* the solution that provides the clearest distinction between the phases.

Phase boundaries and multicritical points can be determined by the behavior of *ψ* (*S*, *η*) and its derivatives at the solutions of Equation (14).

## *3.1. Behavior at Zero Field*

We first offer a short review of the phase diagram at zero field (*H* = 0), shown in Figure 1 and already discussed in [4,5,24].

**Figure 1.** Phase diagram at zero field, in the *T* × Δ plane. The color code shows the behavior of (**a**) the uniaxial order parameter *S*, which is zero in the high-temperature isotropic phase (I), and (**b**) the biaxial order parameter *η*, which is nonzero only in the low-temperature biaxial (B) phase. Notice that at low and intermediate temperatures the value of *S* associated with the smallest equilibrium value of |*η*| changes sign along the line Δ = 1 on which the Landau point (Δ*L*, *TL*) = (1, 1) is located. Solid (dashed) lines indicate continuous (first-order) transitions. The calamitic (U+) and discotic (U-) uniaxial phases are also indicated.

When Δ = 0, nematogens are intrinsically (calamitic) uniaxial, so that no biaxial phase is expected. Therefore, there is only a first-order transition between the uniaxial and the isotropic phases, which occurs at the temperature for which

$$\frac{\partial \boldsymbol{\psi}}{\partial S}\bigg|\_{\left(S\_{0},0\right)} = 0 \quad \text{and} \quad \boldsymbol{\psi}\left(S\_{0},0\right) = \boldsymbol{\psi}\left(0,0\right). \tag{15}$$

As the biaxiality parameter Δ is increased, a biaxial solution appears at a sufficiently small temperature. For fixed Δ, the temperature corresponding to the continuous transitions between the biaxial and the uniaxial phases is obtained from

$$\left. \frac{\partial \psi}{\partial S} \right|\_{\left(\mathcal{S}\_1, \mathcal{O}\right)} = \left. \frac{\partial^2 \psi}{\partial \eta^2} \right|\_{\left(\mathcal{S}\_1, \mathcal{O}\right)} = 0. \tag{16}$$

The lines of first-order uniaxial-isotropic and continuous biaxial-uniaxial transitions meet at the Landau multicritical point (Δ*L*, *TL*), the location of which is determined by the conditions

$$
\left.\frac{\partial^2 \psi}{\partial S^2}\right|\_{\left(0,0\right)} = \left.\frac{\partial^3 \psi}{\partial S^3}\right|\_{\left(0,0\right)} = 0,\tag{17}
$$

derived by noticing that at that point *S*<sup>0</sup> and *S*<sup>1</sup> both tend smoothly to zero. For *H* = 0, the above equations reduce to

$$
\beta A \left(\Delta^2 + 3\right) - 4 = \beta^2 A^2 \left(\Delta^2 - 1\right) = 0,\tag{18}
$$

so that, measuring *T* in units of the energy scale *A*, the solution to Equation (17) is

$$\left(\Delta\_{L}, T\_{L}\right) = \left(1, 1\right). \tag{19}$$

Further increasing the value of Δ, the temperature of the biaxial–uniaxial transition decreases, becoming zero again at Δ = 3, when the nematogen is once more a (discotic) uniaxial object. At the intermediate-temperature uniaxial phases, the uniaxial order parameter is positive (*S*<sup>0</sup> > 0) if Δ < Δ*<sup>L</sup>* and negative (*S*<sup>0</sup> < 0) if Δ > Δ*L*.

We note that this phase diagram qualitatively agrees with the one obtained from computer simulations of a hard biaxial ellipsoid fluid, as a function of the shape anisotropy of the particles [25]. It also reproduces various aspects of phase diagrams of a variety of lyotropic systems [7,26–29], regarding their nematic phases.

## *3.2. Phase Diagrams in a Field*

The effect of a magnetic field is the result of an intricate competition between the magnetic free energy, the entropy and the mutual interactions between nematogens, as dictated by the geometric-mean condition.

The first noticeable effect of a nonzero field is the disappearance of the Landau point, except for the fully symmetric choice *χ*<sup>1</sup> = *χ*<sup>2</sup> = *χ*3, which makes all three principal axes of a nematogen equally likely to lie parallel to the field, having the sole effect of globally shifting all energy levels. This can be checked by fixing *H* and one of *χj*, say *χ*2, and solving for Δ, *βA*, *χ*<sup>1</sup> and *χ*<sup>3</sup> the conditions in Equation (17) supplemented by

$$
\left.\frac{\partial\psi}{\partial S}\right|\_{\left(0,0\right)} = 0 \quad \text{and} \quad \left.\frac{\partial^2\psi}{\partial\eta^2}\right|\_{\left(0,0\right)} = 0,\tag{20}
$$

the latter conditions being necessary for a continuous transition from a biaxial phase to a stable isotropic phase at a candidate Landau point. The only solution ever obtained is Δ = 1, *βA* = 1 and *χ*<sup>1</sup> = *χ*<sup>2</sup> = *χ*3, irrespective of the value of the field.

Nonetheless, for the special choice *χ*<sup>2</sup> < *χ*<sup>1</sup> = *χ*3, which makes the *n*ˆ <sup>2</sup> principal axis of a nematogen less likely to lie perpendicular to the field, while making the *n*ˆ <sup>1</sup> and *n*ˆ <sup>3</sup> principal axes equally likely to lie parallel to the field, the isotropic phase is stable at a point fulfilling the conditions in Equation (17), although the second condition in Equation (20) is not satisfied.

As a typical example of the phase diagram in the presence of a field, let us consider the case *χ*<sup>1</sup> ≡ *χ* = 0 and *χ*<sup>2</sup> = *χ*<sup>3</sup> = 0, which is representative of regime |*χ*1| > |*χ*2| , |*χ*3|. With this choice, the isotropic phase does not minimize the free energy at any finite temperature, being replaced by a paranematic phase in which *η* = 0 but *S* takes a small but nonzero value. This can be rationalized as follows. Assuming *χ* > 0, it can be checked that the magnetic energy of a nematogen is minimized when the principal axes *n*ˆ <sup>1</sup> is parallel to the field direction *z*. From Equations (10) and (11), taking into account that # *n*ˆ*<sup>j</sup>* ⊗ *n*ˆ*<sup>j</sup>* \$*zz* = % *nz j* 2 & and noticing that in a nonzero field the symmetry between ' *nz* 1 2 ( , ' (*n<sup>z</sup>* 2) 2 ( and ' *nz* 3 2 ( is explicitly broken, making ' *nz* 3 2 ( = ' (*n<sup>z</sup>* 2) 2 ( < <sup>1</sup> <sup>3</sup> and ' *nz* 1 2 ( > <sup>1</sup> 3 at high temperatures, it is clear that *S* = *Qzz* > 0. Likewise, assuming *χ* < 0, so that the magnetic energy of a nematogen is minimized when the principal axes *n*ˆ <sup>1</sup> is perpendicular to the field direction *<sup>z</sup>*, at high temperatures we have ' *nz* 3 2 ( = ' (*n<sup>z</sup>* 2) 2 ( > <sup>1</sup> <sup>3</sup> and ' *nz* 1 2 ( < <sup>1</sup> <sup>3</sup> , so that *S* = *Qzz* < 0.

The general aspect of the phase diagram for both *χ* > 0 and *χ* < 0 is shown in Figure 2. In both cases, the Landau point disappears to give rise to a pair of multicritical points: a simple critical point signaling the end of a first-order transition line between the uniaxial and the paranematic phases (which are actually a single phase, in much the same way as the liquid and gaseous phases in simple fluids) and a tricritical point separating regions of first-order and continuous transitions between a biaxial phase and the paranematic phase. The dependence of the uniaxial and biaxial order parameters on Δ and *T* is illustrated in Figure 3, for the case *χ* > 0. Notice the extension of the biaxial phase towards the region occupied, at zero field, by the discotic uniaxial phase.

**Figure 2.** Main plots: Phase diagrams in the presence of a field of intensity *H* = 0.1, with the choice *χ*<sup>1</sup> ≡ *χ* = 0 and *χ*<sup>2</sup> = *χ*<sup>3</sup> = 0, both for *χ* = 1 (**a**) and *χ* = −1 (**b**). The Landau point is replaced by a simple critical point *Pc* and a tricritical point *Pt*. The regions corresponding to the paranematic (PN) and the biaxial (B) phases are also indicated. Insets: Dependence of the values of Δ at *Pc* (orange curve) and *Pt* (blue curve) on the product *χH*2. Apart from small shifts produced by the field, a combined parametric *<sup>T</sup>* <sup>×</sup> <sup>Δ</sup> plot obtained from the values of *<sup>T</sup>* and <sup>Δ</sup> at *Pc* and *Pt* as *<sup>χ</sup>H*<sup>2</sup> is varied would closely follow the first-order lines in Figure 1.

**Figure 3.** Phase diagram at a field *H* = 0.15, with the choice {*χ*1, *χ*2, *χ*3} = {1, 0, 0}, in the *T* × Δ plane. The color code shows the behavior of (**a**) the uniaxial order parameter *S*, which is always nonzero in the paranematic phase, which extends up to high temperatures, and (**b**) the biaxial order parameter *η*, which is nonzero only in the biaxial phase, which now extends to intermediate temperatures when Δ > 1. Notice that both *S* and *η* present a steep but nonsingular behavior in the region corresponding, at zero field, to the continuous transition between the biaxial and the discotic uniaxial phases.

The conditions for locating a first-order transition between the uniaxial and the paranematic phases are

$$\left. \frac{\partial \boldsymbol{\Psi}}{\partial \mathcal{S}} \right|\_{\left(\mathcal{S}\_{0}, 0\right)} = \left. \frac{\partial \boldsymbol{\Psi}}{\partial \mathcal{S}} \right|\_{\left(\mathcal{S}\_{0}', 0\right)} = 0 \quad \text{and} \quad \boldsymbol{\Psi}\left(\mathcal{S}\_{0}, 0\right) = \boldsymbol{\Psi}\left(\mathcal{S}\_{0}', 0\right), \tag{21}$$

and at the associated critical point the solutions *S*<sup>0</sup> and *S* <sup>0</sup> become degenerate, so that

$$\left.\frac{\partial^2 \psi}{\partial S^2}\right|\_{\left(S\_0,0\right)} = \left.\frac{\partial \psi}{\partial S}\right|\_{\left(S\_0,0\right)} = 0.\tag{22}$$

A first-order transition between the biaxial and the paranematic phases occurs for

$$
\left.\frac{\partial\psi}{\partial\eta}\right|\_{\left(\mathcal{S}\_2\eta\_2\right)} = \left.\frac{\partial\psi}{\partial\mathcal{S}}\right|\_{\left(\mathcal{S}\_2\eta\_2\right)} = \left.\frac{\partial\psi}{\partial\mathcal{S}}\right|\_{\left(\mathcal{S}\_0'\mathcal{O}\right)} = 0 \quad \text{and} \quad \psi\left(\mathcal{S}\_2,\eta\_2\right) = \psi\left(\mathcal{S}\_0',0\right),
\tag{23}
$$

while a continuous transition between those same phases satisfies

$$\left.\frac{\partial\psi}{\partial S}\right|\_{\left(\mathcal{S}\_0',0\right)} = \left.\frac{\partial^2\psi}{\partial\eta^2}\right|\_{\left(\mathcal{S}\_0',0\right)} = 0.\tag{24}$$

Finally, the location of the tricritical point can be determined from the conditions

$$\frac{\partial \Psi}{\partial S}\Big|\_{\left(S\_{0}^{\prime},0\right)} = \left.\frac{d^{2}\Psi}{d\eta^{2}}\right|\_{\left(S\_{0}^{\prime},0\right)} = \left.\frac{d^{4}\Psi}{d\eta^{4}}\right|\_{\left(S\_{0}^{\prime},0\right)} = 0.\tag{25}$$

Notice that these last conditions involve total derivatives of the free-energy functional with respect to *η*, which means that, while calculating those derivatives, *S* <sup>0</sup> is implicitly treated as a function of *η*.

The multicritical points originating from the destabilization of the Landau point become increasingly separated as the magnetic field becomes more intense, as the insets in Figure 2 show. In the limit of an infinite field, only a continuous biaxial-paranematic phase transition remains.

The existence of a tricritical point on the line of biaxial to paranematic phase transitions is compatible with the predictions of [16], which assumes a more restrictive pair interaction independent of a biaxiality parameter, as well as with those of [6], in which nematogens are intrinsically uniaxial. In the latter case, the located tricritical point is observed for *χ* < 0 and at the value of the magnetic field necessary for our tricritical point to reach the uniaxial limit Δ = 0; see the inset in Figure 2b. This same qualitative behavior for the appearance of a tricritical point in the uniaxial-paranematic transition was predicted for hard rods or hard plates with *χ* < 0 [30].

Similarly, the critical field associated with a critical end point of the first-order transition line between the uniaxial and the paranematic phases, which has been observed for a fixed biaxiality parameter by various authors [6,11,12,15,31], corresponds to the value of the field, which makes our simple critical point reach the uniaxial limit Δ = 0 for *χ* > 0; see the inset in Figure 2a. Again, the same qualitative behavior for the appearance of a critical end point in the uniaxial-paranematic transition was predicted for hard rods or hard plates with *χ* > 0 [32].

Assuming that the biaxiality parameter of a system shows little variation with temperature or the strength of the magnetic field, which would be reasonable for thermotropic liquid crystals, our results suggest that, starting from a biaxial phase and heating the system under a constant field up to the transition to the paranematic phase, the nature of that transition should change from first-order to continuous (or vice-versa, depending on the sign of *χ*) at a certain value of the field, signaling the tricritical point. Estimating the value of that field would require an estimate of the biaxiality parameter Δ for the system, under the assumption that the geometric mean condition is at least approximately applicable. We expect that lyotropic systems, for which presumably there is a significant dependence of the biaxiality parameter on both the temperature and the composition of the mixture [26,27,29], are better candidates for the observation of a tricritical point.

## **4. Introducing Dilution**

Next we turn our attention to a "nematic lattice gas" model consisting of *V* sites, which can be either empty or occupied by a single object among *N* ≤ *V* nematogens. This possibility can be represented at each site *i* by an occupation variable *ti* ∈ {0, 1}, and if the nematogens interact via the the extension of the Hamiltonian in Equation (4), the canonical partition function of the model at zero field can be written as

$$Z = \sum\_{\{t\}}' \sum\_{\{\tilde{\boldsymbol{\zeta}}\}} \exp\left\{ \frac{\beta \boldsymbol{A}}{2V} \sum\_{a \in \{\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}\}} \left[ \sum\_{i=1}^{V} t\_i \left( q\_i^{aa} + \gamma b\_i^{aa} \right) \right]^2 + \frac{\beta}{2} \sum\_{i=1}^{V} \sum\_{j=1}^{3} t\_i \chi\_j \left( \boldsymbol{\varvarprojlim \boldsymbol{\mathcal{H}}} \boldsymbol{\uppi}\_{j,i} \right)^2 \right\},\tag{26}$$

in which the prime in the summation over the set {*t*} of occupation variables indicates that it should be restricted to those configurations for which

$$\sum\_{i=1}^{V} t\_i = N.\tag{27}$$

This restriction can be relaxed by introducing a chemical potential *μ* (not to be confused with the tensor index in Equation (1)) to enforce that the last equation be satisfied on average. We then have to calculate the grand-partition function

$$\Xi = \sum\_{\{t\}} \sum\_{\{\xi\}} \exp\left\{ \frac{\beta A}{2V} \sum\_{a \in \{\mathbf{x}, y, z\}} \left[ \sum\_{i=1}^{V} t\_i \left( q\_i^{aa} + \gamma b\_i^{aa} \right) \right]^2 + \frac{\beta}{2} \sum\_{i=1}^{V} \sum\_{j=1}^{3} t\_i \chi\_j \left( \vec{H} \cdot \hbar\_{j,i} \right)^2 + \beta \mu \sum\_{i=1}^{V} t\_i \right\}. \tag{28}$$

Following the same strategy as in the previous section and introducing an auxiliary diagonal tensor **Q** = diag *Qxx*, *Qyy*, *Qzz* via three Gaussian identities, which allow us to eliminate the quadratic terms in Ξ, we obtain

$$\begin{split} \Xi &= \left(\frac{\beta AV}{2\pi}\right)^{3/2} \sum\_{\{\mathbf{i}\}} \sum\_{\{\mathbf{j}\}} \int \left[dQ\_{aa}\right] \exp\left\{-\frac{\beta AV}{2} \sum\_{a} Q\_{aa}^2 + \\ &+ \sum\_{a,i} \beta A Q\_{aa} t\_i \left(q\_i^{aa} + \gamma b\_i^{aa}\right) + \frac{\beta}{2} \sum\_{i,j} t\_i \chi\_j \left(\vec{H} \cdot \hbar\_{j,i}\right)^2 + \beta \mu \sum\_{i=1}^{V} t\_i\right\}, \end{split} \tag{29}$$

and performing the sums over the nematogen index *i*, the occupation variables and the states of each nematogen, we can write

$$\Xi = \left(\frac{\beta AV}{2\pi}\right)^{3/2} \int \left[d\mathbf{Q}\_{aa}\right] \exp\left\{-\beta V \Psi'(\mathbf{Q})\right\}\,,\tag{30}$$

with a grand functional

$$\Psi^{\Psi}(\mathbf{Q}) = \frac{A}{2} \sum\_{a} \mathbf{Q}\_{aa}^{2} - T \ln \Phi \left(\mathbf{Q}\right) \,\,\,\,\tag{31}$$

in which

$$\begin{split} \Psi'(\mathbf{Q}) &= 6 + 2\epsilon^{\frac{3}{2}\beta A Q\_{xx} + \frac{1}{4}\beta(\chi\_{2} + \chi\_{3})H^{2} + \beta\mu} \cosh\left[\frac{3\beta A\gamma}{2}\left(Q\_{yy} - Q\_{zz}\right) - \frac{\beta\left(\chi\_{2} - \chi\_{3}\right)H^{2}}{4}\right] \\ &+ 2\epsilon^{\frac{3}{2}\beta A Q\_{yy} + \frac{1}{4}\beta(\chi\_{2} + \chi\_{3})H^{2} + \beta\mu} \cosh\left[\frac{3\beta A\gamma}{2}\left(Q\_{zz} - Q\_{xx}\right) + \frac{\beta\left(\chi\_{2} - \chi\_{3}\right)H^{2}}{4}\right] \\ &+ 2\epsilon^{\frac{3}{2}\beta A Q\_{zz} + \frac{1}{2}\beta\chi\_{1}H^{2} + \beta\mu} \cosh\left[\frac{3\beta A\gamma}{2}\left(Q\_{xx} - Q\_{yy}\right)\right]. \end{split} \tag{32}$$

Adopting again the parametrization in Equation (10) we can write the grand functional in terms of *S* and *η*, obtaining

$$\Psi'(S,\eta) = \frac{3A}{4}S^2 + \frac{A}{4}\eta^2 - T\ln\Phi\left(S,\eta\right),\tag{33}$$

with

$$\Phi\left(S,\eta\right) = 2\varepsilon^{-\frac{3}{4}\beta A(S+\eta) + \frac{1}{4}\beta(\chi\_{2}+\chi\_{3})H^{2} + \beta\mu} \cosh\left[\frac{\beta A\Delta}{4}\left(\eta - 3S\right) - \frac{\beta\left(\chi\_{2} - \chi\_{3}\right)H^{2}}{4}\right]$$

$$+ 2\varepsilon^{-\frac{3}{4}\beta A(S-\eta) + \frac{1}{4}\beta(\chi\_{2}+\chi\_{3})H^{2} + \beta\mu} \cosh\left[\frac{\beta A\Delta}{4}\left(\eta + 3S\right) + \frac{\beta\left(\chi\_{2} - \chi\_{3}\right)H^{2}}{4}\right]$$

$$+ 2\varepsilon^{\frac{3}{4}\beta AS + \frac{1}{2}\beta\chi\_{1}H^{2} + \beta\mu} \cosh\left(\frac{\beta A\Delta}{2}\eta\right) + 6. \tag{34}$$

In the limit *βμ* 1, as expected for a fully occupied lattice, Equation (33) reduces to

$$\Psi\left(\mathcal{S},\boldsymbol{\eta}\right) = \psi\left(\mathcal{S},\boldsymbol{\eta}\right) - \boldsymbol{\mu},\tag{35}$$

in which *ψ* (*S*, *η*) is given by Equation (12).

The equilibrium values of *S* and *η* are obtained by minimizing Ψ (*S*, *η*) with respect to both variables at a fixed chemical potential, yielding the state equations

$$\begin{aligned} \frac{\partial \Psi}{\partial S} &= \frac{3A}{2} S - T \frac{\partial \ln \Psi}{\partial S} = 0\\ \frac{\partial \Psi}{\partial \eta} &= \frac{A}{2} \eta - T \frac{\partial \ln \Psi}{\partial \eta} = 0 \end{aligned} \tag{36}$$

and once again there may exist isotropic, uniaxial and biaxial solutions.

In the presence of a nonzero field, it is straightforward to show that for small dilution we recover the same qualitative results of Section 3 regarding the destabilization of the Landau point, which is replaced by a simple critical point and a tricritical point.

On the other hand, in the limit of zero field (*H* = 0), to which we limit ourselves for the rest of this section, a Landau point occurs when

$$\left.\frac{\partial^2 \Psi}{\partial \mathcal{S}^2}\right|\_{\left(0,0\right)} = \left.\frac{\partial^3 \Psi}{\partial \mathcal{S}^3}\right|\_{\left(0,0\right)} = 0.\tag{37}$$

The above condition on the third derivative yields Δ*<sup>L</sup>* = 1, while the remaining condition leads to

$$\left(\beta A - 1\right)e^{\beta \mu} - 1 = 0,\tag{38}$$

the solution of which defines a line of Landau points corresponding to a varying chemical potential. For *βμ* 1, we recover the result *βLA* = 1 obtained at zero field in Equation (19). As *μ* is reduced, the temperature corresponding to the solution of Equation (38) is also reduced, as it can be checked by calculating the implicit derivative of *β* with respect to *μ* from Equation (38).

Besides having Δ = 1 and fulfilling Equation (38), a true Landau point has to represent a stable isotropic phase. In the present case, this requires that this phase be a minimum rather than a saddle point of the grand functional. It turns out that, at the candidate Landau point, the second derivative of Ψ (*S*, *η*) with respect to *η* is always zero at the isotropic phase (*S*, *η*) = (0, 0), but that the corresponding fourth derivative,

$$\left. \frac{\partial^4 \Psi}{\partial \eta^4} \right|\_{(0,0)} = \frac{27 \beta^3 A^4 \left(1 - e^{-\beta \mu} \right)}{8 \left(1 + e^{-\beta \mu} \right)^2},\tag{39}$$

changes from positive for *μ* > 0 to negative for *μ* < 0, indicating that the Landau point corresponds to a stable phase only if the chemical potential is nonnegative.

In order to obtain phase diagrams involving the concentration *c* of occupied sites rather than the chemical potential *μ*, at every equilibrium point we need to solve for *μ* the equation

*Crystals* **2020**, *10*, 632

$$\mathcal{L} = \left\langle \frac{1}{V} \sum\_{i=1}^{V} t\_i \right\rangle = \frac{1}{\beta V} \frac{\partial}{\partial \mu} \ln \Xi\_{\prime} \tag{40}$$

which is equivalent to solving

$$\mathcal{L} = -\left. \frac{\partial \Psi}{\partial \mu} \right|\_{\left(S, \eta\right)} = \frac{\Phi\left(S, \eta\right) - 6}{\Phi\left(S, \eta\right)},\tag{41}$$

with *S* and *η* assuming their equilibrium values and Φ (*S*, *η*) given by Equation (34).

In particular, at the Landau point the relation between *μ* and the concentration *c* is given by

$$\mathcal{L} = \frac{\mathcal{e}^{\mathcal{G}\mu}}{1 + \mathcal{e}^{\mathcal{G}\mu}} \, \, \, \, \, \tag{42}$$

indicating that the Landau point is stable if *c* > <sup>1</sup> <sup>2</sup> and unstable if *<sup>c</sup>* <sup>&</sup>lt; <sup>1</sup> <sup>2</sup> . Therefore, for small dilution the Landau point is always stable. In this limit, the phase diagram is qualitatively the same as in the fully packed limit, as illustrated in Figure 4.

**Figure 4.** Zero-field phase diagram for the diluted model, when the concentration of nematogens corresponds to *c* = 0.75. Notice the lower temperature at the Landau point and a slightly more asymmetrical uniaxial–biaxial transition line around the line Δ = 1.

## **5. Conclusions**

We use a recently proposed two-tensor formalism, which is based on the physically acceptable pair interactions between nematogenic units, and assume a discrete set of orientations of the microscopic nematic directors, to write a six-state lattice Hamiltonian. The fully-connected version of this system is amenable to detailed statistical mechanics calculations. In the widely used geometric-mean approximation for the model parameters, we regain the well-known phase diagram of nematic biaxial systems, with sequences of biaxial-uniaxial-disordered phase transitions as temperature increases, and the appearance of a well-defined Landau multicritical point. We locate the transition lines and the position of the Landau multicritical point in terms of temperature and a parameter Δ that gauges the degree of biaxiality.

With the addition of an external magnetic field, we show that, except in a quite unusual diamagnetic isotropic case, the topology of this phase diagram is entirely changed. Our calculations indicate that the Landau point can no longer exist. In the phase diagrams, in terms of temperature and the parameter Δ, we show that there appear a pair of multicritical points, and a simple critical point at the end of the first-order border between a uniaxial and a paranematic phase. Also, we show the onset of a tricritical point separating borders of first-order and continuous transitions between a biaxial and a paranematic phase. We expect that lyotropic systems, for which presumably there is a significant dependence of the biaxiality parameter on both the temperature and the composition of the mixture [26,27,29], are good candidates for the observation of both the simple critical point and

the tricritical point. We point out that our predictions for *χ* > 0 are fully compatible with available experimental results for potassium laurate/1-decanol/water mixtures [8], which show, based on measurements performed at three distinct concentrations in the neighborhood of the Landau point, that the application of a magnetic field turns the transition between the rod-like uniaxial and the isotropic phases into a crossover, while preserving the discontinuous nature of the transition between the disk-like uniaxial and the isotropic phases. Locating the simple critical point and the tricritical point in this system would require some extremely fine tuning, but in principle it should be feasible.

On the other hand, in the presence of moderate dilution the zero-field phase diagram is essentially unchanged. For stronger dilution, however, the Landau point also disappears, which is very similar to a phenomenon associated with a dilute ferromagnet.

Finally, we stress that, from a formal point of view, our results are equally applicable to the situation in which a liquid crystal is subject to an electric rather than a magnetic field [33].

**Author Contributions:** Conceptualization, D.D.R., A.P.V. and S.R.S.; methodology, D.D.R., A.P.V. and S.R.S.; software, D.D.R., A.P.V. and S.R.S.; validation, D.D.R., A.P.V. and S.R.S.; formal analysis, D.D.R., A.P.V. and S.R.S.; investigation, D.D.R., A.P.V. and S.R.S.; resources, A.P.V. and S.R.S.; data curation, D.D.R., A.P.V. and S.R.S.; writing—original draft preparation, A.P.V. and S.R.S.; writing—review and editing, A.P.V. and S.R.S.; visualization, D.D.R., A.P.V. and S.R.S.; supervision, A.P.V. and S.R.S.; project administration, A.P.V. and S.R.S.; funding acquisition, A.P.V. and S.R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by INCT/FCx, NAP/FCx, FAPESP (grant 2017/09566-8), and CNPq.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Ordering of Rods near Surfaces: Concentration Effects**

## **Dora Izzo**

Instituto de Física, Universidade Federal do Rio de Janeiro, C.P. 68528, Rio de Janeiro 21941-972, Brazil; izzo@if.ufrj.br

Received: 24 February 2019; Accepted: 16 May 2019; Published: 21 May 2019

**Abstract:** We study the orientation of rods in the neighborhood of a surface. A semi-infinite region in two different situations is considered: (i) the rods are located close to a flat wall and (ii) the rods occupy the space that surrounds a sphere. In a recent paper we investigated a similar problem: the interior of a sphere, with a fixed concentration of rods. Here, we allow for varying concentration, the rods are driven from a reservoir to the neighborhood of the surface by means of a tunable chemical potential. In the planar case, the particle dimensions are irrelevant. In the curved case, we consider cylinders with dimensions comparable to the radius of curvature of the sphere; as they come close to the surface, they have to accommodate to fill the available space, leading to a rich orientational profile. These systems are studied by a mapping onto a three-state Potts model with annealed disorder on a semi-infinite lattice; two order parameters describe the system: the occupancy and the orientation. The Hamiltonian is solved using a mean-field approach producing recurrence relations that are iterated numerically and we obtain various interesting results: the system undergoes a first order transition just as in the bulk case; the profiles do not have a smooth decay but may present a step and we search for the factors that determine their shape. The prediction of such steps may be relevant in the field of self-assembly of colloids and nanotechnology.

**Keywords:** rods; curved surface; Potts

## **1. Introduction**

In the past fifty years, a great deal of attention has been devoted to understanding how liquid crystals are ordered in small cells. Applications in technology have pushed forward the interest in such systems: optical devices with specific properties depend on the anchoring conditions of the mesogens to the surfaces [1]; coating spherical colloidal particles with thin layers of liquid crystals in specific ways provides flexibility for tuning directional interactions between colloids [2]. In that scenario, particles are of the order of a few nanometeres, much smaller than their confining volumes, with sizes in the micrometer range: a description in terms of continuum theories [3] is completely satisfactory.

More recently, larger particles, commonly named macroscopic liquid crystals, have also captured attention. Examples of such structures are liquid crystals formed by viruses [4,5], filamentous biopolymers [6] and carbon nanotubes [7,8], with dimensions of the order of hundreds of nanometers. Not only is confinement a relevant issue, but also self-assembly of colloids mediated by large particles is a new problem to be investigated. From a theoretical point of view, continuum theories break down and one looks to understand how these large particles accommodate finding their optimal conformation when mutual packing and alignment to boundaries compete at the same length scale. Some theoretical descriptions have addressed such problems using Monte Carlo simulations and molecular dynamics [9–12].

Our system consists of a colloidal liquid crystal represented by a set of rod-like particles. They are driven from a reservoir by means of a chemical potential and we are primarily interested in describing how these rods self-organize around a curved surface. The concentration is not fixed and vacancies come into play, leading to *annealed* disorder. Excluded volume interactions are considered and assist organization because the particle dimensions are sizeable with respect to the surface curvature radius, typically on the order of 1/10 (as in the case of filamentous bacteriophage fd-viruses in microchambers); Figure 1 illustrates this situation.

**Figure 1.** Representation of the distribution of rods around a sphere; here, we show a cross section of the three-dimensional system. For large particle concentrations, the rods are forced to align around the sphere; this alignment decreases towards the bulk.

For a given value of a chemical potential (large enough to guarantee minimal occupancy just next to the sphere), the orientation profile can be understood as follows: (i) just outside the sphere there is enough room to accommodate a given number of particles, anchored homeotropically; (ii) beyond this region, in the upper rows, more room is available, which amounts to a decrease in the excluded volume interactions and therefore a reduction in ordering; (iii) for rows further away from the sphere, the "squeezing" effect vanishes and an isotropic configuration sets in.

We model this system in terms of a lattice surrounding a sphere. Particles (rods) occupy the lattice sites, but vacancies may occur. We follow Oliveira and Figueiredo Neto [13] and associate the orientation of each rod to a microscopic state of the three-state Potts model and solve it on a lattice. Nevertheless, the lattice description of the curved three-dimensional system is unfeasible (discussed in the next session) and we are led to solve the two-dimensional version of the problem.

The system described above is the main subject of this paper. Nevertheless, we will also address a related problem: the effect of *quenched* disorder. In this case, it is more realistic to consider the system studied in a previous work [14], in which large rods are located inside a sphere and all lattice sites are occupied. We predict the behavior of this system in terms of a qualitative discussion.

This paper is organized as follows: in Section 2, we define the two-dimensional version of the problem and describe the mean-field approach for the three-state Potts model with annealed disorder; in Section 3, the results are presented, first for a simpler problem, the planar surface (PS) case, and then for the curved surface (CS) case; these results are discussed in Section 4, where we comment on the effect of introducing quenched disorder in a similar system; conclusions are left for Section 5.

## **2. The Problem**

## *2.1. The Two-Dimensional Version*

We model the colloidal liquid crystal as a system of rods in a lattice. A single rod is located on a lattice site but not all sites must be occupied. Because it is not possible to embed a curved object into a cubic lattice, we reduced our problem to two-dimensions redefining the lattice structure and the rods' geometry as follows.

The three-dimensional version of the planar surface case (PS) consists of a cubic lattice bounded by a flat surface. The analogous two-dimensional problem is just a square lattice bounded by a line.

In order to study the two-dimensional analogue of the curved surface case (CS), the sphere is replaced by a disk. We choose a generic cross section of the three-dimensional space that surrounds a sphere with a planar lattice network that represents the symmetry of the problem as illustrated in Figure 2.

**Figure 2.** The network of lattice sites around a disk.

Accordingly, the rods must be replaced by rectangles with linear dimensions *b* and *L* (*b* < *L*). Every geometric center is fixed on a lattice site and rotates on the lattice plane. Because only excluded volume interactions are considered, the alignment of a given rod may be constrained only by its nearest neighbors' orientation. Next, we discuss the geometric constraints and show that, in the CS case, the geometry is responsible for an extra channel of disorder towards the bulk.

Consider a planar lattice such as the one in Figure 2; let *R* be the radius of the disk and *N* the number of lattice sites by row. Accommodating homeotropic anchored rods around the circle requires: *R* > *Nρ*(0)*b*. Towards the bulk, the increase in distance between nearest neighbor rods amounts to a decrease in excluded volume interactions (mutual alignment) and rotation is gradually allowed; at *r*¯ = *Nρ*(*r*¯)*L*, excluded volume interactions finally die out, which characterizes the "bulk". Therefore, the region of interest is *R* < *r* < *r*¯.

## *2.2. The Three-State Potts Model with Annealed Disorder*

We consider a plane (either PS or CS) and use the mapping proposed by Oliveira and Figueiredo Neto [13]. Because we aim to describe occupancy and ordering, we need two order parameters.

The concentration profile is determined by a set of order parameters {*ρ*} associated with the fraction of occupied sites in the system. *ρ* is an average over a microscopic parameter *η* on the th row: *ρ* ≡ *η*0.

The directional order in each row is described by another set of order parameters defined as follows. The orientation of each rod refers to the surface director; in the PS case, the surface is a straight line; in the CS case, it is a circle. We use the three-state Potts model to describe the system: one of the three states corresponds to alignment perpendicular to the surface (homeotropic), whereas the other two are degenerate and refer to alignment in the other orthogonal directions parallel to the surface (planar). The state of ordering is thus described by the set {*q*}, obtained by an average over the microscopic parameters *η* and *s* on the th row: *q* = *η*.*s*0, where the variable *s* describes the microscopic directional state on row . In fact, *ql* specifies the corresponding averaged moment of quadrupole with respect to a homeotropic alignment in the rod frame of reference **Λ** [15]:

$$\mathbf{A}\_{\ell} = \begin{bmatrix} -q\_{\ell}/2 & 0 & 0 \\ 0 & -q\_{\ell}/2 & 0 \\ 0 & 0 & q\_{\ell} \end{bmatrix} \,. \tag{1}$$

In a mean-field description of the annealed three-state Potts model, the expression for the energy per column of a given plane, either in the PS or in the CS cases, is given by [16]:

$$\mathcal{F} = -\mu \sum\_{\ell=0}^{\infty} \rho\_{\ell} - \sum\_{\ell=0}^{\infty} \epsilon\_{\ell} q\_{\ell} q\_{\ell+1} - \sum\_{\ell=0}^{\infty} \epsilon\_{\ell} (q\_{\ell})^2 + \frac{k\_B T}{3} \sum\_{\ell=0}^{\infty} \left[ \ln(\rho\_{\ell} + 2q\_{\ell}) + 2 \ln(\rho\_{\ell} - q\_{\ell}) - 3 \ln 3 (1 - \rho\_{\ell}) \right]. \tag{2}$$

The index labels the rows: = 1 refers to the first row, next to the surface; *ρ*<sup>0</sup> and *q*<sup>0</sup> fix the anchoring conditions as specified below. Each site has two nearest neighbors on the same row and one nearest neighbor on the row just above it. *μ* is a common chemical potential that drives particles from the reservoir and is the interaction energy between a rod on row and its nearest neighbors.

Following the general procedure of the mean-field approach, we minimize expression (2) with respect to *ρ* and *q*, which yeld the recurrence relations:

$$\begin{cases} -\mu + \frac{k\_{\ell}T}{3} [\ln(\rho\_{\ell} + 2q\_{\ell}) + 2\ln(\rho\_{\ell} - q\_{\ell}) - 3\ln 3(1 - \rho\_{\ell})] = 0, \\\\ -\varepsilon\_{\ell}q\_{\ell+1} - \varepsilon\_{\ell-1}q\_{\ell-1} - 2\varepsilon\_{\ell}q\_{\ell} + \frac{2}{3}k\_{B}T\ln\frac{\rho\_{\ell} + 2q\_{\ell}}{\rho\_{\ell} - q\_{\ell}} = 0. \end{cases} \tag{3}$$

This set is solved numerically under constraints in the bulk and on the surface as follows. In the bulk, *<sup>ρ</sup>*<sup>∞</sup> and *<sup>q</sup>*<sup>∞</sup> are obtained from Equation (3), by fixing *<sup>ρ</sup>* = *<sup>ρ</sup>*+<sup>1</sup> = *<sup>ρ</sup>*−<sup>1</sup> = *<sup>ρ</sup>*<sup>∞</sup> and *<sup>q</sup>* = *<sup>q</sup>*+<sup>1</sup> = *<sup>q</sup>*−<sup>1</sup> = *<sup>q</sup>*∞, leading to:

$$\begin{cases} \rho\_{\infty} = \frac{(1 - \rho\_{\infty})^3 \exp 3\mu / k\_B T}{(\rho\_{\infty} - q\_{\infty})^2} - 2q\_{\infty}, \\\\ q\_{\infty} = \frac{k\_B T}{6\kappa\_{\infty}} \ln \frac{\rho\_{\infty} + 2q\_{\infty}}{\rho\_{\infty} - q\_{\infty}}. \end{cases} \tag{4}$$

On the surface, we must specify the anchoring conditions: *ρ*<sup>0</sup> and *q*0. Equation (3) set the limits for stability:

$$\begin{cases} 0 < \rho\_\ell < 1 & \text{and} \\ -1/2 < q\_\ell < \rho\_\ell < 1. \end{cases} \tag{5}$$

The first of the above equations is trivial; the second shows that, for a given occupancy, the upper and lower bounds of *q*<sup>0</sup> correspond to homeotropic and planar anchoring, respectively, while *q*<sup>0</sup> = 0 is associated with a disordered state on the surface.

## **3. Numerical Procedure and Results**

We iterate the equations presented above; three hundred layers are enough to obtain a full picture of the profile. The region far from the surface, but still under the effect of the chemical potential, shal be termed "bulk"; the particle reservoir, free from the chemical potential, is not shown in the profiles below. In what follows, we use reduced non-dimensional parameters: the temperature *θ* ≡ *kBT*/ and the chemical potential *a* ≡ *μ*/.

Here, similarly to what was reported previously [14], the specific values of the homeotropic anchoring conditions *ρ*<sup>0</sup> and *q*<sup>0</sup> affect only the immediate vicinity of the surfaces (not shown here); the profiles are determined primarily by the bulk configuration. Moreover we notice that, for *μ* = 0, the high temperature values of *ρ* and *q* are 0.75 and 0.0, respectively.

We start by presenting the results for the PS case and then turn to the CS case.

## *3.1. The Planar Surface (PS) Case*

The existence of a first order transition is well known in analogous three-dimensional systems. Figure 3 shows the phase diagram describing the orientational ordering transition in the case of non-confinment: regions *I* and *I I* correspond to the ordered (nematic) and isotropic phases, respectively.

In the presence of a surface, the first order phase transition persists; we find that ordering transition. for *μ* = 0, the discontinuous transitions set in for 1.61 < *θ* < 1.72, and occur for *ρ* and *q* simultaneously. This is not the case for *μ* = 0, where *ρ* and *q* have transitions in different windows of temperature. Figure 4 shows their coexistence profiles for *μ* = 0 in the beginning (*θ* = 1.62) and in the end (*θ* = 1.69) of the first order transition window; they show a persistence step in the disordered solution reminiscent from the homeotropic boundary condition on the surface; the steps are wide for *θ* = 1.62, which is closer to the ordered end of the transition, while, for *θ* = 1.69, the steps almost disappear. We emphasize that *θ* is associated with the bulk degree of order, which means that the step width is determined by the bulk configuration.

**Figure 3.** Phase diagram for the bulk ordering transition, described by the non-dimensional chemical potential and temperature. The nematic and isotropic phases correspond to regions *I* and *I I*, respectively. The dashed line represents the discontinuous phase transition.

**Figure 4.** PS case profiles for *a* = 0: (**a**) occupancy and (**b**) orientational order parameter. Ordered solutions for *θ* = 1.62 (continuous lines) and disordered solutions (circles and triangles). Circles: *θ* = 1.69, triangles: *θ* = 1.62, corresponding to the high and low temperature ends of the first order window.

*Crystals* **2019**, *9*, 265

Next, we study the phase coexistence for fixed temperature and varying chemical potential. Similarly to what was reported above, the mean-field solution predicts a window of coexistence, rather than a single point. Figure 5 illustrates this aspect of the transitions: we fix *θ* = 1.69 and vary *a*. In Figure 5a, we notice that *ρ* changes discontinuously; this jump occurs somewhere between *a* 0.20 and *a* 0.82. In Figure 5b, a similar jump is observed in *q*, within exactly the same *a* range.

**Figure 5.** PS case. First order transition as a function of *a* ≡ *μ*/ for fixed *θ* = 1.69. (**a**) occupancy and (**b**) orientational order parameter. Circles: disordered solution, squares: ordered solution. Lines are just guides to the eye.

An important issue to be investigated is the influence of *μ* on the steps shape. Figure 6 illustrates this behavior: the steps go further into the bulk with increasing *μ*. We also verified that, rather than increasing slowly, at *a* = 0.07, the steps disrupt into full plateaus (both for *ρ* and *q*), a manifestation of the first order transition. Concerning the steps' height, the *ρ* profile is more sensitive to the increase in *μ* than the *q* profile.

These results show that the width of the nematic film increases with *a*, lending strong support for complete orientational wetting transition in the limit *a* → *ac*, where *ac* = 0.07. Therefore, we studied the increase of the adsorption Γ for a system of *N* layers, defined as:

$$\Gamma = \frac{1}{N} \sum\_{\ell=1}^{N} \left[ \rho\_{\ell} - \rho\_{\infty} \right] \tag{6}$$

as a function of (*a* − *ac*)/*ac*, shown in Figure 7, confirming the complete orientational wetting scenario.

**Figure 6.** PS case. Effect of increasing chemical potential on profiles step for *θ* = 1.62. (**a**) occupancy and (**b**) orientational order parameter. Dotted line: *a* = 0.00, dashed line: *a* = 0.02 and continuous line *a* = 0.06.

**Figure 7.** PS case. Adsorption Γ as a function of the difference between the reduced chemical potentials below the transition and at the transition (*θ* = 1.62).

We also checked the effect of imposing planar anchoring. Figure 8 shows profiles that illustrate the situation: in this case, the boundary condition on the wall is strongly planar (*q*(0) = −0.49, *ρ*(0) = 0.9), *a* = 0 and we considered two temperatures in the vicinity of the first order transition, *θ* = 1.60 and *θ* = 1.70. For the lower temperature (below the transition), the first layers show planar ordering and eventually reach an condition; for the higher temperature (above the transition), the first layers show also planar ordering that decay into an isotropic state towards the bulk. It seems that no long range

planar ordering can be reached throughout the sample: this is because the mean-field approach for Potts models always favours one direction.

**Figure 8.** PS case. Profile for the orientational order in the case of strongly planar anchoring, *a* = 0.00. Continuous line: *θ* = 1.60, just below the transition; dashed line: *θ* = 1.60, just above the transition.

Next, we turn to the more interesting situation: the curved surface case.

## *3.2. The Curved Surface (CS) Case*

Here, we show the results for decreasing interactions towards the bulk: we fix the interactions values on the surface (*θ<sup>s</sup>* ≡ 1/*kBT*) and in the bulk (*θ<sup>b</sup>* ≡ ∞/*kBT*); we may also impose the interaction decrease rate.

Here, similarly to the last case, first order transitions are obtained within a window of parameters: both for the bulk temperature at fixed chemical potential and for the chemical potential at fixed bulk temperature (not shown here).

In the vicinity of the surface, for sufficiently strong interparticle interactions (low *θs*), an ordered configuration sets in. In order to describe how the system evolves to a less ordered, or even disordered state in the bulk, we vary the interaction decrease rate (*θ* increase rate). Figure 9 shows profiles for ∼, −<sup>1</sup> and ∼ −0.7, for *<sup>θ</sup><sup>s</sup>* = 1.69, *<sup>θ</sup><sup>b</sup>* = 1.80 and *<sup>a</sup>* = 1.0. The occupancy (Figure 9a) decreases to a finite value in the bulk, whereas the orientational order parameter (Figure 9b) decreases to zero (isotropic state). Both order parameter profiles reach the corresponding bulk values presenting a step associated with an ordered region close to the surface; the width of the steps increases as the interaction decay rate becomes smoother.

The dependence of with describes the alignment strength between a rod in row and its neighbors. Let *R* be the circle (surface) radius, *a* the distance between adjacent rows, *d* and *d*<sup>0</sup> the distance between neighboring sites on row and on the surface, respectively. Then, (*R* + .*a*)/*N d*. If <sup>∼</sup> *γ*, where *<sup>γ</sup>* is a generic exponent, then (*d*) <sup>∼</sup> (*d*/*d*<sup>0</sup> <sup>−</sup> <sup>1</sup>)−*γ*. This means that the alignment strength depends on a non-dimensional distance between rods ¯*<sup>d</sup>* like ( ¯*<sup>d</sup>* <sup>−</sup> <sup>1</sup>)−*γ*.

Another aspect regarding steps is that they occur just in situations where the bulk degree of disorder is not too high. As the bulk temperature increases, the step width decreases, up to a situation where the bulk disorder washes out the step, that is, no order is observed throughout the system. Figure 10 illustrates the dependence of the step width on the bulk temperature; the region next to the surface is ordered (*θ<sup>s</sup>* = 1.69), *<sup>a</sup>* = 1.0 and the interaction decays as ∼ −1.

**Figure 9.** CS case. Profiles for: (**a**) occupancy and (**b**) orientational order parameter. Circles: <sup>∼</sup> −1; squares: <sup>∼</sup> −0.7. *<sup>θ</sup><sup>s</sup>* <sup>=</sup> 1.69, *<sup>θ</sup><sup>b</sup>* <sup>=</sup> 1.80 and *<sup>a</sup>* <sup>=</sup> 1.0. Lines are just guides to the eye.

**Figure 10.** CS case. Profiles for: (**a**) occupancy and (**b**) orientational order parameter. <sup>∼</sup> −1; *θ<sup>s</sup>* = 1.69; *a* = 1.0; squares: *θ<sup>b</sup>* = 1.80, circles: *θ<sup>b</sup>* = 1.82, and triangles: *θ<sup>b</sup>* = 1.90.

We also verified the opposite case: a disordered bulk and an ordered surface, such that the bulk temperature is fixed at a high value (*θ<sup>b</sup>* = 1.82) and we keep decreasing the surface temperature, starting from *θ<sup>s</sup>* = 1.65. Appreciable changes in the step height and width are reached only at much lower values of *θs*, of the order of 1.25. With these results (not shown here), we conclude that the steps shape is much more sensitive to the bulk than to the surface temperature.

Concerning the dependence of the step shape on *μ*, we observed that the CS case responds similarly to what was reported for the PS case (Figure 6).

#### **4. Discussion**

## *4.1. Quenched Disorder*

In order to compare the effects of disorder in related systems, we start this section by recalling a problem considered in our previous work [14]. There we had rods confined *inside* a disk-shaped cavity and all lattice sites were occupied, that is, we had a fixed concentration of particles. In order to allow for varying concentration, we introduce here permanent vacant sites: that is, a system with site *quenched* disorder, as opposed to *annealed* disorder studied so far. To our knowledge, there is no complete solution for the Potts model with site disorder in the literature [17], therefore we analyse this problem with an approximate argument.

We recall the expression for the free energy per column:

$$\nabla = \sum\_{\ell=0}^{\infty} \mu\_{\ell} q\_{\ell} - \sum\_{\ell=0}^{\infty} \varepsilon\_{\ell} q\_{\ell} q\_{\ell+1} - \sum\_{\ell=0}^{\infty} \varepsilon\_{\ell} (q\_{\ell})^2 + \frac{k\_B T}{3} \sum\_{\ell=0}^{\infty} \left[ 2(1 - q\_{\ell}) \ln(1 - q\_{\ell}) + (1 + 2q\_{\ell}) \ln(1 + 2q\_{\ell}) - 3 \ln 3 \right], \tag{7}$$

where *μ* is an external potential which vanishes everywhere, except for = 1 (next to the surface).

The corresponding recurrence relations for the order parameter are:

$$q\_{\ell} = \frac{\mu\_{\ell}}{2\varepsilon\_{\ell}} - \frac{1}{2} \left( q\_{\ell+1} + q\_{\ell-1} \right) + \frac{k\_B T}{3\varepsilon\_{\ell}} \ln \frac{1 + 2q\_{\ell}}{1 - q\_{\ell}}.\tag{8}$$

One way to introduce quenched dilution is to claim that the number or nearest neighbor sites is reduced by a factor *h*, where 0 < *h* < 1. In that order, we replace all 's by *h*.'s in Equations (7) and (8), yielding the modified recurrence relations:

$$q\_{\ell} = \frac{\mu\_{\ell}}{2h\varepsilon\_{\ell}} - \frac{1}{2} \left( q\_{\ell+1} + q\_{\ell-1} \right) + \frac{k\_B T}{3h\varepsilon\_{\ell}} \ln \frac{1 + 2q\_{\ell}}{1 - q\_{\ell}}.\tag{9}$$

These last relations show that, in this case of quenched disorder, the introduction of vacancies by means of *h* brings no qualitative novelty: first order transitions are expected as a function of the concentration (or *h*); moreover, the results obtained formerly [14] will be modified by rescaled *T*'s and *μ*'s.

#### *4.2. Annealed Disorder*

The present study, which takes into account annealed disorder, produces a much richer dependence on the concentration, to be discussed next.

A first order transition is observed both in the PS and CS cases, as expected. The profiles depend on the chemical potential (or concentration) and on the strength of the interparticle interactions; nevertheless, the shape of the steps is determined primarily on the degree of ordering (temperature) in the bulk. All results are robust with respect to the strength of the homeotropic anchoring, they do not affect the profiles, except for the immediate neighborhood of the surfaces. Such behavior is not observed in simulations of similar systems [10–12], where profiles are strongly dependent on the anchoring condition. The reason for this discrepancy is that the mean-field approach intrinsically favors a preferential direction.

In the PS case, for vannishing *μ*, we obtained profiles associated with different temperature values, all of them in the first order transition temperature range. The disordered solution in the low temperature end of the coexiste window (unstable solution) has a step, reflecting the presence of the nearby surface; the disordered solution in the high temperature end of coexistence (stable solution) presents just a tiny step. Temperature is associated with the degree of order in the bulk; the step, which arises because of the anchoring at the nearby surface, can be destroyed for sufficiently high bulk disorder. We also characterize the isotropic nematic first order transition in terms as a complete orientational wetting scenario. The results for planar anchoring are presented; nevertheless, they are not very reliable due to the details of the mean-field approach.

In the CS case, the effect of decreasing interactions towards the bulk produces very interesting effects. For an ordered surface and a disordered bulk, the profile may present a step. A first observation regarding the steps is that their widths depend on the slope of the interaction profile, being wider for smoother profiles. We also showed that, similarly to what happens in the PS case, steps are destroyed when interactions in the bulk are weak enough. On the other hand, for a fixed degree of disorder in the bulk, enhancing the interactions' strength on the surfac, produces just a feeble effect on the magnitude of the order parameters and on the step width. These two observations lead us to conclude that the step height depends weakly on the surface temperature and their width is mainly determined by the bulk degree of disorder.

It is possible to compare our results for the bulk densities at the nematic isotropic with those of de Las Heras end Velasco [10]; they obtain a global packing fraction of the order or 0.29, to be compared to ours *ρ* − *ρ*<sup>0</sup> (*ρ*<sup>0</sup> = 0.75), which is approximately 0.15 for *a* = 1.0.

We also studied the effect of varying the chemical potential on the step shape. Increasing *μ* amounts to enhancing the overall concentration and produces wider steps. This behavior has also been observed in previous numerical studies on spherocylinders near flat [9] and curved walls [10].

Our motivation was to describe the distribution of colloids near surfaces, a three-dimensional problem. Solving the two-dimensional problem amounts to disregarding interactions between rods in neighboring planes. In a mean-field description, it consists of considering a smaller number of first neighbors. Having obtained the orientation profile on a given plane, we claim that this simplification reduces the three-dimensional *ρ*'s and *q*'s by a common factor. This affects the corresponding step widths and heights; nevertheless, the qualitative features of the profiles are the same. This is certainly true for flat surfaces and a good approximation for curved surfaces of large radii.

## **5. Conclusions**

Motivated by recent experiments on colloidal particles with dimensions comparable to the confining volumes, we extend a model proposed in a recent work to describe the orientation of rods on planes bounded by straight or curved lines, considering the concentration as an extra parameter: it corresponds to the fraction of occupied sites. The particles are driven from a reservoir to the system by means of a chemical potential, so we allow for annealed disorder. They occupy sites on a semi-infinite lattice above a line or around a disk. We propose a lattice model to describe such systems and map it onto a three-state Potts model with vacancies, which is solved within a mean-field approach. In the curved case, the effect of the curvature is reproduced by layer dependent interparticle interactions. For both types of surfaces, we obtained the occupancy and the orientational profiles, which were studied under various conditions.

The main outcome of our study is to predict the non-smooth behavior of the profiles: the occurrence of steps both in the occupancy and in the orientational order parameters profiles. They may occur because the surfaces impose specific anchoring; nevertheless, their stability is primarily determined by the degree of order in the bulk. Steps like these have been also predicted in the literature, mainly in simulations, albeit with no systematic analysis. Furthermore, we believe that it may be of interest for applications in nanotechnology.

A shortcoming of our study is the weak influence of the anchoring on the overall profile. Moreover, we considered a fixed chemical potential throughout the sample; a more realistic procedure would be to model it as a decreasing function from the surface towards the bulk; this procedure would imply a more complicated interpretation of the results. At this time, an extension of this problem is underway: biaxial particles are being considered.

**Funding:** This research received no external funding.

**Acknowledgments:** We thank M. J. de Oliveira for very useful discussions and comments.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## **References**


c 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Crystals* Editorial Office E-mail: crystals@mdpi.com www.mdpi.com/journal/crystals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18