format: text

text:
# Surface phase transitions in polydisperse hard rod fluids

## I Introduction

A molecular mixture with a finite number of components with different physical properties like characteristic dimension, mass or charge is a well controlled system in the sense that all these properties are known exactly through the study of each component in its pure state. Also the relative proportion of each component is a parameter that can be controlled in experiments. The situation is very different in the industrial colloidal systems because during its manufacture process are involved many no controlled factors at microscopic scales. In this case these physical properties vary in a certain range practically in a continuous way. Then, is more appropriate to treat these system as polydisperse in nature with its properties distributed in some range according to some continuous size distribution function. Although in experiments are known that many systems are polydisperse in nature there is not information about the kind of polydisperse distribution function that follows its physical properties. Example of such systems are colloidal suspensions, polymer blends and emulsions. A common phenomena that have been observed experimentally in many polydisperse systems is that the coexisting phases have different particles sizes distributions Evans  and the simulations also corroborate these results for hard spheres Bolhuis . For anisotropic particles experiments shows that the polydispersity makes the phase diagram much more complex including up to four different coexisting phases Koij . Recent experiments show that the isotropic phase of polydisperse hard disk system under certain conditions can be more dense than the nematic phase Beek .

From theoretical point of view the polydisperse systems composed by anisotropic particles have been in general modeled using a free energy functional for a binary mixture when the polydispersity has a clear bimodal distribution. The fundamental approach is then the mapping of polydisperse system onto a bidisperse one.These bidisperse systems show similar transitions found in experiments with the organophilic polydisperse rod-like colloidal particles (boehmite core-particles) including a three phase isotropic-nematic-nematic region Lekk . Nevertheless the inclusion of bidispersity impose two relevant scales (the two characteristic lengths of the particles) that are absent in the unimodal polydisperse systems.

Recent works include the polydispersity via a continuous size distribution function Speranza . This was motivated by the development of the moment theory for polydisperse systems which makes more accessible the theoretical study because the original infinite number of variables that describe the whole system are reduced in this approximation to some moments of the distribution function through a projection to the space extended by them Sollich . This approximation was applied to the study of melting in polydisperse hard spheres Bartlett . Then, any density functional theory which excess part of free energy depends only on a few moments have a natural extension to the polydisperse case. The fundamental measure theory Rosenfeld  is a clear candidate for applying it to these systems. This functional was used in the study of the local size segregation of hard sphere fluid in a presence of a hard wall Pagona .

Although experimentally was obtained the bulk phase diagram for a few polydisperse systems including also the solid phases Pusey  there is no so many works about inhomogeneous polydisperse systems where the inhomogeneities are created by external fields such as a hard wall or about the study of interfaces between the coexisting phases (there is a recent work Baus  about the study of the surface tension and adsorption properties of the planar interface between coexisting fluid phases).

Recently have been published some works about the study of hard rod monodisperse system in a presence of a hard wall or confined in the slit showing a continuous uniaxial-biaxial nematic transition at density bellow the bulk isotropic-nematic transition followed by the complete wetting by nematic at the wall-isotropic interface Roij . Also was studied in these works the isotropic-nematic interface and the capillary nematization of the slit using the Zwanzig model in the Onsager approximation to describe the hard rod system. The simulations Roij  and experiments Kocevar  confirm in the qualitative level all these results.

In general, the calculation of inhomogeneous density profiles for mixtures of hard core anisotropic bodies using density functional theory is a difficult problem for two reasons. One is the numerical problem caused by the presence of orientational degrees of freedom for each component, besides the spatial one. This complicates functional minimization. The second problem is that there are not so well developed and tested density functionals for mixtures as for one-component systems. Fundamental measure functionals (FMF) are good candidates because of its natural extension to mixtures from one-component systems Rosenfeld . But these functionals which are very good at predicting the freezing of a hard spheres Ros-Schm  or the inhomogeneous density profiles induced by external fields Somoza  have some difficulties if we try to extend them in a classical way to a system of anisotropic particles with continuous rotational degrees of freedom Chamoux . The main problem here is that the Fourier transform of the direct correlation function calculated from the natural extension of the FMF to anisotropic particles is isotropic at $`k=0`$, which is unphysical Chamoux . Nevertheless there is another way to construct a fundamental measure theory (FMT) without using the starting point of the scale particle theory. This is the cavities method Tarazona1 . This method refuses to postulate that the functional should depend only on certain weighted densities which are linear convolutions of the density profiles. The main idea behind it is to ensure the dimensional cross-over of the functional for different configurations of cavities Tarazona1 . The final result is a functional which depends on a nonlinear kernel that can not be decoupled in a set of one particle weighed densities. But this functional is impractical when there are rotational degrees of freedom, because the numerical minimization problem turns out to bee extremely hard.

A simplification of the initial problem is to assume that particles are restricted to have a few discrete orientations. Particles with different orientations are treated as different species and then it is completely justified to use the classical version of FMT for mixtures. The easiest way to implement this method is the Zwanzig model Zwanzig , which describes a system of hard parallelepipeds with three mutual perpendicular orientations. An FMF for parallel hard parallelepipeds which fulfills all dimensional cross-overs Cuesta  has been applied to the study of the nature of the order of nematic-smectic A transition M-R . I will use this functional in the Onsager limit to study the continuous polydisperse system. In Ref. Clarke  it has been presented the general formalism for the study of the bulk phase diagram of the polydisperse Zwanzig model in Onsager’s limit. The bulk phase diagram was obtained for the particular choice of the Schulz size distribution function.

I have two purposes in the study of the polydisperse hard rod system. The first is that as I have already mentioned the knowledge of the specific type of polydispersity inherent to a particular system is commonly very poor. Then is important to study how depends the phase diagram topology on the kind of introduced polydispersity. After have been obtained as general as possible a characterization of the phase behaviour for different families of polydispersity distribution functions (I will suppose that all distributions are unimodal) I will study different interfaces: the wall isotropic, the isotropic-nematic and the nematic polydisperse hard rod system confined between two walls. An interesting questions are how polydispersity affect the surface phase transitions found in Roij  and to study the novel surface behaviour inherent only to polydisperse systems.

This paper is organized as follows. In section II is obtained formally the Onsager limit of the Fundamental Measure Functional for polydisperse hard parallelepipeds and the proper equations to calculate the bulk phase diagram (subsection II.1) and the different interfaces (subsection II.2). In section III are deduced some general results for thermodynamic of polydisperse interfaces, specifically the extension of the interface Gibbs-Duhem equation to polydisperse case. In section IV are presented some important remarks about the validity of the thermodynamic limit in finite polydisperse systems. The results are presented in section V divided in different subsections where I study the possible bulk and surface phase diagrams for rapidly (subsection V.1) and slowly (subsection V.2) decaying distribution functions. The last subsections are devoted to the study of different interfaces: wall-isotropic, isotropic-nematic (subsection V.3) and a polydisperse nematic liquid crystal in the slit (subsection V.4). Finally the conclusions are summarized in section VI.

## II Model

Suppose we have a multicomponent mixture of uniaxial parallelepipeds, all of them with the same cross section width $`\sigma `$ but different lengths $`\{L_\nu \}`$ ($`\nu =1,\mathrm{}c`$, where $`c`$ is the number of species). The orientation of each parallelepiped is restricted to one of the three perpendicular axis $`x`$, $`y`$ or $`z`$ as is sketched in figure 1. This is the well known Zwanzig model Zwanzig . As I want to study planar interfaces I will suppose that the density profiles of different species are inhomogeneous only in one spatial direction and I will label them as follows: $`\rho _\mu (L_\nu ,z)`$, where $`\mu =x,y,z`$ labels the orientation of the species with length $`L_\nu `$.

The excess part of free energy density according to the fundamental measure theory for this model is Cuesta

$`\mathrm{\Phi }_{ex}=n_0\mathrm{ln}(1n_3)+{\displaystyle \frac{𝐧_1𝐧_2}{1n_3}}+{\displaystyle \frac{n_{2x}n_{2y}n_{2z}}{(1n_3)^2}},`$ (1)
where the $`n_\alpha `$’s are weighted densities, i.e. linear convolutions of the density profiles with certain weights:

$`n_\alpha ={\displaystyle \underset{\nu =1}{\overset{c}{}}}{\displaystyle \underset{\mu =x,y,z}{}}\rho _\mu (L_\nu ,z)w_\mu ^{(\alpha )}(L_\nu ,z),`$ (2)
where “$``$” means the convolution with respect to the spatial variables. These weights are characteristic functions whose integrals are fundamental measures of the convex particle Rosenfeld ; Cuesta . The nature of the $`w^{(\alpha )}`$’s can be either scalar or vector. The scalars are

$`\omega _\mu ^{(0)}(L_\nu )=\delta _x^\mu (L_\nu )\delta _y^\mu (L_\nu )\delta _z^\mu (L_\nu ),`$ (3)
$`\omega _\mu ^{(3)}(L_\nu )=\theta _x^\mu (L_\nu )\theta _y^\mu (L_\nu )\theta _z^\mu (L_\nu ),`$ (4)
where for example $`\delta _z^z(L_\nu )=\delta (L_\nu /2|z|)`$ and $`\delta _z^y(L_\nu )=\delta (\sigma /2|z|)`$. These examples allow us to recognize the used notation: the subindex labels the spatial coordinate present in the argument of the Dirac delta function. If the superindex and the subindex coincide then the corresponding parallelepiped edge length in the argument of this function is $`L_\nu `$ otherwise is $`\sigma `$. The meaning of the labels for $`\theta _\sigma ^\mu (L_\nu )`$ is the same as in the preceding case, but $`\theta `$ now represents Hevisaide function $`\mathrm{\Theta }(x)`$.

The vector weighs are

$`𝐰_\mu ^{(i)}(L_\nu )=(w_{\mu x}^{(i)}(L_\nu ),w_{\mu y}^{(i)}(L_\nu ),w_{\mu z}^{(i)}(L_\nu )),`$ (5)
with $`i=1,2`$ and

$`w_{\mu x}^{(1)}(L_\nu )`$ $`=`$ $`\theta _x^\mu (L_\nu )\delta _y^\mu (L_\nu )\delta _z^\mu (L_\nu ),`$ (6)
$`w_{\mu x}^{(2)}(L_\nu )`$ $`=`$ $`\delta _x^\mu (L_\nu )\theta _y^\mu (L_\nu )\theta _z^\mu (L_\nu ).`$ (7)
The component $`w_{\mu x}^{(1)}(L_\nu )`$ is obtained from $`\omega _\mu ^{(0)}(L_\nu )`$ substituting in Eq. (3) $`\delta _x^\mu (L_\nu )`$ by $`\theta _x^\mu (L_\nu )`$ and $`w_{\mu x}^{(2)}(L_\nu )`$ is obtained from $`\omega _\mu ^{(3)}(L_\nu )`$ (see Eq. (4)) doing the opposite substitution. The rest of the components can be obtained using the same procedure.

Now, suppose that we select a natural scale $`L_0`$ for the length of the parallelepipeds, for example the mean length $`L_0=_\mu L_\mu x_\mu `$ (where $`x_\mu `$ is the molar fraction of the specie $`\mu `$). Doing the change of variable $`zz/L_0`$ and $`l_\nu =L_\nu /L_0`$ and taking the limit $`c\mathrm{}`$ (infinite number of components, so the mixture contains parallelepipeds with lengths distributed with certain continuous function) we should replace $`_\nu 𝑑l`$ in any of the equations that define the weighted densities, and the density profiles are now functions of $`l`$ and $`z`$: $`\rho (l,z)`$. We should remember that $`z`$ is in reduced units.

We know from Onsager that the correct variable to describe a nematic phase in the limit of infinite aspect ratio is $`\rho L_0^2\sigma `$, so I will use in the following this scaled density for any specie using the same notation: $`\rho _\mu (l,z)L_0^2\sigma \rho _\mu (l,z)`$. Using the expansions of the convolutions $`\rho _\mu (l,z)\delta _z^\mu `$ and $`\rho _\mu (l,z)\theta _z^\mu `$ ($`\mu =x,y`$) in powers of $`\sigma `$:

$`\rho _\mu (l,z)+{\displaystyle \frac{\rho _\mu ^{\prime \prime }(l,z)}{8}}\sigma ^2+\mathrm{}`$ (8)
$`\text{and}\rho _\mu (l,z)\sigma +{\displaystyle \frac{\rho _\mu ^{\prime \prime }(l,z)}{24}}\sigma ^3+\mathrm{}`$ (9)
respectively (primes denote derivatives w.r.t. $`z`$) and expanding the excess part of the free energy density (1) in powers of $`\lambda ^1=\sigma /L_0`$ is obtained

$`\mathrm{\Phi }_{ex}L_0^2\sigma =\mathrm{\Phi }_0+\lambda ^1\mathrm{\Phi }_1+\mathrm{},`$ (10)
with

$`\mathrm{\Phi }_0`$ $`=`$ $`{\displaystyle \underset{\mu }{}}m_\mu \mathrm{\Psi }_\mu ,`$ (11)
$`\mathrm{\Phi }_1`$ $`=`$ $`m(\rho +\mathrm{\Phi }_0)+{\displaystyle \underset{\mu }{}}\mathrm{\Psi }_\mu \mathrm{{\rm Y}}_\mu +{\displaystyle \underset{\mu =x,y,z}{}}\mathrm{{\rm Y}}_\mu .`$ (12)
where $`\rho =_\mu \rho _\mu `$ and $`m=_\mu m_\mu `$ (was dropped the dependence of $`\rho _\mu `$ and $`m_\mu `$ on $`z`$ for notation simplicity). I have used the following definitions

$`\mathrm{\Psi }_\mu ={\displaystyle \underset{\nu \mu }{}}m_\nu ,\mathrm{{\rm Y}}_\mu ={\displaystyle \underset{\nu \mu }{}}\rho _\nu ,\mu ,\nu =x,y,z`$ (13)
where

$`\rho _\mu (z)`$ $`=`$ $`{\displaystyle 𝑑l\rho _\mu (l,z)},\rho _z(z)={\displaystyle 𝑑l\rho _z(l,z)\delta _l},`$ (14)
$`m_\mu (z)`$ $`=`$ $`{\displaystyle 𝑑ll\rho _\mu (l,z)},m_z(z)={\displaystyle 𝑑l\rho _z(l,z)\theta _l},`$
with $`\mu =x,y`$, $`\delta _l=\delta (l/2|z|)`$ and $`\theta _l=\mathrm{\Theta }(l/2|z|)`$. $`\mathrm{\Phi }_0`$ is then the Onsager limit whereas $`\mathrm{\Phi }_1`$ is the first correction. It is interesting to notice that $`\mathrm{\Phi }_0`$ depends only on the first moments of $`\rho _\mu (l,z)`$ while $`\mathrm{\Phi }_1`$ depends on the zeroth moments as well.

Then, our approximation for the total free energy per unit area when the system is inhomogeneous in one direction is just

$`\beta L_0\sigma /A={\displaystyle 𝑑z\left(\mathrm{\Phi }_{id}(z)+\mathrm{\Phi }_0(z)\right)},`$ (16)
where was done the proper scaling of the ideal part $`\mathrm{\Phi }_{id}\mathrm{\Phi }_{id}L_0^2\sigma `$ which has the form

$`\mathrm{\Phi }_{id}`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle 𝑑l\rho _\mu (l,z)\left[\mathrm{ln}\left(\rho _\mu (l,z)\mathrm{\Lambda }_\mu ^3(l)\right)1\right]},`$ (17)
with $`\mathrm{\Lambda }_\mu ^3(l)`$ the cube of the thermal lengths of specie $`\mu `$.

### II.1 Bulk phases

The bulk phase diagram is calculated using the obvious homogeneities condition $`\rho _\mu (l,z)=\rho _\mu (l)`$. The total density of the rods with length $`l`$ is $`\rho (l)=_\mu \rho _\mu (l)`$, and the mean density of the system $`\rho =𝑑l\rho (l)`$. The fraction of rods with length $`l`$ parallel to the $`\mu `$ axis (calculated from $`P_\mu (l)\rho _\mu (l)/\rho (l)`$) is a measure of the orientational order and so the functional $`\mathrm{\Phi }=\mathrm{\Phi }_{id}+\mathrm{\Phi }_0`$ should be minimized with respect to it with the obvious constraint $`_\mu P_\mu (l)=1`$. After this functional minimization is obtained

$`P_\mu (l)={\displaystyle \frac{e^{2\mathrm{\Psi }_\mu l}}{_\nu e^{2\mathrm{\Psi }_\nu l}}}.`$ (18)
As usual, the chemical potential of rods with length $`l`$ is calculated from $`\beta \mu (l)={\displaystyle \frac{\delta \mathrm{\Phi }}{\delta \rho (l)}}`$. Using this definition and (18) is obtained the result

$`\beta \mu (l)=\mathrm{ln}\left({\displaystyle \frac{\rho (l)}{_\mu e^{2\mathrm{\Psi }_\mu l}}}\right).`$ (19)
Finally, the reduced osmotic pressure of the system is calculated from the thermodynamic relation: $`\beta \mathrm{\Pi }=𝑑l\rho (l)\beta \mu (l)\mathrm{\Phi }`$, with the following result

$`\beta \mathrm{\Pi }=\rho +\mathrm{\Phi }_0.`$ (20)
The coexistence between $`q`$ phases is calculated from the equalities between the chemical potentials and pressures of each type of rod for the $`q`$ different phases. Selecting a reference isotropic phase (referred to as the parent phase and labelled with superindex $`0`$) with some fixed probability length distribution function then the set of coexisting equations are

$`\mu ^{(\alpha )}(l)`$ $`=`$ $`\mu ^{(0)}(l),`$ (21)
$`\mathrm{\Pi }^{(\alpha )}`$ $`=`$ $`\mathrm{\Pi }^{(0)},\alpha =1,\mathrm{},q.`$ (22)
The first set of equations with (19) give us the expressions for the densities of each phase

$`\rho ^{(\alpha )}(l)=e^{\beta \mu ^{(0)}(l)}{\displaystyle \underset{\nu }{}}e^{2\mathrm{\Psi }_\nu ^{(\alpha )}l}.`$ (23)
If we have an initial polydisperse parent phase with density distribution $`\rho ^{(0)}(l)`$ then the $`q`$ coexisting densities $`\rho ^{(\alpha )}(l)`$ should obey the conservation law

$`{\displaystyle \underset{\alpha }{}}x_\alpha \rho ^{(\alpha )}(l)=\rho ^{(0)}(l),`$ (24)
where $`x_\alpha `$ is the volume fraction of the $`\alpha `$ phase ($`_\alpha x_\alpha =1`$).

Using (24) we obtain from (19)and (23):

$`\rho ^{(\alpha )}(l)={\displaystyle \frac{\rho ^{(0)}(l)_\mu e^{2\mathrm{\Psi }_\mu ^{(\alpha )}l}}{_\beta x_\beta _\nu e^{2\mathrm{\Psi }_\nu ^{(\beta )}l}}},`$ (25)
and using the definition $`\rho _\mu ^{(\alpha )}(l)=P_\mu ^{(\alpha )}(l)\rho ^{(\alpha )}(l)`$ we have

$`\rho _\mu ^{(\alpha )}(l)={\displaystyle \frac{\rho ^{(0)}(l)e^{2\mathrm{\Psi }_\mu ^{(\alpha )}l}}{_\beta x_\beta _\nu e^{2\mathrm{\Psi }_\nu ^{(\beta )}l}}}.`$ (26)
The definitions of the first moments, Eq. (LABEL:eq5b) with $`\rho _z(l,z)\theta _l=\rho _z(l)l`$, and Eq. (26) allow us to derive a self consistent set of $`3q`$ equations

$`m_\mu ^{(\alpha )}={\displaystyle 𝑑ll\frac{\rho ^{(0)}(l)e^{2\mathrm{\Psi }_\mu ^{(\alpha )}l}}{_\beta x_\beta _\nu e^{2\mathrm{\Psi }_\nu ^{(\beta )}l}}},\mu =x,y,z,`$ (27)
with $`\alpha =1,\mathrm{}q`$. If $`\rho _0`$ is the mean density of the parent phase $`\rho _0𝑑l\rho ^{(0)}(l)`$ with $`p(l)`$ its length distribution function, $`\rho (l)\rho _0p(l)`$, we have $`4q`$ unknowns: the $`3q`$ moments $`m_\mu ^{(\alpha )}`$, the $`q1`$ independent $`x_\alpha `$ and $`\rho _0`$, and $`4q`$ equations to solve: the set (27) of $`3q`$ equations and the set of $`q`$ equations for the equalities of pressures.

Now I will specify the model for a two phase coexistence between an isotropic and a uniaxial nematic phase. For the nematic phase we have $`m_y^{(N)}=m_z^{(N)}m_{}`$ and $`m_x^{(N)}=m_{}^{(N)}`$, and for the isotropic one $`m_x^{(I)}=m_y^{(I)}=m_z^{(I)}\frac{1}{3}m^{(I)}`$, where $`m^{(I)}`$ is the total moment of the isotropic phase. Defining the new variables

$`\tau =m_{}+m_{}{\displaystyle \frac{2}{3}}m^{(I)},s=m_{}m_{},`$ (28)
and using the shorthand notation

$`E(l)=3x+(1x)ϵ(l),ϵ(l)=e^{2\tau l}\left(e^{2sl}+2\right),`$ (29)
where $`x`$ is the fraction of volume occupied by the isotropic phase, we obtain from (27) the equations for $`\tau `$ and $`s`$:

$`s`$ $`=`$ $`\rho _0\varphi _s,\varphi _s={\displaystyle 𝑑llp(l)\frac{e^{2\tau l}\left(e^{2sl}1\right)}{E(l)}},`$ (30)
$`\tau `$ $`=`$ $`\rho _0\varphi _\tau ,\varphi _\tau ={\displaystyle 𝑑llp(l)\frac{e^{2\tau l}\left(e^{2sl}+1\right)2}{E(l)}}.`$ (31)
The quantities $`\rho ^{(I)}`$, $`\rho ^{(N)}`$ (the total densities of the isotropic and nematic phases), and $`m^{(I)}`$, $`m^{(N)}`$ (the total first moments of the isotropic and nematic phases), are functions of $`\tau `$ and $`s`$:

$`\rho ^{(I)}`$ $`=`$ $`\rho _0\varphi _{\rho ^{(I)}},\varphi _{\rho ^{(I)}}=3{\displaystyle 𝑑lp(l)\frac{1}{E(l)}},`$ (32)
$`\rho ^{(N)}`$ $`=`$ $`\rho _0\varphi _{\rho ^{(N)}},\varphi _{\rho ^{(N)}}={\displaystyle 𝑑lp(l)\frac{ϵ(l)}{E(l)}},`$ (33)
$`m^{(I)}`$ $`=`$ $`\rho _0\varphi _{m^{(I)}},\varphi _{m^{(I)}}=3{\displaystyle 𝑑llp(l)\frac{1}{E(l)}},`$ (34)
$`m^{(N)}`$ $`=`$ $`\rho _0\varphi _{m^{(N)}},\varphi _{m^{(N)}}={\displaystyle 𝑑llp(l)\frac{ϵ(l)}{E(l)}}.`$ (35)
The pressures of the isotropic and nematic phases are related with the new variables $`\tau `$ and $`s`$ via

$`\beta \mathrm{\Pi }^{(I)}`$ $`=`$ $`\rho ^{(I)}+{\displaystyle \frac{2}{3}}m_{}^{(I)}{}_{}{}^{2},`$ (36)
$`\beta \mathrm{\Delta }\mathrm{\Pi }`$ $`=`$ $`\mathrm{\Delta }\rho +{\displaystyle \frac{(\tau s)(3\tau +s)}{2}}+{\displaystyle \frac{2m^{(I)}(3\tau s)}{3}},`$
where $`\mathrm{\Delta }\mathrm{\Pi }=\mathrm{\Pi }^{(N)}\mathrm{\Pi }^{(I)}`$ and $`\mathrm{\Delta }\rho =\rho ^{(N)}\rho ^{(I)}`$. To ensure the coexistence between the isotropic and nematic phases $`\mathrm{\Delta }\mathrm{\Pi }`$ should be equal to zero. Then, from (LABEL:eq22) and (32)-(35) is possible to calculate $`\rho _0`$ as a function of $`\tau `$ and $`s`$:

$`\rho _0={\displaystyle \frac{\varphi _{\rho ^{(N)}}\varphi _{\rho ^{(I)}}}{\frac{1}{2}(\varphi _s\varphi _\tau )(3\varphi _\tau +\varphi _s)\frac{2}{3}\varphi _{m^{(I)}}(3\varphi _\tau \varphi _s)}}.`$ (38)
Finally, to find the coexistence between the isotropic and nematic phases we need to solve the closure set of two equations (30) and (31) (using expression (38)) with two unknowns $`\tau `$ and $`s`$. These are the equations used below to calculate numerically, for a fixed fraction of volume of the isotropic phase, $`x`$, the densities of the isotropic and nematic coexisting phases as a function of the polydispersity. I will use the experimental nomenclature to classify the different coexisting curves. If there is an infinitesimal amount of nematic phase coexisting with an isotropic phase ($`x=1`$) the dependence of the isotropic density with polydispersity is called the cloud isotropic curve, whereas the nematic density as function of polydispersity is called the shadow nematic curve. In the other extreme case where an infinitesimal amount of isotropic phase coexists with a nematic phase ($`x=0`$) we have the nematic cloud and the isotropic shadow curves.

### II.2 Interface

For the inhomogenous system in one spatial direction (the $`z`$ direction) was obtained (Eqs. (16-17) the Helmholtz density functional per unit area. Then the grand potential is, by definition,

$`{\displaystyle \frac{\mathrm{\Omega }L_0\sigma }{A}}`$ $`=`$ $`{\displaystyle \frac{L_0\sigma }{A}}`$ (39)
$`+`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle 𝑑z𝑑l\rho _\mu (l,z)\left(v_\mu (l,z)\mu (l)\right)},`$
where $`v_\mu (l,z)`$ is the external potential acting on a rod of species $`\mu `$ and length $`l`$. A hard wall placed at $`z=0`$ can be described with the following external potential

$`v_z(l,z)=\{\begin{array}{cc}\mathrm{}\hfill & z<\frac{l}{2}\hfill \\ 0\hfill & z\frac{l}{2}\hfill \end{array}`$ (40)
and $`v_x(l,z)=v_y(l,z)=0`$ for all $`z0`$. I will have assumed that at infinite distance from the wall the stable phase is the isotropic one, with mean density $`\rho _0`$, length distribution $`p(l)`$ and mean length fixed to unity. In this case the chemical potential for each species is the same ( $`\mu _x(l)=\mu _y(l)=\mu _z(l)=\mu ^{(0)}(l)`$). From (19) we have

$`\beta \mu ^{(0)}(l)=\mathrm{ln}\left({\displaystyle \frac{\rho _0}{3}}p(l)e^{2l\mathrm{\Psi }^{(0)}}\right).`$ (41)
Inserting (40) in (39) and minimizing $`\mathrm{\Omega }`$ with respect to $`\rho _\mu (l,z)`$ we obtain

$`\rho _\mu (l,z)`$ $`=`$ $`{\displaystyle \frac{\rho _0}{3}}p(l)e^{2l\mathrm{\Delta }\mathrm{\Psi }_\mu (z)},\mu =x,y,`$ (42)
$`\rho _z(l,z)`$ $`=`$ $`{\displaystyle \frac{\rho _0}{3}}p(l)e^{2\mathrm{\Delta }\mathrm{\Psi }_z\theta _l(z)}\mathrm{\Theta }(z{\displaystyle \frac{l}{2}}),`$ (43)
with the definition $`\mathrm{\Delta }\mathrm{\Psi }_\mu (z)=\mathrm{\Psi }_\mu (z)\mathrm{\Psi }_\mu ^{(0)}`$ and $`\mathrm{\Psi }_\mu ^{(0)}=2m_0/3`$ where $`m_0=m_\mu (\mathrm{})`$. A self-consistent set of equations for the first moments can be obtained from (42), (43) and (LABEL:eq5b), with the new variables (28) (identifying $`m_{}`$ with $`m_x`$ and $`m_{}`$ with $`m_y`$ which now depend on $`z`$) and

$`\vartheta (z)=2(m_z(z)m_0),`$ (44)
The result is (I drop the explicit dependence on $`z`$ for notational simplicity)

$`\tau `$ $`=`$ $`{\displaystyle \frac{2}{3}}\rho _0\left({\displaystyle 𝑑llp(l)e^{(\tau +\vartheta )l}\mathrm{cosh}(ls)}1\right),`$ (45)
$`s`$ $`=`$ $`{\displaystyle \frac{2}{3}}\rho _0{\displaystyle 𝑑llp(l)e^{(\tau +\vartheta )l}\mathrm{sinh}(sl)},`$ (46)
$`\vartheta `$ $`=`$ $`{\displaystyle \frac{2}{3}}\rho _0\left({\displaystyle 𝑑lp(l)_{\mathrm{max}(\mathrm{z}\mathrm{l}/2,\mathrm{l}/2)}^{z+l/2}𝑑z^{}e^{2\tau \theta _l(z^{})}}1\right).`$
In terms of these variables the equilibrium density profiles, after finding the solution to (45)-(LABEL:eq30), have the form

$`\rho _x`$ $`=`$ $`{\displaystyle \frac{\rho _0}{3}}{\displaystyle 𝑑lp(l)e^{(\tau +\vartheta s)l}},`$ (48)
$`\rho _y`$ $`=`$ $`{\displaystyle \frac{\rho _0}{3}}{\displaystyle 𝑑lp(l)e^{(\tau +\vartheta +s)l}},`$ (49)
$`\rho _z`$ $`=`$ $`{\displaystyle \frac{\rho _0}{3}}{\displaystyle _0^{2z}}𝑑lp(l)e^{2\tau \theta _l}.`$ (50)
From (48) and (49) we can define an “in plane” order parameter as

$`Q_b={\displaystyle \frac{\rho _x\rho _y}{\rho _x+\rho _y}}={\displaystyle \frac{𝑑lp(l)e^{(\vartheta +\tau )l}\mathrm{sinh}(sl)}{𝑑lp(l)e^{(\vartheta +\tau )l}\mathrm{cosh}(sl)}}.`$ (51)
The model can be further simplify by a locality assumption: is approximated any convolution $`f\theta _l(z)`$ by $`f(z)l`$. Equations (45) and (46) then remain unchanged, but (LABEL:eq30) and (50) change to

$`\vartheta `$ $`=`$ $`{\displaystyle \frac{2}{3}}\rho _0\left({\displaystyle _0^{2z}}𝑑llp(l)e^{2\tau l}1\right),`$ (52)
$`\rho _z`$ $`=`$ $`{\displaystyle \frac{\rho _0}{3}}{\displaystyle _0^{2z}}𝑑lp(l)e^{2\tau l}.`$ (53)
The equations (45)-(46), plus either (LABEL:eq30) or (52), are our starting point for the study of the interface phase behaviour of the polydisperse hard rod system in the presence of a hard wall.

In Roij1  it was studied, theoretically and by simulations, the wetting by nematic at the wall-isotropic interface for the one-component hard rod system, as well as the capillary nematization of the same system confined between two walls. For modeling the system the Zwanzig approximation in the Onsager limit was used. The hard wall generates orientational order because is entropically favored the alignment of rods parallel to the wall. This orientational order enhanced by the wall translates into the density profile in the creation of a nematic layer whose thickness increases continuously from zero to infinity when the bulk density is increased from $`\rho _0^{}`$ (the density at which the second order uniaxial-biaxial nematic transition occurs) to $`\rho _w`$ ( the complete wetting density) Roij1 . An interesting question to study is the effect of polydispersity on the surface phase diagram. Of course this effect depends on the length distribution function of the isotropic bulk phase $`p(l)`$, but as I will show below, the surface phase behaviour for any $`p(l)`$ belongs to one out of two qualitatively very different phase diagrams, which are determined by the behaviour of $`p(l)`$ at long lengths.

Whenever the uniaxial-biaxial nematic transition at the wall is a continuous transition we can calculate the transition density $`\rho _0^{}`$ for any $`p(l)`$ as I will show below. The value of the total density and the first moment at contact can also be obtained analytically.

At the wall, the value of $`\vartheta (0)`$ is $`2\rho _0/3`$ as can be easily seen from Eq. (LABEL:eq30). At infinity we have $`\vartheta (\mathrm{})=0`$ because of the boundary conditions. Then we can solve the set of two equations (45) and (46) for any value of $`\vartheta (z)`$ in the interval $`[2\rho _0/3,0]`$.

Calling $`\kappa =\tau +\vartheta `$, the Eqs. (45)-(46) can be put in a more convenient form to make these calculations:

$`s`$ $`=`$ $`{\displaystyle \frac{2}{3}}\rho _0{\displaystyle 𝑑llp(l)e^{\kappa l}\mathrm{sinh}sl},`$ (54)
$`\kappa `$ $`=`$ $`\vartheta +{\displaystyle \frac{2}{3}}\rho _0{\displaystyle 𝑑llp(l)\left(e^{\kappa l}\mathrm{cosh}sl1\right)}`$ (55)
The variable $`s`$ calculated from the equation $`s=f(\kappa ,s)`$, where $`f`$ is the right hand side of Eq. (54), is a kind of surface order parameter, which is zero below the transition density $`\rho _0^{}`$ and is positive above it (I will assume that the symmetry breaking at the wall occurs in the $`x`$ direction).

If the transition is continuous then $`f_s^{}(\kappa ^{},0)=1`$ should be obey at the transition point. Using this fact and combining Eqs. (54) and (55), we obtain one single equation to calculate the bifurcation value $`\kappa ^{}`$,

$`{\displaystyle 𝑑llp(l)e^{\kappa ^{}l}\left(1\kappa ^{}l\right)}=2,`$ (56)
and from it the bifurcation density,

$`\rho _0^{}={\displaystyle \frac{3}{2}}\left[{\displaystyle 𝑑ll^2p(l)e^{\kappa ^{}l}}\right]^1.`$ (57)
Equations (56) and (57) will be used below for the calculations of the second order uniaxial-biaxial nematic transition at the wall for different length distribution functions.

## III Thermodynamics of polydisperse interfaces

Once the grand potential is minimized we can calculate the surface tension

$`\gamma ={\displaystyle \frac{\mathrm{\Omega }}{A}}+\mathrm{\Pi }_0L,`$ (58)
where $`L`$ is the length of interface and $`\mathrm{\Pi }_0`$ the bulk osmotic pressure. To characterize surface transitions (as wetting) is important to calculate the density adsorption coefficient. For the polydisperse case I will show that, apart from it, the other moment adsorption coefficients of the total density distribution function are also relevant. In the excess part of free energy functional of this model only the first moment appears ,so the quantities we are referring to are

$`\mathrm{\Gamma }^{(0)}(\rho _0)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle _0^{\mathrm{}}}𝑑z\left(\rho _\mu (z){\displaystyle \frac{\rho _0}{3}}\right),`$ (59)
$`\mathrm{\Gamma }^{(1)}(\rho _0)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle _0^{\mathrm{}}}𝑑z\left(m_\mu (z)m_0\right).`$ (60)
In a mixture of a finite number of species the interface Gibbs-Duhem equation relates surface tension changes with changes of the chemical potentials of the different components as in

$`0=d\gamma +{\displaystyle \underset{\nu }{}}\mathrm{\Gamma }_\nu d\mu _\nu ,`$ (61)
where $`\mathrm{\Gamma }_\nu `$ is the adsorption coefficient of component $`\nu `$, and $`\mu _\nu `$ its chemical potential. The natural extension of Eq. (61) to the polydisperse case is

$`d\gamma ={\displaystyle 𝑑l\mathrm{\Gamma }(l)𝑑\mu (l)},`$ (62)
where $`\mathrm{\Gamma }(l)=_0^{\mathrm{}}𝑑z\left(\rho (l,z)\rho ^{(0)}(l,\mathrm{})\right)`$. Suppose the system is composed by a continuous mixture of asymmetric particles with one of its characteristic lengths, say $`l`$ (from the complete set of characteristic lengths that define the geometry of the particle), taking different values according to some fixed length distribution function $`p(l)`$, whereas the other lengths are constant. In general, the chemical potential of the species with length “$`l`$” in the bulk isotropic phase with total density $`\rho _0`$ can be expressed, according to the scale particle theory, as

$`\mu (l)=\mathrm{ln}\left(\rho _0p(l)\right)+{\displaystyle \underset{k=0}{\overset{D^{}}{}}}a_k(\rho _0)l^k,`$ (63)
where the first term comes from the ideal part and the second one from the excess part of the free energy density, which in any scale particle theory is a polynomial on $`l`$ of degree $`D^{}`$ with coefficients depending on $`\rho _0`$. $`D^{}`$ is a number between 1 and $`D`$ (the dimension of the system) which depends on the specific geometry of particles. For example, in three dimensions we have $`D^{}=3`$ for hard spheres, $`D^{}=2`$ for uniaxial oblate parallelepipeds, or $`D^{}=1`$ for uniaxial prolate parallelepipeds.

A change in the bulk density $`\rho _0`$ implies a change in the surface tension, and using Eqs. (62) and (63) we can find an expression for the derivative of the surface tension with respect to the bulk density,

$`{\displaystyle \frac{d\gamma }{d\rho _0}}={\displaystyle \frac{\mathrm{\Gamma }^{(0)}}{\rho _0}}{\displaystyle \underset{k=0}{\overset{D^{}}{}}}a_k^{}(\rho _0)\mathrm{\Gamma }^{(k)},`$ (64)
where

$`\mathrm{\Gamma }^{(k)}={\displaystyle _0^{\mathrm{}}}𝑑z\left(M^{(k)}(z)M^{(k)}(\mathrm{})\right)`$ (65)
is the $`k`$-moment adsorption coefficient. The $`k`$-moment is defined as usual: $`M^{(k)}(z)=𝑑ll^k\rho (l,z)`$.

This important result tells us that the slope of the surface tension depends in general on different moment adsorption coefficients of the density distribution function $`\rho (l,z)`$. It is thus possible to find a qualitatively different behaviour of the surface tension as compared to that of the one-component system, for which the slope depends only on the density adsorption coefficient, as we can see after substituting $`\rho (l,z)=\rho (z)\delta (l_0l)`$ in (64):

$`{\displaystyle \frac{d\gamma }{d\rho _0}}=\left[{\displaystyle \frac{1}{\rho _0}}+{\displaystyle \underset{k=0}{\overset{D^{}}{}}}a_k^{}(\rho _0)l_0^k\right]\mathrm{\Gamma }^{(0)}.`$ (66)
For example, as I will show below, fixing the degree of polydispersity we can find some range of $`\rho _0`$’s for which the density adsorption coefficient $`\mathrm{\Gamma }^{(0)}`$ and the slope $`d\gamma /d\rho _0`$ are both negative, in contrast with the mono-disperse case behaviour (see Eq. (66)), for which $`d\gamma /d\rho _0`$ and $`\mathrm{\Gamma }^{(0)}`$ have always opposite signs. This is due to the fact that the first moment adsorption coefficient $`\mathrm{\Gamma }^{(1)}`$ is still positive as a consequence of the higher depletion effect that the wall induces on the large rods. Then, from Eq. (64) we can see that the final sign depends on the modules of $`\mathrm{\Gamma }^{(0)}/\rho _0`$ and $`a_1^{}(\rho _0)\mathrm{\Gamma }^{(1)}`$ ($`D^{}=1`$ and $`a_0(\rho _0)=0`$ for this model).

For the sake of further illustrating the validity of Eq. (64) I will deduce the surface tension as a function of density $`\rho _0`$ near the complete wetting transition, $`\rho _0\rho _w`$. The adsorption coefficients diverge logarithmically near the transition point, according to

$`\mathrm{\Gamma }^{(k)}\nu _k\mathrm{ln}t+\sigma _k,`$ (67)
where $`t=1\rho _0/\rho _w`$. Considering the Zwanzig model for prolate parallelepipeds in the Onsager approximation, for which all $`a_k`$ are zero except $`a_1=4\rho _0/3`$, and inserting (67) into (64) we have

$`{\displaystyle \frac{d\gamma }{d\rho _0}}({\displaystyle \frac{\nu _0}{\rho _0}}+{\displaystyle \frac{4}{3}}\nu _1)\mathrm{ln}t{\displaystyle \frac{\sigma _0}{\rho _0}}{\displaystyle \frac{4}{3}}\sigma _1.`$ (68)
Integrating (68) with respect to $`\rho _0`$ we obtain

$`\mathrm{\Delta }\gamma =\nu _0(t)\sigma _0\mathrm{ln}(1t)+t\left[\stackrel{~}{\nu }_1(\mathrm{ln}t1)+\stackrel{~}{\sigma }_1\right],`$ (69)
where $`\mathrm{\Delta }\gamma =\gamma \gamma _w`$, with $`\gamma _w`$ the surface tension value at the transition point, $`\stackrel{~}{\nu }_1=4\rho _w\nu _1/3`$, $`\stackrel{~}{\sigma }_1=4\rho _w\sigma _1/3`$, and $`(t)=_0^t𝑑x(1x)^1\mathrm{ln}x`$.

I will use Eq. (69) for testing the validity of (64) with the exact surface tension calculations near the wetting transition.

## IV Polydisperse systems with finite number of particles

In real systems the total number of particles although large is a finite number. For one-component systems is completely justify to take the thermodynamic limit but this is not the case for polydisperse systems. The fundamental reason of this difference is that the finite systems have only a few species with characteristic lengths distributed in the tail (far from the mean value) of the size distribution function. The majority of theoretical calculations on polydisperse systems have in common the hypothesis that for any length there are an infinite number of particles (in the thermodynamic sense) distributed according to the so called parent size distribution function $`p(l)`$ (normally infinite ranged). Although this hypothesis is in general not justified the majority of the obtained results reflect qualitatively the experimental behaviour because normally the decaying law of $`p(l)`$ for large $`l`$’s is so rapidly that as is shown here the particles with large lengths don’t affect too much the properties of the different coexisting phases. For low decaying $`p(l)`$’s like the log-normal distribution the species with large lengths have a strong influence in the phase behaviour of the system. Then, the phase diagrams obtained through theoretical calculations using the preceding hypothesis can differ very much from the real ones. For this reason becomes crucial to estimate the number of particles present in the tail of $`p(l)`$ and if is highly enough to justify the thermodynamic limit.

Suppose we have a polydisperse infinite system described with infinite ranged size distribution function $`p(l)`$ and that we take a sample of $`N`$ particles (with $`N`$ very large) which coincide with our real system. Then, there is a maximum extreme value $`L_{\mathrm{max}}`$ for the lengths of particles which obviously depend on $`N`$ and on the degree of polydispersity. I will obtain this relation using the asymptotic properties of sampling distributions (for more information see the reference Cramer ). Arranging the $`N`$ sample values in order of magnitude and considering the $`N_t`$’th value from the top, for $`N_t=1`$ we obtain the extreme value. If $`l`$ denote the length of the $`N_t`$’th specie from the top, the probability element that among the $`N`$ sample values, $`NN_t`$ are less than $`l`$, and $`N_t1`$ are greater than $`l+dl`$, while the remaining value falls between $`l`$ and $`l+dl`$ is:

$`g_{N_t}(l)dl=NC_{N_t1}^{N1}F(l)^{NN_t}\left(1F(l)\right)^{N1}p(l)dl,`$ (70)
where $`F(l)=_0^l𝑑tp(t)`$ and $`C_{\nu 1}^{n1}`$ is the combinatorial number. Introducing the new variable $`\xi =N(1F(l))`$ we have that the distribution function that follow this variable in the asymptotic limit of large $`N`$’s is $`h_{N_t}(\xi )=\xi ^{N_t1}\mathrm{exp}(\xi )/\mathrm{\Gamma }(N_t)`$, where $`\mathrm{\Gamma }(x)`$ is the Gamma function. The variable $`l`$ depends on $`\xi `$ and $`N`$ by the relation

$`\xi =N{\displaystyle _l^{\mathrm{}}}p(t)𝑑t.`$ (71)
If the asymptotic solution of this equation is $`l=f_N(\xi )+𝒪(u(N))`$ (where $`u(N)`$ is the next term in the asymptotic expansion) we can estimate the mean length (the average is taken over an infinite number of samples with fixed $`N`$) of the $`N_t`$’th specie from the top as

$`L_N(N_t,\mathrm{\Delta })={\displaystyle _0^{\mathrm{}}}𝑑\xi h_{N_t}(\xi )f_N(\xi )+𝒪\left(u(N)\right),`$ (72)
where was shown explicitly the dependence of $`L_N`$ on the degree of polydispersity $`\mathrm{\Delta }`$. Fixing the total number $`N`$, the maximun extreme value is obtained as $`L_{\mathrm{max}}=L_N(1,\mathrm{\Delta })`$. Then, the number of particles with lengths in the interval $`[(1r)L_{\mathrm{max}},L_{\mathrm{max}}]`$ (with $`r<1`$) is the solution of the equation

$`rL_{\mathrm{max}}=L_N(N_t,\mathrm{\Delta }),`$ (73)
with respect to $`N_t`$. Here I will use two different families of length distribution functions. All them have in common that the mean rod length is fixed to one ($`𝑑llp(l)=1`$). Is imposed also that $`lim_{l0}p(l)=0`$, i.e. the probability of finding rods with zero length is zero. As a measure of the degree of polydispersity I select the parameter $`\mathrm{\Delta }=\sqrt{\sigma ^21}`$, where $`\sigma `$ is the mean deviation $`\sigma ^2=𝑑ll^2p(l)`$. The last constraint that I impose to $`p(l)`$ is $`lim_{\mathrm{\Delta }0}p(l)=\delta (l1)`$, so is recovered the one-component system in the limit of small polydispersities.

The first one is defined as

$`p(l)`$ $`=`$ $`{\displaystyle \frac{qc_{\nu ,q}^{\nu +1}}{\mathrm{\Gamma }\left[(\nu +1)q^1\right]}}l^\nu \mathrm{exp}\left[(c_{\nu ,q}l)^q\right],`$ (74)
$`c_{\nu ,q}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }\left[(\nu +2)q^1\right]}{\mathrm{\Gamma }\left[(\nu +1)q^1\right]}}.`$
The number $`q1`$ controls the decaying law at large $`l`$’s, and $`\nu 0`$ controls the degree of polydispersity which maximun value

$`\mathrm{\Delta }_{\mathrm{max}}=\sqrt{{\displaystyle \frac{\mathrm{\Gamma }(q^1)\mathrm{\Gamma }(3q^1)}{\mathrm{\Gamma }^2(2q^1)}}1},`$ (75)
is reached at $`\nu =0`$. When $`q=1`$ (74) coincide with the Schulz’s distribution function for which $`\mathrm{\Delta }_{\mathrm{max}}=1`$.

The second family is

$`p(l)=C_{\alpha ,q}\mathrm{exp}[\alpha |\mathrm{ln}(l/l_0)|^q],`$ (76)
where $`C_{\alpha ,q}`$ and $`l_0`$ depend on $`\alpha `$ and $`q`$ and are calculated from the conditions $`𝑑lp(l)=𝑑llp(l)=1`$. Specially for $`q=2`$ we have the log-normal distribution:

$`p(l)={\displaystyle \frac{1+\mathrm{\Delta }^2}{\sqrt{2\pi \mathrm{ln}(1+\mathrm{\Delta }^2)}}}\mathrm{exp}\left\{{\displaystyle \frac{\mathrm{ln}^2\left[(1+\mathrm{\Delta }^2)^{\frac{3}{2}}l\right]}{2\mathrm{ln}(1+\mathrm{\Delta }^2)}}\right\},`$ (77)
where $`0<\mathrm{\Delta }<\mathrm{}`$. Applying the described formalism to the distribution function family (74) we obtain

$`L_N(N_t,\mathrm{\Delta }){\displaystyle \frac{t^{1/q}}{C_{\nu ,q}}}\left[1d_0\left(1+{\displaystyle \frac{\mu }{t}}+(q1)d_1\right)\right],`$ (78)
where were neglected terms of order $`\mathrm{ln}^{1/q3}(N)`$ and I have defined $`\mu =(\nu +1)/q1`$, $`t=\mathrm{ln}[N/\mathrm{\Gamma }(\mu +1)]`$ and

$`d_k={\displaystyle \frac{1}{qt}}\left({\displaystyle \frac{\gamma _k(N_t)}{2^k}}\mu \mathrm{ln}t\right),`$ (79)
with $`\gamma _k(x)`$ the Polygamma functions of order $`k`$.

For the log-normal $`p(l)`$ (see Eq. 77) we obtain

$`L_N(N_t,\mathrm{\Delta }){\displaystyle \frac{\left(1+\mathrm{\Delta }^2\right)^{\lambda 1/2}\mathrm{\Gamma }(N_t1/\lambda )}{\left(4\pi t\right)^{1/(2\lambda )}\mathrm{\Gamma }(N_t)}},`$ (80)
where were neglected terms of order of unity. Now $`t=\mathrm{ln}N`$ and $`\lambda =\sqrt{2t\mathrm{ln}^1(1+\mathrm{\Delta }^2)}`$. Fixing $`N=10^{25}`$ and $`r=4/7`$ I solved the equation (73) for $`N_t`$ using $`L_N`$ from (78) with $`q=1,2`$ and from (80). The results ($`N_t`$ as a function of $`\mathrm{\Delta }`$) are shown in figure 2. As we can see from the results for low decaying distributions the tails contain less number of particles. This behaviour can be understood if we take into account that fixing $`N`$ the probability to find large particles strongly decrease with $`q`$ and then $`L_{\mathrm{max}}`$ also decrease in so way that the number of particles present in the tail (with length $`rL_{\mathrm{max}}`$) is an increasing function of $`q`$. Another interesting feature is that $`N_t`$ decrease with $`\mathrm{\Delta }`$ (except for $`q=1`$ and $`\mathrm{\Delta }<0.3`$) which can be explained by the same reason because $`L_{\mathrm{max}}`$ is an increasing function of $`\mathrm{\Delta }`$. The main conclusion that we can extract from these results are that the thermodynamic limit for the species belonging to the tail of the length distribution function is less justified from low decaying $`p(l)`$’s and for large polydispersities.

Following the usual theoretical works on polydisperse systems I will fix a parent $`p(l)`$ but taking into account that the system contains a finite number of species will be used a truncated $`p(l)`$. The main purpose is to estimate the number of particles present in the tails. Following the already described procedure I fix now $`L_{\mathrm{max}}`$ and solving equation $`L_N(1,\mathrm{\Delta })=L_{\mathrm{max}}`$ for N I find the dependence $`N(\mathrm{\Delta })`$. Obviously each point in the $`\rho \mathrm{\Delta }`$ phase diagram now represent systems with different number of particles. From equation (73) I calculate $`N_t(\mathrm{\Delta })`$ for the same $`r=4/7`$ which is shown in figure 3 for the log-normal distribution with $`L_{\mathrm{max}}=20,35`$ and $`50`$. Although the relative fraction $`N_t/N`$ is an increasing function of $`\mathrm{\Delta }`$ as should be, $`N_t`$ decrease with polydispersity. Taking into account that the shadow phase fill an infinitesimal amount of the total volume, the number of its particles is less in orders of magnitude compared with the cloud phase. Then the fractionation of $`N_t10^3`$ particles from the parent to the shadow phase can change dramatically the mean length. We can see from figure 3 that until $`\mathrm{\Delta }0.6`$ we have $`N_t10^3`$ so the results presented in the next sections for the log-normal distribution function should be qualitatively correct in this range of polydispersity values.

## V Results

The important quantities to characterize the surface behaviour of our hard rod system interacting with hard wall are the $`k`$-moment adsorption coefficients ($`k=0,1`$) and the surface tension (see Eqs. (58), (59) and (60)). This quantities are calculated after minimization of (39) to obtain the equilibrium density profiles $`\rho _\mu (z)`$ (see Eqs. (48)- (50)). However, we can gain insight into the surface phase diagram for this model without carrying out the full minimization as I will show bellow.

Because the hard wall does not allow a rod with length $`l`$ perpendicular to it to be at a distance less than $`l/2`$, we have $`m_z(0)=0`$, or equivalently $`\vartheta (0)=2\rho _0/3`$. Then Eqs. (45) and (46) can be solved to obtain $`\tau (0)`$ and $`s(0)`$, which yield inserted in (48) and (49) $`\rho _x(0)`$ and $`\rho _y(0)`$. Finally, $`\rho _z^c=𝑑l\rho _z(l,l/2)`$, the total contact density of rods perpendicular to the wall, is calculated using the contact value theorem valid for any mixture of particles interacting via hard core. This theorem for a mixture of hard core particles reads

$`\beta \mathrm{\Pi }={\displaystyle \underset{\mu }{}}\rho _\mu (l_\mu /2,l).`$ (81)
The extension of this theorem to the polydisperse Zwanzig model in the Onsager limit is

$`\beta \mathrm{\Pi }={\displaystyle \underset{\mu =x,y,z}{}}{\displaystyle 𝑑l\rho _\mu (l,l_\mu )}=\rho _x(0)+\rho _y(0)+\rho _z^c,`$ (82)
where $`\beta \mathrm{\Pi }=\rho _0(1+2\rho _0/3)`$ is the bulk pressure and $`l_x=l_y=0`$, $`l_z=l/2`$. Finally, it is necessary to know for which value of the bulk isotropic density $`\rho _0`$ there appears a macroscopically thick layer of the nematic phase at the wall-isotropic interface. For complete wetting this value can be calculated from (30) and (31) with $`x=1`$ (the fraction of volume occupied by the coexisting nematic phase is zero).

Having all this information about the surface structure we can determine the possible scenarios of surface phase transitions occurring before the wetting transition, depending on the form of the length distribution function $`p(l)`$. For example, we can distinguish the nature of the surface transition (second or first order). For the former we obtain the exact value of the transition point $`\rho _0^{}`$ using equation (57), and for the latter we can also calculate the exact value if we use the local model, and produce an estimation for the nonlocal model, as we will show later.

I have found that the surface phase behaviour depends on the decaying law at large lengths of the length distribution function $`p(l)`$. I have found two qualitatively different types of surface phase diagrams, one of them coming from $`p(l)`$’s decaying exponentially or faster, and the other from slowly decaying $`p(l)`$’s.

### V.1 Bulk and surface phase diagrams for rapidly decaying distribution functions

I begin the study with $`p(l)`$ from the first family (74). Taking into account that I pretend to study systems with very large but finite number of particles (see the discussion done in the preceding section) the correct way to do this is to select the length distribution function equal to $`p(l)\mathrm{\Theta }(L_{\mathrm{max}}l)`$, where $`L_{\mathrm{max}}`$ is the maximun extreme value for the lengths of particles and $`\mathrm{\Theta }(x)`$ is the Hevisaide function. With the constraints $`_0^{L_{\mathrm{max}}}𝑑ll^kp(l)=1`$ ($`k=0,1`$) the equation (74) keeps the same but substituting $`\mathrm{\Gamma }(x)`$ by the Incomplete Gamma function $`\mathrm{\Gamma }[(c_{\nu ,q}L_{\mathrm{max}})^q,x]`$ where $`c_{\nu ,q}`$ is now the solution of the equation

$`c_{\nu ,q}={\displaystyle \frac{\mathrm{\Gamma }[(c_{\nu ,q}L_{\mathrm{max}})^q,(\nu +2)q^1]}{\mathrm{\Gamma }[(c_{\nu ,q}L_{\mathrm{max}})^q,(\nu +1)q^1]}}.`$ (83)
Using Eqs. (56) and (57) I have calculated the second order uniaxial-biaxial nematic transition (or isotropic-nematic transition in the plane of the wall), i.e. the bulk density value $`\rho _b`$ at which the solution $`m_x(0)m_y(0)`$ bifurcates, using this family of $`p(l)`$’s for several values of $`q`$. When the bulk density $`\rho _0`$ is greater than this value we find a phase with preferential alignment of rods parallel to the $`x`$ axis, i.e. $`\rho _x>\rho _y,\rho _z`$ with $`\rho _y\rho _z`$ due to the presence of the wall. This is the reason why in Roij1  this phase was called biaxial nematic phase. All results using the truncated length distribution function from the family (74) and $`q1`$ with $`L_{\mathrm{max}}>20`$ are indistinguishable from the infinite ranged $`p(l)`$ which shows that for rapidly decaying $`p(l)`$’s the rods with lengths distributed in the tail of $`p(l)`$ don’t affect the phase behaviour of the polydisperse mixture. In Fig. 4 the values $`\rho _b`$ for which this transition occur are plotted as a function of the polydispersity $`\mathrm{\Delta }`$, for different choices of $`q`$ using the infinite ranged $`p(l)`$. Also are plotted the cloud isotropic curves (the values of $`\rho ^{(I)}`$ when $`x=1`$) as the solution of the coexistence equations (30) and (31). For the complete wetting transition at the wall-isotropic interface the cloud isotropic density $`\rho ^{(I)}`$ coincide with the wetting density value $`\rho _w`$. As we can see from the figure, both curves are practically the same for any distribution function $`p(l)`$ belonging to the family (74) with $`q1`$. The general trend is the decrease of the transition densities with polydispersity. This is due to the fact that the alignment of long rods parallel to the wall is entropically favored. The amount of long rods increases upon increasing polydispersity, therefor this entropic effect increases the total first moment parallel to the wall $`m(0)=m_x(0)+m_y(0)`$ reaching the threshold value for which a second order nematic transition appears at the wall.

Figure 4 also shows the relative stability of the biaxial nematic phase with respect to the wetting transition, measured by the magnitude $`1\rho _b/\rho _w`$. As we can see, the stability does not change too much (from 0.18 to 0.21) and in general is an increasing function of $`\mathrm{\Delta }`$. This means that the bulk transition density $`\rho ^{(I)}=\rho _w`$, is more sensitive to the change in polydispersity than the surface biaxial transition density $`\rho _b`$.

### V.2 Bulk and surface phase diagram for slowly decaying distribution functions.

The other selected family of length distribution functions $`p(l)\mathrm{\Theta }(L_{\mathrm{max}}l)`$ with $`p(l)`$ from (76) and with $`C_{\alpha ,q}`$ and $`l_0`$ calculated from the usual constraints (the equality to one of the zero and first moment). There are experimental works with oil-in-water emulsions where the droplet size distribution is well described by a log-normal distribution function McDonald . Also was proposed a simple statistical model to describe the side distribution of nuclei in films which grow by the absorption of atoms which arrive either by impingement or via surface diffusion (example is adsorption of Ag on amorphous carbon Faure ) The resulting length distribution function is log-normal Granqvist . In some polydisperse homopolymer systems in a single solvent was found untypical topology of the cloud-point curves Konigsveld . Was concluded on the basis of experimental evidence that this behaviour is due to high asymmetry in the molecular weight distribution, conclusion that agree with theoretical calculations using the Flory-Huggins theory for polydisperse polymers with log-normal distribution function Solc . As we can see this type of distribution is not rare in polydisperse systems with particles which dimensions are determined by some coalescence process.

In the set of equations that determine the coexistence between the isotropic and nematic polydisperse phases there are involved integrals of the type $`𝑑ll^kp(l)\mathrm{exp}(\beta l)`$, with $`k=0,1`$. To calculate the first moment or the density of the component parallel to the director in the nematic phase we find integrals of this type with $`\beta >0`$. Then for $`p(l)`$ decaying at large distances slower than an exponential these integrals diverge when $`L_{\mathrm{max}}\mathrm{}`$. This means that in the phase diagram the unique stable phase is the nematic at any density. When $`L_{\mathrm{max}}`$ has a finite value there exist a surface phase diagram very different for the preceding one, as I will show bellow.

The bulk phase diagram of the polydisperse hard rod system with the length distribution function family (76) turns out to be more complex, so I will now describe in more details the whole bulk phase diagram and not only the cloud isotropic branch (or wetting transition line), as in the preceding case.

Carrying out the same calculations for coexistence as in the previous case, but now using (76) with $`q=2`$ and $`L_{\mathrm{max}}=35`$, (i.e. the truncated log-normal distribution function), we find the bulk phase diagram plotted in Fig. 5. In this figure the cloud isotropic-shadow nematic ($`x=1`$) and the cloud nematic-shadow isotropic ($`x=0`$) densities are shown as a function of $`\mathrm{\Delta }`$.

As we can see from the figure there is a polydispersity value ($`\mathrm{\Delta }=0.3834`$) for which in the cloud isotropic-shadow nematic curves appear a triple coexistence between an isotropic an two different nematic phases (a triple point). This behaviour is peculiar of distribution functions decaying slower than an exponential for large $`l`$’s and strongly depends on the cutoff $`L_{\mathrm{max}}`$. I have found that the two coexisting nematics have different distribution functions. The nematic with the highest density and first moment ($`N_2`$) has a bimodal first moment distribution function $`l\rho ^{(N)}(l)/\rho _0`$ with the two values of its maxima having the same order of magnitude, whereas the other nematic ($`N_1`$) has a unimodal distribution (see Fig. 6).

In Fig. 7 I have plotted the total moment $`m`$ of the two coexisting isotropic and nematic phases as a function of polydispersity. The most interesting feature is the prominent maximun exhibited by the shadow nematic curve. Its existence will be explained below.

The origin of the second maximum in the function $`l\rho (l)/\rho _0`$ at large polydispersities is very easy to explain if we take into account that from Eq. (35), the total moment of the coexisting nematic phase for $`x=1`$ is proportional to the integral of a term which, for large $`l`$’s, can be approximated by

$`f(l)=l\mathrm{exp}\left(\alpha \mathrm{ln}^2({\displaystyle \frac{l}{l_0}})+4\mathrm{\Delta }ml\right),`$ (84)
where $`\mathrm{\Delta }m=m_{}^{(I)}m_{}^{(N)}>0`$. The minimum value of this expression is reached at $`l^{}`$, the largest solution to $`f^{}(l^{})=0`$, i.e.

$`l^{}={\displaystyle \frac{2\alpha \mathrm{ln}(l^{}/l_0)1}{4\mathrm{\Delta }m}}.`$ (85)
If $`l^{}<L_{\mathrm{max}}`$ we find that the function $`f(l)`$ has a second maximum at $`l=L_{\mathrm{max}}`$, apart the other one near $`l_0`$. From Eq. (85) we conclude that increasing polydispersity (decreasing $`\alpha `$) $`l^{}`$ decreases, so the value $`f(L_{\mathrm{max}})`$ increases as we can see directly from (84). Also from (84) and the condition $`l^{}<L_{\mathrm{max}}`$ it is clear that $`f(L_{\mathrm{max}})`$ is an increasing function of $`\mathrm{\Delta }m`$, thus taking into account that $`m_{}^{(N_2)}<m_{}^{(N_1)}`$ and $`m_{}^{(I_1)}=m_{}^{(I_2)}`$, i.e. the most ordered nematic has less amount of rods perpendicular to the nematic director, we find that for $`N_2`$ the second maxima value is higher. For the less ordered nematic ($`N_1`$), although $`l\rho ^{(N)}(l)/\rho _0`$ can develop a second maximum at $`L_{\mathrm{max}}`$, its value is less than the first maximum by several orders of magnitude. This is the reason why we call the $`N_1`$ distribution unimodal distribution.

The maximun in the shadow nematic curve (total moment $`m`$ as a function of $`\mathrm{\Delta }`$ (see Fig. 7)) as function of $`\mathrm{\Delta }`$ can be explained again by resorting to Eq. (84). From the general cloud isotropic curve behaviour as a function of polydispersity $`\mathrm{\Delta }m(>0)`$ is a decreasing function of $`\mathrm{\Delta }`$. From (84) the value of $`f(l)`$ for large $`l`$’s depends on the two terms in the exponential. For moderate polydispersities (near the triple point) a rapid decrease in $`\alpha `$ determines the increase of $`𝑑lf(l)`$, whereas for large polydispersities, where $`\alpha `$ changes slowly with $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }m`$ keeps on decreasing, the second term dominates the behaviour of this integral, which decreases with $`\mathrm{\Delta }`$. Thus at intermediate values $`m`$ pass trough a maximun.

When $`\mathrm{\Delta }`$ is larger than its triple point value the stable nematic is the bimodal one. The cloud nematic and the shadow isotropic curves ($`x=0`$) have the usual behaviour with polydispersity, in the sense that they show, for any $`\mathrm{\Delta }`$, a single I-N coexistence, so is expected that for a certain $`x`$ between 1 and 0 this triple coexistence will disappear.

From equation (29) we can estimate the value $`x`$ for which the contributions of the two coexisting phases to the integral have the same order at $`l=L_{\mathrm{max}}`$ (the length value for the second maxima in the $`N_2`$ first moment distribution function). I obtain $`x13\mathrm{exp}(4\mathrm{\Delta }mL_{\mathrm{max}}/3)`$. Taking into account that $`L_{\mathrm{max}}=35`$ and that $`\mathrm{\Delta }m`$ is of order of unity we conclude that the range of $`x_N1x`$ for which take place the preceding described behaviour should be of order $`10^{20}`$. Thus the triple coexistence will extend until very small values of $`x_N`$. To find the triple equilibrium I have solved the set of five equations: $`s_\alpha =\rho _0\varphi _{s_\alpha }`$, $`\tau _\alpha =\rho _0\varphi _{\tau _\alpha }`$ (with $`\alpha =N_1,N_2`$) and $`\mathrm{\Delta }\mathrm{\Pi }(\tau _{N_2},s_{N_2})=0`$ (see equations (30), (31) and (LABEL:eq22)) whereas $`\rho _0`$ is calculated from $`\mathrm{\Delta }\mathrm{\Pi }(\tau _{N_1},s_{N_1})=0`$, i.e. from equation (38) with substitution $`NN_1`$. The unknowns are $`s_\alpha `$, $`\tau _\alpha `$ and $`x_N`$. The definition of $`E(l)`$ is now

$`E(l)=3(1x_N)+x_{N_1}ϵ_1(l)+x_{N_2}ϵ_2(l),`$ (86)
where $`x_N=x_{N_1}+x_{N_2}`$ ($`x_{N_\alpha }=x_N\theta _\alpha `$, with $`\theta _1+\theta _2=1`$) is the total fraction of volume occupied by the nematic phases and $`\theta _i`$ measure the relative fraction of the $`N_i`$ phase with respect to the other. In figure 8 is plotted in the $`x_N\mathrm{\Delta }`$ plane the closed area where the triple coexistence was found. Fixing some $`\mathrm{\Delta }`$ in the range $`[0.3834,0.4058]`$ and increasing $`x_N`$ from zero is reached the value for which appear $`IN_1N_2`$ triple coexistence with infinitesimal amount of the $`N_1`$ phase ($`\theta _1=0`$). Changing the relative fraction of the nematic $`N_1`$ with respect to the $`N_2`$ we move through the triple coexistence line which finish when the nematic $`N_2`$ fill an infinitesimal amount of volume ($`\theta _1=1`$).

The first moment values for the three phases with $`x_N`$ and $`\mathrm{\Delta }`$ inside the triple coexisting region are plotted as function of $`\mathrm{\Delta }`$ in figure 9 where the bottom figure is a zoom of the top figure for values of $`m`$ corresponding to the coexisting nematics near the tricritical point. The shaded area has the same meaning as in figure 8, i.e. represent the allowed values of $`m`$ inside the triple coexisting region. As we can see the gap between the values of $`m`$ corresponding to coexisting nematics decrease monotonically with $`\mathrm{\Delta }`$ from the triple point corresponding to $`x_N=0`$ (the filled circles) until the tricritical point for $`x_N10^{16}`$ (the filled square).

In figure 10 are plotted the first moment distribution function $`l\rho (l)/\rho _0`$ of the three coexisting phases for $`\mathrm{\Delta }=0.395`$ and $`x_{N_2}=x_N`$. As we can see the phase $`N_2`$ has the second maxima at $`l<L_{\mathrm{max}}`$.

I also studied the effect that has the change of the maximun extreme value $`L_{\mathrm{max}}`$ on the bulk phase diagram. In figure 11 I plot the nematic shadow curve for $`L_{\mathrm{max}}=35`$, $`30`$, $`25`$ and $`20`$. As we can see the position of the triple point moves to larger $`\mathrm{\Delta }`$’s when $`L_{\mathrm{max}}`$ decrease, and also there is some $`L_{\mathrm{max}}`$ (between 30 and 25) for which disappear the triple coexistence. Is reasonable to predict this behaviour because the second maxima value of the shadow nematic ($`N_2`$) first moment distribution function is an increasing function of $`L_{\mathrm{max}}`$ and $`\mathrm{\Delta }`$. Then if $`L_{\mathrm{max}}`$ decrease the value $`\mathrm{\Delta }`$ needed for keeping the triple coexistence should increase. We conclude also from figure 11 (the bottom figure) that the coexisting densities for $`\mathrm{\Delta }`$’s greater than the triple point value decrease with polydispersity when $`L_{\mathrm{max}}`$ increase.

The effect of the decaying law of the length distribution function for fixed $`L_{\mathrm{max}}`$ on the bulk phase diagram is similar to the preceding behaviour as we can see from figure 12. In this case was changed the parameter $`q`$ in the family (76). Increasing $`q`$ the triple point position moves to larger $`\mathrm{\Delta }`$’s and disappear at some value between 2.25 and 2.5. The second maxima is higher for low decaying length distribution functions and then is necessary less amount of polydispersity to exhibit the $`N_1N_2`$ transition.

Now I describe the surface phase diagram for the log-normal distribution function. Fixing $`L_{\mathrm{max}}=35`$ the values for the surface parameter $`s=m_x(0)m_y(0)`$ as a function of $`\rho _0`$ (the bulk density at infinite distance from the hard wall) obtained as the solution of the equations (45) and (46) for different values of $`\mathrm{\Delta }`$ are plotted in figure 13.

For small values of polydispersity (less than those shown in the figure) the surface phase behaviour is similar to the characteristic behaviour of the fast decaying distribution functions where at some $`\rho _0`$ (in figure this value is determined by the intersection of the curve with the $`\rho `$ axis) appears a biaxial nematic phase through a second order transition. Then the order parameter $`s`$ increase with density until reaching the complete wetting density value (these transition points are marked with black circles).

When $`\mathrm{\Delta }`$ is larger than some value $`\mathrm{\Delta }^{}0.4`$ the graphics of $`s`$ as a function of $`\rho `$ develop an hysteresis branch. This is an indication of a first order transition between two nematic phases at the wall. This behaviour is common for all truncated length distribution function with slow decay at large $`l`$’s and have not analog in the distribution function family with fast decay as we have already seen. In figure 13 are shown three of these hysteresis branches for different $`\mathrm{\Delta }`$’s. The scenarios is the following: when $`\mathrm{\Delta }>\mathrm{\Delta }^{}`$ after the second order transition from isotropic in plane phase ($`m_x(0)=m_y(0)`$) to the nematic in plane phase (the triangles at figure 13 represent the bulk density value for which this transitions occur) at some point occur a first order surface transition at the wall between two nematics (see the squares in figure 13) and then for larger $`\rho `$ the complete wetting transition (the circles in figure 13). The curves in figure 13 for $`\mathrm{\Delta }=0.4112`$ and $`0.4298`$ represent this behaviour. Increasing $`\mathrm{\Delta }`$ the first order nematic-nematic transition density value move to the second order isotropic-nematic transition value and then for some $`\mathrm{\Delta }^c`$ these two transition points coalesce. For $`\mathrm{\Delta }>\mathrm{\Delta }^c`$ the surface first order transition occur between an isotropic and nematic phase and then the complete wetting transition as usual (see in figure 13 the curve for $`\mathrm{\Delta }=0.4705`$).

The first order transition points can be calculated exactly using the local approximation without carrying out the full minimization of the grand potential with respect to the density profiles. I describe now the followed procedure to find them. The excess surface free energy density for the non local model which integral over “z” is the surface tension ($`\gamma =𝑑z\psi (z)`$) can be expressed as a function of $`\mathrm{\Delta }\rho (z)`$ ($`\mathrm{\Delta }\rho (z)=\rho _0\rho (z)`$), $`\theta (z)`$, $`\tau (z)`$ and $`s(z)`$ as

$`\psi =\mathrm{\Delta }\rho +{\displaystyle \frac{s^2\tau ^2}{2}}{\displaystyle \frac{2}{3}}\rho _0(2\tau +\theta )\theta \tau .`$ (87)
The gradient of $`\psi (z)`$ is

$`\psi `$ $`=`$ $`\tau (\theta +{\displaystyle \frac{2}{3}}\rho _0){\displaystyle \frac{2}{3}}\rho _0p(2z)e^{2_0^{2z}𝑑z^{}\tau (z^{})}`$ (88)
$``$ $`{\displaystyle \frac{2}{3}}\rho _0{\displaystyle _0^{2z}}𝑑lp(l)\tau \delta _l^{}e^{2\tau \theta _l},`$
where I have introduced the notation

$`\tau \delta _l^{}=\tau (z{\displaystyle \frac{l}{2}})\tau (z+{\displaystyle \frac{l}{2}}).`$ (89)
The gradient of $`\tau (z)`$ is related to $`\theta (z)`$ by

$`\tau =\theta l^2_c{\displaystyle \frac{1l^2_c(1r)}{1l^2_c^2(1r^2)}};r={\displaystyle \frac{l^2_s}{l^2_c}},`$ (90)
$`l^2_{s,c}={\displaystyle \frac{2}{3}}\rho _0{\displaystyle 𝑑ll^2p(l)e^{(\tau +\theta )l}\left(\begin{array}{c}\mathrm{sinh}(ls)\\ \mathrm{cosh}(ls)\end{array}\right)},`$
where

$`\theta `$ $`=`$ $`{\displaystyle \frac{2}{3}}\rho _0{\displaystyle _0^{\mathrm{}}}𝑑lp(l)e^{2_z^{z+l}𝑑z^{}\tau (z^{})}`$ (91)
$``$ $`{\displaystyle \frac{2}{3}}\rho _0{\displaystyle _0^z}𝑑lp(l)e^{2_{zl}^z𝑑z^{}\tau (z^{})}.`$
For the local model (see eqs. (52) and (53)) the equation (88) transforms to

$`\psi (z)={\displaystyle \frac{2}{3}}\rho _0p(2z)\mathrm{exp}\left(4\tau (z)z\right).`$ (92)
From (92) we conclude that the excess energy density profile $`\psi (z)`$ for the local model is a monotonic decreasing function of z and that $`\psi (0)=0`$. In the range of $`\rho _0`$ values where there are two solutions (I will call them 1 and 2) of the set (45)-(46) (two nematics or one isotropic and the other nematic) the more ordered phase has a higher value of $`\tau (z)`$ near the wall. Then, from (92) we see that the excess energy density profile of the more ordered branch has a smaller slope. Then the transition point (the value of $`\rho _0`$) can be calculated exactly from $`\psi _1(0)=\psi _2(0)`$ which taking into account the relation $`\psi (z)=\mathrm{\Pi }_0\mathrm{\Pi }(z)=0`$ (where $`\mathrm{\Pi }(z)`$ is the local pressure, i.e. the pressure evaluated at $`\rho _\mu (z)`$)) becomes in $`\mathrm{\Pi }_1(0)=\mathrm{\Pi }_2(0)`$. The total integral of $`\psi (z)`$ should be minimized to obtain the equilibrium profile, then the length of the more ordered phase for higher $`\rho _0`$’s is calculated from the intersection between two branches $`\psi _1(z)`$ and $`\psi _2(z)`$. The length of this layer increase continuously from $`l^{}`$ as a function of $`\rho _0`$. The value $`l^{}`$ depends on the form of $`p(l)`$ for small $`l`$’s. For length distribution functions that have a continuous decay to zero when $`l0`$, $`l^{}=0`$ and the transition from less ordered to more ordered density profile near the wall is continuous, but if $`p(l)`$ is zero for $`l<l_0`$, then $`l^{}=l_0`$ and at the transition point coexist two profiles with different adsorption coefficients, one of them including a layer of length $`l^{}`$ of more ordered phase, so the transition becomes first order.

We conclude that although at the transition point there is a “2D” first order transition between two different phases, in “3D” the order of this transition depends on the form of $`p(l)`$ for small $`l`$’s. The squares in figure 13 were calculated using this procedure and the complete surface phase diagram for the log-normal polydisperse system with $`L_{\mathrm{max}}=35`$ is shown in figures 14 and 15. When the polydispersity value is less than the value corresponding to the square in figure 15 we find after the continuous uniaxial-biaxial nematic transition (the transition density is shown with dotted line) the complete wetting transition (the transition density is shown with solid line). When $`\mathrm{\Delta }`$ is between the square and the triangle there is also a 2D first order N-N transition in between. Finally for $`\mathrm{\Delta }`$’s larger than the position of the triangle this 2D first order transition is between isotropic and nematic phases. Just at the triangle coalesce a second and first order transitions.

The solution of the set of equations (45)-(46) at $`z=0`$ and the condition $`\psi _1(0)=\psi _2(0)`$ is equivalent to find the coexistence between two phases in the two dimensional bulk system of polydisperse hard rods ($`\rho _z(0)=0`$ and $`m_z(0)=0`$) with the excess energy density $`\mathrm{\Phi }^{exc}=2\mathrm{\Phi }_{2D}^{exc}`$, where $`\mathrm{\Phi }_{2D}^{exc}`$ is the excess part of the bulk free energy density of the 2D system. The reason why I have found two coexisting nematics in the bulk phase diagram when $`x_N=0`$ only for $`\mathrm{\Delta }=0.3834`$ (the bulk triple point) whereas for $`\mathrm{\Delta }0.4`$ exists at the wall a whole range of first order N-N transition is the following: The equations (45) and (LABEL:eq30) for $`z=0`$ can be put in the form

$`\tau ={\displaystyle \frac{2}{3}}\rho _0\left({\displaystyle 𝑑llq(l)e^{\tau l}\mathrm{cosh}(sl)}1\right),`$ (93)
$`s={\displaystyle \frac{2}{3}}\rho _0{\displaystyle 𝑑llq(l)e^{\tau l}\mathrm{sinh}(sl)},`$ (94)
which form coincide exactly with the set of equations resulting from the equality of chemical potentials between the cloud isotropic and shadow nematic phases for $`2D`$ system with the following variable definitions: $`\tau =(m_x^{(2D)}+m_y^{(2D)})/2m_0^{(2D)}`$, $`s=(m_x^{(2D)}m_y^{(2D)})/2m_0^{(2D)}`$, $`\rho _0^{(2D)}=4\rho _0/3`$ and $`q(l)=p(l)e^{2\rho _0l/3}`$. Then the surface-parent phase from which separate the different coexisting phases have not a log-normal distribution form. If the bulk density $`\rho _0`$ is large enough then the “surface-parent” length distribution function $`q(l)`$ has a second maxima at $`z=L_{\mathrm{max}}`$ and with a bimodal parent distribution function is possible to find a nematic-nematic coexistence.

For the non local model is very difficult to gain insight about the monotonic behaviour of the density profiles as can be seen from the equations (88)-(91). In principle they can be non monotonic functions of $`z`$. Then, the location of the transition points between different surface phases through a first order transition may differ from the local model results.

### V.3 The wall-isotropic and isotropic-nematic interfaces

I have solved the equations (45)-(LABEL:eq30) using the Schultz’s length distribution function for the bulk isotropic phase to obtain the equilibrium wall-isotropic profiles. I used the Gauss-Laguerre quadrature to calculate the integrals over variable $`l`$ using 10 points and using the fast Fourier transform algorithm to calculate the convolutions $`\rho _z\theta _l`$ for each $`l`$ with a grid spacing $`\mathrm{\Delta }z=l_{\mathrm{min}}/10`$ ($`l_{\mathrm{min}}`$ is the smallest quadrature point). I have compared the results with those obtained by increasing the number of points used in the quadrature until 20 and the density profiles and the surface tension values does not change appreciably for $`0<\mathrm{\Delta }<0.55`$.

For polydispersity value $`\mathrm{\Delta }=0.064`$ and bulk isotropic density $`\rho _0=1.02>\rho _b`$ ($`\rho _b`$ is the density of the isotropic-biaxial nematic transition which is shown as a function of $`\mathrm{\Delta }`$ in the inset of figure 16) I obtain the density profiles $`\rho _x(z)`$, $`\rho _y(z)`$ and $`\rho _z(z)`$ shown in figure 16. In the upper figure are plotted the order in plane parameter $`Q_b=(\rho _x\rho _y)/(\rho _x+\rho _y)`$ for different values of $`\rho _0`$. When $`\rho _0\rho _w`$ ($`\rho _w`$ as a function of $`\mathrm{\Delta }`$ is also plotted in the inset of the upper figure) the point $`z^{}`$ where separates the $`\rho _x`$ and $`\rho _y`$ branches go to infinity continuously as should be in the complete wetting transition, with the final profile composed by macroscopic layer of the nematic phase. I have checked that $`\gamma _{WI}=\gamma _{WN}+\gamma _{NI}`$ for $`\rho _0=\rho _w`$ (with $`\gamma _{WI}`$, $`\gamma _{WN}`$ and $`\gamma _{NI}`$ the wall-isotropic, wall-nematic and nematic-isotropic surface tensions) which should be fulfilled at the complete wetting transition. I also have found that $`\mathrm{\Gamma }^{(0)}`$ and $`\mathrm{\Gamma }^{(1)}`$ diverges logarithmically with coefficients $`\nu _i`$ and $`\sigma _i`$ depending on $`\mathrm{\Delta }`$ as we can see from figure 17. The general trend is that the prefactor ($`\nu _0`$) of the logarithmic divergence of the density adsorption coefficient is a decreasing function (in absolute value) of $`\mathrm{\Delta }`$, whereas the prefactor ($`\nu _1`$) of the logarithmic divergence of the first moment adsorption coefficient is an increasing function (in absolute value) of $`\mathrm{\Delta }`$. This result is due to the fact that the density gap between the wetting density $`\rho _w=\rho ^{(I)}`$ and the density of the coexisting nematic phase $`\rho ^{(N)}`$ decrease with polydispersity as we can see from figure 18, whereas the moment gap is an increasing function of $`\mathrm{\Delta }`$ (see figure 18). The surface tension as a function of $`\rho _0`$ is plotted in figure (19) for three different values of $`\mathrm{\Delta }`$. As I have already shown in section III the behaviour of $`\gamma `$ is governed by certain combination of $`\mathrm{\Gamma }^{(0)}`$ and $`\mathrm{\Gamma }^{(1)}`$ (see equation (64)), thus taking into account that these quantities have a logarithmic divergence can be used the approximate analytic equation (LABEL:eq22) to estimate the difference between the wall-isotropic surface tension $`\gamma _{WI}`$ for $`\rho _0<\rho _w`$ and its wetting value at $`\rho _0=\rho _w`$. For comparison I have plotted in figure 19 the exact value of $`\gamma _{WI}`$ and the approximation (LABEL:eq22) for $`\mathrm{\Delta }=0.064`$. If we take into the account that the comparison between them should be restricted to the narrow neighborhood of $`\rho _w`$ is surprising that they are similar until $`t=1\rho /\rho _w0.1`$, and that the maximun position is well estimated. Also is plotted $`\gamma _{WI}`$ for $`\mathrm{\Delta }=0.33`$ and $`0.53`$.

In figure 20 (the upper one) are shown the isotropic-nematic interfaces for $`\rho _0=\rho _w`$ and $`\mathrm{\Delta }=0.577`$. In particular are plotted the profiles $`\rho (l,z)/\rho _0(l)`$ for three different $`l`$’s. This ratio should tend to 1 at the bulk isotropic side of the interface. From figure is clear that the nematic phase has more amount of large rods than the isotropic one. This effect is more dramatic for higher $`\mathrm{\Delta }`$’s. In the bottom figure I show for fixed $`l=1.91`$ the profiles for two different $`\mathrm{\Delta }`$’s. By increasing $`\mathrm{\Delta }`$ the amount of rods with this length changes its segregation behaviour. For $`\mathrm{\Delta }=0.333`$ they are lightly segregated at the isotropic phase whereas for $`\mathrm{\Delta }=0.577`$ they strongly segregate to the nematic phase.

I have also calculated the wall-isotropic interface in the local approximation using the log-normal length distribution function for the bulk isotropic phase solving the equations (45), (46) and (52). The main reason to use the local approximation is that is impossible to select an effective quadrature to estimate the integral over $`l`$ due to the specific form of the $`\rho (l,z)`$ for large $`l`$’s (which can develop a second maxima).

Far from the wetting transition in the region of $`\rho _0`$’s where the slope of $`\gamma _{WI}`$ with respect to $`\rho _0`$ is practically zero we have calculated the density profiles $`\rho (z)`$ and $`m(z)`$ for $`\mathrm{\Delta }=0.4298`$. In figure 21 are given these profiles for $`\rho _0=0.8397`$. As we can see from the figure whereas the density adsorption coefficient $`\mathrm{\Gamma }^{(0)}`$ is negative, the first moment adsorption coefficient $`\mathrm{\Gamma }^{(1)}`$ is positive. As I have already shown (see eq. (64)) the slope of the surface tension depends on certain combination of $`\mathrm{\Gamma }^{(0)}`$ and $`\mathrm{\Gamma }^{(1)}`$. For this particular case the slope is negative although $`\mathrm{\Gamma }^{(0)}<0`$ in contrast with the result obtained using the interface Gibbs-Duhem equation for one component system.

For some $`\mathrm{\Delta }`$’s greater than $`\mathrm{\Delta }^{}`$ (the critical point value for the first order $`N_1N_2`$ transition at the wall) were calculated the adsorption coefficients $`\mathrm{\Gamma }^{(0)}`$ and $`\mathrm{\Gamma }^{(1)}`$ as a function of $`\rho _0`$. The results are plotted in figure 22. For $`\mathrm{\Delta }=0.4298`$ (the solid line in the figure) increasing $`\rho _0`$ occur first a second order $`IN_1`$ transition and then a first order $`N_1N_2`$ transition in the plane of the wall. Was already shown that the surface length of the $`N_2`$ phase increase from zero resulting in a second order surface transition although in $`2D`$ is first order. Finally the length of the $`N_2`$ surface phase diverge continuously to infinity at the complete wetting transition point. For $`\mathrm{\Delta }=0.5104`$ (the dotted line in the figure) was found a first order $`IN`$ transition in the plane of the wall and then the complete wetting transition.

As we can see from this figure these surface transitions although second order are very hard, i.e. there is not a jump in the graphic of $`\mathrm{\Gamma }^{(0,1)}`$ as a function of $`t`$ but its slopes have a strong discontinuity at the transition points which translates in a discontinuity in the second derivative of $`\gamma _{WI}`$ with respect to $`t`$ as we can see from figure 23. In the figure is also plotted the surface tension for $`\mathrm{\Delta }=0.3837`$ corresponding to the bulk triple point. For this value the transition to the biaxial nematic surface phase is second order (see figure 15) resulting in the continuity of $`d\mathrm{\Gamma }^{(0,1)}/dt`$ ($`k=0,1`$) at this point and then the surface tension curvature is also continuous.

The coefficients $`\nu _{0,1}`$ and $`\sigma _{0,1}`$ of the logarithmic divergence of $`\mathrm{\Gamma }^{(0)}`$ and $`\mathrm{\Gamma }^{(1)}`$ for log-normal $`p(l)`$ have the same qualitative behaviour (see figure 24) as comparing to the Schultz $`p(l)`$, i.e. $`\nu _0`$ ($`\nu _1`$) is a increasing (decreasing) function of $`\mathrm{\Delta }`$. As I have already shown this behaviour is due to the monotonic dependence of the differences $`\rho ^{(N)}\rho ^{(I)}`$ (decrease with $`\mathrm{\Delta }`$) and $`m_Nm_I`$ (increase with delta) corresponding to the coexistence between shadow nematic and cloud isotropic phases. As we can see from figure 7 $`m^{(N)}`$ has a strong maximun at $`\mathrm{\Delta }0.5`$ so is possible that $`\nu _1`$ becomes non monotonic function for higher $`\mathrm{\Delta }`$’s.

### V.4 Enhancement of capillary nematization by polydispersity

Finally I have studied the capillary nematization of the slit as a function of polydispersity. In Roij1  was studied this phenomena for the one-component system. Here I pretend to study the role of polydispersity in the possible surface transitions governed by confinement. The calculations was done in the following way. Fixing the distance between two hard walls and fixing the density distribution $`\rho _0(l)`$ of the reservoir (which is isotropic) in chemical equilibrium with the fluid in the slit, I minimize (39), this time with an external potential which include one more wall at $`z=h`$. For small values of $`\rho _0`$ the density profile is composed approximately by two wall-isotropic interfaces, but there is some $`\rho _0^{}`$ (which is less than the complete wetting density $`\rho _w`$ for a single wall-isotropic interface) for which the system overcome a first order transition to a phase in which the slit is completely filled by the nematic phase. This was called capillary nematization in analogy with the usual capillary condensation Roij1 . There is some value of the distance between the walls $`h^{}`$ (the critical point value) for which the transition becomes second order. The natural question that arise is how depends this phenomena with polydispersity. Increasing polydispersity the entropic ordering effect that has the wall on the fluid particles also increase due to the presence of large rods. Then is reasonable to predict that the capillary nematization will be enhanced by polydispersity, i.e. fixing the distance between the walls the density value $`\rho _0`$ for which occur the transition will decrease with polydispersity, and the critical point should move to higher $`h`$’s as was obtained by numerical minimization.

The effect that has the polydispersity on capillary nematization is shown in figure (25), where are plotted the transition values $`h^{}`$-$`t=1\rho /\rho _w`$. For fixed $`t`$ the width $`h`$ increase with polydispersity.

## VI Conclusions

Using simple Zwanzig model in the Onsager limit for a polydisperse hard rod system I have study how depends the phase diagram topology on the selected length distribution function $`p(l)`$. I have found that the relevant property which affect the bulk phase behaviour is the decaying law of $`p(l)`$ at large $`l`$’s. For distribution functions decaying exponentially or more rapidly was found a single coexistence between an isotropic and nematic phase for any polydispersity. The cloud isotropic and the shadow nematic coexisting densities are decreasing function of $`\mathrm{\Delta }`$, whereas the first moments are increasing functions. For distribution functions decaying slower than an exponential always we can find a maximun extreme value for the lengths of particles ($`L_{\mathrm{max}}`$) for which the bulk phase diagram turns out to be more complex. Now there is some $`\mathrm{\Delta }`$ for which the cloud isotropic and shadow nematic curves ($`x_N=0`$) develop a triple point: one isotropic ($`I_1`$) coexist with two different nematics ($`N_1`$ and $`N_2`$), one of them ($`N_2`$) with a bimodal distribution function. The region of triple coexistence extend until very small values of $`x_N`$ which strongly depend on the decaying law and on the maximun extreme value $`L_{\mathrm{max}}`$. For values of $`x_N`$ and $`\mathrm{\Delta }`$ outside the triple coexistence region the shadow nematic curve becomes continuous keeping its usual behaviour (a decreasing function of $`\mathrm{\Delta }`$). Extending in the appropriate way the interface Gibbs-Duhem equation for the thermodynamic of interfaces to the polydisperse system I have found that the slope of the surface tension depends on certain combination of different moment adsorption coefficients which was confirmed by numerical minimization. I have studied also the surface phase diagram of the hard rod system interacting with a hard wall for two representative families of distribution functions. Each of them produce two qualitatively different surface phase diagrams. With the particular choice of the Schultz’s distribution function from the first family were calculated the equilibrium profiles and different adsorption coefficients. Increasing the bulk density value was found some $`\rho _0^{}`$ for which the system overcomes a second order uniaxial-biaxial nematic transition The thickness of the new biaxial nematic layer increase from zero and diverges continuously at the wetting transition density $`\rho _w`$. Both $`\rho _0^{}`$ and $`\rho _w`$ are decreasing functions of polydispersity. For the log-normal distribution function (a particular choice from the second family) the surface phase diagram changes dramatically. The second order $`IN_1`$ transition at the wall is followed by a two dimensional first order transition between two nematics ($`N_1`$ and $`N_2`$). The thickness of the $`N_2`$ layer increase from zero continuously and diverge at the wetting density. Increasing $`\mathrm{\Delta }`$ this $`N_1N_2`$ transition coalesce with the second order $`IN_1`$ transition and for higher values we find a single first order $`IN_2`$ at the wall followed by the complete wetting. For the log-normal distribution function I have obtained using the local approximation that although very hard (the slopes of the adsorption coefficients at the transition points have strong discontinuities) these transitions are second order in nature. Was studied also the segregation phenomena occurring at the isotropic-nematic interfaces. The numerical calculations show a preferential segregation of long rods to the nematic side of the interface being very sensitive this phenomena to the degree of polydispersity. Finally was shown that polydispersity enhance the capillary nematization of the slit.

###### Acknowledgements.

It is a pleasure to thank José A. Cuesta and Carlos Rascón for stimulating discussions. Is gratefully acknowledged the support from the postdoctoral position at Wuppertal University and from the Dirección General de Investigación, Conserjería de Educación de la Comunidad de Madrid. This work is part of Project No. BFM2000-0004 of the Dirección General de Enseñanza Superior.