Orientational ordering and phase behaviour of a binary mixture of hard spheres and hard spherocylinders

We study structure and fluid-phase behaviour of a binary mixture of hard spheres (HSs) and hard spherocylinders (HSCs) in isotropic and nematic states using the $NP_nAT$ ensemble Monte Carlo (MC) method in which a normal pressure tensor component is fixed in a system confined between two hard walls. The method allows one to estimate the location of the isotropic-nematic phase transition and to observe the asymmetry in the composition between the coexisting phases, with the expected increase of the HSC concentration in the nematic phase. This is in stark contrast with the previously reported MC simulations where a conventional isotropic $NPT$ ensemble was used. We further compare the simulation results with the theoretical predictions of two analytic theories that extend the original Parsons-Lee theory using the one-fluid and the many-fluid approximation [Malijevsk\'y {\it at al} J. Chem. Phys. \textbf{129}, 144504 (2008)]. In the one-fluid version of the theory the properties of the mixture are mapped on an effective one-component HS system while in the many-fluid theory the components of the mixtures are represented as separate effective HS particles. The comparison reveals that both the one- and the many-fluid approaches provide a reasonably accurate quantitative description of the mixture including the predictions of the isotropic-nematic phase boundary and degree of orientational order of the HSC-HS mixtures.

NP n AT ensemble Monte Carlo (MC) method in which a normal pressure tensor component is fixed in a system confined between two hard walls. The method allows one to estimate the location of the isotropic-nematic phase transition and to observe the asymmetry in the composition between the coexisting phases, with the expected increase of the HSC concentration in the nematic phase. This is in stark contrast with the previously reported MC simulations where a conventional isotropic NP T

I. INTRODUCTION
Advanced materials formed by the self-assembly of non-spherical building blocks has experienced an unprecedented growth due to recent advances in experimental techniques to create nano and colloidal particles with almost any imaginable shape. [1][2][3][4] Functional materials can be engineered by tailoring the properties of the individual building blocks. 5,6 Colloidal particles are particularly attractive as building blocks for the design of mesoscale materials, which are difficult to fabricate by using chemical synthesis, as their interactions can also be modulated by modifying both the surface chemistry of the particles as well and the properties of the solvent media. 6 It is possible to tune the interactions of either stericstabilised or charge-stabilised colloidal particles as nearly as hard body-like by matching the index of refraction of both particles and solvent. 7 Moreover, the self-assembly of these systems can also be controlled by the aid of external forces such as magnetic and electric fields, gravity 8 , and even the use of geometrical confinement. [9][10][11][12][13][14][15][16][17] We refer to these processes in general as directed self-assembly. 18 Anisotropic particles can exhibit many fascinating structures in bulk and confinement including crystals, plastic crystals, and liquid crystal (LC) phases. [19][20][21] LC phases in rod-like particles, for example, are observed in many natural and anthropogenic systems. Examples include suspensions of colloidal particles such as vandium pentoxide (V 2 O 5 ) 22 , and Gibbsite (Al(OH) 3 ) 23,24 , carbon nanotubes 25 , and some biological systems such as protein fibers 26 , tobacco mosaic virus 27 , fd -virus [28][29][30][31][32][33][34] , polypeptide solutions 35,36 , and DNA 37 . Over the years, simple but non-trivial hard-core models have been used to study the formation of LC phases. 38 These models have played an important role to understand the behaviour of real systems. In particular, the hard-spherocylinder (HSC) has been used as a standard model to describe the LC behaviour of rod-like colloidal particles. The HSC model consist of a cylinder of length L and diameter D capped at each end by a hemisphere of diameter D, and it is shown in Fig. 1(a). Depending of the aspect ratio of the model, corresponding to the ratio L/D, suspensions of HSCs can exhibit the formation of isotropic, nematic, smectic, and solid phases. This rich phase behaviour has been confirmed by extensive computer simulations. [39][40][41][42][43] Despite our profound knowledge on the phase behaviour of rod-like particles, our understanding of the phase behaviour of mixtures of anisotropic colloidal particles is still limited due to the large parameter space that has to be explored, i.e., different combinations of concentrations, shapes, and sizes of the components, as well as different thermodynamic conditions. Experimental studies for mixtures of rod-like and spherical colloids have been reported. [44][45][46][47][48][49][50][51] From the modelling perspective, binary mixtures of HSC particles have been studied using computer simulations including rod-rod 31,32,52,53 , rod-disc 54 Onsager's method. Parsons 81 proposed an approximation to decouple the orientational and translational degrees of freedom by mapping the properties of the rods to those of a reference HS system. Lee 82,83 approached the problem in a different way by introducing a scaling relation between virial coefficients of anisotropic particles and HS reference. Following two separate routes, Parsons and Lee reached the same expression for the free energy functional which is commonly known as the Parsons-Lee (PL) theory 24,42,74,84,85 . A straightforward extension of PL theory 62 to the mixtures is the one-fluid approximation whereby one maps the mixture on to an effective one-component HS system. A decoupling approximation is used in the PL approach in which the system is represented as the effective hard sphere of the same diameter while any information about the geometry of the LC particles is included in the term of the factorized excluded volumes. In order to improve the PL treatment for mixtures, a many-fluid (MF) approach has been proposed 86  The focus of our current work is the isotropic-nematic phase behaviour of a HSC-HS mixture. Previous reports of the ordering in the HSC-HS binary system have been presented including direct simulation 60,87,88 and theoretical 62,89,90 studies. The work of Cuetos and co-workers is of particular relevance: the one-fluid PL approach 62 was used to study the isotropic-nematic phase diagram of the HSC-HS system characterized by rods of various lengths and diameters; comparisons where made with NP T Monte Carlo simulations 60 employed to investigate the phase diagram and fluid structure of the mixtures. It is worth noting that in the NP T ensemble the system composition remains constant overall, which will lead to an inadequate description of the phase boundary as one enters metastable states which would otherwise phase separate into phases of distinct compositions.
The purpose of our current work is twofold. First, the many-fluid Parsons theory is used to describe the HS-HSC system and comparisons are made with the one-fluid PL approach.
It should be noted that in a one-component case both approaches reduce to the standard PL theory. Second, we present new Monte Carlo simulation results for the mixture. The local density (packing fraction), local composition, and orientational distributions are determined during the simulations to estimate the locations of the isotropic-nematic transitions of the mixture at various compositions in order to make a proper test of the accuracy of the two theories.

II. THEORY OF NEMATIC PHASE IN MIXTURES OF HARD PARTICLES
In this section, the main steps leading to a formulation with both one-fluid and manyfluid theories are briefly recalled; further details can be found in Ref. 86. Consider an ncomponent mixture system of N unaxial (cylindrically symmetrical) hard anisotropic bodies in a volume V at a temperature T . The free energy functional of the system can be expressed as a contribution from an ideal (entropy) term (F id ) and a residual (configurational) part (F res ): where β = 1/(k B T ) and k B is the Boltzmann constant; the temperature plays a trivial role in this case since only hard repulsive interactions between particles are considered. The ideal free energy accounts for the translational and orientational entropy and can be written as where ρ i , and V i is the de Broglie volume of each species, incorporating the translational and rotational kinetic contributions of the ideal isotropic state. With the introduction of single-particle orientational distribution function f i ( ω), the orientational entropy term σ[f i ] can be expressed as an integration over all orientations ω of a single particle: For the residual part, Onsager's original expression 71 is equivalent to truncating the virial expansion at second-virial level. At higher densities, however, the many-body correlations become progressively more and more important. Following the Parsons approach 81 , we can include higher-body contributions in an approximate manner. Assuming a pairwise additive hard interaction u ij (r kl , ω k , ω l ) between particle k of i-th component and particle l of j-th component, the pressure of a fluid mixture of n components can be written in the virial form as 91 where l = k is used to avoid self interactions and · · · represents the ensemble average. In the canonical ensemble, where the configurational partition function Z is defined as and the total configurational energy is given by The canonical pair distribution function is defined as where δ ij is the Kronecker delta. On integrating Equation (5), the expression for pressure can be written in a compact form: Following the Parsons approach, the interparticle separation r 12 is given in terms of the contact distance σ ij ( r 12 , ω 1 , ω 2 ) by defining a scaled distance y ij = r 12 /σ ij ( r 12 , ω 1 , ω 2 ). The scaled distance does not explicitly depend on the orientations of the two particles and y ij = 1 corresponds to contact value. Using the definition of y ij , the pair distribution function (cf. Equation (8)) can be expressed as a function of scaled distance y ij , i.e., g ij = g ij (y), which decouples the positional and orientational dependencies. In this way, a complicated pair potential u ij is mapped onto the spherically symmetrical hard-sphere potential: and the expression for the pressure becomes where the excluded volume between a pair of particles is V exc ij ( ω 1 , ω 2 ) = 1 3 dr 12 σ 3 ij (r 12 , ω 1 , ω 2 ) andr 12 = r 12 /r 12 . The form of hard repulsive pair interaction is a step function (cf. Equation (10)), thus βdu ij /dy ij = − exp(βu ij )δ(y ij − 1) (for example, see Ref. 92). Integrating over the scaled variable y ij and noting that u 1 + (y) = 0 when y = 1, we then obtain where g ij (1 + ) ≈ g HS ij (1 + ) has been approximated as the corresponding hard-sphere contact value of pair distribution function.
The residual free energy can then be obtained from the formal thermodynamic definition (∂F/∂V ) N T = −P by integrating Equation (12) over the volume: where G ij = ρ −1 ρ 0 dρ ′ g HS ij (1 + ). Onsager's second-virial theory can be recovered with G ij = 1 (i.e., g HS ij (1 + ) = 1) corresponding to the low-density limit. In this way, the theory of Parsons for a one-component fluid can be reformulated to describe a n-component mixtures of anisotropic bodies. As shown in Ref. 86, the standard "one-fluid" approach 62 corresponds to G ij = G PL , where G PL = ρ −1 ρ 0 dρ ′ g HS CS (1 + ), given in terms of the Carnahan-Starling form of the radial distribution function at contact 93,94 , with η = n i=1 ρ i V m,i , and V m,i is the volume of the i-th species. The PL residual free energy can then be expressed as Alternatively, in developing the many-fluid theory proposed in Ref. 86 one treats the size (volume) of each species individually, i.e., An expression for the contact value of the distribution function for hard-sphere mixture is given by Boublik 95 : where the moments of the density are defined as ζ α = (π/6) n i=1 ρ i σ α ii , α = 0, 1, 2, 3. Combining Equations (13) and (17) and noting the separate definition of G ij for each i − j pair, one obtains the many-fluid Parsons (MFP) form of the residual free energy F res,MFP .
In the one-component limit the contact value of the radial distribution function of the HS mixture (Equation (17)) reduces to the Carnahan-Starling expression (cf. Equation (14)). Thereby, the MFP approach yields same descriptions as PL treatment for the pure-component systems. In the standard extension of the PL theory to mixtures one therefore adopts a van der Waals one-fluid (VDW1) approximation using an equivalent hard-sphere system with the effective diameter given by the VDW1 mixing rule to represent the anisotropic mixtures. In contrast to the PL approach, each component is represented as a separate effective hard-sphere component, so that the excluded volume between a pair of i-th component and j-th component is weighted by the corresponding contact value of the HS mixture, g HS ij . The equilibrium free energy of the system is determined from a functional variation with respect to the orientational distribution function f i ( ω) of each component which leads to an integral equation for f i ( ω). The set of integral equations are solved numerically using an iterative procedure, details of which can be found in Ref. 86.
In this work, we assess the adequacy of many-body theories such as the MFP for a binary mixture of hard spheres and hard spherocylinders. The models are depicted in Figure 1: the aspect ratio of the HSC is L/D = 5 and the diameter of the HS is taken to be the same as the diameter of the HSC, i.e., σ = D.
The excluded volumes corresponding to the HSC-HSC, HSC-HS and HS-HS interactions are given as where γ = arccos( ω 1 · ω 2 ) is the relative orientation of the two HSC particles. The total where x HS and x HSC are mole fractions of HS and HSC species, respectively. Since the HSC particles are the anisotropic component in the system, f ( ω) is used to describe the orientation distribution of the HSC rods which is related to the nematic order parameter S 2 of the system through In particular, S 2 = 0 corresponds to the isotropic state and S 2 = 1 for a perfectly-aligned nematic phase.

MIXTURES OF HARD SPHERES AND HARD SPHEROCYLINDERS
There are two common approaches to studying fluid-phase separation by molecular sim- The direct molecular simulation of the isotropic-nematic (I-N) phase transition in mixtures of hard spheres and hard spherocylinders is particularly challenging because of the very low interfacial tension between the two phases; for example, the I-N interfacial tension of a hard-core system of for thin disc-like particles has be estimated to be a few tenth of k B T in units of the particle's area 100 . As a consequence there is a very low energetic penalty associated with the deformation of the interface in such systems leading to large interfacial fluctuations; moreover, in the absence of an external field there is no resistance to the translation of a planar interface. The location of the bulk coexistence regions and the determination of the density and compositional profiles becomes a difficult task as a result. In order to break the symmetry of the system and reduce the effect of the interfacial fluctuations one can introduce an external field by placing the system within structureless hard walls; this corresponds to removing the periodicity in one dimension (say the z direction). An issue with this type of approach is that large systems have to be considered in order to study the true bulk phase behaviour and avoid capillary effects. By keeping the separation between the walls large compared to the dimensions of the particles, one can simulate the phase coexistence in the hard-core HS-HSC mixtures with minimal effect from the hard surfaces.
Alternatively, the phase behaviour can be simulated using a popular Gibbs ensemble 96,97 in which the coexisting phases are retained in separate boxes and coupled volume changes and particles exchanges between the boxes are undertaken to meet the requirements of mechanical and chemical equilibria. However, in the case of hard anisometric particles, the acceptance ratio for the insertion of anisotropic particles will be very low, particularly at the high densities of the dense anisotropic phases of interest, requiring an impracticably large number of trial insertions for a proper equilibration of the system 101 . A conventional simulation of the system within a single box will partially solve the problem since trial insertions of the particles are no longer required. There is however a complication with the simulation of bulk phase equilibria of mixtures with a single simulation cell: though the phase transition between the various states can be traced as for a pure component system, the overall composition remains fixed preventing the fractionation of the different species in the various phases.
In view of the aforementioned issues, we employ a less conventional NP n AT ensemble within a single cell where the component P n of the pressure tensor normal to the interface is kept constant, so that the condition of mechanical equilibria is satisfied within the entire simulation cell 92,102,103 . The advantage of simulating the phase separation of mixtures by simultaneously considering the coexisting phases and the interface in a single cell is that this will allow for inhomogeneities in both the density and the composition of the system.
By introducing an external field such as a hard surface one is able to examine both the bulk and interfacial regions of mixtures of hard core particles without constraining the density or composition of the individual bulk phases.
We perform constant normal-pressure Monte Carlo simulation (NP n AT -MC) for a system of N HSC = 1482 HSC particles of the aspect ratio L/D = 5 where the number of hard spheres is varied depending on mole fraction of the binary mixture x HS,tot = N HS /(N HS + N HSC ).
In this system, the intermolecular potential between any two particles is restricted to a pure repulsion. As shown in Figure 2, the simulation cell is a rectangular box of dimension L x = L y = 25D (corresponding to a fixed surface area A in the x-y plane of A = 625D 2 ) and L z varies according to the value set for P n . The parallel hard walls are positioned at z = 0 and z = L z and standard periodic boundary conditions are applied in x and y directions. Since a fixed normal pressure is imposed along the z axis, the system volume in our NP n AT -MC simulation is allowed to fluctuate by scaling the length of the z axis which moves the walls closer together or farther apart, while the system dimensions of the x and y axes and the x-y surface area are kept fixed.
The NP n AT -MC simulation of the HSC-HS mixture is performed for 5 × 10 6 cycles to equilibrate the system and 5 to 8 × and the local composition is then obtained as x HS (z) = ρ HS (z)/(ρ HS (z) + ρ HSC (z)). The packing fraction η b and composition x HS,b of the bulk phase is then determined to be the values corresponding to the uniform regions of the packing fraction and composition profiles in the centre of the simulation cell.
The nematic order parameter profile S 2 is obtained by determining the local nematic order parameter tensor Q(j) in each bin j: where N HSC,j is the number of HSC particles in the j-th bin and I is the unit tensor. On diagonalising the tensor Q(j), three eigenvalues can be obtained and the largest eigenvalue defines the local nematic order parameter S 2 (z) of the j-th bin. Special care is required in calculating the order parameter profile because of finite-size effects. It has been shown by Eppenga  parameter is to examine a larger system; for example, in the case of a system of 14820 rods (an order of magnitude larger than the system studied here) and 200 bins reduces the error in S 2 (z) to ∼ 0.11. However, the simulation of a system of this size is very computationally expensive. As in our current work the focus is the bulk phase behaviour not the interfacial region, we divide the system into 3 large bins: 2 surface regions adjacent to the hard walls and the bulk region (cf. Figure 2). With n bin = 3 system-size error in S 2 (z) decreases to ∼ 0.04 where there are now an average of ∼ 500 rods in each bin. The value of S 2 (z) corresponding to the central region is then taken to represent the nematic order parameter of the bulk phase S 2,b . The reduced normal pressure is defined as P * n = P n D/k B T .

A. Pure hard spherocylinders
Prior to demonstrating our results for mixtures of HS and HSC particles, it is instructive to begin by studying a system of pure HSC particles with aspect ratio of L/D = 5. As is well known, the simple HSC model of a mesogen exhibits isotropic, nematic, smectic-A, and solid phases as the density of the system is increased 39,40,42,43,73,[106][107][108] . As an preliminary assessment we demonstrate that the bulk isotropic-nematic transition for the L/D = 5 HSC system contained between well separated parallel hard surfaces simulated using our NP n AT -MC methodology is essentially unaffected by the external field. The bulk phase behaviour for the homogeneous system obtained using conventional constant pressure NP T -MC with full three dimensional periodic boundary conditions 42 is compared with the corresponding data obtained using the constant normal-pressure NP n AT -MC methodology for the system between parallel hard surfaces in Figure 3 and corresponding simulation data is reported in Table I.
The predictions of the isotropic and nematic branches with the MFP theory 86 is also shown in Figure 3 which reduces to the commonly employed PL theory 81-83 for one- As in other studies of confined liquid-crystalline systems 109-113 , significant structure is also apparent close to the surfaces; a full analysis of the surface effects such as wetting, de-wetting, surface nematization, and adsorption will be left for future work.
In order to estimate the location of the bulk isotropic-nematic bulk phase transition, we examine the density dependence of the nematic order parameter in Figure 5. The order parameter of a finite system S 2,N HSC (z) converges quickly to the limiting bulk thermodynamic value S 2,N HSC →∞ (z) for states with intermediate to high orientational order (S 2,N HSC 0.5).
In Figure 6 we display the order parameter profiles for states of low to moderate orientational order (0.1 S 2,N HSC 0.5) in the close vicinity of the isotopic-nematic transition.
Snapshots of typical configurations of these two states are shown in Figure 7: the low-density configuration corresponds to an bulk isotropic state, the intermediate-density configuration  The reduced normal pressure P * n is set in the simulation and corresponding bulk values of packing fraction η b , nematic order parameter S 2,b , and box length L z are obtained as configurational averages. The isotropic phase is denoted by Iso, the nematic by Nem, and the pre-transitional states by Pre.  to a pre-transitional state, while the high-density configuration has clearly undergone a transition to a bulk nematic phase. The pre-transition states, assumed here to correspond to nematic order parameters in the range 0.3 S 2,N HSC 0.5, are due to system size effects and can also be exacerbated by the confinement and potentially very slow relaxation order processes in the form of slow nucleation kinetics.
The discrepancy between the values of the nematic order parameter obtained for profiles with n bin = 3 and n bin = 20 histogram bins in the z direction is very small so that the size effects are expected to be small in this case. An isotropic bulk phase region is clearly found in the case of the state with a pressure of P * n = 1.00; in this case the orientational order in the bulk region is found to be low corresponding to an nematic order parameter of S 2,b = 0.116. The order parameter profile for the pre-transitional state with a normal pressure of P * n = 1.12 exhibits a curve with no uniform bulk region; the average of the nematic order parameter of S 2,b = 0.456 obtained in the central part of the cell does not therefore represent that of a true bulk phase. In the case of the denser state corresponding to P * n = 1.14, the value of the bulk nematic order parameter S 2,b = 0.553 obtained as an average over n bin = 3 bins is essentially equivalent to the value of S 2,b = 0.556 obtained with n bin = 20 bins in the homogeneous central region of the simulation cell. The difference in the orientational ordering for the states corresponding to normal pressures of P * n = 1.12 and P * n = 1.14 is also apparent from Figure 7 where the orientations of HSC rods have been colour coded to aid the visualization. Small clusters of nematic domains are seen for the pre-transitional states (Fig 6), which lead to slightly larger values of the nematic order Clearly, the equilibrium state at the pressure of P * n = 1.14 corresponding to a bulk density of η b = 0.419 and nematic order parameter of S 2,b = 0.553 can be taken to correspond to the lowest-density nematic state, while the state at the slightly lower pressure of P * n = 1.12 is seen to exhibit some small nematic clusters which would correspond to a pre-transitional state; large region with random orientations corresponding to a bulk isotropic liquid can be seen in the case of the system with P * n = 1.00. Our estimate of the isotropic-nematic transition for the L/D = 5 HSC system from an analysis of this data is η b,iso = 0.386, η b,nem = 0.419 for the bulk coexisting phases which are in good agreement with the corresponding results estimated from conventional NP T -MC simulation with full three-dimensional periodic boundary conditions (η iso = 0.407, η nem = 0.415) 42 ; the coexistence pressure is arithmetic average between the normal pressures corresponding to the highest-density isotropic state and the lowest-density nematic state, P * n = 1.07, which is slightly lower than that estimated for the system with full periodicity (P * n = 1.19) 42 . A determination of the free energy and chemical potential of the system would allow one to get a more precise estimate of the position of the phase transition, but this is beyond the scope of the current study. A slight stabilization of the isotropic-nematic transition (corresponding to a lowering the transition packing fraction by 1 to 2%) is therefore found in our systems of HSC rods placed between two parallel hard walls. One would certainly expect surface induced nematization to sta- Isotropic (P * n = 1.00), pre-transitional (P * n = 1.12), and nematic (P * n = 1.14) states in the vicinity of the isotropic-nematic transition are examined; the thermodynamic state is characterized by the value of the normal pressure tensor, P * n = P n D 3 /(k B T ). In the case of S 2 (z) the profiles are constructed from histograms with n bin = 3 (symbols) and n bin = 20 (curves) to assess possible system size effects.
bilize the transition to a bulk nematic phase. These surface effects are however not the focus of the current study and will be discussed in detail in subsequent studies. When the two confining walls are well separated the effect on the bulk isotropic-nematic transition is inappreciable. However, the effects of confinement will become increasingly more important as the surfaces are brought closer together, i.e., at higher pressures and/or for small system sizes. For the systems studied in our current work the two hard walls in our simulation cells are far enough apart (even for the highest density states) so that the effect of the surfaces on the bulk phase transition is small.

B. HSC-HS mixture
We now turn our attention to a binary mixture of HS and HSC particles characterized by overall HS composition x HS ranging from 0.05 to 0.20. As with the pure component system the mixture is contained between parallel hard structureless surfaces and the isotropicnematic phase transition is simulated using the NP n AT -MC method (cf. Section III). As well as being of different density the coexisting isotropic and nematic phases will exhibit a fractionation of the two components, with an accumulation of the the rod-like particles in the orientationally ordered nematic phase. The fluid phase separation in hard-core system of this type is an entropy driven process. This is not to be confused with the depletion driven phase behaviour exhibited by anisometric colloidal particles on addition of polymer where the polymer induces an effective attractive interaction (depletion force) between the colloids, that would otherwise only interact in a purely repulsive fashion, giving rise to a van der Waals like "vapour-liquid" transition (see for example references [ 58,66,68,69,[114][115][116][117]  order parameter in the bulk isotropic at P * n = 1.10 and nematic at P * n = 1.13 states can be seen in Fig. 8(c); the V-shape nematic order profile of the pre-transitional state at P * n = 1.21 is also clearly apparent. The corresponding data for the mixture with 10% HS is given in Table II. The dependence of the bulk packing fraction η b on the applied normal pressure P * n for the HSC-HS system with an overall composition of x HS,tot = 0.1 is illustrated in Figure   9(a). The theoretical description of isotropic and nematic branches of the equation of state obtained with the PL and MFP approach at the same overall composition are also plotted in Figure 9(a) for comparison. There is only a negligible difference between the PL and MFP results in the isotropic state. This is to be expected, since both the density and the ordering in this region is small and as a consequence of the effective hard-sphere treatment of the excluded volume contributions with a one-or two-fluid approximation should be similar for the isotropic state. One should note that while within the theories (both PL and MFP) the density and composition are input variables for a given system and the pressure is an output, while for our NP n AT -MC simulation the equilibrium pressure is specified and the bulk density and composition are obtained as averages from the central region of the cell. An observable difference between the one-and two-fluid theories can be seen in the vicinity of the isotropic-nematic transition where the branches of the equation of state obtained by simulation data experience an abrupt change in slope indicating the transition to the orientationally ordered state. From the results depicted in Figure 9, one can infer that predictions with the two-fluid MFP approach is marginally superior to that with the one-fluid PL at least for the nematic phase. The isotropic-nematic transition point can be identified from the abrupt change in nematic order parameter as is apparent from Figure 9 by Iso, the nematic by Nem, and the pre-transitional states by Pre. The corresponding results for the HSC-HS system with the highest overall HS composition by Iso, the nematic by Nem, and the pre-transitional states by Pre. x HS,tot = 0.20, are shown in Figure 10 and in Table III. The isotropic-nematic transition lies somewhere between highest-density bulk isotropic state at P * n = 1.24 and the lowest-density bulk nematic state at P * n = 1.32 where a clear change in the density and nematic order parameter is exhibited: the values of the packing fraction and nematic order parameter for these two states are η b = 0.404 and S 2,b = 0.260 for the bulk isotropic phase and η b = 0.430 and S 2,b = 0.619 for the bulk nematic phase. For completeness the less extensive data for overall hard-sphere compositions of 5% and 15% are given in Tables IV and V. The isotropic-nematic phase boundaries estimated for the HSC-HS mixtures for bulk phase compositions ranging from the pure HSC system (x HS,b = 0) to x HS,b ∼ 0.3 is shown in Figure 11. Here, we also compare the predictions of PL and MFP theories with our On increasing the overall composition of the HS particles,the density gap between the isotropic and nematic coexisting states is seen to become larger as found with our simulations. By contrast,the density gap between the phase boundaries obtained from conventional NP T -MC simulations of the fully periodic system 60 appears to shrink slightly.
In fully periodic NP T -MC simulations of this type the system remains essentially homogeneous so that the composition of the state remains fixed. As the pressure is increased the system will undergo a transition from an isotropic to nematic phase but the states will be constrained to have the same composition. As a consequence of lack of fractionation of the species between the two phases with the NP T -MC simulations it is difficult to differentiate metastable states within the binodal region from those corresponding to the coexistence boundaries.
The NP n AT -MC simulation data for the HSC-HS mixture are also summarized reported in Table VI where the slight composition asymmetry between coexisting isotropic and nematic phases is clearly apparent. For the HSC-HS system, the addition of spherical particles is found to destabilize. The formation of a bulk nematic state predominantly composed of HSC rods will cause a reduction in the concentration of the HS particles in the same region, and as a consequence the orientational ordered state which is of higher density but lower bulk HS composition coexists with an isotropic state of lower-density and higher bulk HS composition.
Finally, a phase diagram in the pressure-composition (P * n -x HS,b ) plane is shown in Figure  12. A very narrow region of isotropic-nematic coexistence is obtained with our NP n AT -MC simulation approach. The coexistence region is seen to be at lower pressures that that predicted with the PL and MFP theories or obtained by fully periodic NP T -MC simulation 60 . It should be noted that the results obtained at higher compositions with the fully periodic simulations is in good agreement with the nematic branch predicted with MFP theory. The quality of both PL and MFP is affected by a shift in pressure for the pure HSC system. The predictions with the two-fluid MFP theory are seen to be much better than one-fluid PL theory, particularly for systems of higher composition.

V. CONCLUSIONS
In this paper we have studied phase behaviour of mixtures of purely repulsive rod-like (HSC) and spherical (HS) particles. Using the new NP n AT Monte Carlo simulations we have constructed the isotropic-nematic phase diagrams of the HSC-HS mixture in the pressuredensity, pressure-concentration and density-concentration projections. The comparison of our results with previously reported fully periodic NP T -MC data 60 reveals good agreement at low HS concentrations but some discrepancy occurs for higher HS concentrations especially in the nematic phase. In conventional NP T -MC simulations the system is essentially homogeneous such that the overall composition of the system remains fixed. This means that unless one has a very large system the coexisting isotropic and nematic states would have identical compositions 119,120 . By contrast, using our NP n AT -MC algorithm we have observed compositional asymmetry between the coexisting phases, with a slight increase of the HSC concentration in the nematic phase due to packing effects. This conclusion is also supported by the fact that our NP n AT results for pure HSC particles are completely consistent with the available data obtained from NP T -MC simulations using full three-dimensional boundary conditions 42 . The advantage of our method is that the fluid phase separation can be studied within a single simulation cell which circumvents the problem of inserting anisotropic particles inherent in other techniques such as Gibbs ensemble MC (GEMC). Additionally, particle exchanges are allowed between the surface and bulk regions so that one is able to treat the two bulk states coexisting at different bulk densities (packing fractions) as well as different bulk compositions. The effect of the pair of the auxiliary hard walls put at the end of the box is not appreciable in the bulk region for sufficiently prolonged systems in the direction perpendicular to the walls and considerably helps to stabilise the phase separated states with a low interfacial tension which as typically exhibited by hard-body systems.
We further compared the new simulation data with two theoretical predictions that go beyond the Onsager second-virial theory and extend the well known PL theory to mixtures in two different ways: using the one-fluid approximation whereby one maps the mixture onto an effective HS mixture, and by treating both species separately as an effective HS mixture which we refer to as many-fluid approximation. The comparison reveals that both the one-and many-fluid approaches provide a reasonably accurate quantitative description of the mixture including the isotropic-nematic phase boundary and degree of orientational order of the HSC-HS mixtures. The many-fluid prediction of the coexisting pressure is arguably found to be slightly better to that obtained with the one-fluid method. However, systems with larger aspect ratios should be considered to make a better assessment of the performance of the theories.
Our work can be directly extended in a number of ways. For instance, we have restricted our attention to systems of HSC rods with an aspect ratio of L/D = 5 and HS particles of diameter σ = D. Considering larger aspect ratios would be a more stringent test of the theoretical predictions, as one would expect a more significant improvement of the manyfluid treatment compared to the one-fluid approach. The compositional asymmetry in the coexisting isotropic and nematic phases is expected to be more appreciable. The method for locating the phase boundaries can be improved using thermodynamic integration. Furthermore, our NP n AT -MC approach can be directly applied to describe the surface phenomena due to both fluid-fluid and fluid-wall interfaces. These interfaces plays an important role in liquid crystalline systems 109,121 due of rich surface-induced effects (e.g., nematization and smectization), characteristic in systems comprising anisotropic particles.

LW thanks Department for Business Innovation and Skills, UK and China Scholarship
Council for funding a PhD studentship. MC simulation in this work is performed using  pressure P * n = P n D 3 /(k B T ) is set during the simulation, bulk packing fractions η b and bulk compositions x HS,b of coexisting isotropic (iso) and nematic (nem) states are obtained as averages from the central region of the cell. The bulk nematic order parameter S 2,b of the lowest-density nematic bulk phase is also shown.
x HS,tot P * n,iso