text:
# Breaking of Nanotube Symmetry by Substrate Polarization

## 1 Introduction

Since discovery of carbon nanotubes in 1991, a deep physics of these one–dimensional nanoscale objects has been demonstrated. Fundamental properties of the nanotubes have been studied in view of possible applications in electronics and other devices. A detailed theoretical description for electronic structure of ideal single–wall nanotubes (SWNTs) was obtained as well as for an effect of various defects and disorder on the SWNT electronic properties ( *e.g.*, Ref.). However, the nanotube systems under experimental study often deviate from a model picture. In this letter we investigate one aspect of a real system: the modification of the electronic properties of the SWNT deposited on a substrate. One expects that symmetry of the nanotube at the substrate will be lower than the symmetry of the nanotube itself in vacuum.

A description of a breaking of the symmetry of SWNT bandstructure due to a charge transfer (or charge injection) between the nanotube and the substrate (or contacts) and calculation of a polarization of the substrate and the nanotube, which follows to the charge transfer, are the goals of our study. Effects of splitting, mixing, and/or anti–crossing of the nanotube subbands, that are caused by the depolarization of the electron charge density, have been almost neglected in literature before. We use term “depolarization” for a number of phenomena including a transverse shift of the electron charge density from its equilibrium distribution profile (effects due to an axial/longitudinal depolarization were discussed elsewhere). We will show that the *transverse depolarization* results in qualitative changes of the nanotube density of states (DOS) near van Hove singularities. In particular, we predict the *splitting* of a doublet state to be likely observable as a function of the injected/induced charge density of the SWNT. We will discuss that in Sec.2.1. In Sec.2.2 we will calculate this injected/induced charge density in a selfconsistent way.

The depolarization and *intrasubband* splitting will be studied for a typical experimental situation: a single SWNT lies on a conductive substrate or separated from the conductor by a thin insulating layer representing an oxide on the surface of a metal. We assume that the nanotube is connected to electron reservoirs, which may be the leads or the conductor substrate itself. A transverse external electric field and/or a work function difference between the SWNT and the substrate/contact induce non–zero electron/hole charge density in the nanotube. This extra charge density polarizes the substrate, which breaks the axial symmetry of the nanotube. This effect is much larger than an electronic structure perturbation caused by the lattice distortion which may happen due to a van der Waals attraction to the substrate. We will demonstrate that a direct action of the *uniform* external electric field is of minor importance as compared to the *nonuniform* field of surface charges on the substrate. We will discuss a modification of our theory of the subband splitting for a case of purely insulating substrate in the last section.

## 2 Perturbation Theory for Bandstructure Modification

### 2.1 Splitting of SWNT subband due to transverse depolarization

To calculate the splitting and shift of the electron energy levels one needs to know matrix elements of the perturbation potential between corresponding wave–functions. In our case, the perturbation is a selfconsistent Coulomb potential (operator in Heisenberg representation) which describes the interaction between the probe electron and the extra charge density on the SWNT and the polarization charge density on the substrate surface:

$`\widehat{V}=e{\displaystyle \underset{L/2}{\overset{L/2}{}}}dZ{\displaystyle \underset{0}{\overset{2\pi }{}}}Rd\beta ({\displaystyle \frac{\widehat{\sigma }}{\sqrt{(zZ)^2+(R\mathrm{cos}\alpha R\mathrm{cos}\beta )^2+(R\mathrm{sin}\alpha R\mathrm{sin}\beta )^2}}}`$
(1)
$`+{\displaystyle \frac{\widehat{\sigma }^{}}{\sqrt{(zZ)^2+(R\mathrm{cos}\alpha R\mathrm{cos}\beta )^2+(R\mathrm{sin}\alpha 2hR\mathrm{sin}\beta )^2}}}).`$
Both the probe electron and the nanotube surface charge are taken on a cylinder of a radius $`R`$. Then $`z`$ and $`\alpha `$ are the electron coordinates in the cylindrical coordinate system. $`\sigma `$ is the surface charge density. It does not depend on the coordinate $`Z`$ along the nanotube because we assume the translational invariance of the problem for clarity of derivation. Although, the theory can be easily extended for the case of slow variation of $`\sigma `$ along the axis. We will show later that one can drop dependence of $`\sigma `$ on the angle $`\beta `$ along the circumferential direction in approximation of a *linear response* (in higher orders of perturbation theory a direct transverse polarization must be taken into account). $`\sigma ^{}`$ is an image charge density which is equal to $`\sigma `$ for the metallic substrate.

The first term of Eq.(1) is the interaction with the charge density on the nanotube, which coincides with the Hartree term for the SWNT in vacuum (without charge injection). The second term in Eq.(1) has also a simple physical meaning: this is the energy of interaction of the electron with the image charge. The separation between the SWNT axis and the surface of the conductor is $`h`$. In case of the metallic substrate it is about the nanotube radius, $`R`$, plus the van der Waals distance for graphite: $`hR+0.34`$ nm.

The matrix element of the Coulomb operator (1) is calculated with the wave–functions of a tight–binding (TB) Hamiltonian. We use envelope wave–functions, obtained similarly to Ref.. This approach has been widely used in the literature, so we skip details and give the wave–functions in the one–band scheme ($`\pi `$ electrons only) in the form:

$$|\psi _{m,k,\zeta }=\frac{1}{\sqrt{2}}(|A+\zeta c_{mk}|B)e^{ikz}e^{im\alpha },$$
(2)
here index $`m`$ labels subbands of the SWNT electronic structure, $`k`$ is a longitudinal momentum, these two are good quantum numbers (discrete and continuum, respectively) for an ideal, long enough nanotube; $`\zeta =\pm 1`$ is a pseudospin. (A pseudospinor vector is formed by a two–component wave–function amplitude defined for two atoms in a graphite unit cell, A and B). Coordinate along the tube is $`z`$, and $`\alpha `$ is the angle along the nanotube circumference.

We assume that our potential is smooth at the scale of the single unit cell (0.25 nm). Then, one may neglect transitions with the pseudospin flip (transitions between sublattices). With use of the orthogonality relation between the spinor components, it yields:

$$m|V|n=\frac{8\pi eR\sigma }{|mn|}i^{mn}\left(\frac{R}{2h}\right)^{|mn|},mn$$
(3)
$$m|V|m=4\pi eR\sigma \mathrm{log}\left(\frac{2h}{R}\right),m=n$$
(4)
where $`\sigma `$, the surface charge density, has to be defined later in a selfconsistent way.

The Eqs.(3,4) are obtained by a direct Fourier transformation of (1) and describe the energy level shift when $`m=n`$ and the mixing of different subbands at $`mn`$. The most interesting term with $`n=m`$ is the mixing between the degenerate electron states within the same subband. By solving a secular equation for the intrasubband mixing of the electron doublet we obtain the splitting of the van Hove singularity at the subband edge (Fig.2). The new subband energy separation reads as:

$$\delta E_m=\frac{8\pi eR\sigma }{m}\left(\frac{R}{2h}\right)^{2m}.$$
(5)
Let us now calculate the injected/induced charge density $`\sigma `$ which will allow us a numerical estimation for the $`\delta E_m`$ splitting.

### 2.2 Charge injection due to the Fermi level shift

The Eqs.(35) are written for the given charge density $`\sigma `$ which has to be derived in this section. When the SWNT is placed in a real device, one must consider the work function difference between the nanotube and the contact or the conducting substrate and/or the external potential which may be applied to the SWNT. The potential shifts the Fermi level in the SWNT. As a result of this the positive/negative charge is injected into the nanotube:

$$\sigma =\frac{e}{2\pi R}\underset{0}{\overset{\mu (\sigma )}{}}\nu (E)𝑑E,$$
(6)
here $`\nu (E)`$ is a bare one–dimensional DOS (independent of $`\sigma `$ in a linear response theory); $`\mu =\mathrm{\Delta }We\phi ^{xt}e\phi ^{ind}(\sigma )`$ is the shift of the electrochemical potential of the SWNT (with respect to a charge neutrality level $`E=0`$) which depends on the work function difference, $`\mathrm{\Delta }W`$, on the external potential, $`\phi ^{xt}`$, applied between the nanotube and the reservoir and on the potential $`\phi ^{ind}`$ induced by the charge density of the nanotube, $`\sigma `$. This last term is proportional to the intrasubband term ($`m=n`$) of the Coulomb interaction given by the Eq.(4).

This selfconsistent equation for $`\sigma `$, (6), is readily solved analytically if the electrochemical potential is below the second subband edge. We follow Ref. in derivation of $`\sigma `$: the induced potential is obtained by direct integration of the charge density along the SWNT as in . As it is shown in Figure 1, for a metallic SWNT, the charge is a product of the constant DOS, $`C_Q`$, and the electrochemical potential, $`\mu `$. Then the solution of Eq.(6) is as follows:

$$\sigma _A=\frac{\mathrm{\Delta }We\phi ^{xt}}{2\pi Re\left(2\mathrm{log}(2h/R)+C_Q^1\right)},$$
(7)
and for a semiconductor SWNT, which has the DOS $`C_QE\theta (E\mathrm{\Delta })/\sqrt{E^2\mathrm{\Delta }^2}`$, the charge is as follows:

$$\sigma _Z=\sigma _A\frac{\sqrt{\left(\frac{\mathrm{\Delta }}{\mathrm{\Delta }We\phi ^{xt}}\right)^2\left(4\mathrm{log}^2(2h/R)C_Q^2\right)+C_Q^2}2\mathrm{log}(2h/R)}{2\mathrm{log}(2h/R)C_Q^1}\theta (\mathrm{\Delta }We\phi ^{xt}\mathrm{\Delta }).$$
(8)
Here $`\theta (x)`$ is the Heaviside unit step function and $`\mathrm{\Delta }`$ is $`1/2`$ of the energy gap. We introduced a quantum capacitance of the SWNT following Ref.:

$$C_Q=\frac{8e^2}{3\pi b\gamma }$$
(9)
which is the one–dimensional analog of the quantum capacitance proposed for a two–dimensional electron gas system by Luryi. Here $`b1.4`$Å is the interatomic distance, $`\gamma 2.7`$ eV is the hopping integral for the graphite–like systems. We notice that despite the $`\sigma _Z`$, as given by Eq.(8), comes from a massive subband (in contrast to $`\sigma _A`$, as in Eq.(7) where the lowest subband is massless, see in Fig.1), the linear dependence of $`\sigma _{A/Z}`$ on $`\phi ^{xt}`$ preserves as long as the potential $`\phi ^{xt}`$ is large enough. This reflects the fact that a classical one–dimensional charge density is a linear function of a classical electrostatic potential .

## 3 Results and Discussion

In the last section we obtained the selfconsistent expression for the surface charge density as a function of the external potential and the work function difference which may be considered as a built–in potential. Substituting Eq.(7) into Eq.(5) we obtain the splitting of the degenerate subbands $`|\pm m`$ of the metallic SWNT (when the Fermi level is within the first subband) as follows:

$$\delta E_m=\frac{4\left(\mathrm{\Delta }We\phi ^{xt}\right)}{m\left(2\mathrm{log}(2h/R)+C_Q^1\right)}\left(\frac{R}{2h}\right)^{2m}.$$
(10)
The splitting decreases with $`m`$ exponentially, hence, the effect is likely observable for the lowest degenerate subband. Then, for the following parameters: the SWNT radius $`R6.7`$ Å, the distance to the metal substrate $`h=10.1`$ Å, and the quantum capacitance $`C_Q^10.69`$, we obtain a numerical estimate for the subband splitting $`\delta E_10.15\left(\mathrm{\Delta }We\phi ^{xt}\right)`$. Experimental data for the work function of SWNTs scatters from 4.9 to 5.05 eV . For the SWNT on the gold substrate we use as an estimate $`\mathrm{\Delta }W0.3`$ eV. In absence of the external potential, this work function difference results in $`46`$ meV gap between two split peaks of the density of states (Fig.2), which is larger that $`kT`$ at room temperature. We also calculated the contribution of all other subbands, which is negligible in the splitting, but it shifts the doublet as a whole. As a result, two new peaks in the Fig.2 appear not symmetrical with respect to the original DOS singularity.

The splitting of $`\pm m`$ doublet is an analogue of a degenerate level Stark effect for the nanotube in a multipole potential of the image charge. The lower subband has $`x`$ symmetry and the upper subband has $`y`$ symmetry (with corresponding wave–functions $`|x=1/\sqrt{2}(|+m+|m`$ and $`|y=1/\sqrt{2}(|+m|m`$) because of an attraction energy of the electron to its image charge is lower for the second combination.

We predict a similar effect for the semiconductor nanotube, although, the total external potential causing the charge density injection must be larger than one half of the gap in this case. As we study in this paper only the effect which is linear in the external potential, all high order terms in Eq.(8) have to be discarded.

### 3.1 Dipole polarization correction

The charge injection in the nanotube may be readily achieved by applying external electric field. One may naively argue that the external field itself can break the bandstructure symmetry and result in some level splitting. Although, it is correct statement in general, the direct splitting of the SWNT orbital doublet $`\pm m`$ by the uniform electric field is forbidden by symmetry. These degenerate states do not mix together due to the selection rules of the problem. The matrix element for an intrasubband splitting in the uniform external field $`_{xt}`$ equals zero by parity: $`m|e_{xt}y|m=0`$.

In order to calculate the subband splitting in this case we have to compute the charge injection, which is proportional to the applied field. The potential is equal to $`\phi ^{xt}=_{xt}h`$, where $`h`$ is the distance between the axis of the tube and the metal surface, which has to be substituted into Eq.(10).

In Sec.2.1 we assumed that the charge density $`\sigma `$ has no dependence on the angular coordinate $`\beta `$ along the nanotube circumference. This is an accurate approximation since a dipole (and higher multipole) component of the $`\sigma `$ is small as compared to what is given by Eqs.(7,8). Let us prove this assumption for the dipole polarization of the SWNT.

The non–uniform external potential causes a deviation of the surface density from the uniform equilibrium value, $`\sigma `$, which is given by following expression:

$$\delta \sigma (\beta )=e\underset{ij}{}\frac{(f_if_j)i|V|j}{E_iE_j}j|\beta \beta |i,$$
(11)
where $`f_i`$ are the occupation numbers, the matrix element $`i|V|j`$ is given by Eqs.(3,4), $`E_i`$ are the energies of subbands and $`\beta |i`$ are corresponding wave–functions.

We define the nonuniform dipole part of the charge density of a SWNT as $`\delta \sigma _1\underset{0}{\overset{2\pi }{}}\mathrm{sin}\beta \sigma (\beta )𝑑\beta `$. Then, the dipole component of the surface charge is as follows:

$$\delta \sigma _1=\frac{ie}{8\pi R^2}\underset{i}{}\frac{(f_if_{i\pm 1})i|V|i\pm 1}{E_iE_{i\pm 1}}.$$
(12)
Let us remind that according to Eq.(3) $`i|V|i\pm 1=i8\pi R^2e\sigma /(2h)`$.

In case of the electrochemical potential equals zero (no charge in the nanotube), the transverse polarization includes transitions from the valence to conduction band only $`v|V|c`$ (the details of the calculation are presented elsewhere). Here, we study an extra component of the polarization, which is due to the induced charge density. Thus, we need to consider only transitions from the levels above the charge neutrality level, $`E=0`$, and below the Fermi level, $`E=E_F`$, (which is the shaded area in Fig.1). Hence, the dipole polarization is proportional to the net charge density $`\sigma `$, and the dipole charge density of the armchair SWNT is given by the following expression:

$$\delta \sigma _1=\frac{\sqrt{3}C_Q^2}{32\pi }\frac{\left(2\pi R\sigma _A\right)^2}{e}\frac{R}{h}\mathrm{log}\frac{2h}{R}_{xt}^2$$
(13)
where $`C_Q3.2`$ is the dimensionless quantum capacitance. We single out term $`2\pi R\sigma _A`$, which is the specific one–dimensional charge density of the SWNT and proportional to the external potential and, thus, to the external field.

The Eq.(13) shows that the effect of the transverse polarization on the bandstructure is quadratic in the external field, in good agreement with a plain dielectric response theory, while the effect of the image charge is linear in $`_{xt}`$. Thus, the degenerate level splitting due to the dipole component of the polarization will be less important than the splitting due to a uniform component: $`\sigma _0_0^{2\pi }\sigma (\beta )𝑑\beta `$, at least, in a weak field regime discussed in the article. This proves post factum our assumption of $`\sigma `$ to be independent of $`\beta `$.

### 3.2 Depolarization at the insulator substrate

For the sake of completeness we present here also a modification of our theory to the case of a dielectric substrate. In this case the screening of the charge density in the nanotube is weaker. It results from (i) underscreening of the Coulomb interaction between the nanotube carriers and (ii) lower charge density induced in the substrate. The second factor can be taken into account by substituting an effective image charge density $`\sigma ^{}=\sigma \frac{1\epsilon }{\epsilon +1}`$ in the second term of Eq.(1), where $`\epsilon `$ is the dielectric function of the substrate (in case of highly conductive substrate it equals $`\mathrm{}`$), instead of the bare image charge density $`\sigma `$. This results in substituting $`\sigma ^{}`$ in Eqs.(38) where appropriate.

Now, the fields of the image charge and the charge in the SWNT do not cancel each other in contrast to the case of the metallic substrate. As a result, the underscreening of the Coulomb interaction happens. This modifies the equations for the energy level shift ( *intrasubband* matrix elements as in Eq.(4)) and, thus, the electrochemical potential shift. One must substitute $`\mathrm{log}(2h/R)`$ term everywhere by $`\mathrm{log}(2h/R)+2/(\epsilon 1)\mathrm{log}(L/R)`$ where $`L`$ is the length of the nanotube (or distance between metal leads to it). This expression diverges with the length of the nanotube which reflects the one–dimensional character of the Coulomb interaction. These changes have to be made through Eqs.(710).

The first term of Eq.(1) does not appear in the calculation of the *intersubband* matrix elements as in Eqs.(3,5). Hence, no additional correction is required in the equations (1112) of the last section.

We assumed in this paper that the perturbation theory in a linear approximation in $`\mu `$ (or equivalently in $`\sigma `$) is applicable. Restrictions which may follow from this assumption are as follows. The external potential has to be small. We neglect here the dipole term in the induced charge density (and higher multipoles as well). It is equivalent to a weak intersubband mixing which assumption may not hold for wide nanotubes or strong external fields. The effect of the strong field on the bandstructure is discussed elsewhere. In this paper we used the Eq.(6) for the equilibrium charge density in the SWNT. One may consider transport devices on equal basis, as long as the charge of the nanotube is still given by the quasi–equilibrium charge density. However, for non–zero current flowing through the nanotube, one must use an expression for the charge density which differs from Eq.(6) (to be discussed elsewhere ).

## 4 Conclusions

In summary, we have developed a microscopic quantum mechanical theory for a charge transfer between a SWNT and a conductive substrate (and/or metallic leads). This charge injection results from a natural work function difference between the nanotube and the substrate or/and from an external potential applied between those. A surface charge density of the SWNT is calculated selfconsistently within an envelope function formalism of tight–binding approximation.

We demonstrated for the first time that the influence of this charge transfer on the electronic structure of the SWNT is not negligible for typical material parameters of the problem. Because of the breaking of the axial symmetry of the system, the SWNT DOS changes qualitatively: degenerate subbands $`\pm m`$, where $`m0,n`$, split. It has a simple physical interpretation — the electrons with $`x`$ and $`y`$ polarizations are no longer equivalent as their attraction to the substrate is different. This effect can be related to a degenerate level Stark effect with an appropriate choice of external field of the image charge. The gap between the new $`x`$ and $`y`$ subbands is constant in $`k`$–space (for the external field which is uniform along the tube) so it shows up dramatically at the subband edge. The van Hove singularity splits, and the distance between two peaks of the DOS is about 46 meV for the  armchair SWNT on the gold substrate.

We obtained analytical expressions for matrix elements of the image charge field, which yields the mixing of different subbands and can be used to describe the level anticrossing (to be discussed elsewhere). Same matrix elements enter the expression for the multipole polarizabilites of the SWNT. We estimated a major contribution to the dipole polarizability of the armchair SWNT, which comes from intraband transitions for non–zero charge injection. The analytical expression for the dipole component of the surface charge density is shown to be proportional to the square of the external potential and, hence, appears in the second order of the perturbation theory which corroborates post factum our assumption of uniformity of the induced/injected charge along the SWNT equator.

We show that the modification of our theory to the case of semiconductor substrate is straightforward. The analytical expressions for the van Hove singularity splitting and induced charge density are obtained.

Acknowledgments. S.V.R. acknowledges support of DoE through grant DE-FG02-01ER45932, NSF through grant ECS-0210495 and Beckman Fellowship from the Arnold and Mabel Beckman Foundation. A.G.P. is grateful to the Beckman Institute for hospitality during his work in Urbana. Authors are indebted to Professor K. Hess for valuable discussions.


format: text

text:
# Smectic ordering of parallel hard spherocylinders: An entropy-based Monte Carlo study

## I Introduction

Recently we investigated the formation of liquid-crystalline mesophases in a fluid of hard spherocylinders, i.e., cylindrical segments of length $`L`$ and diameter $`D`$ capped at each end by a hemisphere of the same diameter micali . More specifically, we analyzed the ordering of the homogeneous and isotropic fluid into a nematic, smectic or solid phase within the framework provided by the one-phase entropy-based criterion originally proposed by Giaquinta and Giunta gg . This criterion can be implemented through the well known multi-particle correlation expansion of the configurational entropy ng :

$$s_{\mathrm{ex}}=\underset{n=2}{\overset{\mathrm{}}{}}s_n.$$
(1)
In the above formula $`s_{\mathrm{ex}}`$ is the excess entropy per particle in units of the Boltzmann constant and $`s_n`$ is the $`n`$-body entropy that is obtained upon re-summing spatial correlations between up to $`n`$ particles. In particular, the pair entropy per particle of a homogeneous and isotropic fluid of nonspherical molecules can be written as lapa :

$`s_2(\rho )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\rho }{\mathrm{\Omega }^2}}{\displaystyle }\{g(𝐫,\omega ^2)\mathrm{ln}[g(𝐫,\omega ^2)]`$ (2)
$`g(𝐫,\omega ^2)+1\}d𝐫d\omega ^2,`$
where $`g(𝐫,\omega ^2)`$ is the pair distribution function (PDF) which depends on the relative separation $`𝐫`$ between two molecules and on the set of Euler angles $`\omega ^2[\omega _1,\omega _2]`$ that specifies their absolute orientations in the laboratory reference frame. The quantity $`\mathrm{\Omega }`$ represents the integral over the Euler angles of one molecule, while $`\rho `$ is the particle number density.

In order to identify the ordering threshold of the fluid, we monitor the behavior — as a function of the number density — of the so-called residual multiparticle entropy (RMPE). This quantity is defined as the difference

$$\mathrm{\Delta }ss_{\mathrm{ex}}s_2.$$
(3)
At variance with the pair entropy, the RMPE exhibits a nonmonotonic behavior as a function of $`\rho `$. In particular, it is negative at low densities, and becomes positive as a more ordered phase is approached. The relevance of the condition $`\mathrm{\Delta }s=0`$ as a one-phase ordering criterion has been documented for a variety of phase transitions, both in continuous fluids as well as in lattice-gas model systems spg1 .

Even for systems composed of non-spherical molecules, such as hard spherocylinders, it turns out that the ordering thresholds detected through the zero-RMPE condition systematically correlate with the corresponding phase-transition points, whatever the nature of the higher-density phase coexisting with the isotropic fluid micali . Cuetos and coworkers have successfully applied the RMPE criterion to hard spherocylinders with an attractive square well and to spherocylinders with a soft repulsive core cuetos .

In this paper we intend to analyze the predictions of the RMPE approach in the case of the nematic-smectic transition undergone by parallel hard spherocylinders with aspect ratio $`L/D=5`$. The onset of smectic order out of a nematic phase represents the next step in the process which, upon compression of the isotropic fluid, eventually leads to the formation of the fully crystalline solid. The phase behavior of this model has been investigated with several numerical simulations slf ; verfre ; koda  as well as theoretical studies (see e.g. hol ; tay ; baus  and references contained therein) over the whole $`L/D`$ range. Such studies have ascertained the existence of a second-order phase transition from a low-density nematic state to an intermediate smectic phase. On the other hand, the general features of the higher-density region of the phase diagram has long been debated as far as the stability a columnar phase is concerned verfre .

The paper is organized as follows: In Sections 2 and 3 we describe the theoretical framework and the numerical simulation technique. The results are presented in Section 4 while Section 5 is devoted to concluding remarks.

## II Theory

The formalism developed in micali  for the isotropic fluid needs to be modified in order to describe the nematic and smectic mesophases. Even for uniaxial molecules, the numerical computation of the full pair distribution function $`g(𝐫,\omega ^2)`$ is a formidable task. Costa and coworkers have already shown that, in order to reproduce the phase boundaries of the isotropic fluid, it is enough to take into account the dependence of $`g(𝐫,\omega ^2)`$ on the centers-of-mass separation $`r`$ and on the angle $`\theta `$ formed by the molecular axes of two spherocylinders micali . This angle shows up as the critical parameter which accounts, by itself, for the reduction of orientational states that is ultimately offset by the gain of translational entropy at high densities ons . On the other hand, it is rather obvious that the angle $`\theta `$ conveys no useful information on the structural process that may eventually lead to the formation of a smectic phase out of a nematic phase.

In order to resolve ordering effects associated with the orientational degrees of freedom on one side, and with the modulation of the density along the nematic director on the other side, we shall restrict our analysis to a system composed of parallel spherocylinders. Assuming that the molecules are perfectly aligned does not significantly alter the main aspects of the phenomenology that we want to investigate. As emphasized above, different mechanisms drive the phase transformation of an isotropic or of a nematic fluid into a more ordered phase. Moreover, at the nematic-smectic transition threshold the nematic fluid is typically characterized by a very high degree of relative alignment of the molecules. For example, the nematic order parameter is about 0.85-0.89 for spherocylinders with aspect ratio $`L/D=5`$ jawil . The constraint on the relative orientation of the particles hugely simplifies the expression of the PDF. In fact, in order to specify the spatial configuration of a given spherocylinder relative to that fixed at the origin one just needs two parameters: the distance $`r`$ and the angle $`\vartheta `$ formed by $`𝐫`$ and by the nematic director. The other polar angle $`\varphi `$ is actually averaged out on account of the macroscopic cylindrical symmetry of the model. We emphasize that in this case, at variance with the freely rotating model already investigated in micali , no approximation is needed in order to compute the PDF of the system. Correspondingly, the pair entropy of the fluid can be written as:

$`s_2(\rho )`$ $`=`$ $`\pi \rho {\displaystyle _0^{\mathrm{}}}r^2dr{\displaystyle _0^\pi }\mathrm{sin}\vartheta d\vartheta `$ (4)
$`\{g(r,\vartheta )\mathrm{ln}[g(r,\vartheta )]g(r,\vartheta )+1\}.`$
Using a formalism analogous to that introduced in micali , we can extract from Eq. (4) the excluded-volume contribution to the pair entropy of the fluid. This contribution arises from the space integration carried out over the regions where $`g(r,\vartheta )=0`$:

$$s_2=B_2\rho +s_2^{}.$$
(5)
In Eq. (5) $`B_2`$ is the second virial coefficient of hard parallel spherocylinders that is just four times the volume, $`V_{\mathrm{hsc}}`$, of a spherocylinder:

$$B_2=4[(\pi /4)D^2L+(\pi /6)D^3].$$
(6)
The residual contribution, $`s_2^{}`$, explicitly accounts for the decrease of the pair entropy associated with the onset of interparticle correlations at short and medium-range distances. For parallel spherocylinders, a further separation of this term into a translational and an orientational contribution as done in micali  does not make much sense.

The total excess entropy of the nematic fluid can be evaluated upon integrating the equation of state (EoS):

$$s_{\mathrm{ex}}(\rho )=s_{\mathrm{ex}}(\overline{\rho })_{\overline{\rho }}^\rho \left[\frac{\beta P}{\rho ^{}}1\right]\frac{d\rho ^{}}{\rho ^{}}.$$
(7)
In Eq. (7) $`P`$ is the pressure, $`\beta `$ is the inverse temperature in units of the Boltzmann constant, and $`\overline{\rho }`$ is the density of a suitably chosen reference state. The EoS of the fluid was obtained through a Monte Carlo sampling of the system at different pressures, while the quantity $`s_{\mathrm{ex}}(\overline{\rho })`$ was calculated using the Widom test-particle insertion method widom ; smfr  in order to calculate the excess chemical potential, $`\mu _{\mathrm{ex}}`$, at low enough densities. The excess entropy can then be obtained through the thermodynamic equation:

$$s_{\mathrm{ex}}=\beta \mu _{\mathrm{ex}}+\frac{\beta P}{\rho }1.$$
(8)
## III Simulation

We investigated the phase diagram of a system of parallel hard spherocylinders with elongation $`L/D=5`$, spanning a density range which goes from a rather dilute nematic state up to the smectic transition threshold. We carried out Monte Carlo (MC) simulations at constant pressure as is usual for systems of non-spherical hard-core particles where it may be difficult to calculate the equation of state in a constant-volume simulation through the contact values of the distribution functionssmfr .

The typical sample was composed of $`N=1500`$ particles, aligned along the $`z`$ axis and enclosed in an orthorhombic box with edges $`L_x=L_y=1/3L_z`$. In order to quantify the influence of the size of the system as well as of the shape of the simulation box, we also investigated the behavior of a system composed of $`500`$ particles aligned along the main diagonal of a cubic box and of a fluid of $`768`$ particles enclosed in a orthorhombic cell with the same relative dimensions used for the main sample.

All thermodynamic states at constant pressure were sequentially generated from a translationally disordered low-density configuration upon gradually compressing the nematic fluid. The equilibration period was typically $`10^5`$ MC cycles, a cycle consisting of an attempt to change sequentially the center-of-mass coordinates of each molecule followed by an attempt to modify the volume of the sample. Simulation data were obtained by generating chains consisting of $`5\times 10^520\times 10^5`$ MC cycles, depending on the pressure. Equilibrium averages and standard deviations were computed by dividing chains into independent blocks. During the production runs, we cumulated different histograms of the PDF. In particular, $`g(r,\vartheta )`$ was sampled at intervals $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }\vartheta `$ of $`0.05D`$ and $`1^{}`$, respectively. Different choices of $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }\vartheta `$ were investigated for $`P^{}=1.0`$ and $`P^{}=2.0`$. As for the Widom insertion method, $`100`$ trial insertions per MC cycle turned out to be sufficient to insure a stable statistics for the excess chemical potential of the fluid.

In the presentation of the results, we shall refer to the reduced density $`\rho ^{}=\rho /\rho _{\mathrm{cp}}`$, where $`\rho _{\mathrm{cp}}=2/[\sqrt{2}+(L/D)\sqrt{3}]`$ is the maximum density attained by parallel spherocylinders at close packing, and to the reduced pressure $`P^{}=\beta PV_{\mathrm{hsc}}`$.

## IV Results

The current MC results for the EoS of the model are presented in Fig. 1. The comparison between the data collected for $`1500`$ particles (in an orthorombic box) and those for $`500`$ particles (in a cubic box) shows that the EoS of the fluid is rather weakly affected by the size of the sample only at high densities. The results for $`768`$ particles in an orthorombic box were not reported in the graph since they can be hardly resolved from those pertaining to the largest size investigated. The present EoS is in good agreement with that computed by other authors slf ; koda  over the whole density range explored, apart from a modest deviation, at high densities, from the molecular dynamics results obtained by Veerman and Frenkel verfre . Both these authors as well as Stroobants and coworkers slf  noted a weak change of slope in the EoS at the transition that was estimated to occur for $`\rho ^{}=0.46`$. Other structural properties of the system were evaluated in slf  such as the distribution functions associated with parallel and perpendicular correlations, below and above the observed cusp. Given the absence of either translational order within the layers or hysteresis effects associated with a compression/expansion cycle, Stroobants and coworkers concluded that the nematic phase actually transforms into a smectic phase through a continuous phase transition. Veerman and Frenkel later estimated the transition density by studying the critical slowing down of the intermediate scattering function evaluated along the nematic director verfre . Their result agrees with that reported by Stroobants and coworkers slf .

We verified that, during the simulation, the system spontaneously transforms into a smectic phase, as can be seen through the snapshots of equilibrated configurations that are shown in Fig. 2. Such snapshots were taken for densities across the transition threshold. Traces of smectic order can be hardly detected for $`P^{}=2.00`$, a pressure corresponding to an average density equal to $`0.44`$. The ordering of the fluid is already apparent just beyond the transition threshold for $`P^{}=2.30`$ and $`\rho ^{}0.47`$. Fully developed arrangements of several well-separated smectic layers distributed along the $`z`$ direction are clearly visible for $`P^{}=2.50`$ and $`\rho ^{}0.49`$. Instead, the transversal arrangements of the spherocylinders along the $`x`$-$`y`$ plane look disordered in all states investigated (compare left and right panels in Fig. 2).

We show in Fig. 3 the excess entropy of the fluid plotted as a function of the reduced density. The data were obtained using both the thermodynamic integration of Eq. (7) as well as Widom’s ghost-particle method implemented at moderately low densities. We also show the datum recently reported by Koda and Ikeda koda  which was obtained using a multistage Widom test, based on the gradual insertion of a ghost particle with a variable shape. The modest discrepancy between this finding and the present results is likely due to the multistage method, that is known to work better at high densities. On the other hand, our results for the excess entropy are closely interpolated by the expression obtained upon integrating the five-term virial fit of the EoS reported in slf .

The PDF of the fluid, $`g(r,\vartheta )`$, was plotted in Fig. LABEL:fig:4 as a function of $`\vartheta `$ for a set of interparticle distances at increasing pressures across the transition point ($`P^{}2.2`$). We note, for separations $`r<L+D`$, the existence of a forbidden range around $`\vartheta =0`$ and $`\pi `$ that is due to the overlap of two spherocylinders. This correlation gap decreases with increasing intermolecular separations. The maximum attained by $`g(r,\vartheta )`$ corresponds to the hard-core contact between spherocylinders and its height consistently decreases with increasing interparticle distances. For $`rL+D`$, the entire angular range between $`0`$ and $`\pi `$ can be sampled by a second spherocylinder. Of course, there is a blow up of the contact value of the PDF for $`\vartheta =0,\pi `$, which progressively decreases for increasing distances. An increase of the pressure enhances the overall structure of the PDF; however, no specific signature of the nematic-smectic transition can be detected in the resulting modification of the PDF.

Upon plugging the PDF into Eq. (4), we obtained the pair entropy that was plotted together with the excess entropy in Fig. 5 as a function of the density. The RMPE exhibits a change of sign from negative to positive values for $`\rho ^{}0.453`$. This threshold practically coincides with the currently accepted nematic-smectic transition density ($`\rho ^{}=0.46`$slf ; verfre . The simulations carried out for $`500`$ particles enclosed in a cubic box show a similar behavior of the RMPE whose morphology is not substantially modified with respect to that of $`1500`$ spherocylinders in an orthorombic box. This finding is noteworthy in that it demonstrates the high sensitivity of the RMPE to the structural changes occurring in the fluid, irrespectively of the difficulty to accommodate a well resolved smectic layering in a relatively small cubic box.

Figure 6 shows the quantity that, upon integration, yields the pair entropy in Eq. (4). For very short interparticle distances the $`\vartheta `$ range which can be sampled by two neighboring particles is rather limited. As a result, angular correlations are strong and give rise to the deep well in the integrand function. A steady modulated increase follows up to $`r=L+D`$ where a negative jump witnesses the onset of new strong correlations between spherocylinders lying on top or below the central reference particle. This effect is the integrated counterpart of the behavior that was already discussed for $`g(r,\vartheta )`$ in Fig. LABEL:fig:4.

The pair entropy was finally resolved into an excluded-volume and a correlation term in Fig. 7 (see Eq. 5). It is clear that the second-order virial term cannot account by itself for the crossover between $`s_{\mathrm{ex}}`$ and $`s_2`$ that is indicative of the smectic ordering of the nematic fluid.

As for the numerical reliability of the current results, we note that the statistics cumulated on $`g(r,\vartheta )`$ was such that a smooth integration over the sampled points yielded stable values for $`s_2`$, the dispersion being always lower than a percent of the average value. However, we observed a moderate sensitivity of the angle-dependent quantities in Eq. (4) on the resolution of the angular width $`\mathrm{\Delta }\vartheta `$. In order to gain a better insight into the numerical accuracy of the calculations, we performed a series of test runs for two different pressures ($`P^{}=1.0`$ and $`2.0`$) with several $`r`$ and $`\vartheta `$ grid meshes. The tests were reported in Fig. 8 where the higher sensitivity of $`s_2`$ on the angular grid size — as compared with the radial one — is quite manifest. We estimated, in the limit of $`\mathrm{\Delta }r,\mathrm{\Delta }\vartheta 0`$, an asymptotic shift of $`s_2`$ toward lower values not larger than $`2\%`$. As a result, we expect a comparable shift of the transition threshold signalled by the RMPE to lower densities.

## V Concluding remarks

In this paper we have analyzed the residual multiparticle entropy (RMPE) of parallel hard spherocylinders with aspect ratio $`L/D=5`$ across the nematic-smectic transition. The correspondence between the intrinsic ordering threshold detected through the vanishing of the RMPE ($`\rho ^{}0.453`$) and the independently ascertained phase-transition density ($`\rho ^{}=0.46`$) is quantitative. It also turned out that the indication of the zero-RMPE criterion is not significantly affected by the size and shape of the simulation box, even when the smectic layering of the fluid cannot be easily accommodated in the sample as is the case of a small cubic box. This finding further corroborates the sensitivity of the RMPE to the ordering of the fluid on a local scale.