format: text

text:
# Isotropic-nematic phase behavior of length polydisperse hard rods

## I Introduction

Dispersions of highly anisometric rod- or platelike colloidal particles exceeding a certain concentration are known to undergo a transition from an isotropic state (I), in which the particles are oriented in random directions to an aligned nematic state (N) Zocher ; Langmuir ; Bernalnature . The theoretical basis for this phenomenon has been established by Onsager in the 1940s. In a seminal paper Onsager , he showed that the phase transition can be explained on the basis of purely repulsive interactions between the particles. The basic mechanism behind the transition is a competition between orientional entropy which favors the isotropic state and the entropy effect associated with the orientation-dependent excluded volume of the anisometrical particles which favors the ordered nematic.

A main characteristic of systems of (anisometric) colloidal particles is their inherent polydispersity, i.e. the particles may differ in size and shape Wierenga ; Buiningwater . The issue of polydispersity and its effect on the interpretation of experimental results has already been addressed by Onsager in his original paper. Later on, extensions of the Onsager theory allowing for phase diagram calculations for bidisperse Lekkerkerker84 ; OdijkLekkerkerker ; birshtein ; LekkerVroeg ; vanroij96/2 ; vanroij98  and tridisperse systems LekkerVroegtridisperse  of hard rods as well as binary mixtures of hard platelets wensinkbidikte  revealed a rich variety in behavior, most notably a widening of the coexistence region, a fractionation effect (i.e. segregation of the species among the coexisting phases), a reentrant phenomenon and, most interestingly, the possibility of a demixing of the nematic phase which may give rise to isotropic-nematic-nematic triphasic equilibria.

So far, very few theoretical attempts have been made to study the isotropic-nematic phase behavior of truly polydisperse systems, i.e. systems characterized by a continuous distribution in particle size, within the Onsager treatment. As to rodlike systems, there are some studies in which the effect of polydispersity on the isotropic-nematic transition is accounted for using perturbation theories sluckin ; chen . The major drawback of this approach, in which the effect of polydispersity is considered as a small perturbation to a monodisperse reference system, is its limited applicability; it is only justified for slightly polydisperse systems with very narrow length distributions. Although the perturbation approach qualitatively predicts some generic features such as a broadening of the coexistence area and a fractionation effect, with the longest rod going preferentially into the nematic phase, other interesting phenomena which are expected to occur at much higher polydispersities (in particular polydispersity-induced demixing transitions) are not treated and therefore remain elusive.

Solving the phase equilibrium conditions for systems with arbitrary size distributions is by no means trivial and requires considerable numerical effort. In particular, the presence of almost infinitely many components in a polydisperse system requires an equally large number of coexistence conditions to be solved simultaneously which obviously is a formidable task. Recently, a number of studies have appeared in which the fully polydisperse Onsager model, albeit in simplified form, was subjected to a numerical treatment. Speranza and Sollich sollichonsagerP2 ; sollichonsagerP2fattail  investigated the model using the so-called $`𝒫_2`$-approximation which consists of truncating the angular dependence of the orientation-dependent excluded-volume term after the first nontrivial term in an expansion in terms of even Legendre polynomials $`𝒫_{2n}`$. A remarkable outcome of these calculations is that for rod length distributions with sufficiently fat tails (e.g. log–normal distributions) triphasic isotropic-nematic-nematic equilibria are predicted to occur in a small interval of polydispersities of the parent system. However, the simplified orientation distribution function (ODF) for the nematic phase pertaining to the ‘$`𝒫_2`$-model’ is only valid for the description of very weakly aligned nematic phases. The behavior predicted from this model should therefore be considered with some care, particularly in those regions where the fractionation effect is strong and the phase behavior is dominated by the effect of the longest rods in the system. The presence of very long rods in a nematic state may force the entire system into a strongly aligned nematic configuration, so that a more appropriate form for the nematic ODF is required in these cases.

In this paper we use the Gaussian trial ODF approach to calculate the isotropic-nematic phase behavior of hard rod systems which can be described by either a Schulz or a log-normal length distributions with arbitrary polydispersities. The benefit of using the Gaussian Ansatz is twofold. First, all necessary integrals for the monodisperse Onsager model are analytically tractable so that only numerical integrations over the length distributions need to be considered for the polydisperse case. Second, the Gaussian ODF allows for a qualitatively better description of highly ordered nematic states compared to the $`𝒫_2`$-approximation which makes it a suitable tool for describing polydisperse systems, particularly the ones with a ‘fat-tailed’ length distribution. While work on this subject was still in progress, Speranza and Sollich reported a numerical analysis of the exact Onsager model sollichonsagerexact , i.e. using the numerically exact ODF. Also there, triphasic equilibria were predicted for both Schulz-type parent distributions and ‘fat-tailed’ log-normal forms. However, due to the numerical complexity of the problem only the onset of nematic ordering from an isotropic reference phase was considered there so that no information could be obtained about the properties of the isotropic and nematic phases across the coexistence region. Consequently, no conclusive insight could be gained as to whether the triphasic equilibria constitute a significant part of the phase diagram. Within the Gaussian approach it is possible to access the coexistence region with only limited additional numerical effort. An important consequence is that it also enables us to study the triphasic demixing in more detail and to gain insight in the extent of its stability region. Although we cannot calculate the binodal curves for these equilibria (which locate the precise onset), we are able to localize the spinodal points for the coexisting nematic phase across the two-phase region. The presence of these spinodal points indicates a local instability of the nematic phase which implies that a triphasic demixing transition must take place.

This paper is structured as follows. In Sec. II we describe the Onsager theory generalized to polydisperse systems. The conditions for phase equilibria are outlined in general terms in Sec. III-A, and specified to the onset of isotropic-nematic phase separation and the full phase split situation in Secs. III-B and III-C respectively. In III-D some details of the truncated parent distributions we used in our calculations are given. The numerical results for the different parent distributions will be presented and discussed in detail in Secs. IV and V. Finally, we summarize and conclude in Sec. VI. Technical aspects of the subject are treated in several appendices. In App. A we provide details about the numerical procedure we adopted to solve the coexistence conditions. In App. B we explicitly show that the Gaussian Ansatz yields the exact scaling results for parent distributions with infinite cutoff lengths and in Sec. C we establish a spinodal instability criterion for the nematic phase inside the coexistence region which is used to detect possible isotropic-nematic-nematic triphasic equilibria along the coexistence trajectory.

## II Polydisperse Onsager theory; starting equations

Let us consider a system of hard rodlike cylinders with equal diameters $`D`$ but different lengths $`L`$ , in a macroscopic volume $`V`$. To characterize the rod lengths in our polydisperse system we introduce the relative rod length $`l=L/L_0`$ (with $`L_0`$ some reference rod length $`L_0`$) which is assumed to be a continuous variable. We may then take the limit $`L_0/D\mathrm{}`$ (infinitely thin rods) at constant values for the relative lengths $`l`$. In the Onsager approach, the excess free energy describing the excluded volume interactions between the particles is truncated after the second virial term. This approach can be shown to yield the exact free energy in the limit of infinitely thin rods. A generalization of the original Onsager model to include polydispersity leads to the following expression for the total Helmholtz free energy density $`f`$ (in units $`k_BT\beta ^1`$)

$`f`$ $``$ $`{\displaystyle \frac{b\beta F}{V}}{\displaystyle c(l)[\mathrm{ln}c(l)1]𝑑l}+{\displaystyle c(l)\omega (l)𝑑l}`$ (1)
$`+{\displaystyle c(l)c(l^{})ll^{}\rho (l,l^{})𝑑l𝑑l^{}}.`$
All irrelevant contributions linear in $`c`$ arising from the standard chemical potentials of the particles are omitted since they only depend on the solvent chemical potential and the temperature. The concentrations $`c`$ are rendered dimensionless by relating them to the orientationally averaged excluded volume per particle between two reference rods, $`b=\pi DL_0^2/4`$, via $`c(l)=bN(l)/V`$ where $`N(l)dl`$ is the number of particles with relative length between $`l`$ and $`l+dl.`$ The density distribution over lengths $`c(l)`$ can be decomposed according to $`c(l)=c_0p(l)`$, with $`p(l)`$ the normalized length distribution ($`𝑑lp(l)1`$) and $`c_0`$ the total dimensionless rod concentration.

The free energy Eq. (1) consists of several entropic contributions. The first term represents the exact ideal free energy of the polydisperse system. The second term contains parameter $`\omega `$ as a measure for (the negative) of the orientational entropy Onsager

$$\omega (l)\psi (l,\mathrm{\Omega })\mathrm{ln}[4\pi \psi (l,\mathrm{\Omega })]𝑑\mathrm{\Omega },$$
(2)
where $`\psi (\theta ,l)`$ is the normalized ODF for species $`l`$ describing the distribution of the particles’ solid angle $`\mathrm{\Omega }`$. In the isotropic state, all orientations are equally probable so that $`\psi _{\text{iso}}`$ is simply a constant ($`1/4\pi `$) independent of $`l`$. In the nematic state, however, the ODFs are peaked functions (generally different for each species $`l`$), due to the fact that the rods are aligned along a nematic director. Note that $`\omega `$ (and hence the orientation free energy) attains its minimum ($`\omega =0`$) in the isotropic state, whereas $`\omega >0`$ in the nematic state.

The last term in Eq. (1) describes the excess free energy which accounts for the particle interactions. In a second virial approximation, the interactions between hard particles may be expressed as an excluded volume entropy depending on the excluded volume between two particles. A measure for the average excluded-volume interaction between rods of relative length $`l`$ and $`l^{}`$ is given by the following angular average

$$\rho (l,l^{})\frac{4}{\pi }\left|\mathrm{sin}\gamma (\mathrm{\Omega },\mathrm{\Omega }^{})\right|\psi (l,\mathrm{\Omega })\psi (l^{},\mathrm{\Omega }^{})𝑑\mathrm{\Omega }𝑑\mathrm{\Omega }^{}.$$
(3)
Using the isotropic average $`\left|\mathrm{sin}\gamma (\mathrm{\Omega },\mathrm{\Omega }^{})\right|=\pi /4`$ we obtain $`\rho (l,l^{})1`$ for the isotropic state. This indicates that the excluded volume (or packing) free energy is indeed maximal in the isotropic phase but decreases as soon as the rods align to form a nematic phase.

In this study we use Gaussian trial ODFs with variational parameter $`\alpha (l)`$ to describe the angular distribution of rods with relative length $`l`$ in the nematic state OdijkLekkerkerker . The Gaussian Ansatz consists of supposing

$$\psi (l,\theta )\{\begin{array}{cc}\frac{\alpha (l)}{4\pi }\mathrm{exp}[\frac{1}{2}\alpha (l)\theta ^2]\hfill & 0\theta \frac{\pi }{2}\hfill \\ \frac{\alpha (l)}{4\pi }\mathrm{exp}[\frac{1}{2}\alpha (l)(\pi \theta )^2]\hfill & \frac{\pi }{2}\theta \pi \hfill \end{array},$$
(4)
where $`\alpha `$ is now a function of $`l`$. Note that, due to the uniaxial symmetry of the nematic phase, the ODFs only depends upon the polar angle $`\theta `$ between the particle orientation vector and the nematic director. Inserting Eq. (4) in Eq. (2) and straightforward integration yields for the orientational entropy

$$\omega (l)\mathrm{ln}\alpha (l)1.$$
(5)
For the excluded volume entropy in the nematic phase $`\rho _{\text{nem}}(l,l^{})`$ only the leading order term of its asymptotic expansion for large $`\alpha `$ will be retained

$$\rho _{\text{nem}}(l,l^{})\sqrt{\frac{8}{\pi }\left(\frac{1}{\alpha (l)}+\frac{1}{\alpha (l^{})}\right)}+𝒪[\alpha ^{3/2}(l),\alpha ^{3/2}(l^{})].$$
(6)
Substituting Eqs. (5) and (6) into Eq. (1) and minimizing the free energy density with respect to the non-conserved orientational degrees of freedom by means of a functional differentiation with respect to $`\alpha (l)`$ gives

$`{\displaystyle \frac{\delta f}{\delta \alpha (l)}}`$ $`{\displaystyle \frac{c(l)}{\alpha (l)}}\left({\displaystyle \frac{8}{\pi }}\right)^{1/2}{\displaystyle \frac{lc(l)}{2\alpha ^2(l)}}`$ (7)
$`\times {\displaystyle }l^{}c(l^{})({\displaystyle \frac{1}{\alpha (l)}}+{\displaystyle \frac{1}{\alpha (l^{})}})^{1/2}dl^{}.`$
Applying the stationarity condition $`\delta f/\delta \alpha (l)0`$ and some rearranging leads to the following self-consistency equation

$$\stackrel{~}{\alpha }(l)=2l^2\left\{l^{}p^{(N)}(l^{})\left[1+\frac{\stackrel{~}{\alpha }(l)}{\stackrel{~}{\alpha }(l^{})}\right]^{1/2}𝑑l^{}\right\}^2,$$
(8)
Here, we have factorized the Gaussian variational parameter function $`\alpha (l)`$ into a concentration-dependent part and a contribution $`\stackrel{~}{\alpha }(l)`$ only related, via Eq.(8), to the normalized length distribution in the nematic phase $`p^{(N)}(l)`$. Hence we write

$$\alpha (l)=\stackrel{~}{\alpha }(l)\frac{4c_0^2}{\pi }.$$
(9)
showing that for all $`l`$ the variational parameter $`\alpha `$ depends quadratically on $`c_0`$ just as in the monodisperse case Vroege92 . An approximate analytical solution to Eq. (8) valid for infinitely narrow distributions sluckin  (denoted by subscript $`\delta `$) can be obtained by substituting a delta function $`p^{(N)}(l)=\delta (l1)`$ which gives

$$\stackrel{~}{\alpha }_\delta (l)=\frac{1}{2}\left(\sqrt{8l^2+1}1\right).$$
(10)
This result may be interpreted as a measure for the nematic alignment of a single rod with relative length $`l`$ added to a nematic bulk system of monodisperse rods with reference length $`L_0`$. Eq. (10) shows that $`\stackrel{~}{\alpha }_\delta (l)`$ and hence the order parameter Vroege92  $`P(l)13/\alpha (l)`$ are in general, as we might have anticipated, increasing functions of the relative rod length, i.e. $`\stackrel{~}{\alpha }_\delta (l)l`$ for large $`l`$. Moreover, $`\stackrel{~}{\alpha }_\delta (0)=0`$ which means that there is no ordering for rods of zero length, as formally must be the case. However, it should be pointed out that rods with lengths close to zero must be excluded from our model because the normalization factors for the Gaussian ODFs in Eq. (4) do not allow for a correct description of isotropically distributed or weakly aligned species in the nematic state. For consistency reasons we must therefore introduce a lower limit ($`l_{\mathrm{min}}>0`$) in all length distributions.

## III I-N phase coexistence

### III.1 Equilibrium conditions for polydisperse systems

The conditions for phase equilibrium are that the coexisting isotropic and nematic phases must have equal chemical potential $`\mu (l)`$ for all relative rod lengths $`l`$, as well as equal osmotic pressure $`\mathrm{\Pi }`$. The chemical potential can be derived by functional differentiation of the free energy with respect to the length distribution $`c(l)`$

$$\beta \mu (l)=\frac{\delta f}{\delta c(l)}.$$
(11)
Using Eqs. (5) and (6) together with the isotropic values $`\omega 0`$ and $`\rho 1`$ we obtain

$`\beta \mu _{\text{iso}}(l)`$ $`=`$ $`\mathrm{ln}c^{(I)}(l)+2lc_1^{(I)}`$
$`\beta \mu _{\text{nem}}(l)`$ $`=`$ $`\mathrm{ln}c^{(N)}(l)+\mathrm{ln}\left[{\displaystyle \frac{4}{\pi }}(c_0^{(N)})^2\stackrel{~}{\alpha }(l)\right]1`$ (12)
$`+\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l),`$
where $`c_1`$ denotes the first moment density following from the definition

$$c_k=c_0m_k=𝑑ll^kc(l),k=0,1,2,\mathrm{}$$
(13)
Here, $`m_k`$ denotes the $`k`$-th moment of the (normalized) distribution. The excess chemical potential for the nematic phase $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ is given by

$$\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)=\text{}2^{3/2}l𝑑l^{}p^{(N)}(l^{})l^{}\left(\frac{1}{\stackrel{~}{\alpha }(l)}+\frac{1}{\stackrel{~}{\alpha }(l^{})}\right)^{1/2}.$$
(14)
and is independent of the concentration of the nematic phase. Similarly to Eq. (8) we can straightforwardly obtain an analytical solution for $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ valid for near monodisperse distributions by substituting $`p^{(N)}(l)=\delta (l1)`$, which yields the following scaling result

$$\stackrel{~}{\mu }_{\text{ex,}\delta }^{(N)}(l)l\sqrt{\frac{1}{\stackrel{~}{\alpha }_\delta (l)}+1}.$$
(15)
Using Eq. (10) it follows that $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)l,`$ for very large $`l.`$ This means that the excess chemical potential (i.e. the reversibel work required to insert a single rod of length $`l`$ in a nematic system of reference rods ) increases with length, which is consistent with intuition. The osmotic pressure can be written in terms of the chemical potential and the free energy via

$$b\beta \mathrm{\Pi }f+\beta 𝑑lc(l)\mu (l),$$
(16)
which immediately yields for the isotropic phase

$$b\beta \mathrm{\Pi }_{\text{iso}}c_0^{(I)}+(c_1^{(I)})^2.$$
(17)
For the nematic phase Eq. (16) this formally gives

$$b\beta \mathrm{\Pi }_{\text{nem}}c_0^{(N)}+f_{\text{ex}}^{(N)}.$$
(18)
However, this result can be simplified considerably by noting that the volume fraction of the average excluded volume (per particle) in the nematic phase is a constant, namely

$$N\frac{V_{\text{excl}}_{l,l^{}}}{V}c_0^{(N)}ll^{}\rho (l,l^{})_{l,l^{}}=2.$$
(19)
The brackets denote averages over the normalized length distribution. This result, which is due to Odijk odijkLC , generally holds for both monodisperse and polydisperse systems, independent of their composition. From the free energy Eq. (1) it then follows that $`f_{\text{ex}}^{(N)}=2c_0^{(N)}`$ so that the osmotic pressure of the nematic phase reduces to

$$b\beta \mathrm{\Pi }_{\text{nem}}3c_0^{(N)},$$
(20)
like for a monodisperse system Vroege92 . We can now state the conditions for the coexistence between the isotropic and nematic daughter phases into which a parent phase (henceforth denoted with superscript 0) with length distribution $`c^{(0)}(l)`$ is assumed to have split. From Eq. (12), equality of chemical potentials of both phases is obeyed exactly if the distributions in the phases have the following form

$$c^{(a)}(l)=W(l)\mathrm{exp}[\xi ^{(a)}(l)],a=I,N$$
(21)
where $`W(l)\mathrm{exp}[\beta \mu (l)]`$ must be a function common to both phases, since $`\mu ^{(I)}(l)=\mu ^{(N)}(l)=\mu (l)`$. The functions $`\xi ^{(a)}(l)`$ are given by

$`\xi ^{(I)}(l)`$ $`=`$ $`2lc_1^{(I)}`$
$`\xi ^{(N)}(l)`$ $`=`$ $`\left(1\mathrm{ln}{\displaystyle \frac{4}{\pi }}\right)2\mathrm{ln}c_0^{(N)}\mathrm{ln}\stackrel{~}{\alpha }(l)\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l).`$
Furthermore, conservation of matter requires

$$c^{(0)}(l)=(1\gamma )c^{(I)}(l)+\gamma c^{(N)}(l),$$
(23)
where $`\gamma `$ denotes the fraction of the system volume occupied by the nematic phase. Using Eq. (23), we can express $`W(l)`$ in terms of the parent distribution $`c^{(0)}(l)`$ which gives

$`c^{(a)}(l)`$ $`=`$ $`c^{(0)}(l){\displaystyle \frac{\mathrm{exp}[\xi ^{(a)}(l)]}{(1\gamma )\mathrm{exp}[\xi ^{(I)}(l)]+\gamma \mathrm{exp}[\xi ^{(N)}(l)]}}.`$
(24)
These functions represent the equilibrium rod length distributions for the coexisting phases. The phase equilibria can now, in principle, be obtained by solving a set of self-consistency equations for the moment densities of both phases and for the functions $`\stackrel{~}{\alpha }(l)`$ and $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ pertaining to the nematic phase. These equations will be worked out below for a specific situation, namely at the onset of isotropic-nematic phase separation.

### III.2 The onset of I-N phase separation; cloud and shadow curves

In this section we aim at locating the onset of isotropic-nematic phase separation indicated by so-called cloud and shadow points. A cloud point marks the density where a parent phase starts to split off an infinitesimal amount of a new coexisting phase, called the shadow phase. Accordingly, at the isotropic cloud point only an infinitesimal amount of nematic phase (shadow phase) has emerged and so the distribution of the isotropic phase is only negligibly perturbed away from the parent. Hence, for the isotropic cloud point we may set $`\gamma =0`$ in Eq. (24) so that,

$$c^{(I)}(l)=c^{(0)}(l),$$
(25)
implying that the distribution in the isotropic phase at the cloud point is equal to the parent distribution and hence $`c_0^{(I)}=c_0^{(0)}`$. The (normalized) rod distribution in the nematic shadow phase (with density $`c_0^{(N)}`$) is now given by

$`p^{(N)}(l)`$ $`=`$ $`{\displaystyle \frac{c^{(0)}(l)}{c_0^{(N)}}}\mathrm{exp}\left[\xi ^{(N)}(l)\xi ^{(I)}(l)\right]`$ (26)
$`=`$ $`𝒦_N{\displaystyle \frac{p^{(0)}(l)}{\stackrel{~}{\alpha }(l)}}\mathrm{exp}\left[2c_0^{(0)}l\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)\right],`$
where $`𝒦_N=`$ $`\pi ec_0^{(0)}/4(c_0^{(N)})^3`$ and $`p^{(0)}(l)`$ the normalized parent distribution. Note that Eq. (26) is an implicit expression for $`p^{(N)}(l)`$ because it still depends on the unknown functions for the variational parameter $`\stackrel{~}{\alpha }(l)`$ and the excess chemical potential $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ for each species in the nematic shadow phase. Explicit solutions for these functions can be obtained by substituting Eq. (26) into Eqs. (8) and (14) and numerically solving the resulting self-consistency equations.

The concentrations of the isotropic cloud phase and the coexisting nematic shadow are found by imposing the normalization condition for the distribution in the nematic shadow phase,

$$p^{(N)}(l)1,$$
(27)
and the condition of equal osmotic pressure

$$3c_0^{(N)}=c_0^{(0)}+(c_1^{(0)})^2.$$
(28)
Using this simple equation to eliminate e.g. $`c_0^{(N)},`$ we may conveniently combine Eqs. (27) and (28) into one self-consistency equation for the concentration of the isotropic cloud point, which we can solve in an iterative fashion. However, since $`\stackrel{~}{\alpha }(l)`$ and $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ also depend on $`c_0^{(0)}`$ \[via $`p^{(N)}(l),`$ Eq. (26)\] this equation has to be solved along with the coupled set of self-consistency equations, Eqs. (8) and (14), so that we end up with a set of three coupled nonlinear equations. Obviously, solving this set is not a trivial task but requires some numerical effort. For this reason we have devoted an Appendix A to this issue in which we describe some details of the numerical procedures adopted in this study.

We can now perform a similar analysis to obtain expressions for the nematic cloud point and the associated isotropic shadow point, which locate the onset of I-N equilibrium coming from a dense nematic parent phase. Since the latter now coexists with an infinitesimal amount of an isotropic shadow phase we may set $`\gamma =1`$ in (24) so that $`c^{(N)}(l)=c^{(0)}(l)`$ and $`c_0^{(N)}=c_0^{(0)}`$. The (normalized) rod distribution in the isotropic shadow phase (with density $`c_0^{(I)}`$) is then given by

$`p^{(I)}(l)`$ $`=`$ $`{\displaystyle \frac{c^{(0)}(l)}{c_0^{(I)}}}\mathrm{exp}\left[\xi ^{(I)}(l)\xi ^{(N)}(l)\right]`$ (29)
$`=`$ $`𝒦_Ip^{(0)}(l)\stackrel{~}{\alpha }(l)\mathrm{exp}\left[2c_1^{(I)}l+\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)\right],`$
where $`𝒦_I=4(c_0^{(0)})^3/\pi ec_0^{(I)}.`$ Since the functions $`\stackrel{~}{\alpha }(l)`$ and $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ correspond to the nematic parent phase we may substitute $`p^{(N)}(l)=p^{(0)}(l)`$ in Eqs. (8) and (14). Consequently, as the normalized parent distributions $`p^{(0)}(l)`$ have a predescribed form and do not depend on any concentration, the coupled self-consistency equations need to be solved only once for a given $`p^{(0)}(l)`$ . The forms of the parent distributions will be specified in Sec. III-D.

Once the solutions for $`\stackrel{~}{\alpha }(l)`$ and $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ have been obtained, the concentrations of the nematic cloud and shadow point can be calculated by requiring self-consistency for the zeroth moment (normalization condition) and the first moment of the isotropic shadow distribution, i.e.

$`{\displaystyle p^{(I)}(l)}`$ $`=`$ $`1,`$
$`{\displaystyle lp^{(I)}(l)}`$ $`=`$ $`m_1^{(I)}={\displaystyle \frac{c_1^{(I)}}{c_0^{(I)}}}.`$ (30)
These conditions can, similarly to the previous case, easily be rewritten into a set of consistency relations for the densities of the nematic cloud and shadow points, which we can solve iteratively (see also Appendix A).

### III.3 Inside the coexistence region

We will now focus on the coexistence region between the isotropic and nematic cloud points, where both phases coexist in noninfinitesimal amounts, i.e. $`0<\gamma <1.`$ According to Eq. (24), the equilibrium length distributions in the coexisting phases are then given by

$`c^{(I)}(l)`$ $`=`$ $`{\displaystyle \frac{c^{(0)}(l)}{\gamma \mathrm{exp}\left\{\mathrm{\Delta }\xi [\stackrel{~}{\alpha }(l),\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)]\right\}+(1\gamma )}},`$
$`c^{(N)}(l)`$ $`=`$ $`{\displaystyle \frac{c^{(0)}(l)}{(1\gamma )\mathrm{exp}\left\{\mathrm{\Delta }\xi [\stackrel{~}{\alpha }(l),\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)]\right\}+\gamma }}.`$ (31)
with $`\mathrm{\Delta }\xi (l)=\xi ^{(N)}(l)\xi ^{(I)}(l)`$, given by Eq. (LABEL:xiparameters). Note that both distributions are now different from the parental one. Solving the coexistence problem is done in a similar way to the one described in Sec. III-B for the isotropic cloud and shadow points. From an experimental point of view, we are only interested in results located on so-called physical dilution lines along which the overall system number density ($`c_0^{(0)}`$) is changed (by e.g. adding or evaporating solvent) while the overall composition of the species ($`p^{(0)}(l)`$) remains fixed. The parent distributions will be specified below.

### III.4 Parent distributions

The numerical method described in the previous sections allows us to calculate the isotropic-nematic phase diagram for in principle arbitrary parent distributions. In our study we specify two types of distributions. The first one is the Schulz distribution which has the form

$$p^{(0)}(l)=Nl^z\mathrm{exp}[(z+1)l],$$
(32)
with normalization factor $`N`$. In order to exclude rods with potentially zero length we truncate the distribution at some lower cutoff length $`l_{\text{min}}`$. Henceforth we fix $`l_{\text{min}}=0.01`$. For calculational purposes (see Appendix A) we must also have some finite cutoff length $`l_{\mathrm{max}}`$ at large $`l`$. Of course, introducing finite cutoff lengths is also reasonable from a physical standpoint. The first and second moment (defined as $`m_k=l^k`$, $`k=1,2)`$ of the Schulz distribution are $`m_1=1`$ and $`m_2=(z+2)/(z+1)`$ only for the unbounded case. However for finite cutoff lengths the moments will deviate from these values. Although the corrections are generally small, in particular for large $`l_{\text{max}}`$, they cannot be neglected. Therefore, we choose to calculate all relevant moments of the parent distribution numerically via $`m_k=_{l_{\text{min}}}^{l_{\text{max}}}l^kp^{(0)}𝑑l`$. The polydispersity $`\sigma `$ is defined as

$$\sigma ^2=\frac{l^2}{l^2}1,$$
(33)
and would yield $`\sigma ^2=(1+z)^1`$ for the unbounded Schulz distribution.

The second distribution we consider is the log-normal one. The ‘fat-tailed’ log-normal distribution decays much slower at large $`l`$ than the Schulz one and therefore possesses a significantly larger contribution of long rods. The log-normal distribution reads

$$p^{(0)}(l)=Nl^1\mathrm{exp}\left[\left(\frac{\mathrm{ln}l\mu }{2w^2}\right)^2\right].$$
(34)
For the unbounded log-normal distribution, $`w`$ is directly related to the polydispersity via $`w^2=\mathrm{ln}(1+\sigma ^2)`$ and the parameter $`\mu `$ is chosen such that $`m_1=1`$, giving $`\mu =w/2.`$ The second moment is then given by $`m_2=1+\sigma ^2`$. Also here, truncation of the distribution at some finite values $`l_{\text{min}}`$ and $`l_{\text{max}}`$ leads to deviations for which we correct numerically.

## IV Results for the onset of I-N phase separation

### IV.1 Schulz distributions

In Fig. 1 to 5 we have depicted the results for a Schulz parent distribution with cutoff lengths $`l_{\text{min}}=0.01`$ and $`l_{\text{max}}=100`$. The curves describing the densities of the cloud and shadow phases are shown in Figs. 1 and 2.

A striking broadening of the coexistence gap can be detected, mainly due to a dramatic increase of the concentration of the nematic cloud phase. In Fig. 2 we see that the volume fraction of the nematic cloud phase increases by several orders of magnitude at $`\sigma >0.4`$. Although the nematic shadow curve crosses the corresponding isotropic cloud curve at $`\sigma 0.5`$ in Fig. 1, the volume fraction (and hence the mass density) of the nematic shadow remains higher than that of the isotropic cloud phase throughout the phase diagram as we see in Fig. 2. Fig. 3 shows the extent of fractionation, i.e. the repartitioning of the long and short rods, among the coexisting phases at the onset of phase separation. A marked feature is the rapid increase of the average length in the nematic shadow around some ‘transitional’ polydispersity $`\sigma _t0.5`$. This indicates that the nematic phase becomes preferably populated by the longest rods in the system. Note that there is a similar effect in the isotropic shadow phase around $`\sigma _t0.7`$ where the shortest rods completely dominate the isotropic shadow phase at higher polydispersities. The same effects are reflected somewhat clearer in Fig. 4 showing the evolution of the polydispersity of the shadow phases. At $`\sigma =\sigma _t,`$ the polydispersity of the shadow phases show a kink. The strong decrease at higher $`\sigma `$ is due to the effect that the shadow phases become more and more enriched in either the longest or the shortest rods in the distribution. The dramatic change of the composition across the $`\sigma `$-range is shown explicitly in Fig. 5 where we have depicted the normalized length distributions in the nematic shadow phase for various $`\sigma `$. A similar picture is obtained for the distributions in the isotropic shadow phase (not shown here) but with the peak of the distribution shifting rapidly towards the lower cutoff length $`l_{\text{min}}=0.01`$.

In summary, we can state that there are two fractionation regimes for the onset of phase separation. First, at low polydispersities ($`\sigma <\sigma _t`$) moderate fractionation is observed and the shadow phases are mainly populated by rods with slightly higher (or lower) than average length. Second, at higher polydispersity ($`\sigma >\sigma _t`$) strong fractionation occurs such that the shadow phases are completely dominated by the longest (or shortest) rods in the distribution. In a small interval around $`\sigma =\sigma _t`$ the location of the peak of the length distribution shifts rapidly, upon increasing $`\sigma ,`$ from a value slightly different from one (pertaining to the low-$`\sigma `$ regime) to a value close to the the cutoff length (corresponding to the high-$`\sigma `$ regime).

So far, we have only discussed the results for a single set of cutoff lengths. Although the results for different cutoff lengths can be significantly different, particularly in the ‘cutoff-dependent’ regime sollichonsagerP2fattail  $`\sigma >\sigma _t`$, the global phase behavior remains qualitatively the same. Therefore we conclude that the aforementioned fractionation scenario generally holds for any Schulz parent with sufficiently extreme cutoff lengths ( $`l_{\mathrm{min}}1`$ and $`l_{\text{max}}1`$). An interesting limiting case however is the behavior for Schulz parents with infinitely large cutoff lengths. For this specific case, we could obtain simple scaling relations which describe the global behavior of the nematic shadow in the limit of an unbounded Schulz parent, i.e. $`l_{\text{max}}\mathrm{}`$. The scaling analysis, worked out briefly in Appendix B, is closely related to a more elaborate analysis presented in Ref. sollichonsagerexact  for the exact Onsager model. In particular, we show in Appendix B that the Gaussian Ansatz must yield the exact high-cutoff scaling relations. The reason for this is that our high-cutoff scaling form for $`p^{(N)}(l)`$ (for the nematic shadow phase) is completely analogous to the exact scaling form.

Finally, we remark that we do not observe a real jump in the shadow curves (and a kink in the associated cloud curves), as found in Ref. sollichonsagerexact . The presence of a jump in the shadow curve indicates that, at some point, the isotropic cloud phase coexists with two different nematic shadow phases and that a region of stable triphasic equilibria developes. Therefore we conclude that, within the Gaussian Ansatz, the Schulz form does not give rise to a three-phase separation, at least up to $`l_{\text{max}}=100`$.

### IV.2 Log-normal distributions

The results for the log-normal case are presented in Fig. 6 to 8. The cloud and shadow curves shown in these figures correspond to log-normal parent distributions with the same lower cutoff length $`l_{\text{min}}=0.01`$ but two different higher cutoffs $`l_{\text{max}}=10`$ and $`100`$. We see that the phase behavior is globally the same as for the Schulz case; there is a generic broadening of the biphasic region (Fig. 7) and a very pronounced fractionation effect, particularly for the isotropic cloud and nematic shadow phases, as visible in Fig. 8. At low polydispersities the distribution in the nematic shadow phase is very similar to the parent distribution (albeit with a higher average length). At higher polydispersities, however, we enter a regime characterized by extreme fractionation, i.e. the nematic shadow phase is dominated by the longest rods in the distribution. Remarkably, we do not see a similar transition in the isotropic shadow phase in Fig. 8, as we did in the Schulz case. This implies that the fractionation from the nematic cloud phase is much weaker for log-normal distributions than for Schulz ones.

A crucial difference with the previous results is that the transition between the regimes occurs discontinuously, that is, at the transition polydispersities $`\sigma _t`$ the isotropic cloud curves show a kink and the associated nematic shadow curves exhibit a jump. Precisely at the kink, the isotropic cloud phase coexists with two different nematic phases, one containing mostly rods with slightly higher than average length (denoted by $`N_I`$) while the second one ($`N_{II}`$) is dominated by the longest rods in the distribution. Therefore, this special point marks the beginning of a stability region for $`IN_IN_{II}`$ triphasic equilibria for log-normal distributions. In Fig. 8 we see that the position of the kink (in terms of $`\sigma _t`$) rapidly shifts to lower polydispersities as $`l_{\text{max}}`$ increases. From this, we anticipate that the triphasic equilibrium sets in at almost zero polydispersity (i.e. near monodisperse systems) for very large cutoff lengths.

Like for the Schulz case, we can obtain information about the global phase behavior for parent distributions at infinite cutoff lengths $`l_{\text{max}}`$ from the high-cutoff scaling results, shown in Appendix B. The most important outcome is that the concentrations of the isotropic cloud and nematic shadow phases go to zero for large cutoff length rather than approaching asymptotic forms such as in the Schulz case. Furthermore, it is shown explicitly that the fractionation between the isotropic cloud and nematic shadow phases is stronger than for the Schulz case.

So far, we have only looked at the onset of phase equilibrium by analyzing the properties of the cloud and shadow phases. The next step is to explore the coexistence region in more detail. An intriguing issue is to verify the region of stability for the isotropic-nematic-nematic triphasic region for the log-normal case. This will be dealt with in the next section, where we discuss the phase diagram for a log-normal parent with cutoff length $`l_{\text{max}}=10`$ in more detail.

## V Inside the I-N coexistence region

Across the coexistence region the equilibium length distributions of the coexisting phases, which originate from a parent phase with a prescribed distribution $`p^{(0)}(l)`$), change continuously as the overall density of the parent is $`c_0^{(0)}`$ is changed. In the actual calculations however it is more convenient to impose the fraction $`\gamma `$ occupied by the nematic phase rather than $`c_0^{(0)}`$ and calculate the corresponding densities self-consistently. In Fig. 9 to 12 we show the properties of the coexisting isotropic and nematic phases for a log-normal parent distribution with polydispersity $`\sigma =0.4.`$ Fig. 9 shows that the average length in the nematic daughter phase decreases rapidly in the regime $`c_0^{(0)}3`$ whereas only weak changes are notable at higher $`c_0^{(0)}`$. Furthermore, we see that the volume occupied by the nematic phase is extremely small in this regime. The same feature is observed in the volume fraction representation in Fig. 10.In particular, the non-monotonicity of the nematic branch is reflected somewhat clearer here. The fact that the isotropic branch runs extremely close to the dilution line for $`c_0^{(0)}<3`$ indicates that the fraction of nematic phase formed must indeed be very small. The rather exotic oscillations in the behavior of the polydispersities of the daughter phases in Fig. 11(a) reflect the dramatic change of the shape of the length distribution in the nematic phase in the first part ($`\gamma <10^2`$) of the dilution trajectory as shown in Fig. 11(b). Note that the distribution of the nematic shadow phase is in fact bimodal, with a small peak around $`l=1`$ and a much larger one at $`l=l_{\text{max}}`$. When the overall density is increased the second peak shifts to lower values of $`l`$ and eventually coincides with the first peak. When the overall density has reached about $`c_0^{(0}=4`$ (corresponding to a nematic phase volume fraction of about 10 %) the distribution of the nematic phase resembles the parental one, albeit with a slightly higher average length.

In Fig. 12 we have plotted the evolution of the average length for a parent with polydispersity $`\sigma =0.3`$. A peculiar behavior is observed, which is also reflected in Fig. 13 where the coexistence pressure is plotted versus $`\gamma `$. Clearly, there is a region where the pressure decreases as function of $`\gamma `$ which suggests an instability (or van der Waals) loop indicating a possible destabilization of the nematic phase. In Appendix C we explicitly show that the local extrema in the osmotic pressure in Fig. 13 represent spinodal points for the nematic phase which indicate that the coexisting nematic phase indeed becomes locally (and hence globally) unstable. In the region between the spinodal points (where $`\delta (b\beta \mathrm{\Pi })/\delta \gamma `$ is negative), the coexistence between the isotropic and a single nematic phase also becomes unstable such that a triphasic isotropic-nematic-nematic ($`IN_IN_{II}`$) demixing occurs.

It should be stressed that the actual onset of the three-phase separation is marked by binodal points which we have not located in this study, except for the kink at $`\sigma _t`$. In general, binodal points are located at a lower concentration than the spinodal points so that the demixing usually occurs well before the point where the system becomes locally unstable. This becomes clear in Fig. 14 where we show the details of the phase diagram in the vicinity of the kink including the spinodal curves (in terms of ($`c_0^{(0)}`$)) for the coexisting nematic phase. At the kink $`\sigma _t=0.264`$ the three-phase separation sets in right at the isotropic cloud point but the spinodal points are located at higher concentrations. An important feature in Fig. 14 is the presence of a high-$`\sigma `$ consolute (or critical) point at $`\sigma =0.373\pm 0.001`$ where the spinodal curves meet. This means that the region of stable triphasic equilibria does not extend up to large parent polydispersities but closes off at the consolute point which also constitutes an endpoint of the binodal curves corresponding to the triphasic equilibria.

Consequently, we can distinguish three regimes in the phase diagram depicted in Fig. 14. First, below the kink there is a common isotropic-nematic phase separation involving a moderately fractionated $`N_I`$-phase. Second, in the region between the kink and the consolute point $`0.264\sigma 0.373`$ a strongly fractionated $`N_{II}`$-phase (containing the longest rods) splits off initially at the isotropic cloud point. At higher concentrations (when $`I`$ and $`N_{II}`$ coexist in finite amounts) a second nematic $`N_I`$-phase is formed and an $`IN_IN_{II}`$ triphasic equilibium develops. Upon slightly further concentrating the sample the $`N_{II}`$-phase disappears and a regular $`IN_I`$ biphasic equilibrium is recovered. The third region is located at parent polydispersities above the consolute point $`\sigma >0.373`$. Upon concentrating the isotropic phase a strongly fractionated nematic phase (reminiscent of the $`N_{II}`$ phase) is formed initially but the composition of this phase evolves gradually towards a $`N_I`$-type nematic phase as the biphasic region is crossed.

## VI Discussion and Conclusions

We have numerically investigated Onsager’s second virial theory for polydisperse hard rods within the Gaussian Ansatz. The onset of isotropic-nematic phase separation is obtained from the cloud and shadow curves, which delimit the coexistence region. Within the same numerical framework, we could also explore the properties of the coexisting phases across the coexistence region. In this paper, we focussed on systems of polydisperse hard rods whose lengths can be described by a Schulz or a log-normal distribution. The basic difference between these two forms is that the fat-tailed log-normal one contains a significantly higher fraction of longer rods. For numerical and consistency purposes we truncated the distributions at both ends at sufficiently low and high cutoff lengths. Using truncated distributions is also justifiable from an experimental standpoint. For parent distributions of the Schulz type the phase diagram contains two fractionation regimes. First, at low parent polydispersities moderate fractionation occurs and the average rod lengths in the isotropic and nematic shadow phases are not much different from the average length in the parental cloud phase. Second, at higher parental $`\sigma `$ the fractionation effect is extremely strong and the isotropic and nematic shadow phases are completely dominated by respectively the longest and the shortest rods present in the system.

For the exact Onsager model, Speranza and Sollich sollichonsagerexact  very recently predicted a kink in the isotropic cloud curve (and a jump in the corresponding nematic shadow curve) for Schulz parents with cutoff lengths $`l_{\text{max}}`$ $`>50`$. The presence of this kink indicates a region of stable isotropic-nematic-nematic triphasic equilibria. Here we do not find any indication for such a three-phase separation for Schulz parents at least up to $`l_{\text{max}}=100`$. The discrepancy may be due to the Gaussian Ansatz, which implies that the ODFs are not represented by their correct equilibrium forms. Moreover, the Schulz form might be a borderline case since its tail is too modest to induce a strong demixing but too ‘fat’ to suppress it completely so that the presence of a kink in the isotropic cloud curve depends quite sensitively on the precise representation of the ODF.

Although the Gaussian ODF is not a solution of the exact stationarity condition for the ODF, it does satisfy the exact high-density scaling relation vanroijmulderscaling . This means that the properties of highly ordered nematic states are described very well by the Gaussian form. In fact, the description becomes exact for infinitely aligned states. A manifestation of this is the osmotic pressure for the nematic phase, Eq. (20), which is the exact high-density result vanroijmulderscaling . Consequently, for our polydisperse systems, we expect the Gaussian Ansatz to work increasingly well both for highly concentrated nematic phases and nematics that are dominated by the longest rods. In both cases, the nematic alignment of all species is expected to be very pronounced such that the use of the scaling ODF (for all $`l`$) is justified. To verify this notion, we have plotted the variational parameter as a function of length, for both the nematic shadow phase and the nematic cloud phase corresponding to a Schulz parent in Fig. 15. Since the Gaussian ODF is expected to be the least correct for the shortest rods (which show the weakest alignment), we focus on the interval $`l_{\text{max}}<1`$. In order for the results to be self-consistent, the alignment must be strong enough and hence the variational parameter must be sufficiently large (say $`{}_{}{}^{10}\mathrm{log}\alpha >1`$) for all lengths. Fig. 15 (a) shows that this is not entirely the case; in the regime of low fractionation the shortest rods (with $`l0.4`$) are not sufficiently aligned by the longer rods so that the Gaussian description fails here. In the regime $`\sigma >0.5`$, which is physically the most interesting one, the ordering of the short rods is much higher due to the presence of very long rods in the nematic shadow phase, and the Gaussian Ansatz is fully justified. Similarly for the nematic cloud phase in Fig. 15 (b) we see that the shortest rods are not well represented by the Gaussian ODF at low $`\sigma `$ but much better at $`\sigma >\sigma _t`$ where the variational parameter increases several orders of magnitude due to a dramatic increase of the concentration of the nematic cloud phase (see also Fig. 2). Obviously, for any rod length significantly larger than $`l_{\text{min }}`$ the Gaussian ODF works very well because $`\alpha `$ generally becomes extremely large for any $`\sigma `$. Therefore we conclude that given the fact that the composition of the nematic phases is dominated by the longest rods, particularly in the physically relevant cutoff-dependent regime, the Gaussian description is an appropriate tool in our study.

We now turn to the log-normal case. The fractionation scenario we observe there is qualitatively the same as the one for the Schulz case; weak fractionation occurs at low $`\sigma `$ where the distributions in the shadow phases are reminiscent of the parental one, but dramatic segregation effects take place above some threshold $`\sigma ,`$ particularly between the isotropic cloud and the nematic shadow phase. A crucial difference with the results for the Schulz case however is that the transition between the two fractionation regions shows a discontinuity at some $`\sigma =\sigma _t`$. At this point, the isotropic cloud curve shows a kink which corresponds to a jump in the nematic shadow phase. The jump indicates that an isotropic cloud phase must coexist with two different shadow phases, one containing mostly rods of slightly higher than average length (the $`N_I`$-phase) and the other one predominantly containing very long rods (the $`N_{II}`$-phase). The kink also marks the beginning of a region of stable $`IN_IN_{II}`$ triphasic equilibria. For log-normal distributions with a moderate cutoff-length $`l_{\text{max}}=10`$ we found that the triphasic region opens up at a fairly low polydispersity $`\sigma _t=0.264.`$ This value will decrease for larger cutoff lengths and eventually go to zero when $`l_{\text{max}}`$ approaches infinity. This indicates that adding a very small fraction of long rods to a weakly polydisperse system of much shorter rods can already induce a three-phase demixing. A similar effect is observed in binary mixtures of long and short hard rods with sufficiently large length ratios LekkerVroeg .

By numerically analyzing the spinodal instability criteria for the nematic phase across the coexistence region for the case $`l_{\text{max}}=10`$ we have shown that the triphasic area does not extend up to very large parent polydispersities but terminates at a consolute point located at $`\sigma =0.373\pm 0.001.`$ This means that isotropic-nematic-nematic triphasic equilibria are only expected to occur in a small interval of parent polydispersities, namely $`0.264\sigma 0.373`$. However, it should be mentioned that we have not been able to determine the corresponding binodal points, which mark the actual onset of three-phase separation from the isotropic-nematic biphasic equilibria. Considering the results for the approximate Onsager $`𝒫_2`$-model, as numerically analyzed by Speranza and Sollich sollichonsagerP2fattail , we expect these triphasic equilibria to be limited to a very small density interval across the coexistence region, which makes it very hard to observe the phenomenon in experiment. Another problem is that the fraction of the system volume occupied by the nematic phases is predicted to be at the most 0.1 %, so that it will be very difficult to distinguish (or even detect) the two different nematics. Therefore we must conclude that, although the log-normal distribution contains sufficiently long rods to induce a demixing of the nematic phase, the (mole) fraction of these rods is too small to give rise to an observable fraction of the demixed nematic phases.

Experimentally, phenomena such as a broadening of the biphasic region and a fractionation effect have been observed unequivocally in a number of experimental studies fradenmaret ; donkai93 ; Buining93 ; Itou . Observations of triphasic isotropic-nematic-nematic equilibria were however only reported for systems whose length distributions appear to be more or less bimodal rather than unimodal. These bimodal shapes were either accomplished deliberately by mixing species with different lengths, as done by Itou Itou  with semiflexible schyzophyllan rods or caused by the presence of large aggregates as found by Kajiwara et al. Kajiwara  in systems of rigid imogolite rods. Buining et al. Buining93  observed the formation of second nematic phase in systems of polymer-coated (hard) boehmite rods, albeit a long time after the two-phase isotropic-nematic phase separation had finished. Also there, the triphasic demixing is probably due to the presence of a small number of very long rods or aggregates present in the system causing the length distribution to be (slightly) bimodal Buiningbeharing .

The experimental results therefore seem to indicate that having a three-phase separation with observable fractions of all phases requires some degree of bimodality (with a sufficiently large length ratio between the short and long species) in the parental length distribution. Hence, an intriguing issue left open for future investigation is the question how the triphasic phenomenon predicted for the log-normal distributions change for a parent distribution with a slightly bimodal shape (e.g. with a small second peak just below $`l_{\mathrm{max}}`$). In particular, a bimodal distribution may give rise to enhanced fractionation behavior and more pronounced triphasic equilibria than predicted for the unimodal log-normal form.

## Appendix A: Numerical procedure

The self-consistency equations for $`\stackrel{~}{\alpha }(l)`$ and $`\mu _{\text{ex}}^{(N)}(l),`$ Eqs. (8) and (14), were solved using a numerical grid of lengths. The iterative scheme we used is analogous to the one described by Herzfeld et al. herzfeldgrid  for computing the numerically exact equilibrium ODF of the monodisperse Onsager model. The $`l`$interval $`[l_{\text{min}},l_{\text{max}}]`$ was discretized into $`N`$ (not necessarily equal) parts. The mesh size must be chosen very carefully, particularly for large $`l_{\mathrm{max}}`$, because the distributions become considerably peaked at low polydispersities. Therefore, for parent polydispersities lower than approximately $`\sigma =0.25`$ we chose to divide the integration interval into three regimes. The interval $`[18\sigma ^2,1+8\sigma ^2]`$ in the vicinity of the peak, where the distribution changes rapidly, was discretized into $`\frac{3}{5}N`$ equal parts, and the intervals $`[l_{\mathrm{min}},18\sigma ^2]`$ and $`[1+8\sigma ^2,l_{\text{max}}],`$ where the distribution is generally much smoother, were both discretized into $`\frac{1}{5}N`$ equal parts. For parent polydispersities larger than $`\sigma =0.25`$ the entire interval was discretized into $`N`$ equal parts. It has been proved sufficient to use $`N=150`$ in order to obtain quantitatively reliable results. However, for the calulation of the full phase split (Sec. III-C) a smaller mesh size ($`N=50)`$ was used to limit the computational burden. Increasing the number of mesh points leads to only marginally different results in this case while the calculation time increased dramatically. For small polydispersities we took Eqs. (10) and (15) as initial guesses and the successive iteration was performed until the following convergence criteria were satisfied

$`\underset{n=1,\mathrm{}N}{\mathrm{max}}\left|\stackrel{~}{\alpha }(l_{n+1})\stackrel{~}{\alpha }(l_n)\right|`$ $`<`$ $`10^6,`$
$`\underset{n=1,\mathrm{}N}{\mathrm{max}}\left|\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l_{n+1})\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l_n)\right|`$ $`<`$ $`10^6.`$ (35)
After each iteration step 10% of the new solution had to be added to 90 % of the previous one for the next iteration step to ensure the convergence of the method.

The iteration algorithm we figured out to calculate the phase equilibria, in particular for the isotropic cloud and shadow phases and the full phase split, can be described as follows. First, the corresponding equilibrium forms for $`\stackrel{~}{\alpha }(l)`$ and $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ were calculated for a given set of starting concentrations. These results were then put into the self-consistency equations for the cloud and shadow concentrations to obtain new values. These equations were obtained by recasting the self-consistency conditions for the moments (e.g. Eq. (27)) in an iterative form such that the concentrations could be calculated by a simple fixed-point iteration.

Finally, for the new concentrations, corresponding forms for $`\stackrel{~}{\alpha }(l)`$ and $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ were computed and substituted again into the self-consistency equations. This procedure was repeated until all concentrations had converged to within 10<sup>-5</sup>. To ensure convergence of this iteration loop a damping percentage of 80 % was used, which means that only 20 % of the change was retained at each iteration step.

## Appendix B: High-cutoff scaling results

In this appendix we focus on the properties of the nematic shadow phase in the cutoff-dependent regime ($`\sigma >\sigma _t`$) for systems with infinite cutoff lengths. For the exact Onsager model, Speranza and Sollich sollichonsagerexact  made a detailed analysis of these properties based upon the high-density (and high-cutoff) scaling forms of the exact nematic ODF. Here, we will not reproduce the analysis but merely show that the scaling form for the length distribution in the nematic shadow phase in the limit $`l_{\text{max}}\mathrm{}`$ is analogous to the one obtained in Ref. sollichonsagerexact . Consequently, all scaling properties which follow from the Gaussian approximation must be exactly the same as the ones derived from the exact high-cutoff scaling results.

The first step is to solve the coupled set of consistency equations, Eqs. (8) and (14). In order to obtain analytic solutions for these nonlinear integral equations we exploit the fact that the nematic shadow phase is completely dominated by the longest rods in the system at $`\sigma >\sigma _t`$. When the cutoff length increases the length distribution in the nematic shadow will be more and more peaked at $`l=l_{\text{max}}`$. In the limit of infinite $`l_{\text{max}}`$ it is therefore justified to use the Ansatz $`p^{(N)}(l)=\delta (ll_{\text{max}})`$ which suggests an effectively monodisperse nematic shadow phase only containing the longest rods in the system. Substituting the delta-function in Eqs. (8) and (14) allows us to obtain asymptotic forms for the Gaussian variational parameter $`\stackrel{~}{\alpha }(l)`$ and the excess chemical potential $`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ of the nematic shadow phase. These expressions now read

$`\stackrel{~}{\alpha }(l)`$ $`=`$ $`l_{\text{max}}^4(l/l_{\text{max}}),`$ (36)
$`\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l)`$ $`=`$ $`{\displaystyle \frac{l}{l_{\text{max}}}}2^{3/2}\sqrt{1+^1(l/l_{\text{max}})},`$ (37)
where $``$ is given by

$$(l/l_{\text{max}})=\frac{1}{2}\left[\sqrt{1+8\left(\frac{l}{l_{\text{max}}}\right)^2}1\right].$$
(38)
Note that $`(l/l_{\text{max}})`$ always has a value between zero and unity. A close inspection of Eqs. (36) and (37) reveals that the variational parameter $`\stackrel{~}{\alpha }(l)`$ scales as $`\stackrel{~}{\alpha }(l)l_{\text{max}}^4`$ whereas the excess chemical potential $`\stackrel{~}{\mu }_{\text{ex}}^N(l)`$ remains of the order $`𝒪(1)`$ for all lengths. Using this in Eq. (26) we can write down a scaling expression for the length distribution in the nematic shadow, which in its general form reads

$$p^{(N)}(l)=\text{cst}\frac{c_0^{(0)}}{(c_0^{(N)})^3}l_{\text{max}}^4p^{(0)}(l)\mathrm{exp}\left[2c_0^{(0)}l𝒲(l/l_{\text{max}})\right],$$
(39)
where $`𝒲`$ is a contribution of the order $`𝒪(1)`$:

$$𝒲(l/l_{\text{max}})=\mathrm{ln}(l/l_{\text{max}})+\stackrel{~}{\mu }_{\text{ex}}^{(N)}(l/l_{\text{max}}).$$
(40)
Note that $`𝒲`$ attains its maximum $`𝒲=4`$ for $`l=l_{\text{max}}`$. The scaling solution Eq. (39) is completely analogous to the one found in Ref. sollichonsagerexact , the only differences being the exact form of $`𝒲(l/l_{\text{max}})`$ and the constant cst. However, since these contributions are both subleading in the limit $`l_{\text{max}}\mathrm{}`$ they are irrelevant for the rest of the analysis and hence do not influence the scaling results. The similarity between the exact high density scaling analysis and the Gaussian approximation is also confirmed by the fact that both theories predict the same high-density scaling result for the nematic osmotic pressure, namely $`b\beta \mathrm{\Pi }=3c_0^{(N)}`$.

For the sake of completeness, let us now briefly outline the basic results of the large cutoff scaling analysis. For a comprehensive treatment of this subject the reader is referred to Ref. sollichonsagerexact . For a Schulz parent, we may use Eq. (32) in Eq. (39) to obtain

$$p^{(N)}(l)=𝒦_Nl_{\text{max}}^4l^z\mathrm{exp}\left[\epsilon l𝒲(l/l_{\text{max}})\right],$$
(41)
where $`\epsilon =2c_0^{(0)}(z+1)`$ and $`z=\sigma ^21.`$ For very large $`l`$ the exponent $`\mathrm{exp}[\epsilon l]`$ will be the dominating contribution. At $`\sigma >\sigma _t,`$ the nematic shadow is supposed to be dominated by the longest rods and the distribution $`p^{(N)}(l)`$ should therefore be an increasing function of length. This requires $`\epsilon `$ to be positive and yields the condition $`c_0^{(0)}>\frac{1}{2}(z+1)`$. Since the concentration of the cloud phase appears to decrease with increasing $`l_{\text{max}}`$ this then implies that the isotropic cloud curve (and hence the nematic shadow curve) has a finite lower bound for large cutoff lengths. These limiting solutions, for which $`\epsilon =0`$, therefore read

$`c_0^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2\sigma ^2}},`$
$`c_0^{(N)}`$ $`=`$ $`{\displaystyle \frac{1}{6\sigma ^2}}\left(1+{\displaystyle \frac{1}{2\sigma ^2}}\right),`$ (42)
using Eq.(28), with $`c_1^{(0)}=c_0^{(0)}`$. These results are plotted in Fig. 1. To be consistent, let us now look for a solution for the transition polydispersity $`\sigma _t`$ above which the nematic shadow phase for a Schulz distributed parent is completely dominated by the longest rods. We start with the concentration of the nematic shadow phase which is proportional to the integral over the normalized length distribution, i.e $`c_0^{(N)}p^{(N)}(l)𝑑l`$. From Eq. (41) we thus obtain

$$c_0^{(N)}l_{\text{max}}^4_0^{l_{\text{max}}}l^z\mathrm{exp}\left[\epsilon l𝒲(l/l_{\text{max}})\right]𝑑l.$$
(43)
For the sake of convenience we have set $`l_{\text{min}}`$ equal to zero. Since the integrand is dominated by the exponent $`\mathrm{exp}[\epsilon l]`$ for large $`l`$ we may approximate the integral by bringing all slowly varying contributions in front of the integral sign and evaluating them at $`l=l_{\text{max}}`$, which gives

$`c_0^{(N)}`$ $``$ $`l_{\text{max}}^4\left|l^z\mathrm{exp}[𝒲(l/l_{\text{max}})]\right|_{l=l_{\text{max}}}{\displaystyle _0^{l_{\text{max}}}}𝑑l\mathrm{exp}[\epsilon l]`$
$`c_0^{(N)}`$ $``$ $`l_{\text{max}}^{z3}{\displaystyle \frac{\mathrm{exp}[\epsilon l_{\text{max}}]}{\epsilon l_{\text{max}}}}.`$ (44)
The next step is to recast the latter equation into a scaling relation for $`\epsilon .`$ Taking the logarithm on both sides of Eq. (44) gives

$$\epsilon (3z)\frac{\mathrm{ln}l_{\text{max}}}{l_{\text{max}}}+\frac{\mathrm{ln}\epsilon l_{\text{max}}}{l_{\text{max}}}+𝒪(l_{\text{max}}^1).$$
(45)
From the known limits $`\epsilon l_{\text{max}}\mathrm{}`$ and $`\epsilon 0`$ for $`l_{\text{max}}\mathrm{}`$ we can deduce that $`\epsilon l_{\text{max}}`$ must increase more slowly than linearly with $`l_{\text{max}}`$. Consequently, the second and third terms in Eq. (45) are both subleading contributions so that we retain up to leading order

$$\epsilon (3z)\frac{\mathrm{ln}l_{\text{max}}}{l_{\text{max}}},$$
(46)
which shows that $`\epsilon l_{\text{max}}`$ indeed increases logarithmically rather than linearly as we already anticipated. However, in order to make this result fully self-consistent it is also required that $`z`$ $`<3`$ (and correspondingly $`\sigma >0.5`$) since $`\epsilon `$ must be positive. This means that $`\sigma =0.5`$ is a lower bound for the cutoff-dependent regime in the limit $`l_{\text{max}}\mathrm{}`$. In other words, the transition from the low fractionation regime to the regime where the nematic shadow is completely dominated by the longest rods occurs exactly at $`\sigma _t=0.5`$ for Schulz parents with infinitely high cutoff lengths.

We now turn to the average length in the nematic shadow phase which is related to the first moment density $`c_1^{(N)}lp^{(N)}(l)𝑑l`$. Analogously to Eq. (44) it follows that $`c_1^{(N)}`$ $`l_{\text{max}}^{z2}\mathrm{exp}[\epsilon l_{\text{max}}]/\epsilon l_{\text{max}}`$ and that the average length hence scales as $`l`$ $`c_1^{(N)}/c_0^{(N)}l_{\text{max}}`$. Since the distribution in the nematic shadow phase is dominated by $`\mathrm{exp}[\epsilon l]`$ we expect that only rods whose lengths are of the order $`𝒪(1/\epsilon )`$ smaller than $`l_{\text{max}}`$ contribute to the average length. We therefore can write

$$l=l_{\text{max}}\left[1𝒪\left(\frac{1}{\epsilon l_{\text{max}}}\right)\right],$$
(47)
and from Eq. (46)

$$l=l_{\text{max}}\left[1𝒪\left(\frac{1}{\mathrm{ln}l_{\text{max}}}\right)\right].$$
(48)
This result shows that the average length in the nematic shadow in principle diverges for $`l_{\text{max}}\mathrm{}`$ but the logarithmic correction causes the actual $`l`$ to be significantly lower than $`l_{\text{max}}`$.

A similar treatment can be given for a log-normal parent distribution. However, the analysis for the log-normal case is even more involved and we will only present the basic results and refer to Ref. sollichonsagerexact  for details. First, the concentration of the isotropic cloud phase appears to have the following $`l_{\text{max}}`$-dependence

$$c_0^{(0)}=\frac{\mathrm{ln}^2l_{\text{max}}}{4\mathrm{ln}(1+\sigma ^2)l_{\text{max}}}+𝒪\left(\frac{\mathrm{ln}l_{\text{max}}}{l_{\text{max}}}\right),$$
(49)
which is crucially different from the Schulz case because the concentration of the cloud and shadow phases now tend to zero rather than approaching boundary values as in the Schulz case. Second, the average length can be shown to behave as

$$l=l_{\text{max}}\left[1𝒪\left(\frac{1}{\mathrm{ln}^2l_{\text{max}}}\right)\right].$$
(50)
Similarly to the Schulz distribution, the average length scales as $`ll_{\text{max}}`$ but the correction term is now considerably smaller which implies that the fractionation effect is much more pronounced for log-normal distributions at $`\sigma >\sigma _t`$, as we already noticed by comparing Figs. 3 and 8.

## Appendix C: Polydisperse spinodal instability criterion for the nematic phase

In this appendix we show that the anomalous behavior of the coexistence pressure in Fig. (13) can be related to the local instability of the coexisting nematic phase. In particular, it is established that the (local) extrema in the coexistence pressure are directly related to the spinodal points where the nematic phase becomes locally (and hence globally) unstable. Let us first consider the Gibbs-Duhem equation, which in its general form is given by

$$\left(\frac{P}{\rho _j}\right)_{T,\rho _{ij}}=\underset{i}{}\rho _i\left(\frac{\mu _i}{\rho _j}\right)_{T,\rho _{ij}},$$
(51)
with $`\rho _i=N_i/V`$. For a continuous density distribution $`c(l)`$ this relation can be expressed in terms of functional derivatives

$$\frac{\delta (b\beta \mathrm{\Pi })}{\delta c(l^{})}=𝑑lc(l)\frac{\delta \beta \mu (l)}{\delta c(l^{})}.$$
(52)
We now focus on the evolution of the length distribution of the nematic phase (in coexistence with an isotropic phase) across the coexistence region and denote this distribution by $`c^{}(l)`$. The shape of $`c^{}(l)`$ depends uniquely on the nematic fraction $`\gamma `$, so that the curves in Figs. 12 and 13 represent trajectories parametrized by $`\gamma `$. As to the osmotic pressure, let us now consider an infinitesimal change of the coexistence pressure $`\delta (b\beta \mathrm{\Pi })^{}`$ corresponding to an infinitesimally small dispacement $`\delta c^{}(l)`$ on the trajectory in Fig. (13). Using the Gibbs-Duhem relation Eq. (52), it then follows that

$`\delta (b\beta \mathrm{\Pi })^{}`$ $`=`$ $`{\displaystyle 𝑑l^{}\frac{\delta (b\beta \mathrm{\Pi })}{\delta c(l^{})}}|_{c^{}(l^{})}\delta c^{}(l^{})`$ (53)
$`=`$ $`{\displaystyle 𝑑lc^{}(l)𝑑l^{}\frac{\delta \beta \mu (l)}{\delta c(l^{})}}|_{c^{}(l^{})}\delta c^{}(l^{}).`$
Since $`c^{}(l)>0`$ for all $`l`$ the infinitesimal pressure change is only zero when

$$𝑑l^{}\frac{\delta ^2f}{\delta c(l)\delta c(l^{})}|_{c^{}(l)}\delta c^{}(l)=0,$$
(54)
which is precisely the spinodal instability criterion for polydisperse systems sollichwarrenmomenttheory , albeit restricted to composition fluctuations $`\delta c^{}(l)`$ which lie on the trajectories parametrized by $`\gamma `$. This result therefore shows that the local extrema in the coexisting pressure where $`\delta (b\beta \mathrm{\Pi })^{}=0`$ represent spinodal points indicating local instability of the coexisting nematic phase. The instability direction $`\delta c^{}(l)`$ is defined as the change of the nematic length distribution $`c(l)`$ corresponding to an infinitesimal displacement $`\delta \gamma `$ on the trajectory in Fig. 13. The presence of the spinodal points for the nematic phase imply the existence of a isotropic-nematic-nematic triphasic region.