{'fraction_non_alphanumeric': 0.031739263681061704, 'fraction_numerical': 0.011126064304654293, 'mean_word_length': 4.36148021145878, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simplified Onsager theory for isotropic-nematic phase equilibria of length polydisperse hard rods

Polydispersity is believed to have important effects on the formation of liquid crystal phases in suspensions of rod-like particles. To understand such effects, we analyse the phase behaviour of thin hard rods with length polydispersity. Our treatment is based on a simplified Onsager theory, obtained by truncating the series expansion of the angular dependence of the excluded volume. We describe the model and give the full phase equilibrium equations; these are then solved numerically using the moment free energy method which reduces the problem from one with an infinite number of conserved densities to one with a finite number of effective densities that are moments of the full density distribution. The method yields exactly the onset of nematic ordering. Beyond this, results are approximate but we show that they can be made essentially arbitrarily precise by adding adaptively chosen extra moments, while still avoiding the numerical complications of a direct solution of the full phase equilibrium conditions. We investigate in detail the phase behaviour of systems with three different length distributions: a (unimodal) Schulz distribution, a bidisperse distribution and a bimodal mixture of two Schulz distributions which interpolates between these two cases. A three-phase isotropic-nematic-nematic coexistence region is shown to exist for the bimodal and bidisperse length distributions if the ratio of long and short rod lengths is sufficiently large, but not for the unimodal one. We systematically explore the topology of the phase diagram as a function of the width of the length distribution and of the rod length ratio in the bidisperse and bimodal cases.


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 [1][2][3][4] and biological 2,5-8 systems.
Two main theoretical approaches have been used to explain the I-N phase transition. The Maier-Saupe theory 9 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 10 . The Onsager theory 11 , 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 (θ) of the angles θ 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 (θ). 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 (θ). With two different approximation techniques Lekkerkerker et al. 12 and Lasher 13 obtained more accurate results; the numerically exact solution has also been obtained 12,14 . Isihara 15 , 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 16 , Di Marzio 17 and Alben 18 , 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 16,19 . With a similar lattice model Flory 20 predicted an analogous order-disorder phase transition for semiflexible macromolecules. On the experimental side, Robinson 21 , Hermans 22 and Nakajima et al. 3 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 23,24 which affect the densities of coexisting isotropic and nematic phases and other characteristics of the phase transition 25,26 . These effects can be partially accounted for by adding an attractive potential to the excluded volume interaction 27,28 , which leads to a temperature dependence of e.g. the orientational order parameter.
More importantly for our purposes, polydispersitythe presence of an (effectively) continuous spread of rod lengths and/or diameters -has been shown to have important effects on the phase diagram 29,30 . Onsager 11 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 29,30 and also predicted theoretically, at least within an expansion for narrow length distributions 31,32 . Coexistence of several nematic phases (N-N), occasionally with a third isotropic phase (I-N-N), has also been seen experimentally 29 and been predicted by theory for bi-and tridisperse systems (comprising rods of two and three different lengths, respectively) 12, [33][34][35] .
In a series of papers by Flory et al. [36][37][38][39] , the lengthpolydisperse 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 19 , 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 39 and exponential 38 type. Moscicki et al. applied the same theory to a system with a Gaussian length distribution 40 , 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 [41][42][43] . 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 "P 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 P 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 12,35 , and we find below that there the P 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 P 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 P 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 [44][45][46] (or moment method for short), with its extension to systems with non-conserved densities introduced in the context of the polydisperse Zwanzig model 43 .
This paper is structured as follows. In Sec. II we describe our truncated "P 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 P 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 reentrance 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 0 → 0 at fixed values of l.
The state of a single phase of a system is described by the density distribution ρ(l, Ω), defined so that ρ(l, Ω) dl dΩ/4π gives the density of rods with (normalized) lengths in an interval dl around l, and orientations Ω in a solid angle dΩ. The distribution ρ(l, Ω) is the natural extension for a polydisperse system of the usual orientational distribution P (Ω) used for monodisperse systems of rods 47 . The rod orientation Ω can be parameterized in terms of the angle θ with the nematic axis and an azimuthal angle ϕ; due to the cylindrical symmetry of the nematic phase ρ(l, Ω) is independent of ϕ. Using dΩ = d cos θ dϕ the density distribution as a function of rod length, obtained by integrating over orientations, is therefore Here and below we use the shorthand for the angular integral of an arbitrary function F (θ). The orientational distribution of rods can be factored out from ρ(l, θ) as 43 where P l (θ) represents the probability of finding a rod of given length l in orientation Ω = (θ, ϕ) and is normalized to 1: In the isotropic phase, one has P l (θ) ≡ 1 and ρ(l, θ) = ρ(l).

A. Polydisperse Onsager theory
As mentioned above, the present model is a truncation of the full Onsager theory (see e.g. Refs. 32,47 for reviews), which we now describe. In the Onsager limit D/L 0 → 0, the second order virial approximation becomes exact and gives for the excess free energy in units of Here B 2 is the second virial coefficient which for hard rods is half the excluded volume and given by up to negligible terms of order D/L 0 , with γ the angle between the two rods. Using the cylindrical symmetry of the nematic phase one can perform the integration over ϕ and ϕ ′ in Eq. (5). If we also make the densitiesf and ρ(l, θ) non-dimensional by multiplying by a unit volume (π/4)DL 2 0 we obtaiñ where the angular part, can be expressed as a bilinear expansion in Legendre polynomials 14 P 2n (cos θ) with k 0 = π/4, k 2 = −5π/32 etc. and P 0 (cos θ) = 1, P 2 (cos θ) = 1 2 (3 cos 2 θ − 1). The excess free energy is therefore an infinite sum of moments of the density distribution ρ(l, θ), with weight functions lP 2n (cos θ): B. The P2 Onsager model From Eq. (9) we see that the full Onsager theory does not give a truncatable excess free energy 45 , i.e. one that is just a function of a finite number of moments ρ i of the density distribution ρ(l, θ). To obtain a more manageable theory, to which the efficient moment free energy method [44][45][46] can be applied, we now truncate the excess free energy after the second Legendre polynomial. We will thus consider the "P 2 Onsager model" defined by the excess free energỹ where ρ 1 = dl dθ l ρ(l, θ) and ρ 2 = dl dθ lP 2 (cos θ) ρ(l, θ) are moments of the density distribution with weight functions and the numerical constants are c 1 = (8/π)k 0 = 2, c 2 = −(8/π)k 2 = 5/4. In our units ρ 1 is the rescaled volume fraction φ of rods, ρ 1 = (L 0 /D)φ. We will denote below by ρ 0 the zeroth moment of the density distribution, ρ 0 = dl ρ(l), with weight function w 0 (l, θ) = 1; this is the number density of rods. The ratio ρ 1 /ρ 0 = l then gives the rod length averaged over the normalized length distribution P (l) = ρ(l)/ρ 0 , and ρ 2 /ρ 1 ≤ 1 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 12 for the monodisperse case has investigated the convergence of the densities of coexisting isotropic and nematic phases when the expansion (9) is truncated at higher and higher orders. Already for the truncation after P 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 P 2 Onsager model could come from the fact that its excess free energy (10) is very similar to that of the polydisperse Maier-Saupe model for thermotropics (see e.g. Refs. 32). 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. (10), 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 ρ 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. 32). 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. 32,44,43) Using the decomposition (3) and the relation (1), the total free energy (density) of our model is therefore: withf given by (10). For a given density distribution ρ(l), P l (θ) is obtained by minimization of the free energy with respect of P l (θ), subject to the constraint (4). Introducing corresponding Lagrange multipliers κ(l) one has the condition Solving for P l (θ) gives (here and in the following P 2 ≡ P 2 (cos θ)) P l (θ) = exp(lc 2 ρ 2 P 2 ) dθ exp(lc 2 ρ 2 P 2 ) The orientational distributions are therefore determined only by ρ 2 , which from (12) obeys the self-consistency equation If there are several solutions ρ 2 for a given ρ(l), the one with the smaller free energy (16) is the physical one.

C. Phase coexistence equations
At this point we need the expressions for the chemical potential µ(l) and the (osmotic) pressure Π in order to derive the phase equilibrium conditions. The chemical potential is obtained by functional differentiation of the free energy (16) with respect to ρ(l); the variations of the P l (θ) with ρ(l) need not be considered since P l (θ) is chosen to minimize f . After a little simplification one thus finds The pressure follows from the Gibbs-Duhem relation as Imposing equality of chemical potentials in a set of P coexisting phases, labelled by a = 1 . . . P , we can express the density distributions as where R(l) is a function of l common to all phases. Using Eq. (19), the full density distributions over rod lengths and orientations follow as The function R(l) can be found from Eq. (23) and the requirement of particle conservation: if ρ (0) (l) is the overall or "parent" density distribution and v (a) is the fraction of the system volume occupied by phase a, then This leads to so that, from Eqs. (19,23), the density distributions in the coexisting phases are In principle we thus have the conditions for coexistence between P phases: In each phase, ρ 2 obeys the selfconsistency condition (20) and ρ 1 a similar equation obtained by inserting (27) into (11); for the P phase volume fractions v (a) we have the equality of pressure in all phases (P − 1 equations) plus the normalization a v (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. 45 for details of its properties, and Ref. 43 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 ρ(l) = dθ ρ(l, θ) can be added to the free energy (16) without affecting the phase equilibria, since they merely add constants to the chemical potentials µ(l). We are therefore free to replace f by with r(l) an arbitrary function of l. Guided by the intuition that it is the moment densities ρ 1 and ρ 2 appearing in the free energy that drive phase separation, we now allow violations of the lever rule (25) as long as they do not affect these moments. All other details of the density distribution ρ(l, θ) are then found by minimizing the free energy (28) for the given values of ρ 1 and ρ 2 as defined in Eqs. (11,12). Inserting a Lagrange multiplier λ i for each moment ρ i (i = 1, 2), the minimum value of the free energy is then where ρ 0 = dl ρ(l) is the number density of rods as defined previously. The density distribution at which this minimum free energy is obtained is where the w i (l, θ) are the weight functions (13,14) defining the moments; the λ i are determined implicitly by the requirement that the density distribution (30) gives the correct values for the ρ i . Eq. (29) defines the moment free energy; rather than being a functional of the density distribution ρ(l), i.e. of an infinite number of conserved densities, it is a function of the densities ρ 1 and ρ 2 . To find phase equilibria using the moment free energy, one proceeds as follows.
The "moment chemical potentials" are defined as µ i = ∂f mom /∂ρ i and can be written as 45 µ i = λ i + ∂f /∂ρ i , giving in our case Since ρ 1 is conserved while ρ 2 , which contains the rod orientations, is not, µ 1 has to be identical in a set of coexisting phases while µ 2 has to vanish in all the phases.
The pressure calculated from f mom , also has to be identical in all phases. Finally, the lever rule has to be satisfied for the conserved moment ρ 1 , i.e. ρ 1 . 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 µ 1 and vanishing µ 2 imply from Eq. (30) that the coexisting phases calculated from the moment free energy can be written as This is of exactly the form (24) 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 (33) actually coincides with the one derived from the original free energy of the P 2 Onsager model, Eq. (22). 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 ρ 1 but not, in general, for the whole of ρ(l). However, this final requirement will clearly also be fulfilled if all phases in an exactly calculated phase split of the P 2 Onsager model are of the form (30). 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 ρ (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. (27). Explicitly one gets from Eq. (27), with obvious labels for the phases, ρ I (l, θ) = ρ (0) (l) and ; as claimed, these are both of the form (30) with r(l) = ρ (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 (30). Violations of the lever rule can, however, be minimized by enlarging the family (30) 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 (30) we could in principle use moments defined by any weight function w i (l, θ).
Since the corresponding moments ρ i do not appear in the excess free energy, their moment chemical potentials are The values of the extra Lagrange multipliers must therefore be the same in all phases, and from Eq. (30) we can write the density distributions predicted by the moment method as Comparing with Eq. (27), one sees that the second exponential, in which the extra weight functions appear, ought to provide a good approximation to the denominator of Eq. (26), which is not present in Eq. (30). 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 ρ(l) = const.× ρ (0) (l) is included in the family (30), 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" 43 . 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 ρ tot (l) = a v (a) ρ (a) (l) which in general differs from the parent ρ (0) (l) (though agreeing in the value of the conserved moment ρ 1 ). A convenient measure of this violation of the lever rule is the "log-ratio" ln ρ tot (l)/ρ (0) (l). Since from Eq. (35) 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 ρ tot (l) closer to ρ (0) (l); it extends the family (30) 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 43 . 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 43 and we could easily reach values of δ < 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 P 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. 43, 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 P 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. 43). We will find strong differences compared to the results for the bi-and tridisperse (full) Onsager model. By contrasting with the predictions of the P 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 where the normalized length distribution is of the Schulz form: 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 σ, which from now on we will refer to simply as "polydispersity": The limit z → ∞ therefore corresponds to the monodisperse case, where the Schulz distribution degenerates to a δ-peak. Decreasing z will give increasing polydispersity; we limit ourselves to z ≥ 0 and thus σ ≤ 1 since otherwise the Schulz distribution (36) has a power-law divergence for l → 0. In Fig. 1 we show, for each value of The isotropic cloud curve has the lower density of the two cloud curves and joins the isotropic shadow at σ = 0, where the system becomes monodisperse. The nematic cloud and shadow curves similarly join in the limit. σ, 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 43 , is the crossing between the cloud curves and the corresponding shadow curves; for values of σ 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 ρ 1 rather than the number densities ρ 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 ρ 1 /ρ 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.  . 3. Plot of the average rod lengths in the isotropic (solid) and nematic (dashed) shadow phases. The average rod length in the cloud phases, being identical to that of the parent, is always one. Thus the nematic phase always has a larger average length than its coexisting isotropic phase. Fractionation disappears when the system becomes monodisperse, σ → 0, as expected.

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 35,48 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 P 2 Onsager model. To ascertain which of these applies, we now study the phase diagram of a bidisperse system of rods within the P 2 Onsager model. The normalized parent length distribution is then 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 The polydispersity is given by and σ 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).  . 4. Phase diagram of a bidisperse system for rod length ratio q = 2.5. We only show the limits of the the I-N coexistence region, which are given by the isotropic and nematic cloud curves, and omit the shadow curves. All densities shown are therefore those of the parent phase.
In Fig. 4 we show the phase diagram 49 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 threephase 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 29,35 , it will also be useful to show our phase diagrams in the representation employed in Refs. 29,35, which uses the variables φ 1 = (L 1 /D)φ 1 andφ 2 = (L 1 /D)φ 2 instead of our x and ρ 0 . Here φ 1 and φ 2 are the volume fractions of short and long rods, and L 1 = L 0 l 1 is the unnormalized length of the short rods. With N 1 and N 2 the number of short and long rods, we havẽ The term in brackets is just our dimensionless rod number density ρ 0 , so that 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. 35; we vary l 1 and l 2 according to Eq. (39) to maintain a constant average length in the parent, while in Ref. 35 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 (42,43), the values of theφ 1 andφ 2 remain unaffected. We show theφ 1 ,φ 2 -representation for the case q = 2.5 in Fig. 5. By including Legendre polynomials up to P 14 , and determining the orientational distribution functions P l (θ) numerically, Lekkerkerker et al. 12 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 29 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 (θ), Odijk and Lekkerkerker 33 showed that re-entrance is generically to be expected for large q. Vroege and Lekkerkerker 35 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 ρ 0 , corresponding to horizontal lines in a (ρ 0 ,x) representation of the phase diagram or to lines through the origin in a (φ 1 ,φ 2 ) plot. Van Roij and Mulder 34 later showed that these statements are exactly true (only) in the limit of large ρ 0 ; this means in particular that the N-N coexistence region is not closed off by an N-N critical point. Birshtein et al. 19 , using a trial function of the type used by Onsager 11 , had earlier come to the opposite conclusion. However, they emphasized (as did Abe and Flory 36 ) 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 P 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 (ρ 0 , x andφ 1 ,φ 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 50 . The N-N region is closed off by an N-N critical point, where the moments ρ 1 and ρ 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 condi- tions 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 threephase region shifted to higher q. The only qualitative difference caused by the fact that we are using the truncated P 2 Onsager model is that, in contrast to the full Onsager theory 34 , 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 dis-  . 8. Phase diagram of a bidisperse system with q = 12, in terms of the parent density ρ0 and the long rod number fraction x. As shown in detail in the inset, a three-phase I-N-N region appears, bordered by a region of N-N coexistence; at larger values of x, the nematic phase is re-entrant. tinction between the P 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 ρ 0 where nematic phases will be strongly ordered. In the full Onsager theory, it can be shown 34,35 that the excess free energy density is simplyf = 2ρ 0 for large ρ 0 so that the excluded volume per particle isf /ρ 0 = 2. This simple scaling with ρ 0 arises because the typical rod angles θ are ∼ 1/ρ 0 for large ρ 0 ; the angular part K(θ, θ ′ ) of the excluded volume interaction, as given by Eq. (8), scales linearly with the typical values of θ and θ ′ and is thus also proportional to 1/ρ 0 . The high density phase behaviour thus arises solely from a competition between the two components of the ideal part of the free energy (16), 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 P 2 truncation, on the other hand, K is approximated as K(θ, θ ′ ) ∼ c 1 − c 2 P 2 (cos θ)P 2 (cos θ ′ ) which becomes ∼ a + b[θ 2 + (θ ′ ) 2 ] for small angles, with constants a and b. For large ρ 0 , when the typical values of θ and θ ′ 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 θ → 0 we have ρ 2 = ρ 1 and the excess free energy density becomes 1 2 (c 1 − c 2 )ρ 2 1 . 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 P 2 Onsager model in the limit of large ρ 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(θ, θ ′ ) ∼ a + b[θ 2 + (θ ′ ) 2 ] of the orientational part of the excluded volume for small angles; the correct linear scaling with θ and θ ′ 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 ρ 2 , but these would all saturate to values proportional to ρ 1 at high density, giving again the quadratic scalingf ∼ ρ 2 1 . 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 34 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(θ, θ ′ ) 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 P 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 reentrant 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 Here S(l) denotes the Schulz distribution (36) and l 1 and l 2 are again given by Eq. (39) to constrain the average rod length to unity. Let us call σ S the normalized width (standard deviation) of the Schulz distribution, to distinguish it from the normalized overall standard deviation σ of the parent distribution; the distribution (44) is then characterized by the three parameters q, x and σ S . For σ S → 0, we recover the bidisperse case; on the other hand, as σ 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 σ 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 σ 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 (σ S = 0) we have a three-phase region and re-entrance in the nematic phase. We then increase σ S and monitor the corresponding variation of the (ρ 0 , x) phase diagram. When σ 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 σ 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 FIG. 10. Examples of our bimodal parent length distribution. In the limit where the width of each of the two peaks tends to zero, the bidisperse case is recovered. As the width is increased, on the other hand, the distribution eventually becomes unimodal; in the inset, we show the case where the minimum and the second maximum have just merged into a turning point.
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 σ 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 σ 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 σ 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 σ 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 σ S , on the other hand, do have a significant effect. If the width of the peaks is increased to σ 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 σ 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 ex-  11. Phase diagram for a bimodal length distribution with rod length ratio q = 12 and normalized width of the two peaks of σS = 0.6. Compared to the bidisperse case with the same q (Fig. 8), the three-phase (I-N-N) region as well as the associated N-N region have disappeared. The re-entrance of the nematic phase is still present, but is now visibly less pronounced than in the bidisperse case or for small σS, e.g. σS = 0.1. Towards x = 0 and x = 1 the coexistence region is now also rather broader, since for the limiting values of x we now recover a system with a Schulz distribution of rod lengths, rather than a monodisperse system. ist 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 σ S , up to the maximum value σ 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 σ S = 0, and then find the value of σ 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 σ S above a certain value. This threshold value of σ S is smaller for larger x, showing that in the (ρ 0 ,x) phase diagram the top corner of the three-phase region shrinks towards smaller values of x as σ S is increased. The x-value marking the bottom end of the three-phase region will also depend on σ S . One would suppose that it increases with σ 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 σ S below the curve. Conversely, horizontal sections through the graph give us the extent in x of the three-phase region for given σ 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 σ 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 σ. Eq. (40) showed that for a bidisperse system the overall polydispersity σ 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 σ 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 σ cannot be universally true. In fact, it is easy to see this from Fig. 14. When we increase the width σ S of the two peaks in a bimodal distribution the overall polydispersity σ 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 σ 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 51,52 . We will present our analysis of this case in a separate publication 53 .

V. CONCLUSION
We have analysed the phase equilibria of the P 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 P 2 (cos θ). 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 (ρ 1 and ρ 2 ) of the density distribution -and used the moment method, from which numerical solutions for phase equilibria can be obtained with wellcontrolled accuracy; the onset of nematic ordering from the isotropic side is found exactly by construction of the moment free energy. FIG. 14. Shown is, for a mixture of two Schulz distributions with length ratio q = 12, the value of σS at which we lose the bimodality of the parent (solid), together with the maximum value of σS for which we have a three-phase region for the given value of x (dashed). The parent is bimodal in the region below the solid curve. Note that, from theoretical considerations, the dashed curve has to bend over as x decreases and eventually meet the horizontal axis at some small but positive value of x.
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 P 2 Onsager model, we next studied the bidisperse case within the P 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 P 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 34 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 σ S peaked at different rod lengths; the bidisperse system is recovered for σ S → 0. For small σ S , i.e., small "broadening" of the bidisperse "peaks", the phase diagrams remained essentially as in the bidisperse limit, while larger values of σ 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 σ 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 σ 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 53 ; 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 σ 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 σ of the system is maximal; see Eq. (40). This is also the value of x where the volume fractions occupied by the short and long rods become equal; see Eq. (A10). We have not been able to find a simple phys-ical explanation for this observation.
In summary, our study of the P 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 P 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 P 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 P 2 Onsager model. As discussed above, the theoretical work of van Roij and Mulder 34 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 19 appear to be artefacts of numerical approximations. However, the results of Ref. 34 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 P 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 q ≈ 3.17. It would be worth revisiting the bidisperse Onsager model to confirm this scenario.

APPENDIX A: THE STRONG NEMATIC ORDERING APPROXIMATION
In this appendix we describe an approximation that is useful for analysing the phase diagram of the P 2 Onsager model at high densities. This strong nematic ordering (SNO) approximation becomes exact in the limit ρ 0 → ∞ and allows one to show formally that there can be no N-N coexistence in this limit. For finite ρ 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 θ) that the P 2 Onsager model and the full Onsager theory differ most strongly, so the analysis here focusses on the properties that are particular to the P 2 Onsager model.
The SNO approximation is based on the assumption that all rod angles θ in any nematic phases are reasonably small, so that we can approximate cos 2 θ ≃ 1 − θ 2 and hence From Eq. (19), the orientational distributions are then P l (θ) ∼ exp[− 3 2 lc 2 ρ 2 θ 2 ] and so the SNO assumption of small angles θ corresponds to ρ 2 being large enough to have lρ 2 ≫ 1 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 where the symmetry of the integrand under θ → π − θ has been exploited and corrections that are exponentially small in lρ 2 have been neglected. The chemical potentials (21) thus become µ N (l) ≃ ln ρ(l) + l(c 1 ρ 1 − c 2 ρ 2 ) + ln(3lc 2 ρ 2 ) (A3) for a nematic phase, while for an isotropic phase (ρ 2 = 0) the angular integral (A2) equals unity and the chemical potentials are µ I (l) = ln ρ(l) + lc 1 ρ 1 (A4) Using the SNO assumption (A1) we can also evaluate the integrals appearing in the self-consistency equation (12) for ρ 2 ; the latter can then be solved explicitly to give, after some algebra, In the limit of large density ρ 0 , since ρ 1 = l ρ 0 ∼ ρ 0 , the square root approaches unity and thus ρ 2 ≃ ρ 1 . This is as expected, since the weight functions (13,14) for ρ 1 and ρ 2 only differ by a factor of P 2 (cos θ) which becomes equal to one for maximal nematic order, i.e. for θ → 0. The equality ρ 2 = ρ 1 implies that ρ 2 becomes large as the density increases; for ρ 0 → ∞ 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 (33) becomes Π = ρ 0 + 1 2 (c 1 − c 2 )ρ 2 1 ≃ 1 2 (c 1 − c 2 )ρ 2 1 . Coexisting nematics would therefore need to have the same ρ 1 , hence also the same ρ 2 . But then the excess parts of the chemical potentials µ(l), Eq. (21), would also be the same in the two phases. Chemical potential equality then also forces equality of the density distributions ρ(l), so that the supposed coexisting nematics turn out to be the same phase.
Using the SNO approximation (A3,A5) 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, ρ(l) = ρ 0 P (0) (l), with ρ 1 = ρ 0 since the parent distribution has unit average rod length; its value of ρ 2 is given by (A5). We need to find the value of the density ρ 0 where this phase first begins to coexist with a new isotropic or nematic phase as ρ 0 is lowered. For I-N coexistence, equality of the chemical potentials (A3,A4) implies that the isotropic phase has density distribution For given β and ρ 0 , ρ I 1 is given by ρ I 1 = dl lρ I (l) and the onset of I-N phase coexistence can thus be determined simply by solving numerically for the values of β and ρ 0 that fulfill Eq. (A7) and give equal osmotic pressure (33) in the two phases.
For N-N coexistence one similarly finds that the new nematic phase has density distribution and one can again solve numerically, e.g. for β, ρ 0 , ρ N 0 and ρ N 1 . (The four conditions are then Eq. (A9), ρ N 0 = dl ρ N (l), ρ N 1 = dl lρ N (l), and the osmotic pressure equality; ρ N 2 is given by Eq. (A5) in terms of ρ N 0 and ρ N 1 .) 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 reentrance 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 σ 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, σ S = 1. We were not able to check this prediction by calculating exact phase diagrams for σ S = 1, but it does suggest that the reentrance 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 β, Eq. (A9), and the coincidence of all moments ρ 0 , ρ 1 , ρ 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 (ρ 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 → ∞ 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 ξ = 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 → ∞. Intuitively, ξ 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/(1 − x), so that the volume fractions have the ratio which tends to ξ for x = ξ/q → 0. The numerical solution of the resulting simplified system of equations shows that there is always a stable region of N-N demixing for q → ∞; interestingly, this is true even for bimodal distributions with the largest peak width σ S = 1 that we consider. As explained in Sec. IV C, this is presumably due to the fact that in the limit q → ∞ (with x = ξ/q and ξ held constant) the parent always has two separate peaks, whatever the value of σ S . The N-N phase Phase diagram as derived from the SNO approximation, for a bidisperse system with q → ∞; the scaling variable on the vertical axis is ξ = xq. The straight line near the bottom left corner of the plot indicates the onset of N-N coexistence, which is stable -since it occurs at higher density than I-N phase separation -at sufficiently small ξ. The same is still true if we move to a bimodal distribution, even for the maximal peak widths σS = 1 that we consider.
boundary moves towards larger densities as q increases, and in fact becomes remarkably simple as q → ∞: the N-N critical point moves towards ξ = 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 ξ = 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. The N-N phase boundary for bidisperse rod length distributions with a range of rod length ratios q, as derived from the SNO approximation, with the scaling variable ξ = xq on the vertical axis. The N-N critical point moves towards ξ = 0 for q → ∞, while the N-N boundary as a whole moves towards higher densities in the same limit.