text:
# Phase diagrams of Zwanzig models: The effect of polydispersity

## I Introduction

Although the Zwanzig model was introduced long ago as a very simplified model to study the isotropic-nematic transition in liquid crystals Zwanzig , the determination of its phase diagram, including all inhomogeneous bulk phases, has never been carried out so far. With the help of the scaled-particle theory (equivalent to a $`y3`$ expansion), developed for a variety of hard particle fluids Barboy  including the hard parallelepipeds with restricted orientations which define the Zwanzig model Moore , the phase behavior of the homogeneous bulk phases has been investigated. Also the the isotropic-nematic interface has been accessible to calculations through a smoothed-density approximation Moore , which consists in evaluating the bulk free energy at some smoothed density, following the same recipe of Tarazona for hard spheres Tarazona .

The usual second virial approximation (exact for freely rotating models but only approximate for restricted orientation models) has recently been applied to the Zwanzig model in order to elucidate the bulk phase behavior of a polydisperse mixture of hard rods Clarke . With an inhomogeneous version of this approximation different interfaces of a monodisperse hard rod fluid Roij , as well as its polydisperse counterpart Martinez-Raton , have been studied.

Unlike for hard objects with free orientations—like spherocylinders, which have been extensively simulated Bolhuis1 , there is only one simulation of the Zwanzig model in a square lattice Casey , for parallelepipeds with dimensions $`5\times 1\times D`$ (with $`D=`$1, rods, or 5, plates) Casey . The main reason for this lack of simulation data is that any Monte Carlo movement which includes reorientation of parallelepipeds to any of the three directions has a very low acceptance ratio due to the huge overlap between particles as the length-to-breadth ratio and the density increase.

The second virial approximation of the free energy of the Zwanzig model has also been employed to investigate the phase diagram of symmetric mixtures of rods and plates Van-Roij1 . Stimulated by theoretical calculations made in the early 70’s Alben  which show that a binary mixture of rods and plates can stabilize a biaxial nematic phase (in which the symmetry axes of particles of different types point along mutually perpendicular directions), van Roij and Mulder studied the relative stability of this phase against a nematic-nematic phase separation Van-Roij1 . The authors calculated the phase diagram for different length-to-breadth ratios $`\kappa L/D`$. They showed that for $`\kappa >8.8`$ the biaxial phase is stable in a relative small window above a multicritical point. How the topology of these phase diagrams is modified when a scaled-particle approximation (which includes higher order virial coefficients) is used instead of the second virial is one of the open questions which we will try to answer here. Another one is the influence of polydispersity in the phase behavior of the mono and bidisperse Zwanzig model. This is most relevant since with a recently derived a fundamental measure functional (FMF) for hard parallelepipeds Cuesta  (equivalent to the scaled-particle theory for homogeneous phases), we have shown that polydispersity may enhance the thermodynamic stability of a biaxial nematic phase Martinez-Raton1 .

There are recent experiments on true mixtures of hard rods and plates which show a very rich phase behavior, including triple coexistence between an isotropic and two nematics—each rich in one of the species—as well as inhomogeneous phases, like a columnar phase Kooij . This phase behavior has been quantitatively accounted for using the Parson approximation of the free energy functional of binary mixture of rods and plates in the Onsager limit Wensink . Nevertheless, the long searched biaxial phase has not been found in these experiments. An unavoidable ingredient of experiments is polydispersity. In Ref. Martinez-Raton1  we have shown that this otherwise undesirable element may cause the stabilization of the biaxial nematic phase. It is very important, though, that the system maintains rod-plate symmetry, in contrast with what happens in the existing experiments Kooij . An open question in this work is the topology of the phase diagram at high packing fractions. Addressing this problem requires to study the possible phase transitions to inhomogeneous phases—like smectic, columnar and solid phases. The calculation of phase equilibria in polydisperse liquid crystals when one of the bulk phase is inhomogeneous is thus far a theoretical challenge. There is only one such work, which makes use of an approximate functional for length-polydisperse parallel hard cylinders to make two important predictions: the existence of a terminal polydispersity beyond which the smectic phase is no more stable, and the enhancement of the columnar phase stability Bohle . These results have been confirmed in simulations of freely rotating hard spherocylinders in the Onsager limit Bates .

The paper is organized as follows. In Sec. II.1 we describe both the model and the fundamental measure functional we will use to describe its free energy (A); the equations for phase equilibria between phases of a polydisperse system (B); and the formalism to determine spinodal instabilities in mono, bi or polydisperse systems (C). Section III describes the phase diagrams of the Zwanzig model without (A) and with polydispersity (B), as well as the bidisperse rod-plate model without (C) and with polydispersity (D). Finally we conclude in Sec. IV.

## II Theory

### II.1 Model

Let us consider a length-to-breadth polydisperse mixture of uniaxial oblate and prolate parallelepipeds, with their symmetry axes pointing along one of the three coordinate axes. Let us fix the volume of any particle to 1; thus if $`\lambda =L/D`$ is the length-to-breadth ratio (with $`L`$ the length and $`D`$ the breadth) of a parallelepiped, then

$$L=\lambda ^{2/3},D=\lambda ^{1/3}.$$
(1)
Let us define $`\rho ^\nu (\lambda )`$ to be the density distribution function of the species parallel to the $`\nu `$ ($`=x,y,z`$) axis, and let the total number density to be

$$\rho =_0^{\mathrm{}}𝑑\lambda \rho (\lambda ),\rho (\lambda )=\underset{\nu }{}\rho ^\nu (\lambda ).$$
(2)
The (temperature reduced) free energy density of a multicomponent mixture is given by $`\beta FV^1\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{id}}+\mathrm{\Phi }_{\mathrm{exc}}`$, where $`\mathrm{\Phi }_{\mathrm{id}}`$ is the ideal part, whose exact form is

$$\mathrm{\Phi }_{\mathrm{id}}=\underset{\nu ,i}{}\rho _i^\nu \left(\mathrm{ln}\rho _i^\nu 1\right),$$
(3)
$`i`$ labeling different components, and the excess part, $`\mathrm{\Phi }_{\mathrm{exc}}`$, can be approximated by the FMF for hard parallelepipeds Cuesta , namely

$$\mathrm{\Phi }_{\mathrm{exc}}=n_0\mathrm{ln}(1n_3)+\frac{_\nu n_1^\nu n_2^\nu }{1n_3}+\frac{_\nu n_2^\nu }{(1n_3)^2},$$
(4)
where the weighted densities $`\{n_\alpha \}`$ have the form

$`n_0={\displaystyle \underset{\nu ,i}{}}\rho _i^\nu \rho ,n_1^\nu ={\displaystyle \underset{i}{}}\left[\rho D_i+\rho _i^\nu \left(L_iD_i\right)\right],`$ (5)
$`n_3={\displaystyle \underset{\nu ,i}{}}\rho _i^\nu L_iD_i^2,n_2^\nu ={\displaystyle \underset{i}{}}\left[\rho L_i\rho _i^\nu \left(L_iD_i\right)\right]D_i.`$ (6)
Specializing the above expressions for our polydisperse mixture and using Eq. (1) we obtain

$`\mathrm{\Phi }_{\mathrm{id}}`$ $`=`$ $`{\displaystyle \underset{\nu }{}}{\displaystyle _0^{\mathrm{}}}𝑑\lambda \rho ^\nu (\lambda )\left[\mathrm{ln}\rho ^\nu (\lambda )1\right],`$ (7)
$`\mathrm{\Phi }_{\mathrm{exc}}`$ $`=`$ $`\rho \mathrm{ln}(1\rho )+{\displaystyle \frac{_\nu \xi _{}^\nu \xi _+^\nu }{1\rho }}+{\displaystyle \frac{_\nu \xi _+^\nu }{(1\rho )^2}},`$ (8)
$`\xi _\pm ^\nu `$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}\lambda ^{\pm 1/3}\left[\left(\lambda ^11\right)\rho ^\nu (\lambda )+\rho (\lambda )\right].`$ (9)
The pressure can be obtained from its definition

$`\beta \mathrm{\Pi }`$ $`=`$ $`{\displaystyle \underset{\nu }{}}{\displaystyle _0^{\mathrm{}}}𝑑\lambda \left[\rho ^\nu (\lambda )\mathrm{\Phi }^\nu (\lambda )\right]\mathrm{\Phi },`$ (10)
$`\mathrm{\Phi }^\nu (\lambda )`$ $``$ $`{\displaystyle \frac{\delta \mathrm{\Phi }}{\delta \rho ^\nu (\lambda )}},`$ (11)
which in this particular case yields

$$\beta \mathrm{\Pi }=\frac{\rho }{1\rho }+\frac{_\nu \xi _{}^\nu \xi _+^\nu }{(1\rho )^2}+\frac{2_\nu \xi _+^\nu }{(1\rho )^3}.$$
(12)
### II.2 Phase equilibria between homogeneous phases

To obtain the phase equilibria we follow the general procedure already reported elsewhere Clarke ; Sollich . Suppose that among the $`n`$ coexisting phases there are $`n_I`$ isotropic, $`n_N`$ nematic and $`n_B`$ biaxial phases ($`n=n_I+n_N+n_B`$). The global density distribution (the parent distribution) is fixed to be

$$P(\lambda )=P_0h(\lambda ),𝑑\lambda h(\lambda )=1,$$
(13)
so total mass conservation of each species is expressed by the lever rule

$$P(\lambda )=\underset{a=1}{\overset{n}{}}\gamma _a\rho _a(\lambda ),$$
(14)
with $`\rho _a(\lambda )`$ the total density distribution of phase $`a`$ and $`\gamma _a`$ the fraction of the total volume it occupies (there is the obvious constraint $`_a\gamma _a=1`$). Minimizing the free energy density $`\mathrm{\Phi }`$ with respect to the fraction of particles with length-to-breadth ratio $`\lambda `$ oriented along the $`\nu `$ axis, i.e. $`p_a^\nu (\lambda )\rho _a^\nu (\lambda )/\rho _a(\lambda )`$, and using Eq. (14), the equality of the chemical potentials of each species in different phases

$$\beta \mu _a(\lambda )=\underset{\nu }{}\rho _a^\nu (\lambda )\mathrm{\Phi }_a^\nu (\lambda ),$$
(15)
leads to the following expressions for the coexisting densities

$$\rho _a^\nu (\lambda )=P_0h(\lambda )\frac{e^{\mathrm{\Phi }_a^\nu (\lambda )}}{_b\gamma _b_\tau e^{\mathrm{\Phi }_b^\tau (\lambda )}}$$
(16)
The parent distribution we are going to use throughout the paper is characterized by

$`h(\lambda )=\lambda ^1\left[\zeta f(\lambda /\kappa )+(1\zeta )f(\lambda \kappa )\right],`$ (17)
$`f(z)=K_0(\alpha )^1\mathrm{exp}[(\alpha /2)\left(z^2+z^2\right)],`$ (18)
where $`K_\nu (\alpha )`$ $`(\alpha >0)`$ is the $`\nu `$th-order modified Bessel function, and $`\kappa >1`$. This choice is motivated by its rod-plate symmetry for $`\zeta =1/2`$. The best way to appreciate this symmetry is changing the variable to $`\mathrm{ln}\lambda `$: $`\stackrel{~}{h}(\mathrm{ln}\lambda )=\lambda h(\lambda )`$ has two identical humps (wider the smaller $`\alpha `$) centered at $`\mathrm{ln}\kappa `$ (rods of “typical” aspect ratio $`\kappa `$) and $`\mathrm{ln}\kappa `$ (plates of typical aspect ratio $`\kappa ^1`$). The parameter $`0\zeta 1`$ allows one to tune the overall composition of the mixture, since the molar fraction of the rods is given by $`x_r(\zeta )=_1^{\mathrm{}}h(\lambda )𝑑\lambda `$ \[and that of plates by $`x_p(\zeta )=1x_r(\zeta )`$, of course\]. Thus we can select polydisperse rods (plates) by setting $`\zeta =1`$ ($`\zeta =0`$). The moments of this distribution are given by $`\lambda ^m=K_{m/2}(\alpha )K_0(\alpha )^1[\zeta \kappa ^m+(1\zeta )\kappa ^m]`$, explicitly showing the symmetry of the mixture. A quantitative characterization of the polydispersity can be given if we determine the dispersion in $`L`$ and $`D`$ as obtained from $`h(\lambda )`$ for $`\zeta =0`$ or 1. This yields

$$\mathrm{\Delta }_{L,D}\sqrt{\frac{\lambda ^{2\nu /3}}{\lambda ^{\nu /3}^2}1}=\sqrt{\frac{K_{\nu /3}(\alpha )K_0(\alpha )}{K_{\nu /6}(\alpha )^2}1},$$
(19)
where $`\nu =2`$ for $`\mathrm{\Delta }_L`$ and $`\nu =1`$ for $`\mathrm{\Delta }_D`$.

The number of independent moments in the set $`\{\rho ,\xi _\pm ^\nu \}`$ which completely determine each coexisting phase is 3, 5 and 7 for the isotropic, nematic and biaxial phases respectively (which amounts to a total of $`3n_I+5n_N+7n_B`$ unknowns). They can be obtained through the definitions (2) and (9), together with the distribution functions (16). The remaining unknowns—the $`n1`$ independent $`\gamma _\alpha `$’s—are calculated from the equality of pressures in every phase. This leaves the global dilution, $`P_0`$, as the control parameter. Alternatively, we can fix an external pressure, $`\mathrm{\Pi }_0`$, and eliminate $`P_0`$ in terms of this new control parameter. As $`P_0`$ increases, the fractions of volume of each phase change; thus, for practical purposes, it is computationally simpler to use one of the $`\gamma _\alpha `$’s as control parameter and obtain $`P_0`$ as a function of it.

A particularly important case is the two phase coexistence between a phase (say $`\beta `$) which fills the whole volume (cloud phase) and an incipient new phase (say $`\sigma `$) which fills an infinitesimally amount of volume (shadow phase). In a cloud-shadow coexistence the parent $`P(\lambda )`$ coincides with the distribution function of the cloud phase. Then $`\rho `$ and $`_\nu \xi _\pm ^\nu `$ are fixed for this phase, so the number of unknowns reduces by 3 (the number of constraints).

### II.3 Spinodal instabilities with respect to inhomogeneous phases

In order to be as general as possible in developing the formalism let us consider the problem of finding the conditions for the stability of an arbitrary multicomponent system against spatial modulations of the densities $`\rho _i`$ ($`i`$ label the species). This way we can replace at the end $`_i_\nu `$ for the monodisperse Zwanzig model and $`_i_\nu 𝑑\lambda `$ for the polydisperse one, to obtain the respective spinodal instabilities. The reason for so proceeding is that, as we will show later, the spinodal equations can be expressed in two different ways, each of which suits one of the two cases.

The condition of equality of chemical potentials between two phases is equivalent to the minimization of the grand potential

$$\mathrm{\Omega }=𝑑𝐫\left\{k_\mathrm{B}T\mathrm{\Phi }(\{n_\alpha \})\underset{i}{}\mu _i\rho _i(𝐫)\right\},$$
(20)
at fixed chemical potentials of each species $`\mu _i`$, with respect to the density profile of each species $`\rho _i(𝐫)`$. We specify this grand potential for a general model whose free energy density $`\mathrm{\Phi }`$ has an excess part depending on the densities only through certain weighted densities

$$n_\alpha (𝐫)=\underset{i}{}\left[\rho _i\omega _\alpha ^i\right](𝐫),$$
(21)
as it is the case of Rosenfeld’s fundamental measure theory Rosenfeld . In the definition of $`n_\alpha `$ the symbol “$``$” stands for the convolution of two functions, i.e. $`fg(𝐫)𝑑𝐫^{}f(𝐫^{})g(𝐫𝐫^{})`$.

We are concerned with the case when one of the coexisting phases is homogeneous and the other one is inhomogeneous, with its density profile consisting in a small perturbation around the homogeneous phase with density distribution $`\rho _i`$, i.e.

$$\rho _i(𝐫)=\rho _i\left[1+ϵ_i(𝐫)\right],ϵ_i(𝐫)1.$$
(22)
The result of minimizing (20) with respect to $`\rho _i(𝐫)`$ can be cast as

$$\rho _i(𝐫)=\rho _i\mathrm{exp}\left\{\underset{\alpha }{}\left[\mathrm{\Delta }\varphi _\alpha \omega _\alpha ^i\right](𝐫)\right\},$$
(23)
where $`\mathrm{\Delta }\varphi _\alpha (𝐫)=\varphi _\alpha (𝐫)\varphi _\alpha `$ and $`\varphi _\alpha =\varphi /n_\alpha `$. We are implicitly using magnitudes without spatial variable arguments to denote homogeneous phase quantities. Expanding $`\mathrm{\Delta }\varphi _\alpha (𝐫)`$ to first order in $`ϵ_i(𝐫)`$ yields

$$\mathrm{\Delta }\varphi _\alpha (𝐫)=\underset{\beta }{}\varphi _{\alpha \beta }\underset{j}{}\rho _j\left[ϵ_j\omega _\beta ^j\right](𝐫)+\mathrm{},$$
(24)
where $`\varphi _{\alpha \beta }=^2\mathrm{\Phi }/(n_\alpha n_\beta )`$. Inserting (24) in (23) and expanding the exponential again to first order in $`ϵ_i(𝐫)`$ we arrive at the integral equation

$$ϵ_i(𝐫)=\underset{\alpha ,\beta }{}\varphi _{\alpha \beta }\underset{j}{}\rho _j\left[ϵ_j\omega _\alpha ^i\omega _\beta ^j\right](𝐫),$$
(25)
which is handier to write in Fourier space

$$\epsilon _i(𝐪)+\underset{\alpha ,\beta }{}\varphi _{\alpha \beta }\mathrm{\Omega }_\alpha ^i(𝐪)\underset{j}{}\rho _j\epsilon _j(𝐪)\mathrm{\Omega }_\beta ^j(𝐪)=0,$$
(26)
where $`\epsilon `$ and $`\mathrm{\Omega }`$ are respectively the Fourier transforms of $`ϵ`$ and $`\omega `$ \[with the usual definition $`f(𝐪)=𝑑𝐫e^{i\mathrm{𝐪𝐫}}f(𝐫)`$\].

Equation (26) gets clearer if written in matrix form. Defining the functions

$$m_{ij}(𝐪)\underset{\alpha ,\beta }{}\varphi _{\alpha \beta }\mathrm{\Omega }_\alpha ^i(𝐪)\mathrm{\Omega }_\beta ^j(𝐪),$$
(27)
the matrices $`\widehat{\mathrm{M}}(𝐪)\left(m_{ij}(𝐪)\right)`$ and $`\widehat{\mathrm{P}}\left(\rho _j\delta _{ij}\right)`$, and the vector $`𝜺(𝐪)\left(\epsilon _i(𝐪)\right)`$, Eq. (26) becomes

$$\left(\mathrm{I}+\widehat{\mathrm{M}}(𝐪)\widehat{\mathrm{P}}\right)𝜺(𝐪)=0,$$
(28)
where I is the identity matrix. This linear system has nontrivial solutions if and only if

$$D(\rho ,𝐪)\mathrm{det}\left[\mathrm{I}+\widehat{\mathrm{M}}(𝐪)\widehat{\mathrm{P}}\right]=0.$$
(29)
This equation is suitable for the cases in which there is a small number of components, like the monodisperse Zwanzig model with only three components: $`\nu =x,y,z`$.

Alternatively we can multiply (26) by $`\rho _i\mathrm{\Omega }_\gamma ^i(𝐪)`$, sum over $`i`$ and define the functions

$$u_\gamma (𝐪)\underset{i}{}\rho _i\epsilon _i(𝐪)\mathrm{\Omega }_\gamma ^i(𝐪);$$
(30)
then Eq. (26) becomes

$$u_\gamma (𝐪)+\underset{\alpha ,\beta }{}\varphi _{\alpha \beta }n_{\alpha \gamma }(𝐪)u_\beta (𝐪)=0,$$
(31)
where we have introduced new functions

$$n_{\alpha \gamma }(𝐪)=\underset{i}{}\rho _i\mathrm{\Omega }_\alpha ^i(𝐪)\mathrm{\Omega }_\gamma ^i(𝐪).$$
(32)
Equation (31) can also be written in matrix form, but now the indices run over the set of weights, not the components. Hence, defining the matrices $`\widehat{\mathrm{\Phi }}\left(\varphi _{\alpha \beta }\right)`$ and $`\widehat{\mathrm{N}}(𝐪)\left(n_{\alpha \gamma }(𝐪)\right)`$, and the vector $`𝐮(𝐪)\left(u_\alpha (𝐪)\right)`$, Eq. (31) becomes

$$\left[\mathrm{I}+\widehat{\mathrm{N}}\widehat{\mathrm{\Phi }}\right]𝐮=\mathrm{𝟎}.$$
(33)
Again this system has nontrivial solutions if and only if

$$D(\rho ,𝐪)\mathrm{det}\left[\mathrm{I}+\widehat{\mathrm{N}}\widehat{\mathrm{\Phi }}\right]=0.$$
(34)
This alternative characterization of the spinodal is suitable for polydisperse systems (actually, it is a generalization of the formalism developed in Ref. CuestaEPL  for homogeneous phases). The reason is that it replaces integral operators by finite matrices, in which the number of components is limited by the number of weights of the theory. For instance, for the polydisperse Zwanzig model with a free energy functional given by the FMF derived in Ref. Cuesta , functions $`n_{\alpha \beta }`$ take the form

$$n_{\alpha \beta }(𝐪)=\underset{\nu }{}𝑑\lambda \rho ^\nu (\lambda )\mathrm{\Omega }_\alpha ^\nu (𝐪,\lambda )\mathrm{\Omega }_\beta ^\nu (𝐪,\lambda ),$$
(35)
where $`\alpha ,\beta =0,(1x,1y,1z),(2x,2y,2z),3`$ and

$`\mathrm{\Omega }_0^\nu (𝐪,\lambda )`$ $`=`$ $`{\displaystyle \underset{k}{}}w_0(q_k\mathrm{\Lambda }_{\nu k}),`$ (36)
$`\mathrm{\Omega }_3^\nu (𝐪,\lambda )`$ $`=`$ $`{\displaystyle \underset{k}{}}\mathrm{\Lambda }_{\nu k}w_3(q_k\mathrm{\Lambda }_{\nu k}),`$ (37)
$`\mathrm{\Omega }_{1,k}^\nu (𝐪,\lambda )`$ $`=`$ $`\mathrm{\Lambda }_{\nu k}w_1(q_k\mathrm{\Lambda }_{\nu k})\mathrm{\Omega }_0^\nu (𝐪,\lambda ),`$ (38)
$`\mathrm{\Omega }_{2,k}^\nu (𝐪,\lambda )`$ $`=`$ $`[\mathrm{\Lambda }_{\nu k}w_1(q_k\mathrm{\Lambda }_{\nu k})]^1\mathrm{\Omega }_3^\nu (𝐪,\lambda )`$ (39)
($`k=x,y,z`$). Here $`q_k`$ are the components of the vector $`𝐪`$, and

$`\mathrm{\Lambda }_{\nu k}`$ $`=`$ $`\lambda ^{1/3}+(\lambda ^{2/3}\lambda ^{1/3})\delta _{\nu k},`$ (40)
$`w_0(x)`$ $`=`$ $`\mathrm{cos}(x/2),`$ (41)
$`w_3(x)`$ $`=`$ $`2\mathrm{sin}(x/2)/x,`$ (42)
$`w_1(x)`$ $`=`$ $`w_3(x)/w_0(x).`$ (43)
As there are 8 different weights, $`D(\rho ,𝐪)`$ is in this case the determinant of an $`8\times 8`$ matrix.

Either way we express the equation, (29) or (34), the spinodal instability condition of a homogeneous phase with density distributions $`\rho ^\nu (\lambda )`$ (or simply $`\rho ^\nu `$ in the monodisperse case) with respect to the transition to an inhomogeneous phase occurs at the values of the total density $`\rho `$ and wave vector $`𝐪`$ for which the function $`D(\rho ,𝐪)`$ first vanishes. More precisely, this function has an oscillatory behavior as a function of $`𝐪`$, and is positive as long as the homogeneous phase is stable, so the spinodal corresponds to the smallest density $`\rho `$ for which the absolute minimum of $`D`$ with respect to $`𝐪`$ equals zero. This amounts to finding, for a given $`\rho `$, the solutions to

$$\mathbf{}D(\rho ,𝐪)=0.$$
(44)
When we later consider the rod-plate mixture we will have to locate the nematic-biaxial phase transition. As it is a continuous transition it can also be found as a solution to Eq. (34) with $`𝐪=0`$ (both phases are homogeneous), using

$$\rho ^\nu (\lambda )=P_0h(\lambda )\frac{e^{\mathrm{\Phi }^\nu (\lambda )}}{_\tau e^{\mathrm{\Phi }^\tau (\lambda )}}$$
(45)
as the distribution functions of the nematic phases.

## III Results

### III.1 Monodisperse Zwanzig model

As a first step towards the calculation of the whole phase diagram of the monodisperse (pure rods or pure plates) Zwanzig model using the FMF of Ref. Cuesta  we will determine the isotropic-nematic (I–N<sup>±</sup>) coexistence curves as well as the spinodals for the nematic-smectic (N<sup>±</sup>–S<sup>±</sup>), nematic-columnar (N<sup>±</sup>–C<sup>±</sup>) and isotropic-plastic (orientationally disordered) solid (I–PS) transitions, as a function of the length-to-breadth ratio $`\kappa `$. The $`+`$ and $``$ superscripts label prolate (rods) and oblate (plates) parallelepipeds respectively. If the transitions are first order, the location of the coexistence curves will differ from that of the corresponding spinodals. However, as the main concern of this paper is the effect of polydispersity, we defer the calculations of coexistence with and between inhomogeneous phases to a forthcoming publication.

We have obtained the spinodals for this model by solving Eqs. (29) and (44) using the monodisperse version of (45). If we choose the nematic director parallel to the $`z`$ axis, then the N<sup>±</sup>–S<sup>±</sup> spinodals can be calculated setting $`𝐪=(0,0,q)`$, whereas the N<sup>±</sup>–C<sup>±</sup> ones follow from taking $`𝐪=(q,0,0)`$ or $`𝐪=(0,q,0)`$ (both are equivalent due to the nematic symmetry). For the I–PS spinodal all three previous vectors give the same result.

Figure 1 shows the results of these calculations. We observe in this figure that there is I–N<sup>±</sup> transition only for $`\kappa <0.210`$ and $`\kappa >5.02`$. At those two limiting values of $`\kappa `$ the N<sup>+</sup> phase is destabilized by a S<sup>+</sup> at a packing fraction $`\rho =0.280`$, and the N<sup>-</sup> by a C<sup>-</sup> at a packing fraction $`\rho =0.261`$ (remember the particle volume has been set to 1). Also, the N<sup>+</sup>–S<sup>+</sup> spinodal is always below the N<sup>+</sup>–C<sup>+</sup>, and the N<sup>-</sup>–C<sup>-</sup> spinodal is always below the N<sup>-</sup>–S<sup>-</sup> one. This is what intuition tells us but, as we will see later, it is a peculiarity of the monodisperse system. Despite this, all these spinodal curves converge asymptotically to the same packing fraction, $`\rho =0.314`$, as $`\kappa \mathrm{}`$ or $`\kappa 0`$. This is precisely the packing fraction at which the continuous freezing transition occurs in a system of parallel hard cubes ($`\kappa =1`$) yuri1 . The reason for this is that upon increasing $`\kappa `$ the number of rods with orientation perpendicular to the director becomes vanishing small, and then the system is, after rescaling the $`z`$ direction, almost equivalent to a system of parallel cubes. The same holds for plates upon decreasing $`\kappa `$. What this means is that most likely the N<sup>+</sup>–S<sup>+</sup> and N<sup>-</sup>–C<sup>-</sup> transitions will be metastable with respect to freezing in a large portion of the phase diagram. For $`0.182<\kappa <4.93`$ we have found an I–PS spinodal instability. The peculiarity of this curve is that it exhibits strong oscillations as the aspect ratio changes. These oscillations reflect the packing efficiency of randomly oriented parallelepipeds as a function of their size (the better the packing the lower the curve).

The available simulation results for freely rotating hard spherocylinders show that the I–S<sup>+</sup> and N<sup>+</sup>–S<sup>+</sup> transitions begin at $`\kappa =4.1`$ and $`4.5`$ respectively Bolhuis1  (notice that for hard spherocylinders the length-to-breadth ratio is $`\kappa =\frac{L+D}{D}`$). On the other hand, simulations of hard cut spheres show that for $`\kappa =0.2`$ there is an I–C<sup>-</sup> transition (the isotropic phase might instead be a peculiar ‘cubatic’ phase) and for $`\kappa =0.1`$ a N<sup>-</sup>–C<sup>-</sup> one Veerman . We can see that despite the different particle geometry and the restricted orientations of the Zwanzig model, the agreement with the threshold $`\kappa `$’s at which spatial instabilities destabilize the homogeneous phases is rather good. Also the qualitative picture is similar: elongated rods form smectics, while flat disks form columnars, and more symmetric particles form solids instead of smectics or columnar. Thus, this simple model seems to capture the essence of the entropy driven phase transitions between phases with different symmetries and their relation to particle anisotropy.

The only existing simulation of the Zwanzig model has been performed on a lattice model of parallelepipeds with length-to-breadth ratios $`15:3`$ ($`\kappa =5`$) and $`3:15`$ ($`\kappa =0.2`$) Casey . For these $`\kappa `$’s the authors find I–C<sup>-</sup> for the discs and and I–DS<sup>+</sup> for the rods, where DS<sup>+</sup> stands for a novel phase called ‘discotic smectic’ by the authors (in this phase the axes of the particles point perpendicular to the normal of the smectic layers while there is no orientational order within the layer). The packing fractions they obtain for the transitions are higher than those of the corresponding spinodals of Fig 1, but they should decrease upon decreasing the lattice spacing, as the results for the freezing of parallel hard cubes on a lattice (occurring at $`\rho =0.568`$ for edge-length 2 lattice spacings, at $`\rho =0.402`$, for edge-length 6 lattice spacings, and at $`\rho =0.314`$ for the continuum) illustrate Lafuente . All this is again in agreement with the phase diagram of Fig. 1.

### III.2 Polydisperse Zwanzig model: Unimodal distribution

Let us study now the effect polydispersity has on the phase behavior of the Zwanzig model discussed above. We have introduced polydispersity through the unimodal parent probability density resulting when setting $`\zeta =1`$ (rod-like parallelepipeds) or $`\zeta =0`$ (plate-like parallelepipeds) in function (17). Before reporting the results we have obtained, a few words on polydisperse phase diagram plots are on purpose.

Generally speaking, plotting polydisperse phase diagrams would require an infinite dimensional space, for as we vary the fraction of total volume occupied by the coexisting phases they change their composition. A phase diagram like that of Fig. 1 plots densities vs. aspect ratios, thus carrying no information whatsoever about compositions. So we must restrict ourselves to a given composition, and this is the parent one, since all pure phases have this composition. This means that only cloud lines can be plotted, delimiting regions where we can only know that the system is decomposed into two or more coexisting phases. The result resembles the coexisting lines of monodisperse phase diagram, but the meaning is completely different. Although we could connect cloud points in the two lines delimiting the coexistence region, they are definitely no coexisting states.

Similar considerations hold for spinodal lines. We can represent the spinodal of a system in a pure phase having the parent distribution as composition, but if the system reaches a cloud point an incipient shadow phase with a totally different composition—hence a totally different spinodal line—coexists with it. If the density of this shadow phase is above its own spinodal line, the phase transition will not be stable.

Having all this in mind we have solved the coexistence equations of the isotropic and nematic phases and determined the spinodal lines of inhomogeneous instabilities for a system with the parent distribution (17). The results for the choice $`\alpha =1`$ (corresponding to polydispersity parameters $`\mathrm{\Delta }_L=0.288`$ and $`\mathrm{\Delta }_D=0.143`$) are shown in Fig. 2.

Apart from the usual broadening of the isotropic-nematic transition already observed for this model in Ref. Clarke, , we find a clear enhancement of the stability of the homogeneous phases i.e. the spinodal instabilities occur at higher densities. This effect is more pronounced for the N<sup>±</sup>–S<sup>±</sup> or N<sup>±</sup>–C<sup>±</sup> transitions. Another remarkable feature of Fig. 2 is that the columnar phases are by far more stable than the smectic ones for any aspect ratio.

This general behavior is in agreement with the results obtained from both density functional theory for length-polydisperse parallel hard cylinders Bohle  and computer simulations for length-polydisperse freely rotating hard spherocylinders in the Onsager limit Bates . Although in Ref. Bohle  some strong approximations were made, such as the decoupling of the density profile as $`\rho (𝐫,L)=g(L)\rho (𝐫)`$, and the effect of fractionation between the nematic and columnar phases was ignored, the results agree qualitatively with the simulations: both show a terminal polydispersity beyond which the smectic phase is no more stable, being replaced by a columnar phase (which tolerates a higher degree of polydispersity). The nematic-columnar transition is first order. While fractionation is indeed negligible in the nematic-smectic transition, the smectic-columnar transition clearly segregates long rods into the columnar phase.

It is interesting to notice that in the limits $`\kappa \mathrm{}`$ or 0 none of the spinodal lines tend to the density of the freezing of parallel hard cubes ($`\kappa =1`$). The reason is that no trivial scaling can be applied to a polydisperse system of perfectly aligned parallelepipeds to transform it into a system of parallel hard cubes. On the other hand the S<sup>+</sup> and S<sup>-</sup> spinodals tend to the same limiting density when $`\kappa 0,\mathrm{}`$, and so do the C<sup>+</sup> and C<sup>-</sup> spinodals. The second virial term of the free energy is dominant for large (small) $`\kappa `$’s; this and the rod-plate symmetry of the chosen parent distributions explains this behavior.

### III.3 Binary mixture of Zwanzig rods and plates

We have calculated the phase diagram (coexistence of homogeneous phases as well as their spinodal instabilities to inhomogeneous phases) for the binary mixture of rods and plates with $`\kappa =5`$ and 15. The results for $`\kappa =5`$ are shown in Fig. 3. Apart from the usual I–N<sup>±</sup> transitions and N<sup>-</sup>–N<sup>+</sup> phase separations which coalesce in a multicritical point we see that the nematic phases increase their stability with respect to the inhomogeneous ones as the relative composition of the fluid tends to the equimolar mixture. Another important feature is that the biaxial nematic is not stable.

All these characteristics are also present in the phase diagram obtained using a second virial approximation Van-Roij1 , but there are important differences as well. The most prominent one is that the whole phase diagram is shifted down in density, with all transitions and spinodals appearing below $`\rho =0.4`$ (those obtained with a second virial approximation occur above that value Van-Roij1 ). This result should not surprise if one takes into account that in a $`y3`$ expansion the excess free energy is written in terms of the variable $`y=\rho /(1\rho )`$, which grows faster than $`\rho `$. Another striking difference with respect to the second virial approach is the rod-plate asymmetry of the phase diagram, a consequence of the inclusion of higher virial coefficients—which do not share the symmetry of the second one.

We pass now to describe the loss of stability of the nematics with respect to inhomogeneous phases. The intersection between the N<sup>+</sup> line of the I–N<sup>+</sup> coexistence and the N<sup>+</sup>–C<sup>+</sup> spinodal in the rod-rich part of the phase diagram (marked in the figure with a full circle; the open circle corresponds to its coexisting isotropic phase) is consistent with the existence of a first order phase transition between the I phase and an inhomogeneous (presumably columnar) phase. (Recall that simulations on the lattice find a first order I-DS<sup>+</sup> transition for parallelepipeds with $`\kappa =5`$ Casey ). Nevertheless a definitive answer about the relative stability between these inhomogeneous phases—including also the solid—can only be given after carrying out coexistence calculations.

Increasing the fraction of plates (which in the C<sup>+</sup> phase have their principal axes oriented perpendicular to the director) seems to favor the smectic alignment of rods. This results in the intersection of the N<sup>+</sup>–C<sup>+</sup> and the N<sup>+</sup>–S<sup>+</sup> spinodals (marked with a full square on the right of Fig. 3). In the second virial approach the N<sup>+</sup>–S<sup>+</sup> spinodal is always below the N<sup>+</sup>–C<sup>+</sup> one. At the plate-rich part of the phase diagram the nematic exhibits a spinodal to the C<sup>-</sup> phase, as usual in discotic fluids Veerman . Increasing the fraction of rods beyond a threshold value (indicated with a full square on the left part of the phase diagram) makes the smectic more stable than the columnar, which agrees qualitatively with the results from the second virial approach Van-Roij1 .

The phase diagram for $`\kappa =15`$ appears in Fig. 4. Its main difference with the previous one is the presence of a thermodynamically stable biaxial nematic phase in an inverted triangular right above the multicritical point. The window is limited to the left and to the right by continuous phase transitions to uniaxial nematic phases, and to the top by a B–N<sup>+</sup> coexistence which becomes a N<sup>-</sup>–N<sup>+</sup> phase separation as density increases. As discussed in Ref. Van-Roij1 , the driving mechanism which determines the preference of this system for phase separation instead of biaxial ordering is the larger exclude volume between unlike particles compared to that of like particles (the rod-plate excluded volume divided by the rod-rod one scales as $`\kappa ^{2/3}`$ for large $`\kappa `$). When the gain in free volume compensates the loss in mixing entropy (and this strongly depends on concentration and composition) phase separation occurs.

The slight asymmetry of the phase diagram is again due to the presence of higher virial terms in the free energy. The second virial approximation predicts a small region of triple coexistence between the N<sup>-</sup>, N<sup>+</sup> and B phases right on top of the biaxial phase region. The absence of this triple coexistence region in Fig. 4 is due to the asymmetry of the phase diagram. Simulations of a mixture of prolate and oblate ellipsoids Camp  show a similar asymmetric phase diagram, only tilted in the opposite direction (there is a B–N<sup>-</sup> demixing on top of the biaxial nematic, instead of the B–N<sup>+</sup> of the present model), probably due to the different particle geometry. Despite the asymmetry, we should admit that for $`\kappa =15`$ there is good overall agreement between the second virial approach and ours.

We have not studied the changes an asymmetry in the mixture (a difference in the volume of rods and plates, for instance) would induce in the phase diagram. Intuitively, asymmetry acts against the stability of the biaxial phase, something that is confirmed by a recent work Wensink1  which studies a binary rod-plate mixture of freely rotating oblate and prolate cylinders in the Onsager approximation. This work explicitly shows that the region of stability of the biaxial phase decreases with increasing the asymmetry of the mixture (although, according to this work, this phase can be found in highly asymmetry mixtures).

With respect to the inhomogeneous phases, Fig. 4 shows a remarkable behavior: although pure hard rods in the N<sup>+</sup> phase undergo a spinodal instability to a smectic one, by adding a tiny fraction of plates the instability takes place to a columnar phase. This phenomenon has already been observed in a different system, namely, in simulations of binary mixtures of parallel hard spherocylinders Stroobants . Although pure parallel hard spherocylinders exhibit a continuous nematic-smectic transition, a binary mixture of them with length-to-breadth ratios 2 and 2.9 show a first order nematic-columnar transition instead. This result was explained by the poorer packing of rods of different length in the smectic phase as compared to that of rods of the same length. Borrowing the argument, in our rod-plate model the plates can fit into the interlayer spacing with their axes parallel to the smectic director as long as there are few of them; but upon increasing their molar fraction some of them are forced to get into the smectic layers with their axes perpendicular to the smectic director, thereby destabilizing the smectic phase.

Finally, at the other extreme of phase diagram the N<sup>-</sup> phase losses its stability to a columnar phase of plates at densities lower than those of the N<sup>+</sup>–C<sup>+</sup> spinodal.

### III.4 Polydisperse Zwanzig model: Bimodal distribution

In this section we are going to study the effect of polydispersity on the phase behavior of the previous rod-plate mixture. In Ref. Martinez-Raton1  we have shown how polydispersity can stabilize the biaxial phase even for mixtures with relative small aspect ratios—like $`\kappa =5`$, for which Fig. 3 shows that the biaxial phase is absent.

Apart from the increase in mixing entropy that polydispersity carries, the other mechanism behind the enhancement of stability of the biaxial phase is the decrease of the ratio of the average exclude volumes of like and unlike particles. The excluded volume between two parallel rods with lengths $`L_i`$ and breadth $`D_i`$ ($`i=1,2`$) is

$$v_{rr}=\left(L_1+L_2\right)\left(D_1+D_2\right)^2,$$
(46)
and the excluded volume between a rod with length $`L_1`$ and breadth $`D_1`$ and a plate with length $`L_2`$ and breadth $`D_2`$, with axes perpendicular to each other, is

$`v_{rp}=\left(L_1+D_2\right)\left(L_2+D_1\right)\left(D_1+D_2\right).`$ (47)
According to (1), if the volume of all particles is set to one, in terms of the length-to-breadth ratios $`\lambda _i`$, $`L_i=\lambda _i^{2/3}`$ and $`D_i=\lambda _i^{1/3}`$. Taking a double average of $`v_{rr}`$ over $`_i\left[h(\lambda _i)\mathrm{\Theta }(\lambda _i1)\right]`$ and of $`v_{rp}`$ over $`\left[_ih(\lambda _i)\right]\mathrm{\Theta }(\lambda _11)\mathrm{\Theta }(1\lambda _2)`$ \[$`h(\lambda )`$ is the parent distribution (17)\], and fixing $`\zeta =1/2`$ (the equimolar composition) we arrive at

$$\frac{v_{rp}}{v_{rr}}=\frac{\frac{1}{2}+m_1^2+m_1^2+2(m_1m_2+m_1m_2)}{\frac{1}{2}+4m_1m_1+2m_2m_2},$$
(48)
where $`m_\beta _1^{\mathrm{}}𝑑\lambda h(\lambda )\lambda ^{\beta /3}`$. It can be shown numerically that, for a given $`\kappa `$, the ratio (48) decreases with polydispersity (i.e. with decreasing $`\alpha `$). Analytic expressions for this ratio can be obtained in the limit of high particle anisotropy ($`\kappa 1`$) and high polydispersity ($`\alpha 1`$) with the constraint $`\alpha \kappa ^21`$ (in terms of the parameters $`\mathrm{\Delta }_\nu `$, $`\nu =L,D`$, this constraint is equivalent to $`1\mathrm{\Delta }_\nu \sqrt{\mathrm{ln}\kappa }`$, implying that the fraction of cubic-like particles is vanishing small). Using the asymptotic expressions

$$m_\beta \frac{\kappa ^{\beta /3}}{\alpha ^{\beta /6}(\mathrm{ln}\alpha )},\mathrm{\Delta }_\nu \left(\mathrm{ln}\alpha \right)^{1/2},$$
(49)
we obtain

$$\frac{v_{rp}}{v_{rr}}\kappa ^{2/3}e^{c_\nu \mathrm{\Delta }_\nu ^2},\nu =L,D$$
(50)
with $`c_\nu `$ a positive constant. This asymptotic relationship explicitly shows both the scaling of this ratio discussed for the pure rod-plate mixture Van-Roij1  and the exponential attenuation of this scaling with polydispersity.

We have estimated the threshold value of $`\kappa `$ (denoted $`\kappa ^{}`$) beyond which the biaxial phase begins to be thermodynamically stable, as a function of polydispersity. The results are summarized in Table 1. Without polydispersity this value is $`\kappa ^{}=6.9`$, smaller than the second virial estimation Van-Roij1  $`\kappa ^{}=8.8`$. This is an indication that three (and higher) body correlations increase the stability of the biaxial phase. As we can see in Table 1, $`\kappa ^{}`$ decreases upon increasing polydispersity, an illustration of the enhancement of stability of the biaxial ordering induced by polydispersity.

In Fig. 5 we plot the phase diagram of this system for $`\kappa =5`$ and $`\alpha =0.1`$ (corresponding to $`\mathrm{\Delta }_L=0.610`$ and $`\mathrm{\Delta }_D=0.302`$). According to Table 1 for this parameter choice the biaxial phase is thermodynamically stable, and this is clearly shown in Fig. 5(a), where we can see a small window of biaxial nematic just above the multicritical point. Figure 5(b) is a zoom of the upper border of this window. There is a narrow region of N<sup>-</sup>–B coexistence limited from above by a N<sup>-</sup>–N<sup>+</sup> coexistence. The dividing line (the dotted line in the diagram) represents a continuous transition between the B and the N<sup>+</sup> phases which coexist with the N<sup>-</sup>. This line has no analog in the binary rod-plate mixture. To obtain it one has to solve the equations for two phase coexistence at all values of $`\gamma `$ in (16) (the volume occupied by one of the coexisting phases) and then find the one for which the B–N<sup>+</sup> transition occurs.

The high density part of the phase diagram of Fig. 5(a) was not determined in Ref. Martinez-Raton1 . The most remarkable feature it shows is that pure phases are hardly stable: the diagram is dominated by coexistence regions. The reason is that polydispersity is so high that the mixture has, beside rods and plates, a significant amount of cubic-like particles. This favors entropic phase separation in the three different phases I, N<sup>-</sup> and N<sup>+</sup>.

In the upper part of Fig. 5(a) there is a region of triple coexistence I–N<sup>-</sup>–N<sup>+</sup>. The lines limiting this region were obtained by solving the coexistence equations for the case when one of the coexisting phases is a shadow phase (i.e. it occupies a vanishing fraction of the total volume). As an illustrative example, if we take $`\gamma _\mathrm{I}=0`$, $`\gamma _\mathrm{N}^{}=\gamma `$ and $`\gamma _\mathrm{N}^{}=1\gamma `$ and solve the coexistence equations for the $`N=15`$ unknowns (see the discussion about phase equilibria in polydisperse systems in Section II.2) we obtain the curve joining the two points marked with full squares in Fig. 5(a). The other two curves limiting the triple coexistence region were calculated selecting $`\gamma _{\mathrm{N}^+}=0`$ (the left one) and $`\gamma _\mathrm{N}^{}=0`$ (the right one).

We have checked that the nematics below the triple region are always stable with respect to spatial density modulations of any kind. In order to estimate the packing fraction at which these instabilities occur we solved Eqs. (34) and (44) using the density distributions $`\rho ^\nu (\lambda )`$ of all the coexisting phases along the borders of the triple coexistence region. The left and right open circles correspond respectively to the packing fractions at which the coexisting N<sup>+</sup> and N<sup>-</sup> (the shadow phase in both cases) lose their stability to the C<sup>+</sup> and C<sup>-</sup> phases. Above these points the phase diagram will possibly include coexistences between four phases: I, N<sup>-</sup>, N<sup>+</sup> and C<sup>-</sup> or C<sup>+</sup>) (the experiments of Ref. Kooij  find a similar scenario).

Figure 6 shows the density distribution functions of the phases coexisting at the points marked with full squares in Fig. 5(a). As already discussed, the curves illustrate the high fraction of cubic-like particles partly responsible for the strong demixing this system exhibits. The figure also illustrates the strong fractionation which takes place between the coexisting phases: the isotropic phase is rich in cubic-like particles, the N<sup>-</sup> in plates and the N<sup>+</sup> in rods.

In contrast, the phase diagram of the polydisperse rod-plate mixture for $`\kappa =15`$ does not change qualitatively compared to the binary mixture (cf. Figs. 7 and 4): the size and shape of the biaxial regions are similar and the general topology of the phase diagram is basically the same. The only important difference (minor in terms of the size of the portion of phase diagram it involves) appears just above the biaxial: one can observe both B–N<sup>+</sup> and B–N<sup>-</sup> coexistence (the former occupying a larger region), as well as a triple coexistence zone separating them \[see Fig. 7(c)\]. As for $`\kappa =5`$, the B–N<sup>±</sup> coexistences become N<sup>+</sup>–N<sup>-</sup> coexistence through a second order phase transition of the B phase to the N<sup>±</sup> \[dotted lines in Fig. 7(b) and (c)\].

The N<sup>-</sup> and N<sup>+</sup> lose their stability to columnar phases at relative high packing fractions \[see the dotted lines in the upper part of Fig. 7(a)\]. In both cases it is the N<sup>-</sup> phase which first becomes unstable with respect to C<sup>-</sup> at the N<sup>+</sup>–N<sup>-</sup> coexistence (the N<sup>-</sup> cloud on the left, marked with a full circle in Fig. 7(a), and the N<sup>-</sup> shadow on the right, marked with an empty circle in this figure).

Figure 8 shows the density distribution functions corresponding to the point marked by a full square in Fig. 7(c). The distributions of the N<sup>+</sup> and B phases are very similar, except for the fact that the N<sup>+</sup> has a slightly higher proportion of rods (and correspondingly less of plates) than the $`B`$ phase. It is clear from the figure that the proportion of cubic-like particles is negligible, so this mixture can be regarded as a true (polydisperse) rod-plate mixture.

## IV Conclusions

We have studied the effect of polydispersity on the phase diagrams of several variants of the Zwanzig model for liquid crystals. We have first determined the spinodal instabilities of the homogeneous phases to the smectic, columnar or plastic solid phases for one component fluids of rods or plates. There is a general qualitative agreement with phase diagrams of oblate and prolate particles obtained from Monte Carlo simulations Bolhuis1 ; Veerman . We have also shown that even a small degree of polydispersity has dramatic effects on the transitions to highly ordered phases (like the smectic or the columnar). In particular, we have shown that upon increasing polydispersity the transition to the columnar phase of rods preempts the N<sup>+</sup>–S<sup>+</sup> transition, in agreement with simulations Polson .

Next we have obtained the phase diagrams of symmetric rod-plate binary mixtures with length-to-breadth ratios $`\kappa =5`$ and 15. While for large $`\kappa `$’s the differences between these phase diagrams and those calculated with a second virial approximation Van-Roij1  are only quantitative, for small $`\kappa `$’s both methods differ qualitatively. For instance, all the transitions occur at packing fractions well below those obtained with the second virial approach, and sometimes the relative stability between different inhomogeneous phases changes.

We have shown that the introduction of polydispersity in a rod-plate mixture in a very symmetric way enhances the stability of the biaxial phase even for as small an aspect ratio as $`\kappa =5`$. For this case the amount of polydispersity needed to stabilize the biaxial phase is so high that the mixture becomes very unstable with respect to phase separation at practically any composition. There are large regions of three phase and possibly four phase coexistence (I–N<sup>-</sup>–N<sup>+</sup> and I–N<sup>-</sup>–N<sup>+</sup>–C<sup>±</sup>). These results agree with what is observed in experiments with rod-plate mixtures with a high degree of polydispersity Kooij .

###### Acknowledgements.

This work is part of the research Project No. BFM2000-0004 (DIGI) of the Ministerio de Ciencia y Tecnología (Spain). Y. M-R. is supported by a postdoctoral grant of the Consejería de Educación de la Comunidad de Madrid.


format: text

text:
# Topological Q-Solitons

## 1 Introduction

Topological solitons (stable, localized, topologically non-trivial solutions in a field theory) have long been of great interest, both for their mathematical properties and for their applications in many areas of physics. In any system admitting such solitons, the non-trivial topology is not sufficient to ensure that the solitons are stable. One obvious reason for this is that topology is, by definition, a non-metric structure, and so it cannot determine the size of the solitons; for that, one needs to balance the forces acting on the soliton in such a way that it has a preferred size. Recall, for example, the O($`n`$) sigma model in $`D+1`$ dimensions (with trivial boundary condition at spatial infinity), which is the subject of this article. This system admits topological configurations (textures) whenever the homotopy group $`\pi _D(S^{n1})`$ is non-trivial; in particular, for $`(D,n)`$ equal to $`(2,3)`$, $`(3,3)`$ or $`(3,4)`$. If $`D=3`$, then solitons tend to shrink — in the pure sigma model, there are no static solutions. In the $`D=2`$ case, there are static solutions, for example the Belavin-Polyakov solitons  in the O(3) system; but these are unstable .

A soliton can always be prevented (or rather discouraged) from spreading out by the addition, if necessary, of a potential (a term involving only the field, and not its gradient). In order to stabilize the soliton size, we also need to introduce something which prevents it from shrinking. There are several possibilities for such an anti-shrinking mechanism: for example, a Skyrme term involving four (or more) powers of the field gradient; or a gauge field suitably coupled to the sigma-field; or periodic time-dependence (rotation in an internal space). This third possibility also underlies non-topological solitons (Q-balls). In this paper, we investigate to what extent the Q-ball mechanism is effective at stabilizing topological sigma-model solitons. We shall see that stationary topological Q-solitons exist in $`D=2`$ spatial dimensions, but not for $`D=3`$. This result is analogous to that for the Landau-Lifshitz equation (Heisenberg model of ferromagnetism), as one might have surmised since the static Landau-Lifshitz system is identical to the static sigma model.

## 2 Q-Solitons in the O(3) Sigma Model

The O(3) sigma model involves a scalar field taking values on $`S^2`$; this field can be represented as a unit 3-vector $`\stackrel{}{\varphi }=(\varphi _1,\varphi _2,\varphi _3)`$ with $`\stackrel{}{\varphi }\stackrel{}{\varphi }=1`$. The Lagrangian is

$$={\scriptscriptstyle \frac{1}{2}}(_\mu \stackrel{}{\varphi })(^\mu \stackrel{}{\varphi })V(\varphi _3),$$
(1)
where $`V`$ is some potential function (which, for simplicity, we take to depend only on $`\varphi _3`$). The space-time coordinates are $`x^\mu =(t,x^j)`$, with $`j=1,\mathrm{},D`$. The system has a global SO(3) symmetry which is broken, by the potential term, to SO(2). This SO(2) acts only on $`\varphi _1`$ and $`\varphi _2`$, namely by changing the phase of $`\varphi :=\varphi _1+\mathrm{i}\varphi _2`$. The corresponding conserved quantity is

$$Q=\mathrm{Im}(\dot{\varphi }\overline{\varphi })d^Dx.$$
(2)
Minimizing the energy of a configuration subject to $`Q`$ being fixed implies , in particular, that $`\stackrel{}{\varphi }`$ has the form

$$\varphi (t,x^j)=\psi (x^j)\mathrm{e}^{\mathrm{i}\omega t},\varphi _3=\psi _3(x^j),$$
(3)
with $`|\psi |^2+(\psi _3)^2=1`$. Without loss of generality, we shall assume that $`\omega 0`$. Note that $`Q=\omega I`$, where $`I=|\psi |^2d^Dx`$. The energy of a configuration of the form (3) is $`E=E_d+E_k+E_p`$, where

$`E_d`$ $`=`$ $`{\scriptscriptstyle \frac{1}{2}}{\displaystyle \left[|_j\psi |^2+(_j\psi _3)^2\right]d^Dx},`$
$`E_k`$ $`=`$ $`{\scriptscriptstyle \frac{1}{2}}I\omega ^2={\scriptscriptstyle \frac{1}{2}}Q^2/I,`$
$`E_p`$ $`=`$ $`{\displaystyle V(\psi _3)d^Dx}.`$
The boundary condition is $`\psi _31`$ as $`r\mathrm{}`$; so we need $`V(1)=0`$.

A stationary Q-lump is a critical point of the energy functional $`E[\stackrel{}{\psi }]`$, subject to Q having some fixed value. Such a Q-lump is (classically) stable if this critical point is a local minimum of $`E`$. The usual (Derrick) scaling argument shows that any stationary Q-lump must satisfy

$$(2D)E_dDE_p+DE_k=0.$$
(4)
Let the positive constant $`m`$ be defined by $`V^{}(1)=m^2`$; in other words, $`V(\psi _3)m^2(1\psi _3){\scriptscriptstyle \frac{1}{2}}m^2|\psi |^2`$ for $`\psi _31`$. Then, near spatial infinity, the Euler-Lagrange equations corresponding to $`E`$ imply that

$$^2\psi (m^2\omega ^2)\psi =0.$$

So in order to satisfy the boundary condition $`\psi 0`$ as $`r\mathrm{}`$, we need $`\omega m`$. The solitons are exponentially localized if $`\omega <m`$, but less-localized solitons with $`\omega =m`$ may also exist.

The parameter $`m`$ is part of the specification of the system, and the parameter $`Q`$ is set by the initial data. Each of these two parameters has dimensions; the combination $`Qm^{D1}`$ is dimensionless, whereas the combination $`(Q/m)^{1/D}`$ has dimension of length, and determines the size of the soliton. Configurations in this system are classified topologically by their topological charge $`N`$ (an integer); if $`D=2`$, then $`N`$ is the winding number, while if $`D=3`$, then $`N`$ is the Hopf number. Let $`E(N,Q)`$ denote the energy of a configuration (or rather, of data $`\stackrel{}{\varphi },\stackrel{}{\varphi }_t`$) with topological charge $`N`$ and Noether charge $`Q`$.

Let $`V_0`$ be the normalized potential function $`V_0:=2V/m^2`$. Note that $`V_0(\psi _3)/|\psi |^21`$ as $`\psi _31`$. It is clear from (4) that if $`V_0|\psi |^2`$ everywhere, with $`V_0>|\psi |^2`$ somewhere, then there can be no solution. So the constant $`K`$ defined by $`K=\mathrm{min}\left[V_0(\psi _3)/|\psi |^2\right]`$ should satisfy $`K1`$. It then follows that

$$E_p=m^2V_0m^2KI=(Km^2/\omega ^2)E_k(Km^2/\omega ^2)E_p,$$
(5)
where the final inequality comes from (4). As a consequence, we have

$$Km^2\omega ^2m^2.$$
(6)
In the $`D=3`$ case, the first inequality is strict: $`Km^2<\omega ^2m^2`$.

In two spatial dimensions, it is possible to have $`K=1`$, which corresponds to the choice $`V(\psi _3)={\scriptscriptstyle \frac{1}{2}}m^2(1\psi _3^2)`$. So here $`\omega =m`$. This system , and generalizations in which the target space is some other Kähler manifold, arise naturally by dimensional reduction from ‘pure’ sigma models in one dimension higher . The energy satisfies a Bogomolny bound $`E(N,Q)4\pi N+mQ`$, which (for $`N2`$) can be saturated: for each value of $`N2`$ and $`Q`$, there is an explicit family $``$ of stationary multi-soliton solutions such that $`E(N,Q)=4\pi N+mQ`$. There is no force between the individual solitons: in particular, the total energy has the additive property $`E(N_1,Q_1)+E(N_2,Q_2)=E(N_1+N_2,Q_1+Q_2)`$. One may use moduli-space methods (as was done in  for other sigma-model systems) to investigate the scattering of moving solitons ; this involves finite-dimensional mechanics on $``$. The dynamics turns out to be rather exotic (as is also the case  for non-topological Q-balls). The solitons are only polynomially localized, and the non-existence of an $`N=1`$ soliton is related to this; an $`N=1`$ configuration tends to shrink in size, and there is no stationary $`N=1`$ solution.

On the other hand, if $`K<1`$, then 1-solitons can exist. Different choices of $`V`$ (having $`K<1`$) seem to lead to similar behaviour, but this has yet to be fully investigated; in what follows, we take $`V_0={\scriptscriptstyle \frac{1}{2}}(1\psi _3^4)`$, so $`K={\scriptscriptstyle \frac{1}{2}}`$. Let us consider, first, the thin-wall limit , where $`Qm^{D1}1`$. In this limit, the (bulk) contributions $`E_p`$ and $`E_k`$ to the energy are very much greater than the (surface) contribution $`E_d`$. So the energy is approximately $`E{\scriptscriptstyle \frac{1}{4}}m^2(1\psi _3^4)+{\scriptscriptstyle \frac{1}{2}}Q^2/(1\psi _3^2)`$. Without loss of generality, we may assume that $`\psi _3=1`$ outside of some compact set. So space is partitioned into three regions: one (with infinite volume) where $`\psi _3=1`$, the second (with volume $`A`$) where $`|\psi _3|1`$, and the third (with volume $`B`$) where $`\psi _3=1`$. Note that $`E`$ depends on $`A`$ (and on the value that $`\psi _3`$ takes on $`A`$), but not on $`B`$, and that

$$\frac{\delta E}{\delta \psi _3}=(\omega ^2m^2\psi _3^2)\psi _3.$$

So for fixed $`A`$, the function $`E`$ has a minimum for $`\psi _3=0`$ ($`\psi _3=\pm \omega /m`$ are local maxima). Hence we should set $`\psi _3=0`$ on $`A`$, and the energy becomes

$$E={\scriptscriptstyle \frac{1}{2}}Q^2/A+{\scriptscriptstyle \frac{1}{4}}m^2A.$$

Thus for a given value of $`Q`$, the energy has the minimum value $`E_{\mathrm{min}}=mQ/\sqrt{2}`$ when $`A=Q\sqrt{2}/m`$. Note that $`\omega =m/\sqrt{2}`$, at the lower end of its allowed range (6).

To make further progress, we need to include the effect of surface tension, in other words include the term $`E_d`$. Let us consider, first, the planar case $`D=2`$. For simplicity, we assume rotational symmetry about a point in the plane: the field is taken to have the form $`\psi =\mathrm{sin}(f)\mathrm{exp}(\mathrm{i}N\theta )`$ and $`\psi _3=\mathrm{cos}(f)`$, where $`f=f(r)`$. The boundary conditions are $`f(0)=\pi `$ and $`f(\mathrm{})=0`$, and $`N`$ is the topological charge. The energy functional $`E=E_d+E_k+E_p`$ was minimized numerically, for various values of $`m`$, $`Q`$ and $`N`$. The term $`E_k`$ was used in the form $`E_k=Q^2/I`$, so that one can minimize while keeping $`Q`$ fixed; the quantity $`\omega `$ does not enter explicitly, but can be derived (via the formula $`\omega =Q/I`$) once the minimum has been found. In each case that was investigated, a smooth minimum was reached. For $`N=1`$ and $`Qm=3200\pi `$ (close to the thin-wall limit), the profile function $`f(r)`$ is plotted in Fig 1. We see that there is a region around $`r=0`$ where $`f(r)`$ drops rapidly from $`\pi `$ to $`\pi /2`$ (the term $`E_d`$ prevents this region from shrinking further in size); and then a region (corresponding to $`A`$ in the argument above) where $`f=\pi /2\psi _3=0`$. Outside of this region, the field takes on its asymptotic value $`f=0\psi _3=1`$.

Fig 2 displays results for $`N=1`$, $`m=1`$ and a range of values of $`Q`$. We see that $`E`$ is very close to being linear in $`Q`$ (recall that in the Bogomolny case, it is exactly linear): to a very good approximation, we have $`E=4\pi +3Q/4`$. As for $`\omega `$, we know from (6) that $`\omega `$ has to be in the range $`1/\sqrt{2}<\omega <1`$, and we see from the figure that this is so; furthermore, $`\omega 1/\sqrt{2}`$ as $`Q\mathrm{}`$ (the thin-wall limit ), while $`\omega 1`$ as $`Q0`$ (the thick-wall limit ).

In the Bogomolny case  mentioned previously, the energy $`E(N,Q)`$ of a stationary soliton has the feature that, for a given $`Q`$, the quantity $`E_N:=(4\pi N)^1E(N,NQ)`$ is independent of $`N`$: in fact, $`E_N=1+mQ/4\pi `$. This corresponds to the fact that in this Bogomolny-type system, there is no force between stationary solitons. For the potential $`V={\scriptscriptstyle \frac{1}{4}}m^2(1\psi _3^4)`$, however, there are such forces. This can be seen by examining $`E_N`$ for fixed $`Q=4\pi `$ and for various values of $`N`$. Table 1 shows the results for $`1N7`$. The energy density of the 1-soliton is peaked at the point $`r=0`$, whereas that for the rotationally-symmetric $`N`$-soliton is peaked on a ring. Note that $`E_N`$ is a decreasing function of $`N`$, which suggests that the this $`N`$-soliton is stable against breakup into solitons of lower topological charge. But this remains to be checked; in particular, one should investigate the vibrational modes about these rotationally-symmetric solutions.

Finally, let us turn to the case of $`D=3`$ spatial dimensions. Configurations with nonzero Hopf number $`N`$ look like closed loops, which may be linked or knotted. Static knot-solitons are have been studied in the Faddeev-Skyrme system, where a Skyrme term is added to the Lagrangian: this extra term stabilizes the solitons, which would otherwise shrink . The question here is whether there exist stationary Hopf solitons which are stabilized by internal rotation rather than by a Skyrme term. The answer to this question appears to be negative; what happens is as follows. Consider an $`N=1`$ configuration. Typically, $`\psi _3=1`$ at spatial infinity and on a curve which extends to infinity; let us visualize this curve as the $`z`$-axis. Secondly, $`\psi _3=1`$ on a closed loop $`L`$ around the $`z`$-axis. Finally, $`\psi _3=0`$ on a torus around the $`z`$-axis, with the loop $`L`$ in its interior. So the region $`A`$ in the previous thin-wall analysis is a thickened torus (with the $`B`$-region in its interior) resembling a closed string. Numerical experiments indicate that, roughly speaking, the Q-effect supports the thickness of the string, but not its length; the string has a tension which causes its length to shrink. So the configuration collapses, and there is no stationary Hopf soliton in this system. In principle it remains a possibility that for some potential function $`V`$, some value of $`Q`$, and some nonzero value of the Hopf number $`N`$, there might exist a stationary solution; but this seems rather unlikely.

Another way of viewing the situation is as follows. The Q-mechanism provides a lower bound on the quantity $`I`$ (since $`\omega `$ is bounded above, and $`Q=I\omega `$ is fixed); this in turn means that the volume of the soliton is bounded below. But surface tension then acts to make the soliton spherical. So we are led to the following conjecture: any stationary Q-ball (whether topological or not) with $`\omega <m`$, in $`D`$ spatial dimensions, has O($`D`$) symmetry. (In the Bogomolny case, where $`\omega =m`$, rotational symmetry is not essential .) The instability of Hopf Q-solitons is an immediate consequence of this conjecture, since O(3) symmetry implies that the Hopf number is zero.

## 3 The O(4) Sigma Model in 3+1 Dimensions

In this section, we investigate the analogous problem for the O(4) sigma model in three space dimensions. The details of the system are similar to those of the previous section. The field takes values on $`S^3`$, and is represented as a unit 4-vector $`\stackrel{}{\varphi }=(\varphi _0,\varphi _1,\varphi _2,\varphi _3)`$. The Lagrangian is (1), as before, but with the potential being allowed to depend on $`\varphi _0`$ and $`\varphi _3`$: $`V=V(\varphi _0,\varphi _3)`$. In general, this breaks the global O(4) symmetry to O(2), the subgroup which rotates $`\varphi _1`$ and $`\varphi _2`$, leaving $`\varphi _0`$ and $`\varphi _3`$ fixed. The expressions for $`E_d`$, $`E_k`$ and $`E_p`$ are the same as before, with $`D=3`$, except that in $`E_d`$ there is an extra term involving $`(_j\psi _0)^2`$.

We take the boundary condition to be $`\psi _01`$ as $`r\mathrm{}`$ in $`𝐑^3`$. The mass $`m`$ is defined by $`V(\psi _0,\psi _3)m^2(1\psi _0){\scriptscriptstyle \frac{1}{2}}m^2(|\psi |^2+\psi _3^2)`$ for $`\psi _01`$; the function $`V_0`$ and the constant $`K`$ are defined as before. The virial relation from (4) with $`D=3`$ holds as before, as does the inequality $`Km^2<\omega ^2m^2`$.

In , this system was studied, with the potential $`V_0=2(1\psi _0)`$. Since in that case we have $`K=1`$, no soliton solution can exist. The authors of  reach this conclusion for the topologically trivial case $`N=0`$; they report numerical evidence for a non-trivial solution with $`N=1`$, but this cannot be correct.

In order to allow the possibility of non-trivial solutions, we need a potential $`V_0`$ which has $`K<1`$; for what follows, we shall take $`V_0={\scriptscriptstyle \frac{1}{2}}(1\psi _0^4)`$ as in the previous section. Then there are solutions, but it appears that they all have trivial topology ($`N=0`$). One way to see what happens is to consider the thin-wall limit, where $`m^2Q1`$. So the energy is approximately $`E{\scriptscriptstyle \frac{1}{4}}m^2(1\psi _0^4)+{\scriptscriptstyle \frac{1}{2}}Q^2/|\psi |^2`$; and the corresponding variational equations are

$$\xi \psi _0m^2\psi _0^3=0,\xi \psi \omega ^2\psi =0,\xi \psi _3=0,$$
(7)
where $`\xi :=m^2\psi _0^4+\omega ^2|\psi |^2`$ (the $`\xi `$-term arises from enforcing the constraint $`\stackrel{}{\varphi }\stackrel{}{\varphi }=1`$). These equations (7) have a number of solutions, namely:

* $`\psi _0=0`$, $`\psi =0`$, $`\psi _3=\pm 1`$;
* $`\psi =0`$, $`\psi _3=0`$, $`\psi _0=\pm 1`$;
* $`\psi _0=\psi _3=0`$, $`|\psi |^2=1`$;
* $`\psi _3=0`$, $`\psi _0=\pm \omega /m`$, $`|\psi |^2=1\omega ^2/m^2`$.

So to construct a mimimum-energy configuration, we must partition space $`𝐑^3`$ into regions (separated by infinitesimally-thin walls), on each of which one of these relations holds. It is clear that regions on which $`\psi _30`$ contribute only to $`E_p`$, and that we can reduce the total energy by instead setting $`\psi _3=0`$, $`\psi _0=1`$ on these regions. In other words, $`\psi _3`$ ‘collapses’ to zero, and is replaced by $`\psi _0`$.

This is exactly what one sees in numerical simulations. For example, we may start with the O(3)-symmetric ‘hedgehog’ ansatz

$$\psi _0+\mathrm{i}\psi _j\sigma _j=\mathrm{exp}\left[\mathrm{i}f(r)x^j\sigma _j/r\right];$$
(8)
here $`\sigma _j`$ denotes the Pauli matrices, and the profile function $`f(r)`$ satisfies the usual boundary conditions $`f(0)=\pi `$, $`f(\mathrm{})=0`$. The winding number is $`N=1`$. Note that $`\psi =\psi _1+\mathrm{i}\psi _2`$ vanishes on the $`x^3`$-axis, and that $`\psi _0(0)=1`$. If we now relax the configuration by flowing down the energy gradient, then $`\psi _3`$ approaches zero everywhere except at the single point $`r=0`$; in other words, there is no continuous minimum in this topological class. By contrast, there is a smooth minimum which has $`\psi _30`$: this is topologically trivial, and is essentially a standard (nontopological) Q-ball.

## 4 Concluding Remarks

We begin with a few remarks on the similarities with stationary topological soliton solutions of the Landau-Lifshitz equation

$$\frac{\stackrel{}{\varphi }}{t}=\stackrel{}{\varphi }\times \frac{\delta E}{\delta \stackrel{}{\varphi }}.$$
(9)
Here $`\stackrel{}{\varphi }`$ is a unit 3-vector representing the local orientation of magnetization, and the energy $`E`$ is given by

$$E=\left[{\scriptscriptstyle \frac{1}{2}}(_j\stackrel{}{\varphi })(_j\stackrel{}{\varphi })+U(\varphi _3)\right]d^Dx.$$
(10)
A typical choice for the function $`U`$ is $`U=A(1\varphi _3^2)`$, where $`A`$ is a constant; this corresponds to an easy-axis anisotropy. The boundary condition is $`\varphi _31`$ as $`r\mathrm{}`$. The total magnetization

$$M=(1\varphi _3)d^Dx$$
(11)
is a conserved quantity. It is clear from scaling that the only static solutions are the Belavin-Polyakov solitons  in spatial dimension $`D=2`$, with $`U0`$. But by allowing time dependence, more solutions are possible. In particular, we may allow periodic time dependence, and look for stationary solutions such that

$$\varphi _1+\mathrm{i}\varphi _2=\psi (x^j)\mathrm{e}^{\mathrm{i}\nu t},\varphi _3=\psi _3(x^j).$$
(12)
With this ansatz, the Landau-Lifshitz equation (9) is equivalent to

$$\frac{\delta }{\delta \stackrel{}{\varphi }}\left(E{\scriptscriptstyle \frac{1}{2}}\nu M\right)=0.$$
(13)
We may think of the solutions as critical points of $`E`$, subject to the constraint that $`M`$ has a given value . The simplest example occurs if $`U={\scriptscriptstyle \frac{1}{2}}\nu (1\psi _3)`$, for then the functional appearing in (13) consists of only the gradient term, and the Belavin-Polyakov solitons are solutions (in $`D=2`$); these correspond, of course, to the Bogomolny-type Q-lump solitons . The analysis of stationary topological Landau-Lifshitz solitons leads to rather similar results as for Q-solitons (although the dynamics of moving solitons is quite different). In $`D=2`$, there are topological solutions (called magnetic bubbles — see  for a review); a single soliton is pinned in space, and cannot move . In $`D=3`$, on the other hand, there are no stationary Hopf solitons ; however, such solitons can be stabilized by allowing them to move at constant velocity .

Returning to sigma-model dynamics, we have seen that the Q-mechanism stabilizes topological solitons in $`D=2`$ spatial dimensions, but not in $`D=3`$. Stabilizing vortex rings (Hopf textures) in $`D=3`$ is particularly difficult, since there are two length-scales (the length of the loop and its width), each of which has to be fixed. The Q-effect can stabilize the latter, but not the former. One does get stable loops in systems with a Skyrme term , and also in systems with a magnetic field sufficiently strongly coupled to the scalar field (minimal coupling is not enough) . But in the basic versions of each of these systems, there is only one length-scale; and so the length of the loop is of the same order as (and only slightly greater than) its thickness. It remains an open question as to whether there is a system admitting a stable Hopf soliton in which the two length-scales are significantly different.

A sigma-model soliton in $`D=2`$ can be thought of as a (straight) sigma-model string in three spatial dimensions. So, for example, the Q-stabilized solitons discussed in this paper may find application as cosmic strings. Given an appropriate potential $`V`$, long strings with internal rotational energy will be stable, although closed loops will eventually shrink and decay. In this connection, it is worth recalling that, on a cosmological scale, both the width and the length of sigma-model strings are stabilized by cosmological expansion; but stabilizing the length requires a greater rate of expansion than stabilizing the width .


text:
# Diffusion of Asymmetric Swimmers

## Abstract

Particles moving along curved trajectories will diffuse if the curvature fluctuates sufficiently in either magnitude or orientation. We consider particles moving at a constant speed with either a fixed or with a Gaussian distributed curvature magnitude. At small speeds the diffusivity is independent of the speed. At larger particle speeds, the diffusivity depends on the speed through a novel exponent. We apply our results to intracellular transport of vesicles. In sharp contrast to thermal diffusion, the effective diffusivity increases with vesicle size and so may provide an effective means of intracellular transport.

The thermal Stokes-Einstein diffusivity of a sphere decreases as the particle radius $`R`$ increases Berg93 . For this reason, while diffusive transport is used for individual molecules within living cells Alberts2002 , larger objects such as vesicles and pathogens often use active means of transport. While many intracellular vesicles appear to be transported by molecular motors directed along existing cytoskeletal tracks motorvesicle ; Alberts2002 , undirected actin-polymerization mediated vesicle transport has been reported in some endosomes, lysosomes, other endogenous vesicles, and phagosomes vesicles ; reviews . Active transport is also observed in the actin-polymerization-ratchet motility of certain bacteria reviews ; rickettsiae  and virus particles vaccinia  within host cells. It is important to characterize the transport properties of vesicles that are not moving along pre-existing cytoskeletal tracks.

Existing discussions of the motion of actively propelled microscopic particles, or “swimmers”, assumes that in the absence of thermal fluctuations particles would move in straight trajectories Lovely75 ; Berg93 . Thermal rotational diffusion will then randomly re-orient the trajectory Berg93 , so that over long times diffusive transport will be observed. However in actin-polymerization based motility, particles appear to be attached to their long actin tails Gerbal2000  which in turn are embedded in the cytoskeleton Theriot92 . While thermal fluctuations will thereby be severely reduced, the actin-polymerization itself is a stochastic process with its own fluctuations Alberts2002 ; Mogilner96 . These intrinsic fluctuations can explain the observed curved trajectories, as well as the variation of the curvature over time Rutenberg2001 . The diffusivity of such asymmetrically moving particles has not been previously explored.

In this letter, we study asymmetric swimmers that would move in perfect circles in the absence of fluctuations. We examine both a “broken swimmer” with a fixed curvature magnitude and an axis of curvature that is re-oriented by fluctuations (rotating curvature, RC), and a “microscopic swimmer” with a normally-distributed curvature that is spontaneously generated by fluctuations (Gaussian curvature, GC). In both of these systems, fluctuations lead to diffusion at long times. We use computer simulations to measure the diffusivity of these systems as a function of the root-mean-squared curvature $`K_0`$, the particle speed $`v`$, and the timescale characterizing the curvature dynamics $`\tau `$.

We obtain some exact results from polymer systems, where each polymer configuration represents a possible particle trajectory. Indeed, a broken swimmer with a fixed curvature magnitude in $`d=3`$ is exactly analogous to the hindered jointed chain discussed by Flory Flory89  and we thereby recover the entire scaling function exactly. In that case, the diffusivity is independent of particle speed $`v`$. For Gaussian curvatures and for systems in restricted geometries ($`d=2`$), the polymer analogy gives us the diffusivity only in the limit of slow speeds. At larger speeds, our simulations over $`5`$ decades of speed show that diffusivity depends on particle speed with a non-trivial exponent $`\lambda `$. The diffusivity appears to be dominated by the occasional long straight segments of trajectory that occur when the curvature is small. Scaling arguments based on this observation are consistent with the measured exponent $`\lambda _{2d}=0.98\pm 0.02`$ in $`d=2`$, but do not recover our measured exponent $`\lambda _{3d}=0.71\pm 0.01`$ in $`d=3`$.

A curved path has a curvature magnitude $`K1/R`$, where $`R`$ is the instantaneous radius of curvature. If we describe a particle trajectory by a position $`𝐫(t)`$, then the vector curvature is defined by the cross-product $`𝐊𝐯\times \dot{𝐯}/v^3`$, where $`v=|𝐯|`$ is the speed and the dot $`\dot{}`$ indicates a time-derivative. For uniform motion around a circle, $`R`$ is the radius of the circle, and $`𝐊`$ is oriented perpendicular to the circle along the axis. We consider particles moving at a constant speed and with an instantaneous curvature $`𝐊`$, so that $`\dot{𝐫}=𝐯`$ and $`\dot{𝐯}=v𝐯\times 𝐊`$. For “rotating curvature” dynamics (RC) we fix the curvature magnitude $`|𝐊|=K_0`$ but allow the curvature to randomly rotate around the direction of motion:

$$\dot{𝐊}_{RC}=\xi \widehat{𝐯}\times 𝐊$$
(1)
where the unit-vector $`\widehat{𝐯}=𝐯/v`$, the Gaussian noise $`\xi `$ has zero mean, and $`\xi (t)\xi (t^{})=2\delta (tt^{})/\tau `$ with a characteristic timescale $`\tau `$. This represents the simplest description of a mesoscopic swimmer that has a “locked-in” curvature due to, e.g., an asymmetric shape. For “Gaussian curvature” dynamics (GC) the curvature magnitude changes as well:

$$\dot{𝐊}_{GC}=𝐊/\tau +𝝃$$
(2)
where the noise $`𝝃`$ is perpendicular to $`𝐯`$ with zero mean and $`𝝃(t)𝝃(t^{})=\delta (tt^{})K_0^2/\tau `$, such that $`𝐊^2=K_0^2`$. This represents the simplest description of a microscopic swimmer “trying” to swim in a straight line subject to intrinsic fluctuations in the motion. The resulting curvatures are Gaussian distributed in each component. For particles restricted to two-dimensions with either RC or GC dynamics, we only use the normal ($`\widehat{z}`$) component of the vector-curvature to update the velocity within the plane, i.e. $`\dot{𝐯}=v𝐯\times \widehat{𝐳}K_z`$ in $`d=2`$.

There are two natural timescales. We explicitly introduce $`\tau `$, which controls the noise correlation and so sets the timescale over which the curvature changes. There is also the inverse of the angular rotation rate, $`t_c1/(vK_0)`$. Diffusion will only be observed for elapsed times $`t`$ much greater than any other timescale in the system, i.e. $`tt_c`$ and $`t\tau `$. The diffusivity of a particle is given by $`Dr^2/(2dt)`$ in the limit as the elapsed time $`t\mathrm{}`$, in spatial dimension $`d`$.

A polymer chain with fixed bond lengths ($`\mathrm{}`$) and angles ($`\theta _f`$), and with independent bond rotation potentials ($`V(\varphi _f)`$) Flory89  is statistically identical to the continuous RC trajectory in $`3d`$ if for a discrete time-step $`\mathrm{\Delta }t`$ we take $`\mathrm{}=v\mathrm{\Delta }t`$. The end-to-end distance for a long $`n`$-bond polymer is $`r^2=n\mathrm{}^2C_n`$. The correspondence is complete as the elapsed time $`t=n\mathrm{\Delta }t\mathrm{}`$. The bond and dihedral angles determine $`C_{\mathrm{}}=(1+\mathrm{cos}\theta _f)(1+\mathrm{cos}\varphi _f)/[(1\mathrm{cos}\theta _f)(1\mathrm{cos}\varphi _f)]`$ Flory89 . Swimmers follow continuous paths, so we take the limit of small $`\mathrm{\Delta }t`$ and fix the polymer rotation angle from the curvature in that limit by $`\theta _f=K_0v\mathrm{\Delta }t`$, and rotate the curvature by $`<\varphi _f^2>=2\mathrm{\Delta }t/\tau `$ in agreement with Eqn. 1. In the limit $`\mathrm{\Delta }t0`$ we recover the exact result $`D=1/(3K_0^2\tau )`$ in $`d=3`$. Remarkably, $`D`$ is independent of $`v`$.

For a GC trajectory in $`d=3`$, there is no obvious polymer analogy since the curvature magnitude evolves with time. In the limit of $`\tau 0`$ however, the curvature is independently Gaussian distributed at every point along the trajectory and the diffusivity can be extracted from the “worm-like chain” polymer model originally solved by Kratky and Porod Doi . The result is $`D=1/(3K_0^2\tau )`$ in the limit of small $`\tau `$. Note that the diffusivity diverges as $`1/\tau `$ so this is the leading asymptotic dependence for small $`\tau `$. We obtain the same diffusivity for both RC and GC dynamics for small $`\tau `$.

We use these exact results to define natural dimensionless scaling functions for the diffusivity of microscopic swimmers:

$$\stackrel{~}{D}_{Rd,Gd}(\stackrel{~}{v})DK_0^2\tau $$
(3)
where the index $`Rd`$ or $`Gd`$ indicates both the dynamics (RC or GC) and the spatial dimensionality $`d`$, and

$$\stackrel{~}{v}vK_0\tau $$
(4)
is a dimensionless speed. In terms of these scaling functions we have $`\stackrel{~}{D}_{R3}(\stackrel{~}{v})=\stackrel{~}{D}_{G3}(0)=1/3`$. The same Kratky-Porod approach in $`d=2`$ gives $`\stackrel{~}{D}_{R2}(0)=\stackrel{~}{D}_{G2}(0)=1`$.

We have simulated the trajectories of large-numbers of independent particles with RC and with GC dynamics. For fixed $`v`$ and $`K_0`$, we varied $`\tau `$ to explore the scaled velocity $`\stackrel{~}{v}vK_0\tau `$ over $`5`$ orders of magnitude. For each $`\stackrel{~}{v}`$, we averaged over the trajectories of at least $`1000`$ particles. We explicitly integrated the dynamical equations using a simple Euler update with a small timestep $`\mathrm{\Delta }t`$. In all cases $`t\tau \mathrm{\Delta }t`$ and $`tt_c\mathrm{\Delta }t`$, with separation of timescales by factors of $`10100`$. Systematic errors due to $`\mathrm{\Delta }t`$ and $`t`$ are below our noise levels, and statistical errors (when not shown) are smaller than the size of our plotted points. Consistently, our numerical results agree with all exact results from the polymer analogy. We illustrate the trajectories that we observe in $`d=2`$ in Fig. 1, with both small and large scaled speeds $`\stackrel{~}{v}`$. In both cases the curvature $`K_0=1`$, but particles only complete loops at large $`\stackrel{~}{v}`$. Qualitatively similar trajectories are seen in $`d=3`$ with GC curvature dynamics.

In $`d=2`$, shown in Fig. 2, both rotating curvature (open circles) and Gaussian curvature (filled circles) approach their asymptotic value of $`\stackrel{~}{D}_{R2}(0)=\stackrel{~}{D}_{G2}(0)=1`$ at small $`\stackrel{~}{v}`$. At $`\stackrel{~}{v}1`$ there is a sharp cross-over to a large-$`\stackrel{~}{v}`$ power-law regime, characterized by an exponent $`\lambda _{2d}`$ where $`\stackrel{~}{D}_{R2}\stackrel{~}{D}_{G2}\stackrel{~}{v}^{\lambda _{2d}}`$ for large $`\stackrel{~}{v}`$. We show the effective exponents $`\lambda _{eff}\mathrm{\Delta }\mathrm{log}(DK_0^2\tau )/\mathrm{\Delta }\mathrm{log}\stackrel{~}{v}`$ between consecutive points in the inset of Fig. 2, as well as the best fit exponent $`\lambda _{2d}=0.98\pm 0.02`$. We fit $`\lambda _{2d}`$ from the large-$`\stackrel{~}{v}`$ GC data only, due to the systematic cross-over remaining in the RC data even at large $`\stackrel{~}{v}`$.

Simulations in $`d=3`$ with rotating curvature (RC) dynamics leads to a diffusivity in excellent agreement with the exact result from polymer physics, $`\stackrel{~}{D}_{R3}=1/3`$, as shown by open circles in Fig. 3. Gaussian curvature dynamics ($`\stackrel{~}{D}_{G3}`$, filled circles) has the same behavior for small $`\stackrel{~}{v}`$, but exhibits a sharp crossover at $`\stackrel{~}{v}1`$ to a power-law regime $`\stackrel{~}{D}_{G3}\stackrel{~}{v}^{\lambda _{3d}}`$ for large $`\stackrel{~}{v}`$. We find the best-fit exponent is $`\lambda _{3d}=0.71\pm 0.01`$, as shown by the solid line in the inset of Fig. 3. Because $`\lambda _{3d}<1`$, this scaling curve may be used to uniquely identify the dynamical timescale $`\tau `$ if $`D`$, $`K_0`$, and $`v`$ are measured experimentally.

Is there a simple way of understanding the asymptotic behavior of $`\stackrel{~}{D}`$? For RC dynamics in $`d=3`$ the instantaneous curvature does not change in magnitude even while the curvature axis wanders. The particle will go in a circular trajectory, not contributing to diffusivity, until the curvature axis wanders significantly. The result is a random walk with step size given by the radius of curvature $`\mathrm{\Delta }r1/K_0`$ and an interval between steps of $`\tau `$, leading to $`D1/(K_0^2\tau )`$. This qualitatively explains why the exact result $`\stackrel{~}{D}_{R3}=1/3`$ is independent of $`\stackrel{~}{v}`$.

It is more difficult to understand the $`\stackrel{~}{D}\stackrel{~}{v}^\lambda `$ behavior for large $`\stackrel{~}{v}`$ in the other systems. We start with a simple scaling argument based on the assumption that the relatively straight segments shown in Fig. 1 b) dominate the diffusivity. The interval between periods of small curvature should be on order the autocorrelation time $`\tau `$. The length $`\mathrm{\Delta }r`$ of the straight segments are determined by how long the interval of small curvature lasts, $`\mathrm{\Delta }t`$, since $`\mathrm{\Delta }rv\mathrm{\Delta }t`$. For the segment to be straight, the curvature must be less than the inverse length, i.e. $`K_{max}1/\mathrm{\Delta }r`$. The fraction of the time we have small curvature below $`K_{max}`$ in magnitude should be proportional to the probability of having curvature below $`K_{max}`$. In $`d=2`$ only the normal component of curvature affects the dynamics, so that $`P(K)\mathrm{const}`$ for $`KK_0`$. This applies both to GC and RC. We therefore expect $`\mathrm{\Delta }t\tau K_{max}/K_0`$. We maximize $`K_{max}`$ to maximize the contribution to $`D\mathrm{\Delta }r^2/\tau `$ and find $`\stackrel{~}{D}_{G2}\stackrel{~}{D}_{R2}\stackrel{~}{v}`$ as $`\stackrel{~}{v}\mathrm{}`$. This indicates that $`\lambda _{2d}=1`$, which is consistent with our best-fit value $`\lambda _{2d}=0.98\pm 0.02`$. However, in $`d=3`$ for GC dynamics the same argument leads to $`\lambda _{3d}=2/3`$ since two Gaussian distributed components of the curvature gives $`P_<(K_{max})=_0^{K_{max}}𝑑KP(K)K_{max}^2/K_0^2`$ for $`K_{max}K_0`$. This is inconsistent with our measured value of $`\lambda _{3d}=0.71\pm 0.01`$, with a significant $`4\sigma `$ variation.

At what radius $`R_c`$ does a small spherical particle achieve a higher diffusivity by actively swimming, as compared to passive thermal diffusion characterized by $`D_T=k_BT/(6\pi \eta R)`$ Berg93 ? We can answer this question within the context of actin-polymerization based motility of small intracellular particles, since the size dependence of $`K_0`$, $`v`$, and $`\tau `$ is known, at least approximately. With the approximation that $`n`$ propulsive actin filaments are randomly distributed over a particle of size $`R`$, the curvature of the trajectory will be $`K_01/(R\sqrt{n})`$ Rutenberg2001 . With a size-independent surface-density of filaments we obtain $`K_0A/R^2`$, with a constant of proportionality $`A`$. By observations of Listeria monocytogenes we estimate $`A0.1\mu m`$ Rutenberg2001 . We also conservatively assume size-independent values for cytoplasmic viscosity $`\eta 3Pas`$ Berg93 ; Rutenberg2001 , speed $`v0.1\mu m/s`$, and autocorrelation decay time $`\tau 100s`$ Sechi97 . We find that the micron-scale bacterium L. monocytogenes has $`\stackrel{~}{v}1`$, so that smaller particles will have $`\stackrel{~}{v}>1`$. Using the large $`\stackrel{~}{v}`$ asymptotic behavior of $`\stackrel{~}{D}_{G3}`$ shown in Fig. 3, $`D0.41\stackrel{~}{v}^{\lambda _{3d}}/(K_0^2\tau )`$, and the size-dependence $`K_0A/R^2`$, we obtain

$$D_{G3}R^{42\lambda _{3d}}v^{\lambda _{3d}}/(A^{2\lambda _{3d}}\tau ^{1\lambda _{3d}}),$$
(5)
with a measured $`\lambda _{3d}=0.71\pm 0.01`$. In dramatic contrast to thermal diffusion, $`D`$ increases with increasing particle size. Comparing with $`D_T`$ we find that for all sizes above $`R_c80nm`$ a particle will have a higher diffusivity by actively swimming by the actin-polymerization mechanism than by passive thermal diffusion. Provocatively, this is in the middle of the vesicle-size distribution seen in neural systems vesiclesize .

Our treatment of microscopic swimmers has ignored thermal fluctuations. A “rocket” traveling straight at speed $`v`$ that is re-oriented only by thermal effects will have $`D_u=4\pi \eta R^3v^2/(3k_BT)`$ Berg93 . In comparison with our results for $`D`$, we find that $`D<D_u`$ for particles larger than $`R_u0.07nm`$. For actin-polymerization based motility, intrinsic fluctuations appear to dominate thermal fluctuations at the particle sizes where active transport is advantageous.

In summary, we find that diffusivities of asymmetric microscopic swimmers depends on whether the swimmers are restricted to $`2d`$ or $`3d`$, and whether they have fixed asymmetries (RC) or the asymmetries are spontaneously generated (GC). Diffusivities are independent of particle speed at low speeds, in agreement with analogous polymer systems. At higher speeds an anomalously large diffusivity is observed that depends on the particle speed by $`\stackrel{~}{v}^\lambda `$ where $`\lambda _{2d}=0.98\pm 0.02`$, in agreement with a scaling argument for $`\lambda _{2d}=1`$. However $`\lambda _{3d}=0.71\pm 0.01`$, which significantly differs from our simple scaling result in $`d=3`$. We apply our results to intracellular bacteria, virus particles, and vesicles that move via actin-polymerization. We find that diffusivities due to asymmetric swimming exceed thermal diffusivities for particles larger than approximately $`80nm`$. As a result asymmetric swimming may provide a viable intracellular transport mechanism even for vesicle-sized particles. We find that for the relevant dynamics (GC in $`d=3`$), diffusivities should increase with particle size, speed, and filament turnover rate, and also with smaller curvatures for a given size. It is interesting that the bacterium Rickettsiae rickettsii exhibits actin-polymerization intracellular motility with smaller intra-cellular speeds but straighter trajectories rickettsiae ; Rutenberg2001  — raising the question of whether maximal diffusivity is selected for in this or other biological systems.

This work was supported financially by an NSERC discovery grant. C. Montgomery would also like to acknowledge support from an NSERC USRA.