Structure and interfacial tension of a hard-rod fluid in planar confinement

The structural properties and interfacial tension of a fluid of rodlike hard-spherocylinder particles in contact with hard structureless flat walls are studied by means of Monte Carlo simulation. The calculated surface tension between the rod fluid and the substrate is characterized by a nonmonotonic trend as a function of the bulk concentration (density) over the range of isotropic bulk concentrations. As suggested by earlier theoretical studies, a surface-ordering scenario is confirmed by our simulations: the local orientational order close to the wall changes from uniaxial to biaxial nematic when the bulk concentration reaches about 85% of the value at the onset of the isotropic-nematic phase transition. The surface ordering coincides with a wetting transition whereby the hard wall is wetted by a nematic film. Accurate values of the fluid-solid surface tension, the adsorption, and the average particle-wall contact distance are reported (over a broad range of densities into the dense nematic region for the first time), which can serve as a useful benchmark for future theoretical and experimental studies on confined rod fluids. The simulation data are supplemented with predictions from second-virial density functional theory, which are in good qualitative agreement with the simulation results.

The structural properties and interfacial tension of a fluid of hard-spherocylinder rod-like particles in contact with hard structureless flat walls are studied by means of Monte Carlo simulation. The calculated surface tension between the rod fluid and the substrate is characterized by a non-monotonic trend as a function of bulk concentration (density) over the range of isotropic bulk concentrations. As suggested by earlier theoretical studies, a surface-ordering scenario can be confirmed from our simulations: the local orientational order close to the wall changes from uniaxial to biaxial nematic when the bulk concentration reaches about 85% of the value at the onset of the isotropic-nematic phase transition. The surface ordering coincides with a wetting transition whereby the hard wall is wetted by a nematic film. Accurate values of the fluid-solid surface tension, the adsorption, and the average particle-wall contact distance are reported (over a broad range of densities into the dense nematic region for the first time), which may serve as a useful benchmark for future theoretical and experimental studies on confined rod fluids. The simulation data are supplemented with predictions from a second-virial density functional theory, which are in good qualitative agreement with the simulation results. The interaction between a solid surface and a dense nematic liquid-crystalline fluid can be characterized by the interfacial tension, a challenging property to quantify for fluids of anisotropic molecules. A detailed understanding of the interfacial properties would provide substantial insight into how these complex fluids are influenced by confinement 1 -insight that would be invaluable, for example, in the context of liquid-crystalline displays, a technology that has facilitated a revolution in portable display devices. As a result the surface thermodynamics of anisotropic molecules remains an important area of ongoing research.
Confinement -or the presence of coexisting phasesgives rise to inhomogeneities, increasing the richness of the structure and phase behaviour with an additional surface tension contribution to the thermodynamics of the system associated with the corresponding interfaces. This interfacial tension, or surface free energy, is responsible for many important and fascinating processes. However, considering the challenges involved in its quantification, it is perhaps not surprising that there is still a lack of understanding of the interfacial tension of liquidcrystalline systems in contact with solid surfaces, a deficiency that is particularly striking for higher-density states.
The purpose of our current work is to improve the overall thermodynamic understanding of such systems a) Electronic mail: p.brumby@keio.jp b) Electronic mail: rik.wensink@u-psud.fr c) Electronic mail: a.haslam@imperial.ac.uk d) Electronic mail: g.jackson@imperial.ac.uk and quantify the interfacial tension and surface structural properties of a simple model of a confined liquidcrystalline system, across a broad range of concentrations. Hard spherocylinders represent a suitable simple model of purely repulsive rod-like particles to study the rich morphology of liquid-crystalline matter. The hard-spherocylinder system is known to exhibit the prototypical nematic (orientationally ordered) and smectic (orientationally ordered with one degree of positional order in layers) liquid-crystalline states, as well as the more usual isotropic (orientationally disordered) fluid and solid (orientationally and positionally fully ordered) phases [2][3][4][5][6] . Despite the apparent simplicity of the model, our work will nevertheless provide a benchmark for more complicated effective molecular topologies 7 , interaction potentials 8 , and/or confined structures 9 .
Prior simulation studies on confined liquid-crystalline systems have provided a detailed appreciation of the influence of surface-liquid interactions on the behaviour of these complex fluids. In addition to extensive work on colloidal systems, Dijkstra and co-workers [10][11][12] have simulated hard-spherocylinder fluids confined between two hard planar walls. In agreement with subsequent experimental observations for rod-like particles at interfaces 13 , Dijkstra and co-workers found that the bulk isotropic phase will partially wet the hard planar wall with increased nematic ordering, biaxiality, and planar anchoring. The density profiles of these lower-density states reveals a distinctive peak at a distance of half the length of the particle from the wall surface. In addition it is found that the separation distance between the confining walls must exceed two particle lengths before a first-order isotropic-nematic phase transition can be observed in the 2 "bulk" region of the system 14 . In close proximity to their surfaces, hard confining walls may induce nematic ordering, even if the molecules are not sufficiently elongated to naturally form nematic phases in the completely homogenous system. 15 As well as simulations of fluids of hard spherocylinders confined between hard walls, hard-ellipsoid 16,17 and hard Gaussian overlap [18][19][20] particles have been studied in contact with soft walls. Depending on the nature and degree of the wall-particle interaction, one finds that either planar or homeotropic anchoring of molecules relative to the wall surfaces can be observed. Further modification of the wall-particle interaction can even encourage surface anchoring with tilted orientations 21 . Recent experimental work confirms that similar variations in particleinterface contact angles are possible in real systems 22 .
Introducing surface roughness, as in the simulation studies of Cheung and Schmid 23 with soft ellipsoids, has interesting consequences. The degree of roughness acts to both inhibit the formation of nematic phases at the wall and to shift the bulk isotropic-nematic phase transition towards higher pressures (densities). The presence of particle-particle and/or particle-surface attractive interactions, can also have a substantial effect on the phase transitions and wetting behaviour of the confined fluid [24][25][26] . Other comprehensive studies of surface adsorption and density profiles of confined fluid from the wall surface have involved polymers 27,28 , flexible rods 29 , Gay-Berne attractive particles 8,30-32 , platelets 33,34 , hard rectangles 35 , and Lennard-Jones molecules in contact with grooved surfaces 36,37 to name but a few representative examples.
Our detailed understanding of the effect of confinement on the structure, phase behaviour, and wetting transitions of fluids is, however, in stark contrast to the current knowledge of the solid-fluid interfacial tension of anisotropic particles. In addition to some unresolved fundamental issues and challenges in the analysis of the surface thermodynamics of non-spherical particles, the fact that the interfacial tension is a thermodynamic derivative property leads to practical difficulties; obtaining the solid-fluid tension to an acceptable level of accuracy requires sampling over a very large number of configurations. Unless efficient techniques are employed, the computational expense can, depending on the type and size of system, be insurmountable. Notwithstanding, a variety of different methods are now available for the computation of the interfacial tension in molecular simulations.
The most common methodology for the computation of the interfacial tension is based on the mechanical definition of the restoring force due to a deformation in the area of the system 38,39 : the interfacial tension is calculated from the direct evaluation of the components of the pressure tensor (the negative of the stress tensor) in terms of the virial -specifically, the forces acting in directions normal and tangential to the interface. This is the method of choice in molecular-dynamics simulations, and there are many such examples to be found in the litera-ture [40][41][42][43][44] . Although novel approaches based on the mechanical route are still being developed 45 , the technique is unsuitable when the particle-wall and/or particle-particle interaction potentials are discontinuities or when the interface is not planar. As a consequence of these shortcomings a number of alternatives to the mechanical/virial route have been developed.
The thermodynamic approach of Bennett 46 is amongst the earliest alternative methodologies for the determination of the interfacial tension (surface free energy) where a multi-stage method is employed to sample the freeenergy difference of selected systems. While the technique is well suited for high-density states, the implementation is not straightforward, requiring multiple simulations to obtain the value of the interfacial tension at each state. Salomons and Mareschal 47 have developed the method of Bennett to arrive at an expression that can be readily employed in molecular-dynamics simulation. As with the test-area approach 39 , which will be discussed in more detail later in this section, the implementation of the Bennett method is typically based on the assumption that a decrease in the interfacial area leads an equivalent free-energy change to that brought about by the equivalent increase in surface area. While this holds true for molecules interacting through continuous spherically symmetrical potentials, it is, unfortunately, generally invalid for systems comprising non-spherical molecules.
The finite-size scaling method of Binder 48 and its variants offer an attractive route to the determination of the interfacial tension. Examples of its application include fluid interfaces comprising Lennard-Jones 49 and squarewell 50 molecules. The method relies on the computation of the probability of states in terms of the number of molecules in the system from simulations performed in the grand-canonical ensemble. Accurate sampling of high-density states is therefore difficult to achieve, particularly in the case of non-spherical hard-body systems for which the probabilities of successful particle insertions are acutely small 51,52 . The use of the Binder finite-size scaling method with grand-canonical transition-matrix Monte Carlo and histogram-reweighting can alleviate the problem to a degree 53 . The treatment of dense fluid systems characterized by hard-body interactions nevertheless remains challenging with approaches involving particle insertions. Recent work using finite-size scaling, as part of the ensemble switching method, have further extended the applicability of the Binder approach 54-56 . In doing so the interfacial tension can be computed reliably for vapour-liquid systems and, with a certain degree of success, for solid-liquid interfaces.
The interfacial tension of a fluid phase of anisotropic particles in contact with a surface can be obtained by integrating the Gibbs adsorption equation as described by Mao et al. 57 . This type of thermodynamic integration approach is, however, again implemented within the grand-canonical ensemble and thus suffers from the same drawbacks inherent with the method of Binder of poor insertion statistics when applied to higher-density 3 states. The related Gibbs-Cahn thermodynamic integration technique 58,59 has also been used to determine the wall-fluid interfacial tension of confined systems comprising hard-sphere [59][60][61] and Lennard-Jones 62 particles, although the method requires careful parameterization.
Studies with other approaches which are specific to confined solid-fluid systems have been reported including: the phantom-wall method of Leroy and Müller-Plathe 63 ; the thermodynamic-integration approaches of Hamada et al. 64 and of Das and Binder 65 ; the interface potential analysis method developed by Errington and co-workers [66][67][68] ; and the ensemble mixing method of Deb et al. 69 , wherein one gradually inserts a wall potential in a system without walls. All of these methods provide a route to the interfacial tension -and in some cases the contact angle -of confined fluids.
In a similar manner to the multi-stage simulation scheme of Bennett, obtaining the interfacial tension via the expanded ensemble method 70-73 requires the simulation of multiple connected sub-ensembles. Ideally, these connected systems should have identical properties (i.e., number of particles N , volume V , and temperature T ) excepting an incremental difference in interfacial area between one system and the next. The expanded ensemble thermodynamic method proceeds according to the standard canonical (N V T ) ensemble Monte Carlo algorithm 74 , but with extra trial moves where a jump from one system to another is attempted. The difference in free energy between the systems with different interfacial area is then used to compute the interfacial tension. The related wandering-interface method of MacDowell and Bryk 75 involves a slightly different procedure. In this case, the domain lengths of the system are allowed to fluctuate freely, while maintaining a constant system volume. An analysis of histograms of the probability distribution of the domain lengths allows one to extract the surface tension as the logarithm of the distribution, as it tends towards zero. This thermodynamic approach has the important advantage of not being constrained to the grandcanonical ensemble. In addition, it is equally applicable to both vapour-liquid and solid-liquid systems comprising non-spherical particles. The perturbative wanderinginterface method has been used by Blas et al. 76 to determine the interfacial tension for systems of freely rotating chain molecules, including subsequent work on the effect of long-range interactions on the interfacial tension [77][78][79] .
In our current work we employ a free-energy perturbation methodology to determine the solid-fluid interfacial tension of hard spherocylinders in contact with a hard wall. The approach is related to the aforementioned perturbative test-area method 39 which is not restricted to the grand-canonical ensemble precluding problems associated with high-density states. As its name implies the test-area method is a thermodynamic approach involving the computation of the change in free energy accompanying a vanishingly small perturbation in the surface area of the system. The versatility of the test-area method for the determination of the interfacial tension is apparent from its broad implementation to systems involving diatomic molecules 80 , chain molecules 76,81 , liquid drops 82,83 , mixtures 84,85 , and fluids confined within slitpore 86,87 and cylindrical 88 geometries. It is now a popular alternative for the simulation of the vapour-liquid interfacial tension of fluid systems comprising molecules characterized by continuous interactions. There are a number of inadequacies associated with the direct implementation of the test-area method for systems involving discontinuous potentials or hard non-spherical particles; in this case it is more appropriate to determine the interfacial tension from the components of the pressure tensor obtained from a thermodynamic route involving test-volume perturbations.
The test-volume perturbation technique was first adopted by Eppenga and Frenkel 89 for the determination of the bulk (macroscopic) pressure of the isotropic and nematic phases of systems of purely repulsive hard discs; owing to the discontinuous nature of the potential the pressure can be obtained directly by examining the probability of configurations with overlapping particles resulting from vanishingly small (isotropic) volume perturbations (compressions). Extensions of the methodology to systems with attractive interactions such as Lennard-Jones 90 , square-well 91 , and anisotropic Gay-Berne 92 fluids have been reported. The test-volume approach has also been used to determine the bulk pressure of other systems comprising anisotropic particles such as hard-Gaussian overlap molecules 93 and hard spherocylinders 94 , and has become the method of choice for purely repulsive particles of arbitrary shape 95,96 .
As mentioned earlier, the interfacial tension can be determined from knowledge of the normal and tangential components of the pressure tensor. Test-volume perturbation methods for the calculation of the macroscopic components of the pressure tensor are well suited for the simulation of the interfacial tension of systems with planar interfaces. This type of perturbative approach is non-invasive and is applicable over a wide range of densities, as it is not anchored to a specific simulation ensemble. Test-volume perturbations have been applied for the determination of the pressure tensor of hard spherical 93 and non-spherical 94 particles; it is important to emphasize that there are important subtleties to consider in the implementation of the approach to nonspherical particles owing to the lack of equivalence between the compressive and expansive contributions to the pressure tensor associated with anisotropic volume changes 94 . Jiménez-Serratos et al. 97 further extended its utility and successfully applied the method to general non-spherical systems with discontinuous potentials, such as chain molecules and spherocylinders with squarewell and square-shoulder interactions. For an excellent discussion and comparison of the various techniques for the simulation of the pressure tensor of confined systems characterized by discontinuous interactions the reader is directed to the review by Deb et al. 69 .
It may now have become apparent that volume-4 perturbation methods offer great promise for the simulation of the solid-fluid interfacial tension, especially in the case of high-density states of molecules interacting through discontinuous potentials -states which are inaccessible with the majority of other approaches. The key goal of our paper is to determine, for the first time, the interfacial tension of confined systems of hard spherocylinders confined between structureless hard walls, for states ranging from low-density isotropic phases, through the (bulk) isotropic-nematic phase transition, and deep into the high-density nematic region of the phase diagram. In addition, a detailed analysis of the effect of confinement on the thermodynamic, structural and orientational properties of hard-spherocylinder fluid is made. The remainder of this paper is organized as follows: the simulation methodology is detailed in Sec. I; we examine the phase diagram of the unconfined bulk system in Sec. II; profiles of the density and nematic-order parameter for the confined systems are presented in Sec. III; the test-volume perturbation method is employed to determine normal and tangential components of the pressure tensor and thus the solid-fluid interfacial tension in Sec. IV; an Onsager density functional theory for hard spherocylinders confined in slit-pore geometry and its predictions are compared with the simulation results in Sec. V; finally, some concluding remarks are made in Sec. VI. Details of the perturbation method used in our current work are assigned to the appendices, wherein our simulated interfacial-tension data are also tabulated.

I. SIMULATION METHODOLOGY
In all of the simulations described in our work, a collection of N hard spherocylinders (each formed from a hard cylinder of length L capped by two hard hemispheres of diameter D where the aspect ratio is defined as L/D) are placed in a rectangular box of volume V and dimensions x , y and z . Standard periodic boundary conditions are used in all three directions in the case of the bulk systems. Planar confinement is considered by placing structureless hard walls at positions z = 0 and z = z of the simulation box while maintaining the periodicity of the other system boundaries; the closest distance of approach of a particle from the hard wall is therefore D/2. We explore the phase behaviour of these systems at various concentrations (densities) of the rod particles. Sampling is performed using Wood's adaption 98,99 of the Metropolis Monte Carlo method 74 for the isobaric-isothermal (constant N pT ) ensemble in the case of bulk (unconfined) systems, and using the standard Metropolis method 74 in the canonical (constant N V T ) ensemble for the confined systems (see Refs. 5 and 94 for further details). Trial states are created by performing a series of cycles: N attempts to randomly translate or rotate one of the hard spherocylinders (which are also selected at random for each attempt) are made within each cycle, maintaining the condition of detailed balance; where appropriate a single system-volume expansion or compression is attempted at the end of each cycle in the isobaric-isothermal ensemble. This process is carried out until equilibrium states are attained, and the density and nematic-order parameter are then calculated as configurational averages. The amplitudes of the trial moves are chosen on an empirical basis so as to achieve an overall acceptance ratio of ∼ 30% for the translational and rotational displacements, and system volume changes.
The system density ρ = N/V can be conveniently quantified in terms of the volume (packing) fraction as where for our model v m = πD 3 /6 + πD 2 L/4 is the volume of the hard-spherocylinder particle. A dimensionless concentration is also often employed in theoretical studies to characterize the density of the system, defined in terms of the so-called Onsager limit for the second-virial coefficient of infinitely long rod-like particles: Both measures of density are employed interchangeably in the ensuing discussion, the dimensionless concentration being the more appropriate to facilitate comparison with previous theoretical studies. With our choice of Gibbs dividing surface at z = 0 and z = z , two inaccessible layers of thickness D/2 in the proximity of the walls has to be accounted for; this has implications in the determination of the overall system volume, density, wall-fluid surface tension and surface adsorption (cf. Section IV). As a consequence of the external potential due to the presence of the confining walls, the density is not homogenous throughout the system but will depend on the distance z from the wall. Local density profiles can be obtained by dividing the space between the walls into a number of equal-sized bins. The concentration c(z) of hard spherocylinders within each bin, at a distance z, is determined from Eq. 2, where V and N now represent the volume of the bin and the number of particles with their centre of mass within the bin, respectively. In the case of a sufficiently large separation distance between the two walls, the density in the central (bulk) part of the system ceases to be a function of z; the corresponding values of the packing fraction and concentration of the bulk phase are obtained as averages of the density profiles in the central region of the simulation box and are denoted by η b and c b , respectively. The pressure p of the system is expressed in dimensionless form throughout our work. A convenient measure for the athermal hard-spherocylinder system is where k B is the Boltzmann constant. Equivalent dimensionless expressions for averages of the normal p N and 5 tangential p T components of the pressure tensor p αα are employed in our subsequent analysis of the interfacial tension. In the case of the inhomogenous system confined between two parallel planar walls (in the xy plane) at mechanical equilibrium, the normal component p N = p zz of the pressure tensor is constant and corresponds to the bulk scalar pressure p; for this planar interfacial geometry the tangential component p T (z) = (p xx (z)+p yy (z))/2 of the pressure tensor is cylindrically symmetrical and is a function of the position z from the surface of the wall (only the macroscopic average of p T over the box dimension z needs to be evaluated in order to determine the interfacial tension of the system). 38 The average local orientational order of the spherocylinders can be quantified from the second-rank traceless, symmetric tensor 5,89 : whereû represents the orientational unit vector along the main axis of the hard spherocylinder, ⊗ denotes the dyadic product, I the unit second-rank tensor, and · the canonical ensemble average. The nematic-order parameter S, defined as the largest eigenvalue S = λ + of the tensor Q, is a key parameter to distinguish between nematic (S > 0) and isotropic order (S = 0). The corresponding eigenvectorn corresponds to the nematic director, which for the confined systems indicates the principal direction of alignment with respect to the wall normal. In the case of uniaxial order, the two remaining eigenvalues (λ 0 and λ − ) are equal and there is no preferential orientational order in the plane perpendicular to the nematic director. The presence of a hard wall, however, will break the uniaxial nematic symmetry of the fluid and induce local biaxial nematic order: a non-zero difference ∆ = λ 0 − λ − between the two smallest eigenvalues of Q can be used as a measure to quantify the degree of biaxial order across the range of distances z. As for the concentration profile c(z), the spatial resolution of the orientational-order parameters between the walls can be obtained in a straightforward manner by discretizing z into a finite number of slabs of equal volume (see Ref. 100).

II. BULK PROPERTIES OF THE HARD-SPHEROCYLINDER FLUID
Prior to studying the system under confinement, we first characterize the unconfined fluid of hard spherocylinders, paying particular attention to the degree of orientational order over a range of concentrations. This allows for a better appreciation of the (sometimes subtle) effects of confinement on the liquid-crystalline stats formed by the system, as described later in our paper. For the bulk systems, we study N = 3080 hard spherocylinders, each with an aspect ratio of L/D = 10. Simulations are performed in the isobaric-isothermal (N pT ) ensemble, starting from an initially perfect face-centredcubic crystal-lattice arrangement. The pressure of the system is reduced until an equilibrium low-density disordered isotropic phase is achieved. The system is then compressed from this state by increasing the pressure in incremental steps. A minimum of 2 × 10 6 Monte Carlo cycles are required to equilibrate these systems at each pressure. The average values of the equilibrium density and order parameters are then determined over a further 2 × 10 6 cycles. The values obtained for the compression runs are shown (as the red circles) in Fig. 1. In a similar manner, expansion runs (involving the incremental reduction in the pressure) are performed, starting with a stable nematic state at a density well above the isotropic-nematic phase transition pressure. From the equation of state depicted in Fig. 1 one can estimate the isotropic-nematic coexistence pressure as p * ∼ 1.88 corresponding to the sharp discontinuity in the nematicorder parameter. For the isotropic branch this pressure corresponds to a concentration of c I ∼ 2.29 (η I ∼ 0.24) and nematic-order parameter of S I ∼ 0.04, while for the nematic branch the corresponding values are c N ∼ 2.50 (η N ∼ 0.27) and S N ∼ 0.71.

III. STRUCTURE AND ORIENTATION OF THE HARD-SPHEROCYLINDER FLUID CONFINED BY HARD WALLS
Having simulated systems of hard spherocylinders with periodic boundaries along all three dimensions, we now focus on the systems confined between two hard walls, highlighting the influence of confinement on the isotropicnematic phase transitions (relative to that of the bulk unconfined system). We examine the local properties (density and orientational order) in the regions where the solid walls interact with the hard spherocylinders. As mentioned earlier, the confinement implemented in our current work involve a pair of flat, featureless, and impenetrable walls positioned parallel to each other a distance z apart at the boundaries of the z Cartesian axis; this arrangement is often referred to as a slit-pore geometry.
In our case, we require z to be large enough for a clearly distinct bulk phase to form in the centre of the system. In order to maintain a fixed inter-wall separation distance, these systems are simulated in the canonical (N V T ) ensemble. The initial configuration is created from a high-density state of N = 3000 non-overlapping hard spherocylinders arranged in a perfect face-centredcubic crystal lattice. The specific box length dimensions For the calculation of density profiles in terms of the distance from the surface of one of the walls, the system volume is divided into n bin = 200 bins, positioned adjacent to one another along the z axis. A minimum of 1 × 10 6 Monte Carlo cycles are performed for each state, and averages of the density for each bin are determined. Due to the sensitivity of the calculation of the local nematic-order parameter to the number of particles sampled in each bin, it is necessary to increase the system size. In order to simulate a larger number of particles, the systems are replicated in the x and y axes, such that there are at least N = 25, 000 hard spherocylinders in the simulation cell. The extra periodicity created by this replication process is removed by performing an additional series of equilibration cycles until new equilibrium states are generated. The considerably larger systems are then simulated for a further 1 × 10 5 cycles to determine equilibrium profiles of the nematic-order and biaxiality parameters. To further increase the average number of particles per bin, the overall system volume is divided into a smaller number n bin = 20 of bins along the z axis (reducing uncertainties relating to finite-size effects of the sample, at the cost of a lower resolution).

Isotropic
Nematic (a) A detailed analysis is first made of the local density c(z), uniaxial nematic order S(z), and biaxial order ∆(z) for equilibrium states with stable bulk isotropic phases. The profiles of these properties as functions of the distance from one of the walls are presented in Fig. 2 (left), and representative snapshots of configurations of these systems are displayed in Fig. 3 (a) and (b). At low bulk concentration, a particle dewetting of the wall surface is observed as the rods are depleted from the impenetra-  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 7 ble wall. The nematic orientational order near the wall is anti-nematic lim z→0 S(z) = S w ∼ −1/2 and uniaxial lim z→0 ∆(z) = ∆ w ∼ 0 indicating planar anchoring with rods oriented perpendicular to the wall normal. Finitesize effects are clearly observable for the biaxiality at the wall, as a result of the relatively low average number of hard spherocylinders per bin; this is particularly noticeable for the systems with lower concentrations.
As the bulk density is increased there is a marked change in the surface structure. A sharp surface phase transition is observed at c b ∼ 2 where the orientational order at the surface changes abruptly from uniaxial to biaxial. This is more-clearly reflected in the behaviour of the profiles of the nematic order S(z) in close proximity to the wall surface of Fig. 4 (left). The steep increase in the biaxiality ∆ w of the particles in contact with the wall at concentrations of c b ∼ 2 coincides with a marked increase in the uniaxial nematic-order parameter S w at the wall. This is a direct consequence of a change of nematic director which defines the reference frame for the orientational-order parameters.
At low density, the nematic director is oriented randomly with respect to the wall normal while, at the surface phase transition, the director rotates perpendicular to the wall normal. It is evident that this sudden reorientation signals the formation of a nematic wetting layer close to the wall. The extent of this layer can be gleaned from the peaks in the density profiles c(z) in Fig. 2. The transition from uniaxial to biaxial nematic surface order for a rod fluid in contact with a hard wall has been predicted from Onsager theory. A careful stability analysis reveals 101 that one would expect a wetting transition at c b ∼ 2.79, well before the bulk isotropic-nematic transition which occurs at c b ∼ 3.29. From this one can infer that a surface ordering transition will occur universally when the bulk concentration reaches about 85-87 % of the coexistence value at the isotropic-nematic transition. The precise value will depend only weakly on the aspect ratio provided the rods are sufficiently anisometric, i.e., L/D 1. From purely thermodynamic arguments one would expect that hard-core mesogens confined between parallel hard walls in slit-pore geometry would exhibit a transition to a bulk nematic phase at densities below the isotropic-nematic transition of the unconfined system. Indeed the occurrence of a first-order transition from an isotropic phase (with a biaxial nematic film at each wall) to a condensed nematic phase that fills the slit pore has already been reported for confined hard rods 10,11,14 . The slit-pore geometry essentially stabilizes the nematic phase relative to the isotropic shifting the isotropicnematic transition to lower densities (and higher temperatures in the case of mesogens with attractive interactions such as the Gay-Berne fluid 24 ). This stabilization of the nematic phase due to confinement is now commonly referred to as capillary nematization 10,11,14 . We also observe a capillary nematization transition in our systems of confined hard spherocylinders at a density of c b ∼ 2 which is below the onset of the isotropic-nematic transition of the unconfined bulk system, c b ∼ 2.29.  When the bulk phase of the confined systems is nematic, the trends in the local density and particle ordering differ from those for systems with bulk isotropic phases, but smoothly follow the characteristics of the nematic wetting layer formed at concentrations below the bulk isotropic-nematic phase transition. The concentration c(z), nematic-order S(z), and biaxial-order ∆(z) profiles for these systems are shown in Fig. 2 (right); representative configurations of two systems with different bulk densities corresponding to the nematic state are shown in Fig. 3 ((c) and (d). As the density is increased, the biaxial nature of the nematic phase gradually decreases since the rods progressively align perpendicular to the wall normal, thereby reducing the orientational symmetry breaking imposed by the hard wall; the biaxiality of the high-density states is dominated by that of the bulk nematic region. The overall scenario for the surface ordering of rod-like fluids borne out by our simulations is in line with the findings of earlier simulations and Zwanzig-type density functional theory for hard rods of finite anisotropy 10,14,102,103 as well as the predictions from Onsager density functional theory for infinitely slender rods 101,104 .

IV. INTERFACIAL TENSION, SURFACE ADSORPTION, AND AVERAGE ROD-WALL CONTACT DISTANCE
In this section, we provide a detailed analysis of the interfacial thermodynamic properties pertaining to the evolution of the structure of the rod fluid at the wall surface upon increasing the bulk concentration. Of particular interest is the effect of the onset of local biaxial nematic order and the formation of a nematic wetting layer on the fluid-wall surface tension, a quantity of great practical interest in understanding the behaviour of liquid crystals at surfaces. A test-volume free-energy perturbation method 94 (outlined in Appendix A) is applied to determine the fluid-wall surface tension of our system of hard spherocylinders in planar slit-pore confinement from canonical Monte Carlo (MC-N V T ) simulation. This approach requires the evaluation of averages of the normal p N = p zz and tangential p T = (p xx + p yy ) /2 components of the pressure tensor of the system using the test-volume method with appropriate anisotropic volume perturbations. The fluid-wall surface tension can be expressed in dimensionless form as where the average components of the pressure tensor are expressed in dimensionless form, p * αα = p αα v m /(k B T ), and the factor of a half accounts for the presence of two surfaces in our system. We should note that the bulk packing fraction is included in the denominator in this definition of the surface tension to aid direct comparison with exiting theoretical estimates. The normal and tangential components of the pressure tensor can be expressed as a sum of the ideal contribution (which in our case simply corresponds to the overall packing fraction η of the system) and the corresponding excess contributions determined from expansive p * + αα and compressive p * + αα perturbations: For the slit-pore geometry considered here, p * + N = p * + zz , p * − N = p * − zz , p * + T = p * + xx + p * + yy /2 and p * − T = p * − xx + p * − yy /2. The excess contributions are determined by performing expansive and compressive test-volume perturbations using the relations and where ∆V α i→j+ and ∆V α i→j− represent changes in system volume brought about by the test (non-permanent) anisotropic affine deformations, achieved by increasing and reducing the length of the α axis, respectively. The methodology essentially requires the calculation of the probabilities P + nov and P − nov that the aforementioned perturbations produce configurations without overlaps between particles and/or between particles and either of the two walls.
Following the methodology described in previous studies 93,94 , 20 evenly spaced values of the expansive ∆V α i→j+ and compressive ∆V α i→j− perturbations are considered for each state, selected so that the resulting values of P + nov and P − nov are evenly spread between zero and unity. The values of P + nov and P − nov in the limit of infinitesimal volume changes are obtained using a linear least-squares fit. The anisotropic test-volume perturbations are undertaken every 20 Monte Carlo cycles in order to achieve efficient convergence (see Ref. 105 ). The averages are evaluated over 4.4 × 10 7 cycles for the low-density systems (0 < c b ≤ 2.1046), while the averages for the higherdensity (c b > 2.1046) systems -which require additional sampling due to slow re-orientation of the nematic director -are taken over 3.6 × 10 8 cycles.
The wall-fluid interfacial tension of the hard spherocylinders confined it slit-pore geometry calculated with the test-volume approach across a range of bulk concentrations, from the low-density isotropic to the highdensity nematic state, is displayed in Fig. 5 (a). The variation of the surface tension with concentration in the isotropic bulk phase is non-monotonic. The maximum in the isotropic fluid-wall tension is reached at a bulk concentration c b ∼ 2 which corresponds to the onset of the wetting transition and surface biaxial order. The interfacial tension can then be seen to decrease as the density is further increased into the nematic region and appears to level off to a plateau value at high concentrations close to the bulk nematic-smectic transition, which is estimated to occur at c b ∼ 5 for the aspect ratio of L/D = 10 considered here 6 . The concentration dependence of a dimensionless interfacial tension defined as γ = γD 2 /(k B T ) = (1/2) z D 2 /v m (p * N − p * T ) in the dilute isotropic phase phase is highlighted in the inset of Fig. 5 (a) to enable direct comparison with the corresponding values determined by Mao et al. 57 by integrating the Gibbs adsorption equation using grand canonical simulation data. Good agreement between the testvolume and Gibbs adsorption wall-fluid interfacial tension is found for the isotropic state; it becomes increasingly difficult to employ grand canonical simulations as the density is increased limiting the applicability of the latter approach.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 10 In order to broaden our insight of the interaction of the rod fluid with the hard walls two additional properties are assessed. The surface adsorption Γ which can be determined (in dimensionless units) from 106 where the number of particles N b associated with the bulk phase is computed from the bulk packing fraction as Fig. 5 (b)); it is important to note that in the case of our confined system V denotes the overall volume of the system including the two regions inaccessible to the centre-of-mass of the hard spherocylinders in the layers of thickness D/2 close to the two walls. The average of the shortest distance Z between the wall and the centres-of-mass of the hard spherocylinders in contact with the wall is also determined (see Fig. 5 (c)). This quantity can be established by recording the subset of overlapping hard spherocylinders {1 . . . N ov } generated by each trial displacement (be it a change of rod position, orientation, or system volume), and performing a canonical average · of the centre-of-mass distance with respect to the wall. The contact distance Z can be expressed (in units of L) as withẑ denoting the direction of the wall normal, and u i the orientation of the principal unit vector of the rod i overlapping with the wall. To ascertain whether our simulation reproduces the correct zero-density limit we recall that for an ideal gas of hard spherocylinders the surface tension and average rod-wall contact distance are given by 104,107 lim which corresponds to γ * = Z * = 0.3 for hard spherocylinders with L/D = 10 in accordance with the limiting values apparent from Fig. 5 (a) and (c). In the limit of a perfectly ordered nematic state in planar confinement u i ·ẑ = 0 corresponding to Z * = 0.05 for our system of hard spherocylinders.

V. DENSITY FUNCTIONAL THEORY FOR A CONFINED ROD FLUID
In this section, we compare the results obtained from the Monte Carlo simulations with predictions from a density functional approach based on the seminal secondvirial theory of Onsager 108 . Although the theory is known to give quantitative results only for bulk systems in the limit of infinitely elongated rods with L/D → ∞, we expect the theory to provide a useful qualitative guide for our confined systems of hard spherocylinders with a finite aspect ratio of L/D = 10.
Within the framework of density functional theory 109,110 the grand potential Ω of a fluid of N hard spherocylinders in the presence of an external potential U ext can be expressed formally in terms of the (unknown) single-particle density ρ(s): where s = {r,û} collectively denotes both the position r and orientation vectorû characterizing the phase space of each rod. Furthermore, µ is the chemical potential, and V denotes the total thermal volume of the rods incorporating the translational and orientational contributions to the kinetic energy 111 ; the value of V will prove to be immaterial in the subsequent analysis. The first term in Eq. (13) is exact and represents the ideal translational and rotational entropy of an ensemble of rods, while the second term F ex is the excess Helmholtz free energy describing the interactions between the rods. In the spirit of second-virial theory of Onsager one can write: (14) where Φ = exp(−U r /k B T ) − 1 is the Mayer function of the pair potential U r (s, s ) between the rod particles. In the case of hard spherocylinders, Φ adopts a simple step function signalling overlap of the spherocylindrical hard cores: The slit-pore confinement of parallel hard walls considered in our simulation corresponds to an external potential which encodes the impenetrability of the hard walls. For a hard spherocylinder of length L and diameter D, the particle-wall potential takes the following form: whereẑ denotes the normal direction to the wall, and the contact distance σ between the centre-of-mass of the rod and the wall is Eq. (16) needs to be supplemented with its mirrored form U ext ( z − r ·ẑ,û) to reflect the second wall located at a distance z from the first.
Since the systems considered in our current study correspond to either isotropic or nematic phases in the bulk,  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 11 the imposition of a slit-pore geometry will break the translational invariance of the fluid along a direction normal to the wall. The particle-wall distance is denoted by z = r ·ẑ. The equilibrium profile for the single-particle density ρ eq (z,û) (with 0 < z < z ) then follows from a formal minimization of the grand potential Eq. (13) via δΩ[ρ(s)]/δρ(s)| eq = 0. The equilibrium density profile associated with the minimum of the grand potential is unique for any given set {µ, L/D, z } and can be expressed in the following self-consistent form: with β = 1/(k B T ) and ∆z = z − z . The key parameter here is the overlap area A(|∆z|,û,û ) which corresponds to the two-dimensional excluded volume swept by two spherocylinders at fixed mutual orientation (û,û ) and relative lateral distance ∆z. Noting that the Mayer function Φ depends only on the interparticle centre-of-mass distance ∆r = r − r and orientations (û,û ), the overlap area can be formally derived as where δ(∆r ·ẑ − ∆z) is the appropriate Dirac delta function.
In the case of bulk homogeneous phases, there are no inhomogeneities along the z direction and the main interaction kernel corresponds to the excluded volume v excl (û,û ) for a pair of particles with fixed mutual orientations. The excluded volume can therefore be obtained from A via a simple integration over ∆z. For hard spherocylinders this quantity is well known and can be expressed as v excl (û,û ) = − d∆z A(|∆z|,û,û ) Whilst the expression for the excluded volume is easily obtained by geometric inspection, a similar analytical derivation for A proves to be highly non-trivial and closed expression are known only up to O(LD) 104 . The full result, which includes all contributions including those of O(D 2 ) arising from the overlap of the hemispherical end caps, can be readily obtained numerically by means of a simple two-dimensional integration scheme by employing an efficient overlap routine for hard spherocylinders, similar to the one used in the simulations. It is useful to locate the bulk phase transition from the isotropic to the nematic fluid before embarking on a calculation of the density profiles of the particles from the walls. In the absence of an external field (U ext = 0) the single-particle density reduces to ρ(z,û) = ρf (û) reflecting a spatially homogeneous fluid with number density ρ and orientational distribution function f (û); the latter is normalized according to dû f (û) = 1, which may either describe an isotropic phase with f = 1/(4π), or a nematic phase with f (û ·n) peaked around a common nematic directorn. The numerical procedure has been described in detail elsewhere 112 . The coexistence densities obtained from Onsager's second-virial theory for the bulk isotropic-nematic transition of a fluid of hard spherocylinders with an aspect ratio of L/D = 10 are obtained as c I = 3.39 (η I = 0.36) for the isotropic state and c N = 3.82 (η N = 0.41) for the nematic state with a nematic-order parameter of S = 0.70; the value of the chemical potential at coexistence is µ * = βµ = 10.89 and the bulk pressure is p * = 5.853 (corresponding to the dimensionless pressure of βp(πDL 2 /4) = 54.87 employed in the theoretical studies).
Once the equilibrium single-particle density is established using Eq. (19), the fluid-wall surface tension γ can be obtained from the grand potential Ω = −pV + 2γA, where p is the bulk pressure, V is the total volume of the confined system, and A is the surface area of the wall; the factor of 2 again accounts for the presence of two walls. Note that, here, p is the scalar pressure deep in the bulk where, locally, p = p N = p T . Inserting the equilibrium density profile from Eq. (19) into Eq. (13) one can express the minimum of the grand potential as In the particular case of the bulk homogeneous system, the minimum of the grand potential Ω b,eq can be expressed in terms of the equilibrium single-particle orientational distribution function f eq (û) as The difference between the equilibrium grand potential of the inhomogeneous system and the bulk homogeneous system provides a route to calculate the fluid-wall surface tension: γ = (Ω eq − Ω b,eq )/(2A), again expressed in dimensionless form as γ * = γv m /(k B T η b L) as for Eq. (5). The surface adsorption Γ is readily established from a simple spatial average of the difference between the density of the inhomogeneous system and the bulk homogeneous system over the range of the z direction in the slit-pore geometry: which is equivalent to Eq. (10) employed in the analysis of the simulation data. The average rod-wall contact distance Z can be determined from the equilibrium density profile ρ eq (z,û) with expressed in units of L as for Eq. (11). The self-consistency equation Eq. (19) for equilibrium single-particle density is solved by discretising the three-dimensional space {z,û} on an evenly spaced grid. The two-dimensional orientational phase space dû = dθ sin θdϕ = 4π is routinely parameterized in terms of a polar angle 0 < θ < π and an azimuthal angle 0 < ϕ < 2π, each discretized into 60 evenly spaced grid points. In order to enable comparison with the simulation data the wall-to-wall distance is fixed at z = 43 D. Noting that the profiles are characterized by planar symmetry with respect to the walls, i.e., ρ(z,û) = ρ( z −z,û), we need to consider only half of the spatial interval and the relevant range of values of z: 0 < z < z /2 is compartmentalized into 30 equally spaced sections per spherocylinder length L. Finer grids are found to only marginally improve accuracy whilst greatly increasing the computational burden. The results for a non-interacting system, equivalent to the limit of zero bulk density, ρ b ↓ 0, are readily recovered from Eq. (19) by noting that the exact density profile of an ideal gas exposed to an external wall potential, cf. Eq. (16), corresponds to the Boltzmann factor ρ(z,û) ∝ exp[−βU ext (z,û)]. Substituting this expression into Eq. (21), omitting the first contribution of O(ρ 2 ), and performing the spatio-orientational integration we readily obtain the expression for the fluid-wall surface tension and average particle-wall contact distance for an ideal gas in Eq. (12). The theoretically determined fluid-wall surface tension shown in Fig. 6 (a) reveals a non-monotonic trend with the bulk density, similar to what is observed with the simulation data (cf. Fig. 5 (a)). The density functional theory is seen to overestimate the jump in the surface tension for the transition from the isotropic to the nematic branch. This feature is intimately linked to the fact that the predicted density gap between the isotropic and nematic phases (and the transition pressure) obtained with the second-virial theory of Onsager is generally too large in comparison to the essentially exact simulation data. Moreover, the actual values of the tension are somewhat lower than the corresponding simulated values obtained by simulation. This discrepancy could be attributed to the inadequacy of the second-virial coefficient in describing correlations in confined high-density fluids of rod-like particles where the local density may become very high (cf. the peaks in the density profiles of Fig. 2). Despite these differences, the overall qualitative trends for the dependence of the surface tension, adsorption, and average rod-wall contact distance on the bulk density are in good agreement with the findings from simulation. The sharp drop in the fluid-wall surface tension and the concomitant increase of the adsorption indicate the formation of a nematic wetting layer predicted theoretically at c b ∼ 2.6. Recalling that the onset of the bulk isotropic-nematic transition is at c b = 3.39, the wall-induced nematiza-tion occurs when the concentration exceeds a threshold of 77 % of the value at the bulk phase transition. Close to the bulk isotropic-nematic transition the theory predicts that the surface tension becomes negative, which implies that the free energy of the system actually decreases upon increasing the surface area of the wall in contact with the rod fluid. The negative values of the wall-fluid interfacial tension can be qualitatively understood as a consequence of the formation of a nematic wetting layer at the wall -it becomes thermodynamically favorable to create a nematic-wall (and nematic-isotropic) interface. The tension determined with the second-virial theory in the vicinity of the isotropic-nematic transition is that of the nematic wetting film in contact with the wall rather than that between an isotropic fluid and the wall. The negative values obtained for the nematic film-wall tension (cf. Fig. 6 (a)) indicate the propensity of the particles to spread on the surface. The negative tension obtained with the theory may be due to a finite-size effect in the numerical treatment; it is apparent from the density profiles depicted in Fig. 7 (a) that the bulk isotropic density is not fully reached at large particle-wall distances, particularly at wetting conditions where there is a fastgrowing nematic layer forming at the wall. It is also possible that the effect is an artefact of the second-virial approximation in the theory where higher-body interactions between the rods are neglected. It is important to realize that the theory is highly approximate and the treatment clearly becomes inadequate close to the wetting transition.
The location of the wetting transition can be estimated more precisely from the behaviour of the uniaxial nematic and biaxial nematic order parameters at the wall surface and at the centre of the slab, as depicted in Fig. 8. This suggests that the onset of biaxial order in the vicinity of the wall coincides with a capillary nematization (the transition from a bulk isotropic phase, which wets confining surfaces with a biaxial film, to a condensed nematic 10,11,14 ). This behaviour is somewhat different from what is observed in the so-called Zwanzig model in which rod orientations are constrained to point along one of the three Cartesian axes 10 . In the case of the Zwanzig model the uniaxial-biaxial surface ordering transition is found to pre-empt capillary nematization. In line with the simulation results, marked optima are observed in both the fluid-wall surface tension and adsorption at bulk concentrations below the wetting transition. The extrema do not, however, coincide since the minimum in the adsorption occurs at c b ∼ 1.7 whereas the surface tension exhibits a maximum at c b ∼ 1.9. We remark that as long as z L the surface structure is expected to be fairly insensitive to a variation of the wall-to-wall distance. We have verified this by repeating the theoretical calculations for a wall-wall separation of z = 80D, and indeed find that the profiles and surface thermodynamic properties deviate negligibly from those for z = 43D, indicating that the two walls influence the properties of the fluid virtually independently from each other at these 14 slit widths.
Although the overall trends of the surface ordering and wetting scenario found with the simulation are appropriately captured by the simple second-virial approach adopted here, some important discrepancies are apparent. The main issue pertains to the microstructure of the wetting layer. Whilst the simulated density profiles (cf. Fig. 2(a)) are seen to exhibit a sharp peak in density close to the wall, hinting at the formation of a coherently structured nematic monolayer, the theoretically determined density profiles (cf. Fig. 7(a)) appear much more diffuse at these high isotropic bulk densities. More importantly, the sharp drop in the surface tension which eventually becomes negative close to the bulk isotropicnematic transition is not in accordance with the findings of the simulations. This could point to a fundamental deficiency of the simple second-virial theory in capturing the highly anisometric many-body correlations that prevail in the high-density confined fluid of rod-like particles. Clearly, more realiable approximations are required in the development of accurate density functionals 113,114 for confined fluid of anisometric particles in order to resolve this issue allowing for more quantitatively reliable theoretical predictions of the fluid-wall surface tension and related properties.

VI. SUMMARY AND CONCLUDING REMARKS
We present a detailed study of the surface behaviour and interfacial thermodynamics of a fluid of hard spherocylinders with aspect ratio of L/D = 10 confined in a slit-pore geometry by means of Monte Carlo computer simulation. Anisotropic volume perturbations are employed to determine the normal and tangential components of the pressure tensor of the inhomogeneous fluid as a function of density; the full data is reported in Table I. Accurate knowledge of the components of the pressure tensor and density profiles allows for an efficient calculation of the fluid-wall surface tension and related quantities such as the surface adsorption and average rod-wall contact distance, over a wide interval of bulk densities ranging from the low-density isotropic fluid to dense nematic states (for the first time). The results unveil a marked non-monotonic trend in the surface tension as a function of bulk density, with a pronounced maximum just below c b ∼ 2. The latter density indicates a surface phase transition corresponding to a change in the local orientational order of the rod particles in the vicinity of the wall from uniaxial nematic to biaxial nematic. This behaviour also signals the onset of the formation of a nematic wetting film at the wall surface. The methodology could also be employed to determine the interfacial tension between the wall and bulk smectic phases at higher densities.
The novel data for the fluid-wall surface tension, reported in Table II, serves as valuable benchmarks for testing more sophisticated theoretical predictions from den-sity functional theories such as those based on modifiedweighted density functionals 115 or fundamental measure theory 116 . The simplest non-trivial functional, based on a second-virial theory of Onsager (which is accurate at low density) is assessed here and found to give a reasonable qualitative account of the evolution of the surface microstructure such as the symmetry-breaking of uniaxial nematic order at the wall surface, the capillary nematization, and the non-monotonic trend of the surface tension with bulk density for the isotropic bulk system.
The negative (tensile) values of expansive contributions to the pressure tensor (p * + N and p * + T ) of the hard-rod fluid seen in Table I have a subtle, but potentially far reaching, implication: when the particles are anisotropic, the system as a whole possesses an expansive tensile strength, even in the complete absence of attractive interactions. In turn, this implies that a unidirectional stretching deformation of a liquidcrystalline fluid (either along or perpendicular to the nematic-ordering director) gives rise to an asymmetrical elastic response. It appears that this tensile strength correlates with the orientation of the nematic director. This provides a possible explanation for the significantly larger error bars found for the tangential components of the pressure tensor (p * + T and p * − T ) of the nematic phase when compared to the corresponding values of the normal components (p * + T and p * − T ). To address this issue and obtain an equivalent level of precision, the number of configurations generated must be large enough for the nematic-ordering director to rotate such that it is uniformly sampled in the plane parallel to the interfaces. We note that the errors associated with the overall components of the pressure tensor p * N and p * T , in Table II) are significantly smaller than those of their constituent components (p * + N , p * − N , p * + T and p * − T ). When summed together for each specific state, the positive and negative volume perturbations uniquely capture the total residual resistance to the volume deformation. The asymmetry in the expansive and compressive deformations will also be reflected in the elastic constants (twist, bend, and splay) of the system pertaining to a deformation of the director field (at constant volume).

APPENDIX A: PERTURBATION METHOD FOR THE DETERMINATION OF THE FLUID-WALL INTERFACIAL TENSION
Inhomogeneities in confined fluids cause imbalances between the diagonal elements of the pressure tensor, creating tension at the interfaces between unlike phases. For fluids at equilibrium (with no shear forces present) confined in slit-pore geometry p T = (p xx (z) + p yy (z))/2, where (p xx (z) + p yy (z), and p N = p zz is constant and independent of the distance from the confining walls 38,92 . The expression originally presented by Kirkwood and Buff 117 can be used to used to determine the surface tension: where F and A are the Helmholtz free energy and interfacial area, respectively. For systems with planar interfaces 92 dz p αβ (z) = z p αβ . Conveniently, this allows one to express the wall-fluid surface tension of our systems in terms of the anisotropy of the normal and tangential components of the pressure tensor: The separate components of the pressure tensor are calculated using the test-volume perturbation method presented in Ref. 94 . This involves undertaking a series of anisotropic volume deformations during the course of a standard Monte Carlo simulation. At periodic intervals, non-permanent perturbations are made to the volume V i of configuration i by scaling one axis only; separate perturbations are made to reduce and increase the system volume anisotropically. Evenly spaced values are chosen for ∆V α i→j+ and ∆V α i→j− , the expansive and compressive volume changes, respectively. The average values of p N and p T are each calculated during the course of the simulation from [94]: Here P + nov is the probability that no particle-particle and/or particle-wall overlaps occur upon a volume expansion of ∆V α i→j+ , and likewise P − nov is the probability that no overlaps occur for a volume compression of ∆V α i→j− . Evaluating the limits for infinitesimal volume changes using a mean-squares correlation gives us the excess contributions to the pressure tensor.

16
TABLE I. Results for systems of N hard spherocylinders with an aspect ratio of L/D = 10 confined between two hard walls with a separation distance z = 43D obtained from canonical Monte Carlo (MC-N V T ) simulation. The packing fraction η b , concentration c b , and nematic-order parameter S b in the bulk region of the system are reported. The dimensionless components of the pressure tensor are represented by p * αα , where αα is the axis of the appropriate component. The excess contributions to the normal and tangential components of the pressure tensor (denoted by the subscripts N and T , respectively) relative to the wall surface are calculated from anisotropic test-volume perturbations in the xx, yy, and zz axes. p * + αα and p * − αα are the excess contributions from expansive and compressive changes, respectively, where p * + N = p * + zz , p * − N = p * − zz , p * + T = p * + xx + p * + yy /2 and p * − T = p * − xx + p * − yy /2. TABLE II. Results for systems of N hard spherocylinders with an aspect ratio of L/D = 10 confined between two hard walls with a separation distance z = 43D obtained from canonical Monte Carlo (MC-N V T ) simulation. The packing fraction η b , the concentration c b , the nematic-order parameter S b in the bulk region of the system, the surface adsorption Γ * , the average particle-wall contact distance Z * , the normal p * N and tangential p * T components of the pressure tensor, and the fluid-wall surface tension γ * are reported for a range of densities from the low density isotropic state to the high density nematic.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60