Mixing soft materials with appropriate nanoparticles (NPs) can yield effective materials exhibiting new or anomalously enhanced properties (Balazs et al., 1979; Hamley, 2003). If soft materials exhibit a kind of orientational or translational order, the richness of the resulting qualitatively different effective materials explodes. Liquid crystalline (LC) materials are an excellent example. They combine a unique combination of liquid character, order, softness, and optical transparency (de Gennes, 1995; Kleman and Lavrentovich, 2003; Oswald and Pieranski, 2019). The liquid character enables relatively simple preparation of mixtures. Order in LC matrices can give rise to long-range forces among immersed NPs (Poulin et al., 1979; Pires et al., 2007). Softness can enable relatively strong responses in the LC matrix in the presence of adequate NPs (Lelidis et al., 1993; Li et al., 2006; Palffy-Muhoray, 2007). Furthermore, optical transparency enables relatively simple observation of NP-driven changes (Poulin et al., 1979; Pires et al., 2007). It should be noted that long-range forces among NPs could trigger the self-assembling of NPs, which can open different structural pathways (Hegmann et al., 2007; Bisoyi and Kumar, 2011; Lagerwall and Scalia, 2012; Lagerwall and Scalia, 2017).

Nematic (de Gennes, 1995; Kleman and Lavrentovich, 2003; Oswald and Pieranski, 2019) orientational order represents the simplest LC phase structure. It is described by the nematic molecular field. In bulk nematic equilibrium, the molecular field is uniaxial and spatially homogeneously aligned along a single symmetry-breaking orientation. Local uniaxial order is commonly described in terms of the nematic director field n → and the nematic order parameter S. The unit vector points along the local uniaxial direction, where the states ± n → are physically equivalent. The amplitude field S describes the amount of ordering along n → . In thermotropic LCs, one commonly obtains a nematic phase by cooling it from the isotropic (ordinary liquid) phase in which S = 0. Perturbed nematic order could generally exhibit biaxial states, which requires description in terms of the tensor nematic order parameter Q _ (see Supplementary Appendix SA).

Appropriate NPs could influence the nematic order via different generic mechanisms if inserted into nematic LCs. Important controlling parameters are the concentration of NPs, geometrical shape (Bellini et al., 2000; Kyrou et al., 2018; Kyrou et al., 2020; Skarabot et al., 2022), NP surface treatment (Nobili and Durand, 1992; Lavric et al., 2013a; Oswald and Pieranski, 2019; Kyrou et al., 2020), material characteristics of NPs (Mertelj et al., 2013; Kumar, 2014; Moghadas et al., 2015; Poursamad and Hallaji, 2017; Emdadi et al., 2018a; Emdadi et al., 2018b; Drozd-Rzoska et al., 2019; Jahanbakhsh et al., 2019; Poursamad and Emdadi, 2019; Lahiri et al., 2020; Bury et al., 2022; Ma et al., 2022; Sumandra et al., 2022), and LC material properties (de Gennes, 1995; Kleman and Lavrentovich, 2003; Oswald and Pieranski, 2019). Furthermore, NPs could effectively impose qualitatively different disorders on the LC order (Harris et al., 1973; Cleaver et al., 1996; Crawford and Žumer, 1996; Radzihovsky and Toner, 1997; Popa-Nita and Kralj, 2006). Most studies focusing on the impact of the disorder are performed in LC-aerosil mixtures, which can exhibit at least three qualitatively different disorder characteristics (Jin and Finotello, 2001; Bellini et al., 2002; Leon et al., 2004; Rotunno et al., 2005; Buscaglia et al., 2006; Cordoyiannis et al., 2006; Relaix et al., 2011). Furthermore, NPs could enforce LC matrix topological defects (TDs) (Mermin, 1979; Kurik and Lavrentovich, 1988) in the nematic orientational order. TDs correspond to topologically protected, elastically distorted regions in the orientational order. Nematic LCs could host either point or line defects (Schopohl and Sluckin, 1987; Kurik and Lavrentovich, 1988; Lavrentovich, 1998; Kralj and Virga, 2001). Generally, TDs can strongly interact with NPs, opening the door to predetermined and controlled superstructures stabilized by TD–NP interactions (Kikuchi et al., 2002; Yoshida et al., 2009; Karatairi et al., 2010; Liu et al., 2010; Rozic et al., 2011; Coursault et al., 2012; Wang et al., 2012; Lavric et al., 2013a; Lavric et al., 2013b; Liu et al., 2018).

This paper considers different NP-driven generic mechanisms via which NPs impact nematic thermotropic behavior.

Mixtures of nematic LCs and NPs can exhibit complex configurations. In the following, we present some qualitatively different scenarios and discuss generic mechanisms.

Supplementary Appendix SB illustrates a simple mean-field Maier–Saupe analysis that reveals the type of coupling terms in a mixture of two nematic mesogens. It yields the following free energy density expression: f ≈ f 0 + 1 − p a 0 T − 1 − p T * S 2 − b S 3 + c S 4 , + p A 0 T − p T 2 * S 2 2 − B S 2 3 + C S 2 4 − w S S 2 . (1)

Here, S ≡ S 1 and S 2 label the nematic order parameters of the first and second components. Their volume concentrations are given by p 1 = 1 − p and p 2 = p , respectively. Quantities a 0 , b, c, T * ≡ T 1 * , A 0 , B , C, and T 2 * are material constants, in which we neglect temperature dependencies in the temperature windows of our interest. Furthermore, f 0 determines the contributions of the remaining degrees of freedom.

Let us suppose that the first and second components represent LC molecules and NPs, respectively, and the volume concentration of NPs is given by p ≡ p 2 . The presence of NPs introduces an additional contribution equal to p 1 − p a 0 T * S 2 in the LC condensation term in Eq. 1.

This term resembles the structure of the so-called Flory–Huggins free energy contribution (Onsager, 1949; Flory, 1956; Doi, 1981): f F H = p 1 − p χ , (2) where χ stands for the Flory–Huggins constant. For a positive value of χ , this term enforces phase separation if χ exceeds the critical value χ c . Suppose that in the absence of ordering, it holds χ I < χ c , where χ I denotes the Flory–Huggins constant in the isotropic phase (where S = 0). On entering the nematic phase, the effective Flory–Huggins interaction increases: χ N = χ I + a 0 T * S 2 . (3)

In most LCs, the second term, “switched-on” in the nematic, is relatively strong. Consequently, on cooling an isotropic mixture to the nematic phase, phase separation is generally likely to occur (Anderson et al., 2001). Figure 1 illustrates a typical phase separation in a mixture of nematic and spherical NPs (i.e., S 2 = 0 ) on entering the nematic phase. In the phase separation regime, the regions exhibiting essentially strong and low nematic LC order coexist. These regions exhibit relatively low and high concentrations of NPs, respectively.

Let us suppose that the mixture remains homogeneous. A simple Maier–Saupe-type analysis (Humphries et al., 1972) presented in Supplementary Appendix SB suggests that in a mixture of two nematogenic components, a free energy density f c o u p l term coupling ordering of both ingredients arises; see the last term in Eq. 1. It should be noted that despite its symmetric structure (e.g., in Eq. 1), the coupling term reads f c o u p l = − w S 1 S 2 ); such terms have significantly different impacts on the two components. Here, S i determines the amplitude of nematic ordering in the ith component. Furthermore, in real samples, different structures of coupling terms arise depending on the type of NPs and their surface treatment. Figure 2 illustrates the “sergeant–soldier” type of behavior, where we consider the simplest possible term (Holbl et al., 2022), f c o u p l = − w S 1 2 S 2 2 , (4) which does not affect the qualitative order-disorder behavior of individual components. For simplicity, we assume that both components have identical condensation material constants with the exceptions of bulk phase transition temperatures of isolated components and w > 0 (i.e., this term tends to increase the degree of order in both components). We assume that an isolated ith component exhibits an isotropic-nematic phase transition at T i c , where T 2 c > T 1 c . Thus, the coupling term is absent in the regime T > T 2 c , where both order parameters are melted. However, in the temperature regime T < T 2 c , the second component condensates. For this reason, the first mesogen experiences the effective free energy density ordering field contribution f c o u p l = − w e f f S 1 2 , (5) where w e f f = w S 2 2 > 0 . This term renormalizes the phase transition of the first component. It holds T e f f c = T 1 c + w e f f a 1 , (6) where T e f f c determines the enhanced phase transition of the first component. Therefore, increasing the value of the coupling strength raises the phase transition temperature of the first component, while the phase transition of the second component is unaffected. The resulting temperature behavior of the mixture is shown in Figure 2. It should be noted that in the case of bilinear coupling (see the last term in Eq. 1), the first component also exhibits pre-transitional effects in the temperature interval T ∈ T e f f c , T 2 c if the effective coupling strength is below some critical value. If the latter value is exceeded, the degree of order S 1 exhibits gradual (noncritical) increase with decreasing temperature below T 2 c .

Next, we consider a dilute mixture where the second component consists of anisotropic particles. For illustration, we assume that the free energy F consists of LC condensation, LC nematic elastic, and NP–LC interfacial free energy contributions, where details are given in Supplementary Appendix SA. Of our interest is the impact of NPs on the LC order–disorder phase transition. We assume that NPs relatively weakly disrupt the nematic director field, as illustrated in Figure 3. In the nematic phase, n → is essentially homogeneously aligned along a single symmetry-breaking direction. We assume that a nanoparticle locally enforces orientation along the unit vector e → , which is allowed to fluctuate. We also set S to be essentially spatially homogeneous, i.e., S ∼ S ¯ , where … ¯ indicates the spatial averaging over the system’s volume V. Furthermore, we neglect spatial variations in the nematic director field. With this in mind, it follows (see Supplementary Appendix SA) F V ∼ a 0 T − T * S 2 − b S 3 + c S 4 − p a N P w S v N P   S 2 , (7) where S 2 ≡ P ¯ 2 = 3 e → ∙ n → 2 − 1 2 ¯ measures the orientational order of anisotropic NPs. Quantities N N P , a N P , and v N P stand for the number of NPs, NP’s surface area, and NP’s volume, respectively. For weakly interacting NPs, it roughly holds (van der Schoot et al., 2008) S 2 ∼ ∫ P 2 e α x d x ∫ e α x d x , (8) where e → ∙ n → = x ∈ − 1,1 , α = a N P w S / k B T , k _B , and we imposed cylindrical symmetry.

In the limit α ≫ 1 , it holds S 2 ∼ 1 . Eq. 7 yields f ∼ ∼ r   s 2 − 2 s 3 + s 4 − σ s , (9) where r = T − T * T I N − T * is the reduced temperature, f ∼ = f / a 0 T I N − T * S 0 2 , s = S/S ₀ stands for the scaled nematic order parameter, the scaling unit S 0 = b 2 c measures the nematic order at T I N in bulk equilibrium, and σ ∼ p a N P ξ I N 2 v N P d e is the effective dimensionless field conjugated to the order parameter. Quantities ξ I N and d e stand for the nematic correlation length at T I N and the surface extrapolation length (see Eq. (A6)). The latter length describes the effective strength of LC–NP coupling.

The phase behavior of this effective system is as follows. Typical nematic order parameter temperature variations are depicted in Figure 4. For σ < σ c = 0.5, the system exhibits a first-order transition at r c = 1 + σ , where the bulk I–N phase transition corresponds to r I N = 1 . On the other hand, for σ ≥ σ c , the system exhibits a gradual evolution of nematic order with decreasing temperature.

If α ∼ 1 , the following roughly holds: S 2 ≡ P ¯ 2 ∼ a N P w S ¯ k b T . In this case, NPs only renormalize the I–N phase transition temperature. It should be noted that we have neglected the nematic director field distortions in the previously considered estimates.

In the following, we discuss cases where NPs introduce a certain degree of randomness. The I–N phase transition is extremely susceptible to disorder due to the existence of Goldstone modes in the nematic director field in the nematic phase (Kleman and Lavrentovich, 2003; Palffy-Muhoray, 2007). These are the consequences of continuous symmetry breaking, via which bulk nematic equilibrium is established. According to the Imry–Ma theorem (Larkin, 1970; Imry and Ma, 1975), even infinitesimally weak random-field-type disorder is sufficient to break the long-range order (LRO) of the pure bulk phase, which is reached via continuous symmetry breaking. The resulting phase should possess short-range order (SRO), and the resulting domain-type pattern is predicted to be dominated by a single characteristic domain size ξ d ∝ w R F − 2 4 − d . Here, w R F measures the disorder strength, and d is the space-dimensionality. Therefore, in effectively d = 2 or d = 3 systems, one expects ξ d ∝ w R F − 1 and ξ d ∝ w R F − 2 , respectively. However, several later experimental and theoretical systems reveal that the resulting behavior is much more complex. For appropriate conditions, configurations exhibiting quasi-long-range order (QLRO) (Denholm and Sluckin, 1993; Chakrabarti, 1998; Feldman, 2000; Cvetko et al., 2009; Ranjkesh et al., 2014) or even LRO (Chakrabarti, 1998; Bisoyi and Kumar, 2011; Ranjkesh et al., 2014) are reported.

A convenient model to study the impact of NP-imposed randomness on the I–N phase transition is the Lebwohl–Lasher lattice model (Lebwohl and Lasher, 1972). Its key ingredients are presented in Supplementary Appendix SC. Using it, we probe the impact of NP-imposed random-type behavior on the LC order–disorder phase transition in the orientational order. Despite its simplicity, this approach well captures the essential features of LCs. In modeling, we vary the concentration p of sites by imposing a randomly selected orientation with finite disorder strength W. In typical studies focusing on a system’s range of orientational order, one commonly measures or calculates the orientational order correlation function G(r). This correlation function measures how orientational correlations decay with relative distance r. However, several studies (Denholm and Sluckin, 1993; Cvetko et al., 2009; Ranjkesh et al., 2014) reveal that G(r) dependence relatively poorly distinguishes between SRO and QLRO. Numerically, it is computationally more effective to extract the range of order using finite-size analysis (Eppenga and Frenkel, 1984; Ranjkesh et al., 2014). Namely, according to the central limit theorem, one expects that S ≡ P ¯ 2 in the case of SRO exhibits scaling behavior P ¯ 2 ∝ L − γ , (10) where γ = 3 / 2 in d = 3. If QLRO replaces SRO, it follows 0 < γ < 3 / 2 . On the other hand, LRO is fingerprinted by γ = 0 .

Our simulations used this finite-scaling approach to determine the range of order on the varying concentration of NPs, their imposed anchoring strength W, and temperature. Furthermore, in the presence of disorder, one expects history-dependent behavior. For this reason, we probed systems’ behavior for three different histories, which we refer to as i) annealed history (AH), ii) temperature-quenched history (TQH), and (iii) field-quenched history (FQH). In AH, we gradually decreased temperature stepwise and calculated the structure at a given T by using as the initial “seed” structure the fixed-point configuration calculated at the previous T. In QH, we calculated the structure at a given T by originating from the isotropic phase. Finally, in FQH, the “seed” structure was spatially homogeneously aligned along a single symmetry-breaking direction. In Figures 5, 6, we plot γ for varying several control parameters. In Figure 5, we focus on W-driven behavior for different concentrations of NPs using AH deep in the nematic phase. In Figure 6, we present a more detailed impact of W on LC configurations at one concentration of NPs where we vary the history of samples and probe two different temperatures below the bulk I–N phase transition temperature. One sees that SRO (i.e., configurations characterized by γ = 3 / 2 can be reached only for relatively high values of W, especially in diluted samples. The control parameters T and W are given in the dimensionless scaled form. The references are as follows: bulk I–N phase transition is realized at T ∼ 1.1 , and W is measured with respect to LC–LC neighboring molecular interactions. It should be noted that for low enough values of W and p, one observes LRO.

We illustrated several mechanisms via which NPs could quantitatively or qualitatively modify the phase behavior of the bulk nematic LC phase. The simplicity of nematic ordering explodes into a rich palette of behaviors enabled by varying control parameters. In the study, we addressed only a few of them: concentration of NPs, nature and strength of NP–LC interactions, and anisotropy of NPs. The diversity of phenomena becomes even wider if the ferromagnetic, ferroelectric, or multiferroic properties of NPs are included. The ferromagnetic nanoparticles dispersed within the liquid crystals provide a new way to develop magneto-optic devices (Ji et al., 2019) that can also be used for THz filtering and modulation, in which sensing applications are becoming increasingly widespread (Abina et al., 2022).

In this paper, we focused only on phase behavior. It should be noted that structural behavior could be even richer. For instance, recent technological advances enable the stabilization of nematic structures incorporating diverse configurations of TDs, i.e., lattices of disclinations (Culbreath et al., 2011; Ackerman et al., 2012a; Ackerman et al., 2012b; Evans et al., 2013; Murray et al., 2014; Glazar et al., 2015; Guo et al., 2016; Peng et al., 2017; Wang et al., 2017; Yu et al., 2019; Sohn et al., 2020; Guo et al., 2021) or other assemblies (Harkai et al., 2020). Such structures could be further modified by appropriate NPs, which could make them more robust or introduce additional functionalities in the system. For example, conductive NPs assembled within a line defect can form a conductive wire (Coursault et al., 2012). Furthermore, recent experiments demonstrated (Harkai et al., 2020) that one could efficiently switch among competing line-defect structures using external fields. If these well-controllable defects would drag trapped NPs with them, one could construct rewritable networks of NPs, which could open doors to diverse applications.