Chaining of hard disks in nematic needles: particle-based simulation of colloidal interactions in liquid crystals.

Colloidal particles suspended in liquid crystals can exhibit various effective anisotropic interactions that can be tuned and utilized in self-assembly processes. We simulate a two-dimensional system of hard disks suspended in a solution of dense hard needles as a model system for colloids suspended in a nematic lyotropic liquid crystal. The novel event-chain Monte Carlo technique enables us to directly measure colloidal interactions in a microscopic simulation with explicit liquid crystal particles in the dense nematic phase. We find a directional short-range attraction for disks along the director, which triggers chaining parallel to the director and seemingly contradicts the standard liquid crystal field theory result of a quadrupolar attraction with a preferred [Formula: see text] angle. Our results can be explained by a short-range density-dependent depletion interaction, which has been neglected so far. Directionality and strength of the depletion interaction are caused by the weak planar anchoring of hard rods. The depletion attraction robustly dominates over the quadrupolar elastic attraction if disks come close. Self-assembly of many disks proceeds via intermediate chaining, which demonstrates that in lyotropic liquid crystal colloids depletion interactions play an important role in structure formation processes.

Introduction

Colloidal suspensions of nanometer to micrometer-sized particles in a host fluid can form liquid and crystalline phases, but also liquid crystalline mesophases if colloidal particles are, for example, rod-shaped[1, 2]. Phase behavior and stability of a colloidal system can often be explained based on effective interactions between colloidal particles which arise from integrating out microscopic degrees of freedom of the host fluid. The resulting effective interactions govern colloidal stability, coagulation, flocculation, and structures that eventually self-assemble.

Colloidal mixtures containing different colloidal particles, often of different size, feature additional effective interactions if degrees of freedom of one species are integrated out[3]. The most prominent example of additional effective interactions in mixtures are depletion interactions, as they arise, for example, in a mixture of large and small hard spheres[4, 5]: a short-range attractive depletion interaction between large hard spheres emerges because the excluded volume for the small spheres decreases (and, thus, their available phase space increases) if large spheres approach closer than one diameter of a small sphere.

Effective interactions become even more interesting for colloids suspended in anisotropic fluids such as liquid crystals (LCs), which can also be seen as colloidal mixtures of larger colloidal particles suspended in a liquid of small rod-like particles, in particular for lyotropic LCs[6]. Such LC colloids exhibit anisotropic effective interactions between colloidal particles if the LC is in an ordered, e.g., nematic phase[7]. The first studies of LC colloids have been performed on latex spheres suspended in a lyotropic nematic LC (a micellar nematic phase of discoid type)[8, 9]. The nematic LC phase forms typical defect-structures around a spherical inclusion depending on the anchoring conditions, “Saturn-ring” disclination rings or a satellite hedgehog for normal anchoring and boojums for planar anchoring[10-12]. These defect-structures also induce long-range colloidal interactions between spherical inclusions. These have been first explored systematically by Ramaswamy et al.[13], Ruhwandl and Terentjev[14] and by Poulin et al. who also observed a chaining of water droplets inside a LC[15]. The nature of these elastic LC-mediated interactions strongly depends on the details of the interaction between colloidal particle and LC host, i.e., how the LC molecules or rods are anchored on the colloid (normal, conic or planar, weak or strong anchoring), and can be of dipolar, quadrupolar, or even more complicated nature[16-18]. Such systems are promising candidates to realize a controlled self-assembled structure formation by anisotropic interactions[19-24]. The elastic interactions in combination with different shapes and sizes of the colloidal particles can be used to create a multitude of different colloidal interactions[25-27]. LC colloidal assemblies can be engineered, for example, by tuning the surface anchoring or engineering nematic defect structures[24,28, 29].

Here, we consider a colloidal mixture of hard disks and hard needles, which serves as a simple two-dimensional model system of colloidal spheres suspended in a lyotropic LC. We note that systems of had rods and hard disks in three dimensions are qualitatively different as both disks and rods form LC phases in three dimensions[30,31]. We are interested in the situation, where the hard needles are sufficiently dense to form a nematic LC phase such that they can mediate elastic interactions between the hard disks. This system is interesting from a theoretical point of view as it can be interpreted, on the one hand, as a colloidal mixture, where we expect depletion interactions between two hard disks embedded into a fluid of shorter needles. We also expect, on the other hand, to find long-range elastic interactions mediated by the nematic hard needle LC, at least for needles shorter than the disks, such that a coarse-grained continuum description is appropriate. This is the situation we want to address in this paper. The interplay and competition of (at least) two types of effective interactions—short-range depletion and long-range elastic—has to be unraveled and will have interesting consequences for the total effective interaction between the disks.

Depletion interactions are well documented for the dilute isotropic phase of needles[32-37] but much less is known for a nematic host with the notable exceptions of Refs.[38, 39], where the limit of small spheres in long needles has been considered. In the experimental study in Ref.[39] a chaining of small spherical particles in a host of fd virus rods parallel to the director has been found, which was attributed to depletion attraction.

On the other hand, the theory of colloidal LCs should apply with long-range elastic interactions mediated by director field distortions in the nematic hard needles. Depending on anchoring conditions and dimensionality of the system dipolar interactions (falling off as with the sphere separation in three dimensions)[40-42] or quadrupolar interactions (falling off as ) can occur[13, 14, 43–45]. Because hard needles tend to align tangentially at a hard wall, we expect a quadrupolar elastic interaction, which is characteristic for planar anchoring at the colloidal disk[13, 14, 44–46] but also generic in two dimensions[43]. Both for dipolar and quadrupolar interactions the elastic interaction is attractive and chaining of colloidal spheres has been experimentally observed[20, 47, 48]. Whereas dipolar interactions prefer chaining of spheres parallel to the director axis in three dimensions[41, 42], quadrupolar interactions prefer an angle of approximately with respect to the director axis in three dimensions[44, 46] and a angle in two dimensions[45]. For our two-dimensional system of hard disks and hard needles we thus expect the elastic interaction to favor a angle if disks form chains.

Colloidal mixtures are also computationally challenging systems. Effective interactions are essential to characterize stability and potential self-assembly into crystalline phases but hard to access in a microscopic particle-based simulation. The process of integrating out microscopic degrees of freedom corresponds to the numerical evaluation of a potential of mean force between the colloidal species of interest. In order to measure the potential of mean force all degrees of freedom including, for example, relatively slow large spheres in a bath of small particles must be properly equilibrated.

The colloidal mixture of hard disks suspended in a nematic host of hard needles is particularly challenging as the hard needle system must be fairly dense to establish a nematic phase. While particle-based simulations exist for dilute rods in the isotropic phase[49, 50], the regime of a nematic host is fairly unexplored up to now and simulations resorted to coarse-graining approaches[46]. So far, only single inclusions[51] or confining geometries[52] have been investigated by particle-based simulations. In order to calculate the effective interactions between hard disks and their self-assembly we apply a novel Monte Carlo (MC) technique, the rejection-free event-chain sampling technique. The event-chain algorithm has been originally proposed for pure hard sphere systems[53], and recently generalized to dense hard needle systems[54], and is used here to efficiently equilibrate the system, which allows us to directly obtain the potential of mean force and, thus, the effective interaction between two hard spheres in a nematic hard needle host.

The simulation will reveal a surprising and robust tendency for chaining of disks along the director axis which seems to contradict the chaining in a angle with respect to the director axis as predicted by quadrupolar elastic interactions in two dimensions[45]. This turns out to be a result of a dominant short-range depletion interaction, and the analysis and explanation of this phenomenon are the main topic of the present paper.

Results

We use hard needles as a model system for a lyotropic LC. Hard needles can be viewed as two endpoints connected by an infinitely thin, hard line, see Fig. 1a, and we sample in the MC simulation by moving the two endpoints. The needle length is used as unit length for non-dimensionalization in the following (primed quantities are dimensionless). In two dimensions needles can order at sufficiently high area densities in a Kosterlitz-Thouless transition into a quasi-ordered nematic phase[54-56]. In the following, we focus on a needle system with well within this nematic regime. We suspend hard disks of diameter () as colloidal particles into the nematic hard needle system, see Fig. 1b. The interaction between all particles, i.e., disk-disk, disk-needle and needle-needle, is given by hard core potentials, which prohibit overlaps. For fast sampling, we employ the rejection-free event-chain MC algorithm. More details on the simulation are provided in the Methods section. The resulting needle-mediated effective interactions between disks in this lyotropic system is the focus of our investigation.

Figure 1

(a) Simulated particles: hard disks with diameter and hard needles with characteristic length . (b) Simulation snapshot of disks () forming chains along the director (approximately in horizontal direction). On the outer disks and between disks surface depletion zones with low needle density form. (c) Disk distance r and disks angle towards the director are used to describe the potential of mean force as well as other interactions. (d) Needle center of mass density distribution around a disk, relative to the director for different relative disk sizes (scale bars are one needle length ). For small disks, there are distinct depletion zones in front of and behind a disk. Depletion zone extensions are . Overlap of these depletion zones gives rise to the density-dependent depletion interaction. (Figure created using gnuplot 5.2 (http://www.gnuplot.info/), matplotlib 3.2.1 (https://matplotlib.org/), python 3.8.2 (https://www.python.org/), inkscape 0.92.5 (https://inkscape.org/)).

Chaining and depletion

In the MC simulation, we observe a chaining of the disks along the director axis of the nematic needle phase, see Fig. 1b. This chaining indicates a strongly anisotropic directional attraction caused by the hard needles along the director axis. The chaining axis coincides with two elongated zones depleted of needles in the direction of the director on the surface of the disk, as revealed by the needle density shown in Fig. 1d. Depletion zones have an extension of order . Therefore, they are deep relative to the disk diameter and very distinct for small disks (). For larger disks, the depletion zones smear out and their depth diminishes as compared to the disks size.

Chaining along the director axis also governs self-assembly in simulations containing many disks shown in Fig. 2a,b. Here we observe self-assembly of hexagonal disk crystals via intermediate chain formation and subsequent aggregation of parallel chains. The chaining along the director axis is very robust and occurs up to high needle densities , see Fig. 2c,d, where the nematic LC becomes very stiff.

Figure 2

(a) and (b) Evolution of a multi-colloid system, containing 40 disks with diameter in (a) 160 disks with diameter in (b). In both systems the same area is covered by disks, the needle density is . To generate the initial configuration we place the disks overlap-free and add needles where they fit until the desired density is achieved. The series of snapshots shows the self-assembly of the disks first in chains and then into hexagonal clusters. We show partly periodic images, extending the original system indicated by a black square. (c) and (d) Also in high density systems () disks align with the director. (c) and (d) have different initial configurations, shown in the left column with disks at 0 and 45 angles. Right column shows equilibrated state. (Figure created using gnuplot 5.2 (http://www.gnuplot.info/), inkscape 0.92.5 (https://inkscape.org/)).

Effective disk-disk interaction

To quantify the effective interaction between two hard disks suspended in nematic hard needles, we use the potential of mean force . Here, r is the distance between the disks and the angle between the connecting line and the director , see Fig. 1c; is the pair distribution function and the thermal energy. The event-chain MC simulation technique is ideally suited for dense mixed colloidal systems[54, 57] and allows us to obtain the positional and angular dependence of the potential of mean force as shown in Fig. 3.

Figure 3

Contour plots of the measured interaction U, the calculated density-dependent depletion interaction and the difference of the two , showing the quadrupolar behavior. Columns represent the different interactions and rows describe different disk diameters , 3 and 5. The measurement was performed in a square system with system size and needle density . On the right side are exemplary snapshots of configurations, marked with the black dot in the contour plot next to it. (Figure created using matplotlib 3.2.1 (https://matplotlib.org/), python 3.8.2 (https://www.python.org/), inkscape 0.92.5 (https://inkscape.org/)).

For small angles the interaction is attractive and has its minimum at the disk surface at an angle of . This explains the observed parallel chaining of disks. For larger disks, the interaction also develops a repulsive part around . The range of the interaction is decreasing relative to the disk size and of order of .

The effective interaction can be described as the sum of two interactions, which are the depletion interaction of short range resulting from overlapping depletion zones around a hard disk, and a power-law quadrupolar elastic interaction resulting from the elastic energy of nematic hard needles. Because hard needles tend to align tangentially at a hard wall, we expect a quadrupolar elastic interaction falling off as and with a angular characteristic[43, 45, 58]. Therefore, the elastic interaction is maximally attractive for and repulsive both for and . The repulsive interaction at can thus be explained by dominating quadrupolar elastic interactions, which becomes more relevant for larger disks. The interaction minimum at , however, comes as a surprise in view of the quadrupolar interaction of standard LC field theory. This minimum can only be explained by the dominance of the attractive depletion interaction induced by the nematic hard needles, which must be directed along the depletion zones parallel to the director.

Density-dependent depletion interaction

In order to calculate the depletion interaction caused by the complex anisotropic and smeared shapes of the depletion zones in Fig. 1d, we use a density-dependent description of depletion interactions and findwhere is the local needle density and the average needle density. The depletion interaction is proportional to a generalized overlap area of depletion regions of two disks at distance , where the depletion region of a disk at is given by the region with , which can correctly capture the complex shaped depletion zones in Fig. 1d. A detailed derivation of Eq. (1) is presented in the Methods section.

The depletion zones are mainly shaped by the elastic interactions in the nematic phase via anchoring conditions at the disk surface. Hard rods exhibit planar anchoring at a hard disk. In the simulation the anchoring is only parallel on average, and fluctuations weaken the anchoring considerably. We only find a weak entropic planar anchoring (quantified below) such that the needle orientation only deviates little from the director orientation at the disk surface. This explains the elongated depletion zones of extension of order (see Fig. 1d). This also results in an overlap area from standard geometric arguments[5], such that is expected. Moreover, this gives rise to a strongly anisotropic depletion attraction, which is maximal parallel to the director (). We can calculate the depletion interaction numerically using measured density profiles from simulations and (1); the result is shown in Fig. 3.

Revealing the residual elastic interaction

In order to uncover additional effective interactions apart from depletion, we subtract the numerically calculated depletion interaction from the measured total potential of mean force and obtain the residual interaction . All three interactions are shown in Fig. 3. This reveals the presence of another interaction, which can actually be identified as the elastic interaction from LC field theory. Hard needles tend to align tangentially at a hard wall, such that we have planar boundary conditions. For two disks suspended in a nematic LC with planar boundary conditions, LC field theory predicts[43, 58]This result is based on the one-constant approximation, i.e., assuming an isotropic elasticity of the nematic LC with a single elastic constant K. This elastic constant can be calculated for hard needles in two dimensions as by adapting a method from Straley[59]. The quadrupolar interaction (2) is shown in Fig. 4a and explained in more detail in the Methods section.

Figure 4

(a) Quadrupolar interaction (2) of two disks in a liquid crystal with planar boundary condition on the disk’s surface (with , see text). This is a field theoretical result for the far-field behavior, assuming weak distortions of the director and using a one-constant approximation. (b) Numerically calculated quadrupolar interaction between two disks in a system of size in the presence of elastic anisotropy and weak anchoring ( and , see (4), and ). (Figure created using matplotlib 3.2.1 (https://matplotlib.org/), python 3.8.2 (https://www.python.org/), inkscape 0.92.5 (https://inkscape.org/)).

The field-theoretical quadrupole interaction (2) is qualitatively similar to the residual interaction from Fig. 3 but exhibits characteristic deviations. The main differences are: (1) a much weaker interaction strength, (2) a distorted quadrupolar symmetry, and (3) for small disks (), a missing repulsive regime around and . The low interaction strength of the measured interaction (i) is caused by the weak parallel anchoring of hard rods at the disk surface, which we find in the simulation (and which we quantify below). In the theoretical result (2), on the other hand, strict planar boundary conditions, i.e., infinitely strong parallel anchoring is assumed. The distorted quadrupolar symmetry (ii) becomes manifest in differently strong repulsion zones for and . This is caused by the elastic anisotropy of the nematic hard needle phase, i.e., different elastic constants for orientation fluctuations perpendicular () or parallel () to the director, which are not captured in the one-constant approximation () employed in (2). Adapting the method from Ref.[59], we find for hard needles in two dimensions (our above value is actually defined as . From fitting mean nematic orientation field around a disk we confirm a value of below. The weak repulsion for small disks () at and (iii) can be explained by the pronounced density dependence of the elastic constant () via a correlation effect with the depletion interaction. Since the depletion zones are of very low density and pronounced in this direction, especially for small disks, they also weaken the elastic interaction in this direction.

This interpretation is supported by a numerical calculation of the elastic interaction between two disks in the presence of elastic anisotropy and weak parallel anchoring (as quantified below), which is shown in Fig. 4b. Details of the numerical calculation are explained in the Methods section. The numerically calculated elastic interaction only misses the correlation effect (iii) and, indeed, resembles the residual interaction from Fig. 3 more closely in its distorted shape.

Both depletion and quadrupolar interaction are attractive and promote chaining of disks, however, at angle if depletion dominates and around if the elastic quadrupolar interaction dominates. We see a robust -chaining which points at a dominating depletion interaction. In order to explain this dominance we need to quantify the weak anchoring for the hard needle and disk system, which is central both for directing the depletion interaction in direction as well as for weakening the residual elastic interaction.

Weak anchoring strength and elastic anisotropy in quadrupolar interaction

Both anchoring strength and elastic anisotropy can be quantitatively analyzed by investigating the needle orientation field around a single disk, which is suspended in a nematic phase with asymptotically at the system boundary ( is the angle of the local director field with this asymptotic orientation, i.e., ; is defined with respect to the asymptotic director orientation, which determines the x-axis). By comparing numerical MC results for hard needles and field theory calculations we can deduce both anchoring strength and elastic anisotropy.

The free energy of a two-dimensional LC with anchoring on a surface is given by the Frank-Oseen elastic energy[60] and a surface anchoring potential:The first integral is the approximate Frank-Oseen free energy with elastic constants and (see Methods section) and the second integral over the disk’s surface represents the anchoring energy with anchoring strength W. The anchoring depends on the boundary condition at the disks surface, which is parallel anchoring for hard needles at a hard disk (i.e., ). The numerical minimization of this free energy is discussed in the Methods section. Introducing dimensionless quantities (marked with a tilde) by measuring the free energy in units of the elastic constant and distances in units of the disk radius , we find two control parameters, the elastic anisotropy and the relative anchoring strength .

We fit the free energy minimization results to MC simulation data for the mean orientation field of the needles (see Methods section). Fitting in the range and doing so for nematic densities and disk sizes in systems of size we obtain the best fit forIn principle, the relative anchoring strength w can depend on and . We only find a very weak decrease from to 0.45 with density for . The result for the elastic anisotropy is in good agreement with our above finding using the method from Ref.[59]. Both results (4) are also very similar to results from Ref.[61] for granular rods in two dimensions obtained in a cavity geometry. This confirms that the hard needle nematic phase has indeed a pronounced elastic anisotropy and only a weak effective anchoring strength at the hard disk surfaces; this anchoring is purely entropic as reflected by implying .

Depletion dominates quadrupolar interaction for chaining along director

In the chaining of disks the elastic quadrupolar interaction competes with the depletion interaction. For weak anchoring the quadrupolar interaction is [45]. From Eq. (2) we expect a quadrupolar interaction strength at the disk surface. This competes with a depletion interaction of strength . For weak anchoring with the depletion interaction dominates as observed in our simulations. This explains the robustly observed -chaining along the director. At high densities, the quadrupolar interaction becomes more relevant. We performed simulations up to very high needle densities and still observe robust -chaining as shown in Fig. 2.

Discussion

We modeled lyotropic LC colloids as hard disks in a suspension of hard needles in a two-dimensional system. Simulations with an event-chain MC algorithm showed a chaining of disks along the director axis of the nematic needle phase. This chaining is caused by a depletion interaction due to elongated depletion zones behind and in front of the disks parallel to the director. This depletion interaction is only accessible in efficient microscopic simulations with explicit rods. Elongated depletion zones are caused by the weak planar anchoring of hard needles at hard disks. Calculating the depletion interaction with a density-dependent depletion theory and subtracting the depletion interaction we reveal a residual elastic quadrupolar interaction which is mediated by the director distortions around the disks. A quadrupolar elastic interaction is consistent with the planar anchoring. The elastic interaction is weakened because hard needles are only weakly anchored and the angular dependence is characteristically deformed because of a pronounced elastic anisotropy of the nematic needle phase.

Both depletion and quadrupolar interaction are attractive and can give rise to chaining, in principle, but the depletion favors a -angle with the director axis while the quadrupolar interaction would favor . In our simulations for densities up to and disk sizes up to we only find -angle chaining indicating a robustly dominating depletion, which is mainly due to the weak planar anchoring. This type of chaining also governs the intermediate states in self-assembly of many colloidal disks into crystals as shown in simulation in Fig. 2.

A natural continuation of our work will be the investigation of hard spheres suspended in hard rods (such as spherocylinders) in three spatial dimensions. We believe that the two-dimensional simulation can already provide results, which qualitatively apply also to the the three-dimensional case, while giving a significant simulation speedup.

In three spatial dimensions, experiments[47, 48] and numerical field theory[44, 46] found quadrupolar interactions with an interaction minimum and chaining of colloidal spheres at a -angle for planar anchoring, which corresponds to a in two dimensions[45]. The field-theoretic continuum approaches did not capture depletion interactions, however, while experiments in Refs.[47, 48] were not dealing with lyotropic systems.

With respect to applications, our results show that, in lyotropic LC colloids, depletion interactions can play an important role in structure formation. By exploiting their dependence on particle shape they could provide an additional tool to control the structure formation process. Assuming that our results carry over to three-dimensional systems, we expect that the range of the depletion interaction mediated by the lyotropic LC is given by the rod length . In order to develop similar observable consequences for chaining as in the simulations, should not be orders of magnitude smaller than colloidal spheres; this requires fairly long lyotropic rod-like particles. The depletion interaction dominates for weak planar anchoring as realized by hard rods and hard spheres. Experimental investigations for these types of systems are still rare at present. An ideal system to test the predictions experimentally are fd virus rods. In the experimental study in Ref.[39] very small hard spherical particles () in a host of hard fd virus rods (length ) in the nematic phase have been studied and, indeed, an ordering parallel to the director has been found, which was attributed to depletion attraction and is in line with our observations. Our results suggest that future experiments with micron-sized colloidal particles suspended in a nematic fd virus phase could display dominant depletion interaction effects including colloidal self-assembly via parallel chaining.

Methods

Simulation method

Event-chain Monte-Carlo

To simulate hard disks in a suspension of hard needles the event-chain MC algorithm for many-body interactions is used[53, 54, 62]. This is a rejection-free MC method that performs very well for dense systems[57], which makes it an excellent choice for needle systems in the nematic phase, for which it has been adapted in Ref.[54]. The event-chain MC is a Markov-chain MC scheme, which fulfills maximal (rejection-free) global balance rather than detailed balance. Global balance is achieved by introducing lifting moves[62]. For hard spheres or needles a lifting move is the transfer of MC displacement from one particle to another particle. This means in a MC move only one particle at a time is moved along a line until it contacts another object. For the application to hard needles, one of the two endpoints are moved and, upon collision with another needle, the remaining MC move distance is transferred to one of the endpoints of this needle. To which of the endpoints it is transferred depends on the collision point along the needle[54]. In the presence of additional disks, MC displacement is also transferred to disks if a needle collides with the disk and vice versa. Collision detection is the computational bottleneck, and we use a sophisticated neighbor list system to achieve high simulation speeds also in the nematic phase (see Supplementary Information).

Hard needle liquid crystal

We model the molecules of a lyotropic LC as hard needles, which consist of two endpoints connected by an infinitely thin, hard line (see Fig. 1a). For efficient sampling, the distance of the two endpoints, i.e. the length l of the hard needle can fluctuate around its characteristic length in order to allow for independent motion of both endpoints in the MC simulation; the needle length l is restricted by a hard potential with for and infinite else, such that is the effective needle length. Restricting the needle length is also necessary for our neighbor lists, which accelerate the simulation.

The interaction potential between two hard needles is either infinite when they overlap or zero else. This results in a lyotropic transition from isotropic to nematic upon increasing the needle density . In a system with needles the density is defined using the available free area (reduced by the area occupied by disks). The simulation box is a square with edge length L and periodic boundary conditions. In a typical system of size we simulate needles at a needle density .

Liquid crystalline ordering is described by the standard tensorial order parameterwhich is invariant under needle inversion; is the orientation of needle i.

The scalar order parameter is calculated by diagonalizing the matrix, where is the biggest eigenvalue. The corresponding eigenvector is the director . The scalar order parameter S measures the degree of alignment in the system and is zero in a perfectly isotropic phase and unity in a perfectly aligned nematic phase. Figure 5 shows the order parameter S as a function of density for the two-dimensional hard needle system[54]; the system quasi-orders above [55]. In our two-dimensional system, this transition is a Kosterlitz-Thouless transition into a quasi-ordered nematic state[56].

Figure 5

(a) Measurement of the probability distribution of angles with the director (black dots). The probability function (blue line) matches the MC data for . (b) Scalar order parameter S as a function of the needle density . The transition from the isotropic to the nematic phase occurs roughly at and is of Kosterlitz-Thouless type. We consider needle systems with densities as in the nematic phase and use for the potential of mean force calculation and local needle density measurements. (Figure created using matplotlib 3.2.1 (https://matplotlib.org/), python 3.8.2 (https://www.python.org/), inkscape 0.92.5 (https://inkscape.org/)).

Elastic energy, weak anchoring and elastic anisotropy

For a nematic LC the state of perfect alignment is the energetically preferred configuration. But this state is almost never reached because of thermal fluctuations, boundary conditions, or topological defects. Often effective free energies from field theories are used to describe deviations from homogeneous alignment. For a director field in three dimensions, the energy of small distortions can be captured with the Frank-Oseen free energy [63] with three elastic constants , which determine the energy cost of the three different distortions splay (), twist () and bend (). In a two-dimensional system, twist distortions are not possible and, therefore, resulting in[60]with the needle orientation field related to the director via . We assume corresponding to asymptotically at the system boundary. The last approximation is the leading order in an expansion in deviations from and neglects non-linear terms in [61].

The elastic constants and can be calculated by adapting a method described by Straley[59] to two dimensions. Strictly speaking the pronounced logarithmic orientational fluctuations in two dimensions will give rise to a renormalization at large scales[64], which we ignore in our finite size system. We obtain and for the nematic phase with . This result is based on the angular distribution of needle orientations with matching our MC simulations in the nematic phase with (see Fig. 5). In the one-constant approximation, we assume resulting in ; for this value of K the quadrupolar interaction (2) is shown in 4.

In order to describe colloidal disks in a nematic needle phase faithfully, we need too take into account elastic anisotropy and a finite parallel anchoring with an anchoring strength W. In order to quantify the elastic anisotropy and the relative anchoring strength , we consider the needle orientation field around a single disk, which is suspended in a nematic phase with asymptotically at the system boundary and quantitatively compare MC simulation data with minimization of the free energy (3) containing both the anisotropic 2D Frank-Oseen elastic energy and the anchoring potential. Free energy minimization using the linearized approximation in the Frank-Oseen part (5) results in the partial differential equationin dimensionless units . Analytical solutions are only available in the one-constant approximation , where[65](both and are defined with respect to the asymptotic director orientation, which determines the x-axis). For the full problem (6) we have to resort to numerical solutions by finite element methods (using MATHEMATICA).

In Fig. 6 we plot MC simulation data for the director orientation field as a function of (i.e., along circles) for different rescaled and for different system sizes such that remains fixed. Rescaling of lengths with results in good data collapse for different . The characteristic shoulder of as a function of can only be explained by an elastic anisotropy . Fits with Eq. (7) from elastically isotropic one-constant approximation () remain unsatisfactory. We perform fits for specific and overall fits within the whole range with the numerical solution of the PDE (6) in the same geometry as the MC simulation (square box with and director oriented along the diagonal as fixed by Dirichlet boundary conditions) using and w as fit parameter. The fit results for w are weakly r-dependent (dashed blue lines) and the result of the overall fit (solid blue line) is , see Eq. (4).

Figure 6

MC simulation results for as a function of and for for parameters , , and . Data for different collapse fairly well. Blue dashed lines: Least square fits of data for each with numerical solution of the PDE (6) in the same geometry as simulation, for , and with w as fit parameter. Blue solid lines: Analogous least square fit of data for all . Yellow lines: Least square fits with one-constant result (7) () with w as fit parameter. (Figure created using Wolfram Mathematica 12.0, inkscape 0.92.2 (https://inkscape.org/)).

Quadrupolar elastic interaction

If two colloidal disks are embedded, boundary conditions on the disk surface induce deformations in the director field and, thus, induce an elastic interaction given by the energy (5) of the deformation. There are some results on these elastic interactions in the one-constant approximation (). For two disks with perfect planar boundary conditions there is an approximate analytical result which is the quadrupolar long-range interaction from Eq. (2)[43, 58]. Since this is a field theoretical result, the needles are assumed to be infinitely small relative to the colloids. There are no analytical results for the full two-constant free energy (5).

Therefore, we performed numerical simulations based on the finite element method (using MATHEMATICA), i.e., solving Eq. (6) for a system with two disks with distance r and with an angle with respect to the director axis and using the same anisotropy and weak anchoring strength that we determined numerically. We place both disks symmetrically in a system in a square box with and with a director axis oriented parallel to one edge of the simulation box as fixed by Dirichlet conditions at the system boundaries. The total energy of the system is the sum of elastic and anchoring energy (see Eq. (3)). The interaction energy is obtained as the difference in energy between a system containing two disks and the non-interacting system as obtained by the sum of the energies of two systems containing only one of the disks each. This results in Fig. 4b.

We use the results of Biben et al.[66] and generalize them to anisotropic depletants with a rotational degree of freedom to get a density-dependent depletion interaction for disks in a suspension of hard needles. We consider a system of hard disks with positions and hard needles with positions and orientations . Upper case indices refer to disks, lower case indices to needles. The energy of the system is given byThe disk-disk interaction is given by , the needle-needle interaction by and the disk-needle interaction by . By integrating over the needle degrees of freedom one can derive the effective part of the interaction between the disks [66],(). The corresponding force on disk K is given byHere, we used the single particle density of needles with angle at for fixed disk positions:By using , evaluating the sum to a factor and defining the average over the needle angles as , we obtainfor the case of two disks at and . We use the superposition approximation , where is the density distribution around a single disk. For a single disk the needles are distributed according to the direct interaction potential ,This finally leads toThis effective potential is the density-dependent depletion interaction, which we further approximate by employing a factorization approximation for the angular averages, , which is valid for the isotropic phase and the ideal nematic phase. Since we investigate the effective interaction in the nematic phase this should be a good approximation. This gives our final resultMore intermediate steps of the derivation are given in the Supplementary Information.

Supplementary information