text:
# Simplified Onsager theory for isotropic-nematic phase equilibria of length polydisperse hard rods

## I Introduction

Suspensions of rod-like particles can undergo an orientational disorder-order phase transition from a disordered isotropic (I) phase, in which rods point with equal probability in every direction, to an orientationally ordered nematic (N) phase in which particle orientations cluster around a preferred direction. This behaviour has been observed experimentally in both chemical  and biological  systems.

Two main theoretical approaches have been used to explain the I-N phase transition. The Maier-Saupe theory  neglects density variations and is therefore appropriate for so-called thermotropic materials, where the transition is driven by variations in temperature. It focuses on attractive, long-range interactions between particles; in the original theory these attractions were thought of as arising from Van der Waals forces but in fact the Maier-Saupe approach can be considered a lowest order approximation for general orientation-dependent attractive interactions . The Onsager theory , on the other hand, takes into account only hard core repulsions between the particles on contact; temperature then becomes unimportant and so the theory is appropriate for lyotropic materials, where the phase ordering is driven by changes in density. The I-N phase transition arises from the competition between the orientational entropy (corresponding to the tendency of rods to stay orientationally disordered) and the packing entropy (which is higher for aligned rods due to the excluded volume interaction). The key simplification is that, in the “Onsager limit” of infinitely thin rods, the virial expansion for the free energy truncates after the second-order contribution. The state of the system is described by the density of rods and the distribution $`P(\theta )`$ of the angles $`\theta `$ which the rods make with the nematic axis, i.e., the preferred orientation in the nematic phase. A formal minimization of the free energy with respect to this function leads to a self-consistency equation for $`P(\theta )`$. Solving this gives, in principle, the free energy as a function of density, from which the I-N phase coexistence region can then be obtained by a standard double-tangent construction. Onsager gave an approximate solution for the I-N phase transition using a variational ansatz for $`P(\theta )`$. With two different approximation techniques Lekkerkerker et al.  and Lasher  obtained more accurate results; the numerically exact solution has also been obtained . Isihara , basing his calculations essentially on the Onsager theory, had shown earlier that the order-disorder phase transition will occur at similar densities for particles of different shapes, as long as they are strongly asymmetric (e.g. plates, cylinders, spheroids).

We mention briefly an alternative to Onsager theory that was developed by Flory , Di Marzio  and Alben , based on a lattice model. The interaction is again of excluded volume type, but particles are now embedded in a three-dimensional cubic lattice whose sites are occupied either by segments of the rod-like particle or by solvent. Although qualitetively this model gives the same behaviour as the Onsager theory, the I-N phase transition is predicted at a density about twice that found by Onsager . With a similar lattice model Flory  predicted an analogous order-disorder phase transition for semiflexible macromolecules. On the experimental side, Robinson , Hermans  and Nakajima et al.  found results in good agreement with Flory’s theory for hard rod-like molecules.

Returning now to Onsager theory, a direct comparison between the theoretical predictions and experimental data is difficult due to the complexity of experimental systems. For example, there are corrections due to non-hard interactions in real systems  which affect the densities of coexisting isotropic and nematic phases and other characteristics of the phase transition . These effects can be partially accounted for by adding an attractive potential to the excluded volume interaction , which leads to a temperature dependence of e.g. the orientational order parameter.

More importantly for our purposes, polydispersity – the presence of an (effectively) continuous spread of rod lengths and/or diameters – has been shown to have important effects on the phase diagram . Onsager  already gave the extension of his approach to the case of mixtures of rods with different lengths or different diameters, but the obvious complication of having many (and in the fully polydisperse limit, infinitely many) distinguishable particle species prevented a direct numerical solution, or the systematic use of trial functions as in the monodisperse case. We focus in this paper on length polydispersity alone, assuming that rod diameters are still monodisperse. Such length polydispersity can lead to a strong broadening of the coexistence region, which has been observed experimentally  and also predicted theoretically, at least within an expansion for narrow length distributions . Coexistence of several nematic phases (N-N), occasionally with a third isotropic phase (I-N-N), has also been seen experimentally  and been predicted by theory for bi- and tridisperse systems (comprising rods of two and three different lengths, respectively) .

In a series of papers by Flory et al. , the length-polydisperse case was solved within the lattice model for a bidisperse system as well as for cases with many different rod lengths (which in the lattice model are always multiples of a given unit length). The results obtained are qualitatively similar to the ones predicted by the Onsager theory , with a three-phase I-N-N region and an N-N region at high density for the bidisperse system, and with strong fractionation and broadening of the coexistence region for length distributions of e.g. Poisson  and exponential  type. Moscicki et al. applied the same theory to a system with a Gaussian length distribution , observing strong fractionation and a pronounced dependence of the size of the coexistence region on the width of the distribution.

The theoretical studies referred to above which concern themselves directly with the off-lattice situation are restricted to either mixtures of rods of two or three different lengths or to continuous length distributions of small width. The question of the effects of full length polydispersity within the Onsager theory of thin hard rods therefore remains open, and it is this problem that we address in the present paper. It is clear that the introduction of full polydispersity significantly complicates the analysis of phase behaviour. Each rod length now has its own conserved density associated with it, and a naive application of the Gibbs phase rule suggests that one could have arbitrarily many coexisting phases. This problem is compounded by the presence of the varying rod orientations, so that the full density distribution over lengths and orientations has both conserved and non-conserved parts. A direct attack on the Onsager model with full length polydispersity is therefore extremely difficult.

As a first step towards understanding the phase behaviour of polydisperse hard rods, we have previously analysed the Zwanzig model, a discretization of the Onsager theory which allows rods to point only in one of three orthogonal directions . We found that this exhibited some of the features that are typical of polydisperse systems (broadening of the coexistence region, fractionation of rod lengths across different phases). It was not, however, able to predict the observed N-N and I-N-N phase separations. This could be traced back to the restriction on rod orientations: a separation of a single nematic phase into two nematics containing predominantly shorter and longer rods involves a loss of entropy of mixing and a gain in orientational entropy; with only three different rod orientations available, the possible gain in orientational entropy is too limited for such a phase separation to be favourable.

In the present work we analyse a model which can take the complexities of full length polydispersity into account and still predict the more complex phase coexistences involving several nematic phases that were missing in the Zwanzig model. Our “$`𝒫_2`$ Onsager model” is an approximation of the full Onsager theory obtained by truncating the expansion of the angular dependence of the excluded volume of two rods after the lowest nontrivial term, which involves the Legendre polynomial $`𝒫_2`$. The model is thus no longer exact, even in the limit of thin rods, but nevertheless useful as a solvable approximation to the full Onsager theory which can, for example, provide some information about the possible topologies of the phase diagram. The errors introduced by the truncation can be directly assessed for the bidisperse case, where accurate numerical solutions of the full Onsager theory exist , and we find below that there the $`𝒫_2`$ Onsager model reproduces most of the trends in the evolution of the phase diagram with the rod length distribution. The truncation errors could in principle be reduced by extending the angular expansion to increasingly high orders, giving a series of models whose behaviour will approach that of the full Onsager theory in the limit; the $`𝒫_2`$ Onsager model is the first member of this series. Its free energy also bears a formal similarity to that of the polydisperse Maier-Saupe model, and by imposing the constraint of constant density could yield qualitative information about the phase behaviour of polydisperse thermotropics. A final benefit is that the $`𝒫_2`$ Onsager model is, in contrast to the full Onsager theory, truncatable: its excess free energy only depends on a finite number (two) of moments of the density distribution. Phase equilibria can therefore be found efficiently using the moment free energy method  (or moment method for short), with its extension to systems with non-conserved densities introduced in the context of the polydisperse Zwanzig model .

This paper is structured as follows. In Sec. II we describe our truncated “$`𝒫_2`$ Onsager model” and its link to the full Onsager theory, and give the exact phase coexistence equations for the truncated model. In Sec. III, we briefly review the moment method and then apply it to the $`𝒫_2`$ Onsager model. The numerical method used for calculating phase equilibria is described in some detail, including our procedure for ensuring that the moment method produces numerically accurate predictions for phase equilibria even beyond the point where it is exact by construction (which turns out to be the onset of nematic phase ordering). Our main results are detailed in Sec. IV, where we give phase diagrams for three different types of rod length distributions: Schulz distributions, which are unimodal, bidisperse distributions containing only two different rod lengths, and finally mixtures of Schulz distributions that interpolate between these two extremes. One of our guiding questions will be under which conditions on the rod length distribution the phase diagram shows the more “exotic” features of nematic re-entrance and of I-N-N and N-N coexistence regions. We will find that, while these features are absent for the unimodal Schulz distributions with even the largest widths that we consider, the bidisperse and bimodal systems give a much richer phase behaviour. We summarize in Sec. V and discuss directions for future work. App. A contains a description of a simple approximation that we found useful in elucidating certain trends in the evolution of the phase diagram, particularly in the limit of large length ratios between long and short rods for the bimodal and bidisperse distributions.

## II Model details

We consider a system of hard cylinders of length $`L`$ and diameter $`D`$, capped with hemispheres at both ends. We assume that these rod-like particles are polydisperse in their length $`L`$, but all have identical diameter $`D`$. To be able to take the Onsager limit of long thin rods, we introduce a reference length $`L_0`$ and the normalized rod lengths $`l=L/L_0`$; the Onsager limit is then reached for $`D/L_00`$ at fixed values of $`l`$.

The state of a single phase of a system is described by the density distribution $`\rho (l,\mathrm{\Omega })`$, defined so that $`\rho (l,\mathrm{\Omega })dld\mathrm{\Omega }/4\pi `$ gives the density of rods with (normalized) lengths in an interval $`dl`$ around $`l`$, and orientations $`\mathrm{\Omega }`$ in a solid angle $`d\mathrm{\Omega }`$. The distribution $`\rho (l,\mathrm{\Omega })`$ is the natural extension for a polydisperse system of the usual orientational distribution $`P(\mathrm{\Omega })`$ used for monodisperse systems of rods . The rod orientation $`\mathrm{\Omega }`$ can be parameterized in terms of the angle $`\theta `$ with the nematic axis and an azimuthal angle $`\phi `$; due to the cylindrical symmetry of the nematic phase $`\rho (l,\mathrm{\Omega })`$ is independent of $`\phi `$. Using $`d\mathrm{\Omega }=d\mathrm{cos}\theta d\phi `$ the density distribution as a function of rod length, obtained by integrating over orientations, is therefore

$$\rho (l)=\frac{1}{4\pi }𝑑\mathrm{\Omega }\rho (l,\mathrm{\Omega })=\stackrel{~}{d\theta }\rho (l,\theta )$$
(1)
Here and below we use the shorthand

$$\stackrel{~}{d\theta }F(\theta )=\frac{1}{2}_1^1d\mathrm{cos}\theta F(\theta )$$
(2)
for the angular integral of an arbitrary function $`F(\theta )`$. The orientational distribution of rods can be factored out from $`\rho (l,\theta )`$ as

$$\rho (l,\theta )=\rho (l)P_l(\theta )$$
(3)
where $`P_l(\theta )`$ represents the probability of finding a rod of given length $`l`$ in orientation $`\mathrm{\Omega }=(\theta ,\phi )`$ and is normalized to $`1`$:

$$\stackrel{~}{d\theta }P_l(\theta )=1$$
(4)
In the isotropic phase, one has $`P_l(\theta )1`$ and $`\rho (l,\theta )=\rho (l)`$.

### A Polydisperse Onsager theory

As mentioned above, the present model is a truncation of the full Onsager theory (see e.g. Refs.  for reviews), which we now describe. In the Onsager limit $`D/L_00`$, the second order virial approximation becomes exact and gives for the excess free energy in units of $`k_\mathrm{B}T`$

$$\stackrel{~}{f}=𝑑l𝑑l^{}\frac{d\mathrm{\Omega }}{4\pi }\frac{d\mathrm{\Omega }^{}}{4\pi }B_2(lL_0,\mathrm{\Omega },l^{}L_0,\mathrm{\Omega }^{})\rho (l,\theta )\rho (l^{},\theta ^{})$$
(5)
Here $`B_2`$ is the second virial coefficient which for hard rods is half the excluded volume and given by

$$B_2=DL_0^2ll^{}|\mathrm{sin}\gamma |$$
(6)
up to negligible terms of order $`D/L_0`$, with $`\gamma `$ the angle between the two rods. Using the cylindrical symmetry of the nematic phase one can perform the integration over $`\phi `$ and $`\phi ^{}`$ in Eq. (5). If we also make the densities $`\stackrel{~}{f}`$ and $`\rho (l,\theta )`$ non-dimensional by multiplying by a unit volume $`(\pi /4)DL_0^2`$ we obtain

$$\stackrel{~}{f}=\frac{4}{\pi }𝑑l𝑑l^{}\stackrel{~}{d\theta }\stackrel{~}{d\theta ^{}}ll^{}K(\theta ,\theta ^{})\rho (l,\theta )\rho (l^{},\theta ^{})$$
(7)
where the angular part,

$`K(\theta ,\theta ^{})`$ $`=`$ $`{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\phi }{2\pi }}{\displaystyle \frac{d\phi ^{}}{2\pi }}|\mathrm{sin}\gamma |`$ (8)
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}k_{2n}𝒫_{2n}(\mathrm{cos}\theta )𝒫_{2n}(\mathrm{cos}\theta ^{})`$ (9)
can be expressed as a bilinear expansion in Legendre polynomials  $`𝒫_{2n}(\mathrm{cos}\theta )`$ with $`k_0=\pi /4,k_2=5\pi /32`$ etc. and $`𝒫_0(\mathrm{cos}\theta )=1`$, $`𝒫_2(\mathrm{cos}\theta )=\frac{1}{2}(3\mathrm{cos}^2\theta 1)`$. The excess free energy is therefore an infinite sum of moments of the density distribution $`\rho (l,\theta )`$, with weight functions $`l𝒫_{2n}(\mathrm{cos}\theta )`$:

$$\stackrel{~}{f}=\frac{4}{\pi }\underset{n=0}{\overset{\mathrm{}}{}}k_{2n}\left(𝑑l\stackrel{~}{d\theta }l𝒫_{2n}(\mathrm{cos}\theta )\rho (l,\theta )\right)^2$$
(10)
### B The $`𝒫_2`$ Onsager model

From Eq. (10) we see that the full Onsager theory does not give a truncatable excess free energy , i.e. one that is just a function of a finite number of moments $`\rho _i`$ of the density distribution $`\rho (l,\theta )`$. To obtain a more manageable theory, to which the efficient moment free energy method  can be applied, we now truncate the excess free energy after the second Legendre polynomial. We will thus consider the “$`𝒫_2`$ Onsager model” defined by the excess free energy

$$\stackrel{~}{f}=\frac{c_1}{2}\rho _1^2\frac{c_2}{2}\rho _2^2$$
(11)
where

$$\rho _1=𝑑l\stackrel{~}{d\theta }l\rho (l,\theta )$$
(12)
and

$$\rho _2=𝑑l\stackrel{~}{d\theta }l𝒫_2(\mathrm{cos}\theta )\rho (l,\theta )$$
(13)
are moments of the density distribution with weight functions

$`w_1(l,\theta )`$ $`=`$ $`l`$ (14)
$`w_2(l,\theta )`$ $`=`$ $`l𝒫_2(\mathrm{cos}\theta )`$ (15)
and the numerical constants are $`c_1=(8/\pi )k_0=2`$, $`c_2=(8/\pi )k_2=5/4`$. In our units $`\rho _1`$ is the rescaled volume fraction $`\varphi `$ of rods, $`\rho _1=(L_0/D)\varphi `$. We will denote below by $`\rho _0`$ the zeroth moment of the density distribution, $`\rho _0=𝑑l\rho (l)`$, with weight function $`w_0(l,\theta )=1`$; this is the number density of rods. The ratio $`\rho _1/\rho _0=l`$ then gives the rod length averaged over the normalized length distribution $`P(l)=\rho (l)/\rho _0`$, and $`\rho _2/\rho _11`$ expresses the degree of nematic ordering.

It is clear that the truncated model defined above will give approximate results compared to the predictions of the full Onsager theory. However, previous work  for the monodisperse case has investigated the convergence of the densities of coexisting isotropic and nematic phases when the expansion (10) is truncated at higher and higher orders. Already for the truncation after $`𝒫_2`$ the results were qualitatively correct, and this encourages us to study polydispersity effects within the truncated model defined above.

A different rationale for studying the $`𝒫_2`$ Onsager model could come from the fact that its excess free energy (11) is very similar to that of the polydisperse Maier-Saupe model for thermotropics (see e.g. Refs. ). There the interaction is not of the hard core type, but if one assumes that the interactions of rods of lengths $`l`$ and $`l^{}`$ again scales as $`ll^{}`$ the final expression for the excess free energy is identical to Eq. (11), with $`c_1=0`$ and $`c_2=c_2(T)`$ a function of temperature. The actual predictions of the two theories would nevertheless not necessarily be close, since within Maier-Saupe theory the overall rod density $`\rho _0`$ is assumed to be the same in all phases; this implies, for example, that a finite region of I-N coexistence cannot appear for monodisperse rods (see e.g. Ref. ). We therefore leave a study of polydisperse Maier-Saupe theory for future work.

To complete the specification of our model, we need to add to the free energy the ideal part, i.e. the free energy of an ideal mixture of polydisperse particles (see e.g. Refs. )

$$f_{\mathrm{id}}=\stackrel{~}{d\theta }𝑑l\rho (l,\theta )\left[\mathrm{ln}\rho (l,\theta )1\right]$$
(16)
Using the decomposition (3) and the relation (1), the total free energy (density) of our model is therefore:

$$f=𝑑l\rho (l)\left[\mathrm{ln}\rho (l)1\right]+𝑑l\stackrel{~}{d\theta }\rho (l)P_l(\theta )\mathrm{ln}P_l(\theta )+\stackrel{~}{f}$$
(17)
with $`\stackrel{~}{f}`$ given by (11). For a given density distribution $`\rho (l)`$, $`P_l(\theta )`$ is obtained by minimization of the free energy with respect of $`P_l(\theta )`$, subject to the constraint (4). Introducing corresponding Lagrange multipliers $`\kappa (l)`$ one has the condition

$`0`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta P_l(\theta )}}\left(f+{\displaystyle 𝑑l\stackrel{~}{d\theta }\kappa (l)P_l(\theta )}\right)`$ (18)
$`=`$ $`\rho (l)\left[\mathrm{ln}P_l(\theta )+1\right]+l\rho (l)\left[c_1\rho _1c_2\rho _2𝒫_2\right]+\kappa (l)`$ (19)
Solving for $`P_l(\theta )`$ gives (here and in the following $`𝒫_2𝒫_2(\mathrm{cos}\theta )`$)

$$P_l(\theta )=\frac{\mathrm{exp}(lc_2\rho _2𝒫_2)}{\stackrel{~}{d\theta }\mathrm{exp}(lc_2\rho _2𝒫_2)}$$
(20)
The orientational distributions are therefore determined only by $`\rho _2`$, which from (13) obeys the self-consistency equation

$$\rho _2=𝑑l\stackrel{~}{d\theta }\rho (l)\frac{l𝒫_2(\mathrm{cos}\theta )\mathrm{exp}(lc_2\rho _2𝒫_2)}{\stackrel{~}{d\theta }\mathrm{exp}(lc_2\rho _2𝒫_2)}$$
(21)
If there are several solutions $`\rho _2`$ for a given $`\rho (l)`$, the one with the smaller free energy (17) is the physical one.

### C Phase coexistence equations

At this point we need the expressions for the chemical potential $`\mu (l)`$ and the (osmotic) pressure $`\mathrm{\Pi }`$ in order to derive the phase equilibrium conditions. The chemical potential is obtained by functional differentiation of the free energy (17) with respect to $`\rho (l)`$; the variations of the $`P_l(\theta )`$ with $`\rho (l)`$ need not be considered since $`P_l(\theta )`$ is chosen to minimize $`f`$. After a little simplification one thus finds

$$\mu (l)=\mathrm{ln}\rho (l)\mathrm{ln}\stackrel{~}{d\theta }\mathrm{exp}(lc_1\rho _1+lc_2\rho _2𝒫_2)$$
(22)
The pressure follows from the Gibbs-Duhem relation as

$$\mathrm{\Pi }=f+𝑑l\rho (l)\mu (l)=\rho _0+\frac{c_1}{2}\rho _1^2\frac{c_2}{2}\rho _2^2$$
(23)
Imposing equality of chemical potentials in a set of $`P`$ coexisting phases, labelled by $`a=1\mathrm{}P`$, we can express the density distributions as

$$\rho ^{(a)}(l)=R(l)\stackrel{~}{d\theta }\mathrm{exp}(lc_1\rho _1^{(a)}+lc_2\rho _2^{(a)}𝒫_2)$$
(24)
where $`R(l)`$ is a function of $`l`$ common to all phases. Using Eq. (20), the full density distributions over rod lengths and orientations follow as

$$\rho ^{(a)}(l,\theta )=R(l)\mathrm{exp}(lc_1\rho _1^{(a)}+lc_2\rho _2^{(a)}𝒫_2)$$
(25)
The function $`R(l)`$ can be found from Eq. (24) and the requirement of particle conservation: if $`\rho ^{(0)}(l)`$ is the overall or “parent” density distribution and $`v^{(a)}`$ is the fraction of the system volume occupied by phase $`a`$, then

$$\underset{a}{}v^{(a)}\rho ^{(a)}(l)=\rho ^{(0)}(l)$$
(26)
This leads to

$$R(l)=\frac{\rho ^{(0)}(l)}{_av^{(a)}\stackrel{~}{d\theta }\mathrm{exp}(lc_1\rho _1^{(a)}+lc_2\rho _2^{(a)}𝒫_2)}$$
(27)
so that, from Eqs. (20,24), the density distributions in the coexisting phases are

$$\rho ^{(a)}(l,\theta )=\frac{\rho ^{(0)}(l)\mathrm{exp}(lc_1\rho _1^{(a)}+lc_2\rho _2^{(a)}𝒫_2)}{_bv^{(b)}\stackrel{~}{d\theta }\mathrm{exp}(lc_1\rho _1^{(b)}+lc_2\rho _2^{(b)}𝒫_2)}$$
(28)
In principle we thus have the conditions for coexistence between $`P`$ phases: In each phase, $`\rho _2`$ obeys the self-consistency condition (21) and $`\rho _1`$ a similar equation obtained by inserting (28) into (12); for the $`P`$ phase volume fractions $`v^{(a)}`$ we have the equality of pressure in all phases ($`P1`$ equations) plus the normalization $`_av^{(a)}=1`$. However, finding a starting point from which a numerical solution of this strongly coupled system of nonlinear equations converges successfully is very difficult. We use instead the moment free energy method, which gives exact results for some features of the phase diagram and otherwise allows us to approach the solution of the full phase equilibrium conditions with controllable accuracy.

## III The moment method

We only outline the construction of the moment free energy here and refer to Ref.  for details of its properties, and Ref.  for the extension to non-conserved degrees of freedom such as the rod orientations in our case. Terms which are linear in the conserved densities $`\rho (l)=\stackrel{~}{d\theta }\rho (l,\theta )`$ can be added to the free energy (17) without affecting the phase equilibria, since they merely add constants to the chemical potentials $`\mu (l)`$. We are therefore free to replace $`f`$ by

$$f=\stackrel{~}{d\theta }𝑑l\rho (l,\theta )\left[\mathrm{ln}\frac{\rho (l,\theta )}{r(l)}1\right]+\stackrel{~}{f}$$
(29)
with $`r(l)`$ an arbitrary function of $`l`$. Guided by the intuition that it is the moment densities $`\rho _1`$ and $`\rho _2`$ appearing in the free energy that drive phase separation, we now allow violations of the lever rule (26) as long as they do not affect these moments. All other details of the density distribution $`\rho (l,\theta )`$ are then found by minimizing the free energy (29) for the given values of $`\rho _1`$ and $`\rho _2`$ as defined in Eqs. (12,13). Inserting a Lagrange multiplier $`\lambda _i`$ for each moment $`\rho _i`$ ($`i=1,2`$), the minimum value of the free energy is then

$$f_{\mathrm{mom}}=\underset{i}{}\lambda _i\rho _i\rho _0+\stackrel{~}{f}$$
(30)
where $`\rho _0=𝑑l\rho (l)`$ is the number density of rods as defined previously. The density distribution at which this minimum free energy is obtained is

$$\rho (l,\theta )=r(l)\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(l,\theta )\right)$$
(31)
where the $`w_i(l,\theta )`$ are the weight functions (14,15) defining the moments; the $`\lambda _i`$ are determined implicitly by the requirement that the density distribution (31) gives the correct values for the $`\rho _i`$.

Eq. (30) defines the moment free energy; rather than being a functional of the density distribution $`\rho (l)`$, i.e. of an infinite number of conserved densities, it is a function of the densities $`\rho _1`$ and $`\rho _2`$. To find phase equilibria using the moment free energy, one proceeds as follows. The “moment chemical potentials” are defined as $`\mu _i=f_{\mathrm{mom}}/\rho _i`$ and can be written as  $`\mu _i=\lambda _i+\stackrel{~}{f}/\rho _i`$, giving in our case

$$\mu _1=\lambda _1+c_1\rho _1$$
(32)
$$\mu _2=\lambda _2c_2\rho _2$$
(33)
Since $`\rho _1`$ is conserved while $`\rho _2`$, which contains the rod orientations, is not, $`\mu _1`$ has to be identical in a set of coexisting phases while $`\mu _2`$ has to vanish in all the phases. The pressure calculated from $`f_{\mathrm{mom}}`$,

$$\mathrm{\Pi }=f_{\mathrm{mom}}+\rho _1\mu _1+\rho _2\mu _2=\rho _0+\frac{c_1}{2}\rho _1^2\frac{c_2}{2}\rho _2^2$$
(34)
also has to be identical in all phases. Finally, the lever rule has to be satisfied for the conserved moment $`\rho _1`$, i.e. $`\rho _1^{(0)}=_av^{(a)}\rho _1^{(a)}`$.

It is now easy to see that the moment free energy gives the correct phase equilibrium conditions, apart from the violations of the lever rule that it allows. In fact, equality of $`\mu _1`$ and vanishing $`\mu _2`$ imply from Eq. (31) that the coexisting phases calculated from the moment free energy can be written as

$$\rho ^{(a)}(l,\theta )=r(l)\mathrm{exp}\left(\mu _1llc_1\rho _1^{(a)}lc_2\rho _2^{(a)}𝒫_2\right)$$
(35)
This is of exactly the form (25) that we derived earlier from the full phase equilibrium conditions and the minimization of the free energy w.r.t. the rod angle distributions. The requirement of equal pressures in all phases is also captured exactly by the moment free energy, since the expression (34) actually coincides with the one derived from the original free energy of the $`𝒫_2`$ Onsager model, Eq. (23).

We have seen so far that any phase equilibria calculated from the moment free energy obey the exact conditions of equality of chemical potentials and pressure. The lever rule will be satisfied only for $`\rho _1`$ but not, in general, for the whole of $`\rho (l)`$. However, this final requirement will clearly also be fulfilled if all phases in an exactly calculated phase split of the $`𝒫_2`$ Onsager model are of the form (31). We can guarantee that this is the case at least until the onset of nematic order (the so-called “cloud point” of the isotropic phase), by choosing the so far unspecified function $`r(l)`$ to be equal to the parent density distribution $`\rho ^{(0)}(l)`$. This is the choice we make from now on. It works because at the isotropic cloud point the nematic phase still occupies a negligible fraction of the system volume, so that only the isotropic phase contributes to the denominator of Eq. (28). Explicitly one gets from Eq. (28), with obvious labels for the phases, $`\rho ^\mathrm{I}(l,\theta )=\rho ^{(0)}(l)`$ and $`\rho ^\mathrm{N}(l,\theta )=\rho ^{(0)}(l)\mathrm{exp}[lc_1(\rho _1^\mathrm{I}\rho _1^\mathrm{N})+lc_2\rho _2^\mathrm{N}𝒫_2]`$; as claimed, these are both of the form (31) with $`r(l)=\rho ^{(0)}(l)`$.

Beyond the onset of nematic phase coexistence, phase equilibria constructed from the moment free energy will in general not solve the full phase equilibrium conditions, because the lever rule cannot be satisfied with density distributions from the family (31). Violations of the lever rule can, however, be minimized by enlarging the family (31) to include the actual density distributions occurring in the coexisting phases. This can be done by retaining extra moments beyond those appearing in the excess free energy. Obviously, in the limit of retaining infinitely many extra moments, the moment free energy would eventually recover the full free energy of the model, but any computational efficiency would then be lost. The idea is therefore to add only a few extra moments to reduce the approximation caused by using the moment method, keeping their number sufficiently low to avoid making the computation too slow.

In order to enlarge the family (31) we could in principle use moments defined by any weight function $`w_i(l,\theta )`$. Since the corresponding moments $`\rho _i`$ do not appear in the excess free energy, their moment chemical potentials are $`\mu _i=\lambda _i+\stackrel{~}{f}/\rho _i=\lambda _i`$. The values of the extra Lagrange multipliers must therefore be the same in all phases, and from Eq. (31) we can write the density distributions predicted by the moment method as

$`\rho ^{(a)}(l,\theta )`$ $`=`$ $`\rho ^{(0)}(l)\mathrm{exp}\left[{\displaystyle \underset{i=1,2}{}}\lambda _i^{(a)}w_i(l,\theta )\right]`$ (36)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \underset{i1,2}{}}\lambda _iw_i(l,\theta )\right]`$ (37)
Comparing with Eq. (28), one sees that the second exponential, in which the extra weight functions appear, ought to provide a good approximation to the denominator of Eq. (27), which is not present in Eq. (31). Since that denominator is a function of $`l`$ only, the same needs to be true for the extra weight functions. We will always retain the zeroth moment, i.e. the total number density, with weight function $`w_0=1`$. This ensures that the “dilution line” of density distributions $`\rho (l)=\mathrm{const}.\times \rho ^{(0)}(l)`$ is included in the family (31), a useful feature when we want to scan across the phase diagram in the direction of increasing or decreasing density.

In order to choose the other extra weight functions, we apply the recently developed “adaptive method” . The intuition behind this is as follows. A phase split calculated by the moment method for a given point in the phase diagram will give a total density distribution $`\rho _{\mathrm{tot}}(l)=_av^{(a)}\rho ^{(a)}(l)`$ which in general differs from the parent $`\rho ^{(0)}(l)`$ (though agreeing in the value of the conserved moment $`\rho _1`$). A convenient measure of this violation of the lever rule is the “log-ratio” $`\mathrm{ln}\rho _{\mathrm{tot}}(l)/\rho ^{(0)}(l)`$. Since from Eq. (36) extra weight functions appear exponentially in the density distributions, including an extra weight function $`w_3(l)`$ similar in shape to the log-ratio should drive $`\rho _{\mathrm{tot}}(l)`$ closer to $`\rho ^{(0)}(l)`$; it extends the family (31) in the “right direction” in order to contain the exact distribution. The method therefore proceeds by iteratively including extra weight functions that are fitted to the current log-ratio; a proliferation of extra moments is avoided by combining the extra weight functions appropriately . In this way one never requires more than two extra moments in the calculation, but these are continually adapted to approach the solution of the full phase equilibrium equations. The quality of the approximation is measured by the average square log-ratio

$`\delta ={\displaystyle 𝑑l\left(\mathrm{ln}\frac{\rho _{\mathrm{tot}}(l)}{\rho ^{(0)}(l)}\right)^2}`$

and we could easily reach values of $`\delta <10^6`$, implying that our results are essentially identical to those that would result from a direct solution of the full phase equilibrium conditions for the $`𝒫_2`$ Onsager model. We verified this explicitly at some points in the phase diagram: with the solution obtained from our procedure as a starting point, a Newton-Raphson solver for the full phase equilibrium conditions converged and terminated at an essentially indistinguishable solution.

To fit the extra weight functions to the log-ratio in the above procedure, we initially followed the method of Ref. , representing the weight functions as linear combinations of some fixed set of basis functions. For some of the more numerically difficult situations, such as bimodal distributions with strongly different weights in the two peaks, it turned out to be more efficient instead to represent the extra weight functions as spline (piecewise cubic polynomial) fits to the log-ratio.

## IV Phase diagrams

In this section we will present our phase diagrams for the length polydisperse $`𝒫_2`$ Onsager model. We consider first a Schulz distribution of lengths, which has frequently been used to model simple unimodal distributions (see e.g. Ref. ). We will find strong differences compared to the results for the bi- and tridisperse (full) Onsager model. By contrasting with the predictions of the $`𝒫_2`$ Onsager model for bidisperse systems, we show that these differences are not due to our truncation of the model but rather due to qualitative differences in phase behaviour for discrete and continuous length distributions. Finally, we study a bimodal mixture of two Schulz distributions. This allows us to analyse in detail how the topology of the phase diagram changes as we interpolate between the two extremes of bidisperse and unimodal distributions.

### A Unimodal length distribution

We study a parent phase with density distribution $`\rho ^{(0)}(l)=\rho _0^{(0)}P^{(0)}(l)`$ where the normalized length distribution is of the Schulz form:

$$P^{(0)}(l)=\frac{(z+1)^{z+1}}{\mathrm{\Gamma }(z+1)}l^z\mathrm{exp}[(z+1)l]$$
(38)
This gives an average rod length of one for the parent, a convention we will follow throughout since a different average parent rod length could always be absorbed into the reference length $`L_0`$. The parameter $`z`$ controls the width of the distribution, more precisely its normalized standard deviation $`\sigma `$, which from now on we will refer to simply as “polydispersity”:

$$\sigma ^2=\frac{l^2l^2}{l^2}=\frac{1}{z+1}$$
(39)
The limit $`z\mathrm{}`$ therefore corresponds to the monodisperse case, where the Schulz distribution degenerates to a $`\delta `$-peak. Decreasing $`z`$ will give increasing polydispersity; we limit ourselves to $`z0`$ and thus $`\sigma 1`$ since otherwise the Schulz distribution (38) has a power-law divergence for $`l0`$.

In Fig. 1 we show, for each value of $`\sigma `$, the isotropic cloud and nematic shadow – giving the rod number densities in the coexisting isotropic and nematic phases at the point where the nematic first appears – and the nematic cloud and isotropic shadow, which relate to the point where the isotropic phase disappears as density is increased. As anticipated, the coexistence region, which is delimited by the two cloud curves, broadens strongly with increasing polydispersity. The other interesting feature, which was not seen in our earlier study of polydispersity effects in the Zwanzig model , is the crossing between the cloud curves and the corresponding shadow curves; for values of $`\sigma `$ above the crossing points the number density of the isotropic phase is larger than the one of the coexisting nematic phase.

If we switch to a different representation of the phase diagram by plotting the (scaled) rod volume fractions $`\rho _1`$ rather than the number densities $`\rho _0`$ in the various phases (Fig. 2) the crossing between cloud and shadow curves disappears. So the rod volume fraction in the nematic is always larger than in the isotropic phase, even though the isotropic can have the larger number density. Since $`\rho _1/\rho _0`$ is the average rod length in each phase, this is clear evidence of a strong fractionation effect, with the longer rods being found predominantly in the nematic phase. We show this explicitly in Fig. 3. Bearing in mind that the cloud phases have the same length distributions as the parent and therefore an average rod length of unity, one sees that indeed the average rod length in the nematic shadow is much larger than in its corresponding isotropic cloud, and that of the nematic cloud is much larger than that of the corresponding isotropic shadow.

### B Bidisperse rod lengths

From Figures 1, 2 and 3 we saw that there is no I-N-N three-phase region and associated N-N demixing within the present model, at least for a Schulz distribution of lengths. On the other hand, previous work on bi- and tridisperse systems  shows that such a three-phase separation is possible within the full Onsager theory. The absence of a three-phase region observed above could then be due to either the difference in the length distributions investigated (continuous and unimodal versus discrete with two or three different lengths), or to our truncation of the Onsager theory that gave the $`𝒫_2`$ Onsager model. To ascertain which of these applies, we now study the phase diagram of a bidisperse system of rods within the $`𝒫_2`$ Onsager model. The normalized parent length distribution is then

$$P^{(0)}(l)=\left[(1x)\delta (ll_1)+x\delta (ll_2)\right]$$
(40)
where $`l_1<l_2`$ and $`x`$ is the (number) fraction of longer rods. We again fix the average length to unity, so that $`l_1`$ and $`l_2`$ can be written in terms of their ratio $`q=l_2/l_1`$ as

$$l_1=\frac{1}{1x+xq},l_2=ql_1$$
(41)
The polydispersity is given by

$$\sigma ^2=\frac{x(1x)(q1)^2}{(1x+xq)^2}$$
(42)
and $`\sigma `$ becomes zero for $`x=0`$ and $`x=1`$, where we recover a monodisperse system; its maximum for given $`q`$ is reached at $`x=1/(q+1)`$.

In Fig. 4 we show the phase diagram  for rod length ratio $`q=2.5`$. The density of the parent at the transition from single-phase to two-phase regions is plotted against the number fraction of longer rods on the vertical axis; effectively, we are only showing the cloud curves and omitting the shadow curves. This “phase boundaries only” representation will be useful below when three-phase equilibria appear (which are represented on the cloud and shadow curves only in terms of the exceptional points where direct transitions from one to three phases occur). Comparing with the Schulz distribution case in Fig. 1, where the solid lines correspond to those of Fig. 4, we see that again the coexistence region broadens as we move from either of the monodisperse limits $`x=0`$ and $`x=1`$ towards larger polydispersity.

For comparison with previous results on bidisperse systems obtained from the full Onsager theory , it will also be useful to show our phase diagrams in the representation employed in Refs. , which uses the variables $`\stackrel{~}{\varphi }_1=(L_1/D)\varphi _1`$ and $`\stackrel{~}{\varphi }_2=(L_1/D)\varphi _2`$ instead of our $`x`$ and $`\rho _0`$. Here $`\varphi _1`$ and $`\varphi _2`$ are the volume fractions of short and long rods, and $`L_1=L_0l_1`$ is the unnormalized length of the short rods. With $`N_1`$ and $`N_2`$ the number of short and long rods, we have

$$\stackrel{~}{\varphi }_i=\frac{L_1}{D}\varphi _i=\frac{L_1}{D}\frac{\pi }{4}D^2L_i\frac{N_i}{V}=\frac{N_i}{N}l_1l_i\left(\frac{\pi }{4}DL_0^2\frac{N}{V}\right)$$
(43)
The term in brackets is just our dimensionless rod number density $`\rho _0`$, so that

$`\stackrel{~}{\varphi }_1`$ $`=`$ $`(1x)l_1^2\rho _0`$ (44)
$`\stackrel{~}{\varphi }_2`$ $`=`$ $`xl_1l_2\rho _0`$ (45)
This gives us the relation between the two representations of the phase diagram. As an aside, we note that our convention for choosing the normalized rod lengths $`l_1=L_1/L_0`$ and $`l_2=L_2/L_0`$ differs from that of Ref. ; we vary $`l_1`$ and $`l_2`$ according to Eq. (41) to maintain a constant average length in the parent, while in Ref.  constant values of $`l_1`$ and $`l_2`$ are used. However, these assignments can alternatively be thought of as arising from a different choice of the reference length $`L_0`$ for identical values of $`L_1`$ and $`L_2`$. Since $`L_0`$ does not enter the definitions (44,45), the values of the $`\stackrel{~}{\varphi }_1`$ and $`\stackrel{~}{\varphi }_2`$ remain unaffected. We show the $`\stackrel{~}{\varphi }_1`$, $`\stackrel{~}{\varphi }_2`$-representation for the case $`q=2.5`$ in Fig. 5.

By including Legendre polynomials up to $`𝒫_{14}`$, and determining the orientational distribution functions $`P_l(\theta )`$ numerically, Lekkerkerker et al.  obtained good approximations to the phase diagrams for the full Onsager theory of bidisperse rods. They found that at $`q=5`$ there is a re-entrant nematic phase, while at $`q=2`$ there is no such re-entrance. Buining and Lekkerkerker  also found no re-entrance at $`q=2.5`$; the latter finding agrees with our above results. Using a Gaussian variational approximation for the $`P_l(\theta )`$, Odijk and Lekkerkerker  showed that re-entrance is generically to be expected for large $`q`$. Vroege and Lekkerkerker  later found, using the same approach, that in addition to the re-entrance one obtains three-phase I-N-N and two-phase N-N equilibria for $`q`$ above $`3.17`$. Their approximation predicted that the boundaries of the N-N coexistence region should be independent of the rod number density $`\rho _0`$, corresponding to horizontal lines in a ($`\rho _0`$,$`x`$) representation of the phase diagram or to lines through the origin in a ($`\stackrel{~}{\varphi }_1`$,$`\stackrel{~}{\varphi }_2`$) plot. Van Roij and Mulder  later showed that these statements are exactly true (only) in the limit of large $`\rho _0`$; this means in particular that the N-N coexistence region is not closed off by an N-N critical point. Birshtein et al. , using a trial function of the type used by Onsager , had earlier come to the opposite conclusion. However, they emphasized (as did Abe and Flory ) the difficulties of solving the phase equilibrium conditions at large density, which reduced the reliability of their results.

To compare the above findings from treatments of the full Onsager theory with the $`𝒫_2`$ Onsager model, we now investigate how the phase diagram changes with increasing $`q`$.

Figs. 6, 7, 8 and 9 show the phase diagrams for rod length ratios $`q=3.5`$, $`q=7`$ and $`q=12`$ in both representations ($`\rho _0,x`$ and $`\stackrel{~}{\varphi }_1,\stackrel{~}{\varphi }_2`$). At $`q=3.5`$ a re-entrance of the nematic phase has almost appeared; at $`q=7`$ this is fully developed. At $`q=12`$ (Fig. 8), finally, we have in addition a three-phase I-N-N region bordered by a region (N-N) of coexistence of two nematics . The N-N region is closed off by an N-N critical point, where the moments $`\rho _1`$ and $`\rho _2`$ in the two nematics become identical. These features persist at higher $`q`$; for $`q=30`$, for example, the phase diagram (not shown) has the same structure as Fig. 8 but with both the re-entrance and the three-phase region being located at smaller $`x`$ and extending to higher densities. For even larger $`q`$ the numerical evaluation of the phase diagrams becomes more troublesome because the interesting features of the phase diagram shift towards smaller and smaller $`x`$ where even with the moment method the phase equilibrium conditions cannot be reliably solved numerically. However, a simple approximate treatment (see App. A) can be used to make progress and understand how the phase diagram scales for large $`q`$.

In summary of our results for the bidisperse case, an increase in $`q`$ produces similar qualitative changes in the structure of the phase diagram as for the full Onsager theory. The quantitative details of course differ, with the onset of both the nematic re-entrance and the three-phase region shifted to higher $`q`$. The only qualitative difference caused by the fact that we are using the truncated $`𝒫_2`$ Onsager model is that, in contrast to the full Onsager theory , the N-N region does not extend to large densities but is instead closed by a critical point, whose properties we discuss further in App. A.

This distinction between the $`𝒫_2`$ Onsager model and the full Onsager theory can be made plausible by considering the effect of our truncation on the excluded volume term, i.e., the excess free energy, in the limit of large densities $`\rho _0`$ where nematic phases will be strongly ordered. In the full Onsager theory, it can be shown  that the excess free energy density is simply $`\stackrel{~}{f}=2\rho _0`$ for large $`\rho _0`$ so that the excluded volume per particle is $`\stackrel{~}{f}/\rho _0=2`$. This simple scaling with $`\rho _0`$ arises because the typical rod angles $`\theta `$ are $`1/\rho _0`$ for large $`\rho _0`$; the angular part $`K(\theta ,\theta ^{})`$ of the excluded volume interaction, as given by Eq. (8), scales linearly with the typical values of $`\theta `$ and $`\theta ^{}`$ and is thus also proportional to $`1/\rho _0`$. The high density phase behaviour thus arises solely from a competition between the two components of the ideal part of the free energy (17), the entropy of mixing and the orientational entropy; N-N demixing occurs when the resulting gain in orientational entropy dominates the loss of entropy of mixing. In our $`𝒫_2`$ truncation, on the other hand, $`K`$ is approximated as $`K(\theta ,\theta ^{})c_1c_2𝒫_2(\mathrm{cos}\theta )𝒫_2(\mathrm{cos}\theta ^{})`$ which becomes $`a+b[\theta ^2+(\theta ^{})^2]`$ for small angles, with constants $`a`$ and $`b`$. For large $`\rho _0`$, when the typical values of $`\theta `$ and $`\theta ^{}`$ again become small, this approaches a constant rather than decrease to zero as in the full Onsager theory, giving a quadratic density scaling of the excess free energy. The same conclusion can be reached from the moment description: for $`\theta 0`$ we have $`\rho _2=\rho _1`$ and the excess free energy density becomes $`\frac{1}{2}(c_1c_2)\rho _1^2`$. The quadratic scaling of the excess free energy with density, as opposed to the linear scaling in the full Onsager theory, will act as a driving force against phase separation. This is in line with our findings above regarding the absence of N-N demixing at high density. One can in fact show formally (see App. A) that N-N coexistence is not possible in the $`𝒫_2`$ Onsager model in the limit of large $`\rho _0`$.

The above argument can in fact be extended to all possible truncations of the full Onsager theory, since the inclusion of any finite number of higher order Legendre polynomials will still give a behaviour $`K(\theta ,\theta ^{})a+b[\theta ^2+(\theta ^{})^2]`$ of the orientational part of the excluded volume for small angles; the correct linear scaling with $`\theta `$ and $`\theta ^{}`$ of the full theory is obtained only when all terms in the infinite series are retained. Equivalently, in the moment description one would have a finite number of moments analogous to $`\rho _2`$, but these would all saturate to values proportional to $`\rho _1`$ at high density, giving again the quadratic scaling $`\stackrel{~}{f}\rho _1^2`$. We conclude, therefore, that there cannot be any N-N region at sufficiently high density for any truncation of the full Onsager theory; this is consistent with van Roij and Mulder’s explanation  for Birshtein’s result of an N-N critical point in his approximate treatment of bidisperse Onsager theory. Of course, for any finite density we must eventually recover the correct phase behaviour as more and more higher order Legendre polynomials are included to give an increasingly good representation of $`K(\theta ,\theta ^{})`$ for nonzero angles. One would expect, therefore, that successively higher order truncations of Onsager theory would lead to N-N coexistence regions extending to higher and higher densities, and eventually diverging as the order of truncation is taken to infinity.

Before passing to a different length distribution, let us note that in the volume fraction representation (e.g. Fig. 9) of the phase diagrams of our bidisperse system all tielines connecting coexisting phases are necessarily straight. This implies, in particular, that the boundaries of the three phase region, which are formed by two-phase tielines, are straight. The same argument no longer applies in the scenario that we study next: when more than two different rod lengths are present, the two-dimensional phase diagrams are cuts through a higher-dimensional phase diagrams, and generic tielines no longer lie within this cut plane.

### C Bimodal length distribution

Our results for the bidisperse case suggest that the $`𝒫_2`$ Onsager model correctly reproduces most of the qualitative features of the phase diagram of the full Onsager theory, except for N-N coexistence at high densities. This suggests similar qualitative agreement also in the case of unimodal length distributions. We therefore conclude that the differences between unimodal and bidisperse distributions that we found, in particular the absence of re-entrant features and N-N and I-N-N phase equilibria in the unimodal case, are physical and not due to our truncation of the model. To understand in more detail how these features arise in dependence on the shape of the length distribution, we now turn to a family of length distributions that lets us interpolate between the two limits of bidisperse and unimodal distributions. The normalized parent length distributions we consider are mixtures of two Schulz distributions, of the form

$$P^{(0)}(l)=(1x)l_1^1S\left(l/l_1\right)+xl_2^1S\left(l/l_2\right)$$
(46)
Here $`S(l)`$ denotes the Schulz distribution (38) and $`l_1`$ and $`l_2`$ are again given by Eq. (41) to constrain the average rod length to unity. Let us call $`\sigma _\mathrm{S}`$ the normalized width (standard deviation) of the Schulz distribution, to distinguish it from the normalized overall standard deviation $`\sigma `$ of the parent distribution; the distribution (46) is then characterized by the three parameters $`q`$, $`x`$ and $`\sigma _\mathrm{S}`$. For $`\sigma _\mathrm{S}0`$, we recover the bidisperse case; on the other hand, as $`\sigma _\mathrm{S}`$ increases for given $`x`$ and $`q`$, the distribution eventually becomes unimodal as shown in Fig. 10. For given $`q`$ and $`x`$, the value of $`\sigma _\mathrm{S}`$ at which this happens is determined by the condition that there must be a value of $`l`$ for which $`P^{(0)}(l)`$ has vanishing first and second derivative; the resulting set of equations can easily be solved numerically for $`\sigma _\mathrm{S}`$.

Our aim is now to understand how the occurrence of the more exotic features of the phase diagram, such as an I-N-N three-phase region or a re-entrant nematic phase, depends on the rod length distribution; one might, for example, suspect a dependence on bimodality (the presence of two separate peaks in the distribution) or on the overall polydispersity.

We will proceed in two separate ways. First, we fix a value of $`q`$ (=12) for which in the bidisperse case ($`\sigma _\mathrm{S}=0`$) we have a three-phase region and re-entrance in the nematic phase. We then increase $`\sigma _\mathrm{S}`$ and monitor the corresponding variation of the ($`\rho _0,x`$) phase diagram. When $`\sigma _\mathrm{S}`$ becomes sufficiently large, we have a unimodal length distribution and thus expect to see behaviour more akin to that for the Schulz length distribution. For fixed $`q`$, the three-phase region should thus disappear when $`\sigma _\mathrm{S}`$ reaches a sufficiently large value. The process will essentially be the opposite of what we have done so far: in the previous sections we began from a unimodal phase diagram and moved to the bidisperse case with larger and larger rod length ratio $`q`$. Now we start from the largest value of $`q`$ that we used for a bidisperse distribution, and widen the two peaks of the distribution until we return to a unimodal (though not a simple Schulz) distribution.

The second procedure will be to fix values of $`x`$ and $`q`$ for which we have coexistence of three phases in the bidisperse limit $`\sigma _\mathrm{S}0`$, and to start increasing the polydispersity. This lets us check whether the existence of a three-phase region is linked directly to the bimodality of the parent, by comparing the value of $`\sigma _\mathrm{S}`$ for which the region of three-phase equilibrium disappears to that where the parent length distribution changes from bimodal to unimodal.

Adopting the first procedure, we have found the phase diagrams for a small normalized width $`\sigma _\mathrm{S}=0.1`$ of the two peaks in the bimodal length distribution. For $`q=3.5`$ the result is a phase diagram that is practically identical to Fig. 6; the system is still effectively bidisperse. More interestingly, the same conclusion applies when we introduce a small peak width for $`q=12`$: with $`\sigma _\mathrm{S}=0.1`$, we again get phase diagrams that are essentially indistinguishale from the bidisperse case shown in Figs. 8 and 9. In particular, neither the three-phase region nor the re-entrant nematic phase are significantly affected by such a small degree of polydispersity of the two peaks. Larger values of $`\sigma _\mathrm{S}`$, on the other hand, do have a significant effect. If the width of the peaks is increased to $`\sigma _\mathrm{S}=0.6`$, for example (Fig. 11), the phase diagram changes noticeably. In particular, the three-phase region has now disappeared, or at least shrunk to values below $`x=10^3`$ where our numerical calculations become unreliable, and the nematic re-entrance has become much less pronounced. This is in line with the intuition that for increasing $`\sigma _\mathrm{S}`$ behaviour similar to that for the unimodal Schulz distribution should be recovered.

The results so far confirm that, in order for a three-phase region to exist in the phase diagram, the two peaks of the parent length distribution must, in some sense, be sufficiently well “separated”. This notion is supported by an analysis of the large $`q`$ limit in App. A. There we show that, for large $`q`$ and values of $`x`$ of order $`1/q`$, there is an N-N region and a corresponding three-phase region for all $`\sigma _\mathrm{S}`$, up to the maximum value $`\sigma _\mathrm{S}`$ that we consider; correspondingly, one easily sees that in this limit the length distribution always has two separate maxima.

We now proceed to a closer investigation of the link between the presence of three-phase equilibria and bimodality of the length distribution, using our second approach explained above: we fix $`q`$ and $`x`$ so that the three-phase coexistence occurs in the bidisperse limit $`\sigma _\mathrm{S}=0`$, and then find the value of $`\sigma _\mathrm{S}`$ at which the three-phase region disappears, comparing with the value where the parent changes from bimodal to unimodal.

Figs. 12 and 13 show for $`q=12`$ and two different values of $`x`$ that the three-phase region indeed disappears on increasing $`\sigma _\mathrm{S}`$ above a certain value. This threshold value of $`\sigma _\mathrm{S}`$ is smaller for larger $`x`$, showing that in the ($`\rho _0`$,$`x`$) phase diagram the top corner of the three-phase region shrinks towards smaller values of $`x`$ as $`\sigma _\mathrm{S}`$ is increased. The $`x`$-value marking the bottom end of the three-phase region will also depend on $`\sigma _\mathrm{S}`$. One would suppose that it increases with $`\sigma _\mathrm{S}`$ until the $`x`$-values for the top and bottom eventually coincide and the three-phase region disappears. However, since in the phase diagrams we considered the lower end of the three-phase region is located at extremely small $`x`$-values we were not able to confirm this explicitly.

Repeating the above calculations for a range of values of $`x`$, we arrive at the dashed curve in Fig. 14: at a given $`x`$, we have three phase coexistence for values of $`\sigma _\mathrm{S}`$ below the curve. Conversely, horizontal sections through the graph give us the extent in $`x`$ of the three-phase region for given $`\sigma _\mathrm{S}`$. Since the three-phase region cannot extend all the way to $`x=0`$, this shows that the dashed curve must eventually bend down towards the horizontal axis as $`x`$ decreases and meet the axis at some small but nonzero value of $`x`$. We have not been able to detect this part of the curve, however, since the $`x`$-values concerned are too small for our numerics to be accurate.

We can now analyse whether there is a correlation between three-phase coexistence region and bimodality of the parent length distribution. For values of $`x`$ and $`\sigma _\mathrm{S}`$ below the solid curve in Fig. 14, the parent is bimodal; above the curve we have a unimodal length distribution. The comparison of the regions under the dashed and solid curves suggests that bimodality may be required for three-phase coexistence to occur (except possibly in the region of very small $`x`$, where our data on three-phase coexistence are incomplete), but does not necessary entail it.

One might therefore look for additional properties of the length distribution as predictors of three-phase coexistence, e.g. the overall polydispersity $`\sigma `$. Eq. (42) showed that for a bidisperse system the overall polydispersity $`\sigma `$ vanishes at $`x=0`$ and $`x=1`$ and has a maximum for $`x=1/(q+1)`$; the values of $`x`$ where three-phase coexistence occurs in Fig. 14 are smaller but at least of the same order as those where $`\sigma `$ is maximized. Also, the maximum polydispersity in a bidisperse distribution with a given length ratio $`q`$ increases with $`q`$, in line with the general trend that separation into three phases requires sufficiently large $`q`$. Nevertheless, it is clear that the trend of three-phase separation occurring preferably for systems with large $`\sigma `$ cannot be universally true. In fact, it is easy to see this from Fig. 14. When we increase the width $`\sigma _\mathrm{S}`$ of the two peaks in a bimodal distribution the overall polydispersity $`\sigma `$ also increases; but this actually suppresses three-phase coexistence rather than enhance it.

Loosely speaking, our results show that, in order for the system to be able to separate into three phases, we need a sufficiently large disparity between long and short rods in the system, i.e. large $`q`$; but at the same time the concentration of short and long rods must not be too similar. The second condition would account for the fact that the three-phase region develops only at small values of $`x`$, and disappears as $`\sigma _\mathrm{S}`$ is increased. If this interpretation is correct, then we should be able to see three-phase coexistence even for a unimodal rod length distribution, as long as it contains a sufficiently large number of rods that are much longer than the average. A good candidate for this is a log-normal distribution, which has been shown already in the context of polymer solutions to present some interesting features caused by the presence of very large particles . We will present our analysis of this case in a separate publication .

## V Conclusion

We have analysed the phase equilibria of the $`𝒫_2`$ Onsager model of hard rods with length polydispersity. The model is defined by a truncation of the angular dependence of the excluded volume interaction of the Onsager theory after the second Legendre polynomial $`𝒫_2(\mathrm{cos}\theta )`$. Within this model we derived the exact phase equilibrium conditions, but these are still rather difficult to solve numerically. We therefore exploited the fact that the model is truncatable – the excess free energy is a function of only two moments ($`\rho _1`$ and $`\rho _2`$) of the density distribution – and used the moment method, from which numerical solutions for phase equilibria can be obtained with well-controlled accuracy; the onset of nematic ordering from the isotropic side is found exactly by construction of the moment free energy.

For a fully polydisperse case with a unimodal (Schulz) length we found some of the common features of polydisperse systems. In particular, we observed a strong fractionation (Fig. 3) and a pronounced broadening of the coexistence region (Fig. 1) with increasing polydispersity; the latter is defined as the normalized standard deviation of the length distribution. Fractionation effects were strong enough to lead to a crossing of the cloud point curves with the corresponding shadow curves in the number density representation (Fig. 1), but not in the volume fraction representation (Fig. 2). This implies that while the nematic phase always has a larger rod volume fraction than the isotropic, it can actually have a smaller rod number density.

For the unimodal length distribution we did not observe any three-phase I-N-N or two-phase N-N coexistence, or re-entrant behaviour of the nematic phase, all of which had been found in the Onsager model with bi- and tridisperse length distributions. In order to understand whether this was due to the introduction of a continuous and unimodal length distribution, or to the truncation that defined the $`𝒫_2`$ Onsager model, we next studied the bidisperse case within the $`𝒫_2`$ model. As the length ratio $`q`$ increased, we indeed found (Fig. 6 and 7) the development of a re-entrance of the nematic phase, and eventually also I-N-N and N-N coexistence (Fig. 8) for larger $`q`$. Thus all the qualitative features of the bidisperse phase diagrams of the full Onsager theory are recovered. The only exception is the presence in the $`𝒫_2`$ model of a critical point that closes the N-N coexistence region at high density, while for the full Onsager theory it is known  that for large $`q`$ such a critical point does not exist. We explained this difference in terms of the different behaviour of the excluded volume interaction in the limit of strongly ordered rods, and argued that it would persist for any truncation of the full Onsager theory involving only a finite number of Legendre polynomials.

Finally, in order to interpolate between the different behaviours in the unimodal and bidisperse cases, we introduced polydispersity into the latter system, by considering mixtures of two Schulz distributions of width $`\sigma _\mathrm{S}`$ peaked at different rod lengths; the bidisperse system is recovered for $`\sigma _\mathrm{S}0`$. For small $`\sigma _\mathrm{S}`$, i.e., small “broadening” of the bidisperse “peaks”, the phase diagrams remained essentially as in the bidisperse limit, while larger values of $`\sigma _\mathrm{S}`$ lead to a disappearance of the I-N-N and N-N coexistence regions (Fig. 11). It is natural to suppose that this change is related to the fact that the length distribution changes from bimodal to unimodal as $`\sigma _\mathrm{S}`$ is increased. However, a closer investigation of this point showed that bimodality, although correlated with the occurrence of I-N-N and N-N coexistence, is not sufficient to guarantee that such coexistence will occur. Indeed, this follows already from the fact that in the bidisperse system (which is always bimodal) a minimum value of the length ratio $`q`$ is needed for three-phase coexistence. Our results suggest that, also in the polydisperse case, a large ratio between the typical (e.g., mean) values of the two Schulz length distributions is required. However, the result that three-phase coexistence generally only occurs for small values values of the number fraction $`x`$ of long rods suggests that in addition the length distribution must be rather asymmetric, with fewer long rods than short ones. When the width $`\sigma _\mathrm{S}`$ of the Schulz distributions is increased too much, this (here only vaguely defined) “asymmetry” is lost, and the three-phase region disappears. In order to assess whether this intuition is correct, we are currently investigating the phase behaviour for log-normal length distributions ; while unimodal, these contain rods of very dissimilar lengths and are also strongly asymmetric. The intuition developed above would lead one to expect three-phase coexistence.

We did not investigate in detail the conditions for occurrence of the re-entrance in the nematic phase boundary. This phenomenon seemed rather robust compared to I-N-N and N-N coexistence, although in the unimodal (Schulz distribution) case it was again absent. In the bidisperse case the re-entrance appeared at lower values of $`q`$ than the three-phase coexistence, and also did not require the asymmetry $`x`$ between long and short rod number densities to be quite so extreme, extending e.g. almost up to $`x=0.3`$ for $`q=12`$ (Fig. 8). In the case of a mixture of Schulz distributions case it also survived to larger values of the width $`\sigma _\mathrm{S}`$ of the two peaks of the distribution. An intriguing observation is that, for the bidisperse systems, the re-entrance begins (starting from small $`x`$) very close to the value of $`x=1/(q+1)`$, where the polydispersity $`\sigma `$ of the system is maximal; see Eq. (42). This is also the value of $`x`$ where the volume fractions occupied by the short and long rods become equal; see Eq. (A11). We have not been able to find a simple physical explanation for this observation.

In summary, our study of the $`𝒫_2`$ Onsager model shows that the details of the rod length distribution can have profound effects on the phase behaviour. In particular, while bidisperse systems can show a re-entrant nematic phase boundary as well as I-N-N and N-N coexistence, these features are absent from the phase diagram for a unimodal length distribution; they also disappear as one interpolates between the bidisperse and unimodal extremes. Put differently, the more “exotic” aspects of the phase behaviour of mixtures of hard rods are not expected to be observed in systems with only a moderate spread of rod lengths around a mean value, but require sufficiently wide (and, as we argued above, asymmetric) rod length distributions. Our comparison of the predictions of the $`𝒫_2`$ model and the full Onsager theory for the bidisperse case suggests that these conclusions are not artefacts of the truncation used to the construct the $`𝒫_2`$ model. Rather, we expect that they would qualitatively also be found in an analysis of the effects of length polydispersity in the framework of the full Onsager model; such an analysis remains a challenging problem for future work.

We comment finally on the presence of the N-N critical point within the bidisperse $`𝒫_2`$ Onsager model. As discussed above, the theoretical work of van Roij and Mulder  shows that for sufficiently large ratios $`q`$ of rod lengths such a critical point will not occur in the full Onsager theory; earlier reports to the contrary  appear to be artefacts of numerical approximations. However, the results of Ref.  show only that the N-N coexistence region is open, i.e., extends to arbitrarily large densities, for $`q`$ above $`3.17`$. The possibility remains that for a small range of $`q`$-values below 3.17 an N-N coexistence region would already exist, at moderate densities and closed off by a critical point; the topology of the phase diagram for such $`q`$ would be akin to the $`𝒫_2`$ Onsager phase diagrams that contain N-N and I-N-N coexistence regions. As $`q`$ increases, the critical point would then have to move to larger and larger densities so that the N-N region eventually becomes open at $`q3.17`$. It would be worth revisiting the bidisperse Onsager model to confirm this scenario.

Acknowledgement: PS acknowledges financial support through EPSRC grant GR/R52121/01.

## A The strong nematic ordering approximation

In this appendix we describe an approximation that is useful for analysing the phase diagram of the $`𝒫_2`$ Onsager model at high densities. This strong nematic ordering (SNO) approximation becomes exact in the limit $`\rho _0\mathrm{}`$ and allows one to show formally that there can be no N-N coexistence in this limit. For finite $`\rho _0`$ it is approximate but nevertheless useful for understanding some limits of the phase behaviour, e.g. for large values of the rod length ratio $`q`$ in bidisperse and bimodal length distributions. As discussed in Sec. IV B, it is in the regime of strong nematic ordering (small angles $`\theta `$) that the $`𝒫_2`$ Onsager model and the full Onsager theory differ most strongly, so the analysis here focusses on the properties that are particular to the $`𝒫_2`$ Onsager model.

The SNO approximation is based on the assumption that all rod angles $`\theta `$ in any nematic phases are reasonably small, so that we can approximate $`\mathrm{cos}^2\theta 1\theta ^2`$ and hence

$$𝒫_2(\mathrm{cos}\theta )1\frac{3}{2}\theta ^2$$
(A1)
From Eq. (20), the orientational distributions are then $`P_l(\theta )\mathrm{exp}[\frac{3}{2}lc_2\rho _2\theta ^2]`$ and so the SNO assumption of small angles $`\theta `$ corresponds to $`\rho _2`$ being large enough to have $`l\rho _21`$ for all lengths $`l`$ that contribute significantly. Using the approximation (A1), we can get explicit expressions for the length dependent chemical potentials. The relevant angular integral can be evaluated as

$`{\displaystyle \stackrel{~}{d\theta }e^{lc_2\rho _2𝒫_2}}`$ $``$ $`{\displaystyle _0^{\pi /2}}𝑑\theta \mathrm{sin}\theta e^{lc_2\rho _2(3/2)lc_2\rho _2\theta ^2}`$ (A2)
$``$ $`{\displaystyle \frac{\mathrm{exp}\left(lc_2\rho _2\right)}{3lc_2\rho _2}}`$ (A3)
where the symmetry of the integrand under $`\theta \pi \theta `$ has been exploited and corrections that are exponentially small in $`l\rho _2`$ have been neglected. The chemical potentials (22) thus become

$$\mu ^\mathrm{N}(l)\mathrm{ln}\rho (l)+l(c_1\rho _1c_2\rho _2)+\mathrm{ln}(3lc_2\rho _2)$$
(A4)
for a nematic phase, while for an isotropic phase ($`\rho _2=0`$) the angular integral (A2) equals unity and the chemical potentials are

$$\mu ^\mathrm{I}(l)=\mathrm{ln}\rho (l)+lc_1\rho _1$$
(A5)
Using the SNO assumption (A1) we can also evaluate the integrals appearing in the self-consistency equation (13) for $`\rho _2`$; the latter can then be solved explicitly to give, after some algebra,

$$\rho _2\frac{\rho _1}{2}\left(1+\sqrt{1\frac{4\rho _0}{c_2\rho _1^2}}\right)$$
(A6)
In the limit of large density $`\rho _0`$, since $`\rho _1=l\rho _0\rho _0`$, the square root approaches unity and thus $`\rho _2\rho _1`$. This is as expected, since the weight functions (14,15) for $`\rho _1`$ and $`\rho _2`$ only differ by a factor of $`𝒫_2(\mathrm{cos}\theta )`$ which becomes equal to one for maximal nematic order, i.e. for $`\theta 0`$. The equality $`\rho _2=\rho _1`$ implies that $`\rho _2`$ becomes large as the density increases; for $`\rho _0\mathrm{}`$ it diverges and the SNO approximation becomes exact. One can now easily show that in this large density limit there can be no coexistence between two nematic phases. The osmotic pressure (34) becomes $`\mathrm{\Pi }=\rho _0+\frac{1}{2}(c_1c_2)\rho _1^2\frac{1}{2}(c_1c_2)\rho _1^2`$. Coexisting nematics would therefore need to have the same $`\rho _1`$, hence also the same $`\rho _2`$. But then the excess parts of the chemical potentials $`\mu (l)`$, Eq. (22), would also be the same in the two phases. Chemical potential equality then also forces equality of the density distributions $`\rho (l)`$, so that the supposed coexisting nematics turn out to be the same phase.

Using the SNO approximation (A4,A6) we can furthermore get approximate results for the nematic cloud curves, i.e. the phase boundaries at high densities between the single-phase nematic region and any I-N or N-N coexistence regions. The single nematic phase must have the density distribution of the parent, $`\rho (l)=\rho _0P^{(0)}(l)`$, with $`\rho _1=\rho _0`$ since the parent distribution has unit average rod length; its value of $`\rho _2`$ is given by (A6). We need to find the value of the density $`\rho _0`$ where this phase first begins to coexist with a new isotropic or nematic phase as $`\rho _0`$ is lowered. For I-N coexistence, equality of the chemical potentials (A4,A5) implies that the isotropic phase has density distribution

$$\rho ^\mathrm{I}(l)=3lc_2\rho _2\rho _0e^{l\beta }P^{(0)}(l)$$
(A7)
with

$$\beta =c_1(\rho _0\rho _1^\mathrm{I})c_2\rho _2$$
(A8)
For given $`\beta `$ and $`\rho _0`$, $`\rho _1^\mathrm{I}`$ is given by $`\rho _1^\mathrm{I}=𝑑ll\rho ^\mathrm{I}(l)`$ and the onset of I-N phase coexistence can thus be determined simply by solving numerically for the values of $`\beta `$ and $`\rho _0`$ that fulfill Eq. (A8) and give equal osmotic pressure (34) in the two phases.

For N-N coexistence one similarly finds that the new nematic phase has density distribution

$$\rho ^\mathrm{N}(l)=\frac{\rho _2}{\rho _2^\mathrm{N}}\rho _0e^{l\beta }P^{(0)}(l)$$
(A9)
with

$$\beta =c_1(\rho _0\rho _1^\mathrm{N})c_2(\rho _2\rho _2^\mathrm{N})$$
(A10)
and one can again solve numerically, e.g. for $`\beta `$, $`\rho _0`$, $`\rho _0^\mathrm{N}`$ and $`\rho _1^\mathrm{N}`$. (The four conditions are then Eq. (A10), $`\rho _0^\mathrm{N}=𝑑l\rho ^\mathrm{N}(l)`$, $`\rho _1^\mathrm{N}=𝑑ll\rho ^\mathrm{N}(l)`$, and the osmotic pressure equality; $`\rho _2^\mathrm{N}`$ is given by Eq. (A6) in terms of $`\rho _0^\mathrm{N}`$ and $`\rho _1^\mathrm{N}`$.)

We can obtain in this way approximate nematic cloud curves for bidisperse and bimodal parent length distributions; these reproduce qualitatively the features of our numerically exact phase diagrams, including the re-entrance of the I-N phase boundary. Whether I-N or N-N phase separation actually occurs as density is lowered depends on which of the calculated cloud curves is reached first; due to its overestimate of the nematic ordering, it turns out the SNO predicts stable N-N demixing only above a larger threshold value of the length ratio $`q`$ than in the numerically exact phase diagrams. The overall trends are, however, the same as in the numerically exact phase diagrams of Sec. IV B and IV C. Increasing $`q`$ increases the size of the re-entrant region of the I-N phase boundary; the N-N phase boundary also moves to larger densities and becomes stable over a wider range. Increasing the peak width $`\sigma _\mathrm{S}`$, on the other hand, produces the opposite effects. Interestingly, the SNO predicts that the I-N phase boundary remains re-entrant even for the largest values of the peak width, $`\sigma _\mathrm{S}=1`$. We were not able to check this prediction by calculating exact phase diagrams for $`\sigma _\mathrm{S}=1`$, but it does suggest that the re-entrance is mainly due to the overall polydispersity of the parent, rather than the presence of two well separated peaks in the length distribution.

On the N-N phase boundaries that we calculate within the SNO approximation, a critical point occurs, signalled by the vanishing of $`\beta `$, Eq. (A10), and the coincidence of all moments $`\rho _0`$, $`\rho _1`$, $`\rho _2`$ in the two phases. This critical point is not always observable, however, since it can lie in the metastable region inside the I-N phase boundary. This can happen even if other parts of the N-N phase boundary are stable. We find that this situation occurs for a narrow range of $`q`$, for genuinely polydisperse (bimodal rather than bidisperse) distributions. This prediction of the SNO approximation may appear peculiar but must be expected to be a general feature of systems containing more than two types of rods. Our ($`\rho _0,x`$) phase diagrams are then planar cuts through a higher-dimensional phase diagram; tielines beginning in the cut plane will connect to coexisting phases off the plane. For an appropriately chosen plane, all tielines originating from a stable N-N phase boundary in the plane may then point “to one side” of the plane, with none of them having zero length and thus giving a critical point.

Finally, we use the SNO approximation to analyse the phase diagrams for bidisperse and bimodal length distributions in the limit of large rod length ratio $`q`$; this limit would be difficult to access using direct numerical calculations of the phase diagrams that do not employ the SNO assumption. We also expect the SNO approximation to become more accurate for larger $`q`$. This is suggested by the fact that an increase in $`q`$ generally causes the nematic cloud point to shift to higher density; at such higher densities the nematic will be more strongly ordered, and thus the SNO approximation better justified. While such an argument applies directly to the case of the boundary between the nematic region and the I-N region, the same cannot be said for the N-N phase boundary. There the hypothesis of strong ordering is applied to both nematic phases in the SNO approach, and there is no trivial reason why a nematic shadow phase coexisting with a denser nematic parent is necessarily strongly ordered. Numerically, however, we have found this trend to be confirmed.

It turns out that in the limit $`q\mathrm{}`$ the equations that we need to solve to obtain the nematic cloud curves, i.e. the onset of phase separation into I-N or N-N when decreasing the density, become dependent only on the scaling variable $`\xi =xq`$ rather than separately on $`q`$ and the fraction of long rods $`x`$. This simplification occurs because terms that are exponentially small in $`q`$ can be neglected for $`q\mathrm{}`$. Intuitively, $`\xi `$ is for large $`q`$ simply the ratio of the volume fractions of long and short rods. This follows from the fact that the ratio of number densities is $`x/(1x)`$, so that the volume fractions have the ratio

$$\frac{l_2x}{l_1(1x)}=\frac{qx}{1x}$$
(A11)
which tends to $`\xi `$ for $`x=\xi /q0`$. The numerical solution of the resulting simplified system of equations shows that there is always a stable region of N-N demixing for $`q\mathrm{}`$; interestingly, this is true even for bimodal distributions with the largest peak width $`\sigma _\mathrm{S}=1`$ that we consider. As explained in Sec. IV C, this is presumably due to the fact that in the limit $`q\mathrm{}`$ (with $`x=\xi /q`$ and $`\xi `$ held constant) the parent always has two separate peaks, whatever the value of $`\sigma _\mathrm{S}`$.

The N-N phase boundary moves towards larger densities as $`q`$ increases, and in fact becomes remarkably simple as $`q\mathrm{}`$: the N-N critical point moves towards $`\xi =0`$, and the part of the N-N phase boundary at densities below the critical point collapes onto the horizontal axis. A plot of the N-N cloud curve for different $`q`$ with the scaling quantity $`\xi =xq`$ on the vertical axis shows this evolution explicitly (Fig. 16). Note that in the same representation the I-N phase boundary would change rather little for the same range of values of $`q`$, showing that the relative stability of the N-N phase separation (compared to I-N) is enhanced for large $`q`$. This is consistent with the trend observed in the numerically exact phase diagrams on Sec. IV B.