text:
# Third-generation muffin-tin orbitals

## 1 Introduction

Muffin-tin orbitals (MTOs) have been used for a long time in ab initio, e.g. density-functional (DF), calculations of the electronic structure of condensed matter. Over the years, several MTO-based methods have been developed. The ultimate aim is to find a generally applicable electronic-structure method which is intelligible, fast, and accurate.

In order to be intelligible, an electronic-structure method must employ a minimal and flexible basis of short-ranged orbitals. As an example, the method should be able to describe the valence band and the lower part of the conduction band in $`sp^3`$-bonded materials using merely four short-ranged $`s`$\- and $`p`$-orbitals per atom and, for insulating phases, using merely occupied orbitals such as bond orbitals. Another example is materials with strong electronic correlations. For such materials, one must first construct a small, but realistic Hilbert space of many-electron wave functions, and this requires an accurate and flexible single-particle basis of atom-centered short-ranged orbitals. A small basis of short-ranged orbitals is a prerequisite for a method to be intelligible and fast, but it may be a hindrance for its accuracy, because the orbitals of a smaller basis tend to be more complicated than those of a larger basis.

Most other density-functional methods, such as plane-wave pseudopotential, LAPW, PAW, and LCAO methods, aim at simulation, and are therefore primarily accurate and robust. But they are neither fast nor intelligible in the above-mentioned sense, because they employ basis sets with of order hundred functions per atom. With such methods, understanding can therefore only be attempted after the calculation, by means of projections onto e.g. Wannier functions in case of insulators, charge densities, electron-localization functions (ELFs), partial waves, a.s.o..

The so-called 3rd-generation MTO method (Andersen et al 1994, Andersen et al 1998, Andersen et al 2000, Tank and Arcangeli 2000, Andersen and Saha-Dasgupta 2000) should come close to what we have been aiming for. In the present paper we shall explain what 3rd-generation MTOs are and what they achieve. Emphasis will be on the so-called downfolding and energy-mesh features which enable MTO bases to be small, flexible, and accurate, as we shall demonstrate by exposing them to the above-mentioned $`sp^3`$-test. From the result, the idea emerges, that for band insulators, an MTO basis can be designed a priori to span the Hilbert space of the occupied states only. That is, there is one, and only one, such MTO per electron. To get this count right, one may associate each orbital with a nominal electron (or pair), and leave it to the method to shape the orbitals in such a way that the basis set becomes complete for the occupied states. This can be done because MTOs are selective in energy, in the sense that the MTOs of order N (NMTOs) are shaped in such a way that the NMTO basis set solves Schrödinger’s equation exactly for N+1 single-particle energies, which in the present case must be chosen in such a way that they span the valence band. This ability to generate Wannier functions directly in real space, should be useful for ab initio molecular-dynamics simulations. NMTOs may also prove useful for designing many-electron wave functions, which describe correlated electron systems in a realistic way. The description of spin and orbital ordering is a trivial example. Also for the conduction bands of metals, Wannier-like, low-energy MTOs can be designed a priori. This has been demonstrated in several cases (Müller et al 1998, Sarma and Saha-Dasgupta 2000, Korotin 2000, Valenti et al 2001, Dasgupta et al 2002), most recently for the hole-doped cuprate high-temperature superconductors, where the material-dependent trend of the hopping integrals and their correlation with the maximum $`T_c`$ was discovered (Pavarini et al 2001, Dasgupta et al).

For a description of how we expand the charge density locally in such a way that Poisson’s equation can be solved and the total energy and forces can be evaluated fast and accurately, we refer to previous (Andersen et al 2000, Tank and Arcangeli 2000) and coming publications (Arcangeli and Andersen, Savrasov and Andersen). This part of the 3rd-generation computer code is still under construction.

## 2 Screened spherical waves, kinked partial waves, and muffin-tin orbitals

The 3rd-generation MTO formalism is the multiple-scattering – or KKR (Korringa 1947, Kohn and Rostoker 1954)– formalism for finding the solutions, $`\mathrm{\Psi }_i\left(𝐫\right),`$ of Schrödinger’s equation for an electron in a muffin-tin potential, $`V\left(𝐫\right)=_Rv_R\left(r_R\right),`$ with the following three extensions:

1.: The KKR formalism is proved to hold, not only for superpositions of spherically-symmetric, non-overlapping potential wells, $`v_R\left(r_R\right),`$ but also to leading order in the potential-overlap (Andersen et al 1992). Here, and in the following, $`r_R\left|𝐫𝐑\right|,`$ and $`𝐑`$ are the sites which we label by $`R.`$ The potential, $`v_R\left(r\right),`$ is taken to vanish outside a radius, $`s_R,`$ which should not exceed 1.6 times the radius of touching spheres. i.e.: $`s_R+s_R^{}1.6\left|𝐑𝐑^{}\right|`$ for any pair of sites, $`R`$ and $`R^{}.`$

2.: Exact screening transformations of the spherical waves, $`n_l\left(\kappa r_R\right)Y_L\left(\widehat{r}_R\right),`$ are introduced in order to reduce the spatial range and the energy dependence $`\left(\kappa ^2\epsilon \right)`$ of the wave-equation solutions, $`\psi _{RL}(\epsilon ,𝐫)`$ (Andersen and Jepsen 1984, Andersen et al 1992, Zeller et al 1995). Here, and in the following, $`Llm`$ labels the spherical- (or cubic-) harmonic’s character.

3.: Energy-independent MTO basis sets are derived which span the solutions $`\mathrm{\Psi }_i\left(𝐫\right)`$ with energies $`\epsilon _i`$ of Schrödinger’s equation to within errors proportional to $`(\epsilon _iϵ_0)(\epsilon _iϵ_1)..(\epsilon _iϵ_N),`$ where $`ϵ_0,ϵ_1,\mathrm{},ϵ_N`$ is a chosen energy mesh with N+1 points (Andersen et al 2000, Andersen and Saha-Dasgupta 2000). Such an energy-independent set of Nth-order MTOs is called an NMTO set. By virtue of the variational principle, the errors of the energies $`\epsilon _i`$ will be proportional to $`(\epsilon _iϵ_0)^2(\epsilon _iϵ_1)^2..(\epsilon _iϵ_N)^2.`$

At the top of figure 1 we show the LDA energy bands $`\epsilon _i\left(𝐤\right)`$ of Si in the diamond structure, calculated with the basis set of Si-centered $`s`$-, $`p`$-, and $`d`$-MTOs, i.e. with 9 orbitals/atom, for the 3-point energy mesh $`ϵ_0,ϵ_1,ϵ_2`$ indicated on the right-hand side. These bands have meV-accuracy for the MT-potential, which in the present case was the standard all-electron DF-LDA atomic-spheres potential. Since three energy points were used, the MTOs are of order $`N`$=2, that is, they are quadratic MTOs, so-called QMTOs.

The QMTO, $`\chi _{p111}^{\left(2\right)}\left(𝐫\right),`$ pointing along  from one Si to its nearest neighbor, is shown in the (21̄1̄)-plane by the first contour plot. This orbital is localized and smooth, with a few ”orthogonality wiggles” at the nearest neighbor. The remaining three contour plots show major constituents of this $`p_{111}`$-QMTO: The $`p_{111}`$-kinked partial wave (KPW), $`\varphi _{p111}(ϵ,𝐫),`$ at the central site and at the three energies $`ϵ_0,ϵ_1,`$ and $`ϵ_2.`$ In figure 2 we show this $`p_{111}`$-QMTO, together with the $`p_{111}`$-KPW at the three energies, along the line connecting the two nearest neighbors and proceeding into the back-bond.

In general, the members (labelled by $`R^{}L^{})`$ of the NMTO basis set for the energy mesh $`ϵ_0,\mathrm{},ϵ_N`$ are superpositions,

$$\chi _{R^{}L^{}}^{\left(N\right)}\left(𝐫\right)=\underset{n=0}{\overset{N}{}}\underset{RLA}{}\varphi _{RL}(ϵ_n,𝐫)L_{nRL,R^{}L^{}}^{\left(N\right)},$$
(1)
of the kinked partial waves, $`\varphi _{RL}(\epsilon ,𝐫),`$ at the $`N+1`$ points (labelled by $`n)`$ of the energy mesh. In the present case, the $`L`$-summation is over the nine $`s`$-, $`p`$-, and $`d`$-KPWs, and the $`R`$-summation is over all Si sites. Due to the localized nature of the KPWs illustrated in the figures, the latter summation is limited to the neighbors. The $`RL`$-values for which we have MTOs in the basis set, we label active $`\left(A\right),`$ or low. Expression (1) is the energy-quantized form of Lagrange interpolation,

$$\chi ^{\left(N\right)}\left(\epsilon \right)\underset{n=0}{\overset{N}{}}\varphi \left(ϵ_n\right)l_n^{\left(N\right)}\left(\epsilon \right),l_n^{\left(N\right)}\left(\epsilon \right)\underset{m=0,n}{\overset{N}{}}\frac{\epsilon ϵ_m}{ϵ_nϵ_m},$$

of a function of energy, $`\varphi \left(\epsilon \right),`$ by an $`N`$th-degree polynomial, $`\chi ^{\left(N\right)}\left(\epsilon \right):`$ The $`N`$th-degree polynomial, $`l_n^{\left(N\right)}\left(\epsilon \right),`$ is substituted by a matrix with elements, $`L_{nRL,R^{}L^{}}^{\left(N\right)},`$ the function of energy, $`\varphi \left(\epsilon \right),`$ by a Hilbert space with axes, $`\varphi _{RL}(\epsilon ,𝐫),`$ and the interpolating polynomial, $`\chi ^{\left(N\right)}\left(\epsilon \right),`$ by a Hilbert space with axes, $`\chi _{R^{}L^{}}^{\left(N\right)}\left(𝐫\right).`$

As illustrated in figures 1 and 2, a kinked partial wave is basically a partial wave with a tail joined continuously to it with a kink at a central, so-called hard sphere of radius $`a_R`$. This kink is seen most clearly for the lowest energy, $`ϵ_0.`$ As usual, the partial wave is $`\phi _{Rl}(\epsilon ,r_R)Y_L\left(\widehat{𝐫}_R\right),`$ where the function of energy is the regular solution of the radial Schrödinger equation,

$$\left[r\phi _{Rl}(\epsilon ,r)\right]^{\prime \prime }=\left[\epsilon v_R\left(r\right)l\left(l+1\right)/r^2\right]r\phi _{Rl}(\epsilon ,r),$$
(2)
for the potential-well $`v_R\left(r\right).`$ The tail of the kinked partial wave is a so-called screened spherical wave, $`\psi _{RL}(\epsilon ,𝐫),`$ which is essentially the solution with energy $`\epsilon `$ of the wave equation in the interstitial between the hard spheres, $`\mathrm{\Delta }\psi (\epsilon ,𝐫)=\epsilon \psi (\epsilon ,𝐫),`$ with the boundary condition that, independent of the energy, $`\psi _{RL}(\epsilon ,𝐫)`$ go to $`Y_L\left(\widehat{𝐫}_R\right)`$ at the central hard sphere, and to zero (with a kink) at all other hard spheres. It is this latter confinement, easily recognized in the plots, particularly at the highest energy $`ϵ_2,`$ which makes the screened spherical waves, the KPWs, and the MTOs localized when the energy is not too high. At the same time, it makes the KPW have pure $`L`$-character merely at its central sphere, because outside, it is influenced by the hard spheres centered at the neighbors. The default value of the hard-sphere radii, $`a_R,`$ is 90 per cent of the appropriate covalent, atomic, or ionic radius. The kinked partial wave thus has a kink, not only at its own, but also at the neighboring hard spheres, inside which it essentially vanishes. ’Essentially’ because the above-mentioned boundary condition only applies to the active components of the spherical-harmonics expansions of the screened spherical wave on the hard spheres. For the remaining components, in the present case the Si $`f`$\- and higher components, as well as all components on empty (E) spheres, the screened spherical wave equals the corresponding partial-wave solution of Schrödinger’s equation throughout the MT-sphere. The small bump seen in figure 1 in the lowest KPW contour along the -direction is mainly caused by the $`f`$-character on the nearest neighbor, and so is the finite amplitude seen in figure 2 inside the nearest hard sphere.

## 3 Computational steps

The radial Schrödinger (Dirac) equations (2) are integrated numerically from $`r=0`$ to $`s_R.`$ This yields the radial functions, $`\phi _{Rl}(\epsilon ,r),`$ and their phase shifts, $`\eta _{Rl}\left(\epsilon \right),`$ each of which are obtained by matching the logarithmic derivative of $`\phi _{Rl}(\epsilon ,r)`$ at $`r=s_R`$ to that of

$$\phi _{Rl}^o(\epsilon ,r)j_l\left(\kappa r\right)\mathrm{tan}\eta _{Rl}\left(\epsilon \right)n_l\left(\kappa r\right).$$
(3)
The radial integration must be performed for each potential well and for each $`l,`$ increasing until all further phase shifts vanish due to dominance of the centrifugal term in (2).

The screened spherical waves are specified by a Hermitian structure matrix, whose element $`B_{R^{}L^{},RL}\left(\epsilon \right)`$ is essentially the radial logarithmic derivative of the $`L^{}`$-component in the spherical-harmonics expansion at the hard sphere at site $`R^{}`$ of the screened spherical wave $`\psi _{RL}(\epsilon ,𝐫).`$ What is known analytically, is the element

$$B_{R^{}L^{},RL}^0\left(\epsilon \right)\underset{l\mathrm{"}}{}4\pi i^{l+l^{}l^{\prime \prime }}C_{LL^{}l^{\prime \prime }}\kappa n_{l^{\prime \prime }}\left(\kappa \left|𝐑𝐑^{}\right|\right)Y_{L^{\prime \prime }}^{}\left(\widehat{𝐑𝐑^{}}\right),$$

of the bare KKR structure matrix, which specifies how the spherical wave, $`n_l\left(\kappa r_R\right)Y_L\left(\widehat{r}_R\right),`$ at site $`R`$ is expanded around another site, $`R^{},`$ in regular spherical waves, $`j_l^{}\left(\kappa r_R^{}\right)Y_L^{}\left(\widehat{r}_R^{}\right).`$ Here $`\kappa ^2\epsilon ,`$ and the on-site terms of the bare structure matrix are defined to vanish. Screening of the structure matrix,

$$\left[B\left(\epsilon \right)^1\right]_{RL,R^{}L^{}}\left[B^0\left(\epsilon \right)^1\right]_{RL,R^{}L^{}}+\kappa ^1\mathrm{tan}\alpha _{RL}\left(\epsilon \right)\delta _{RR^{}}\delta _{LL^{}},$$

requires inversion of the matrix $`B_{RL,R^{}L^{}}^0\left(\epsilon \right)+\kappa \mathrm{cot}\alpha _{RL}\left(\epsilon \right)\delta _{RR^{}}\delta _{LL^{}}.`$ This can be done by fixing $`R`$ and limiting $`R^{}`$ to the 10-50 nearest sites. $`\alpha _{RL}\left(\epsilon \right)`$ are the hard-sphere phase shifts for the active channels,

$$\mathrm{tan}\alpha _{RL}\left(\epsilon \right)j_l\left(\kappa a_R\right)/n_l\left(\kappa a_R\right),$$

and, for the remaining channels, $`\alpha _{RL}\left(\epsilon \right)`$ are the proper phase shifts, $`\eta _{RL}\left(\epsilon \right).`$ The latter channels, which will not have KPWs and MTOs associated with them, are said to be downfolded. With appropriate division into active and downfolded channels, the screened structure matrix will have short spatial range and no poles in the energy-range of the occupied states.

A kinked partial wave is defined as:

$$\varphi _{RL}(\epsilon ,𝐫)=\left[\phi _{Rl}(\epsilon ,r_R)\phi _{Rl}^o(\epsilon ,r_R)\right]Y_L\left(\widehat{r}_R\right)+\psi _{RL}(\epsilon ,𝐫),$$
(4)
where $`\phi (\epsilon ,r)`$ is the radial solution for the central well from 0 to $`s,`$ and $`\phi ^o(\epsilon ,r)`$ is the phase-shifted wave (3) proceeding smoothly inwards from $`s`$ to the central $`a`$-sphere, where it is matched with a kink to the screened spherical wave $`\psi (\epsilon ,𝐫).`$ The kinks of the KPW set are then given by the kink matrix,

$$K_{RL,R^{}L^{}}\left(\epsilon \right)\frac{B_{RL,R^{}L^{}}\left(\epsilon \right)+\kappa \mathrm{cot}\eta _{RL}^\alpha \left(\epsilon \right)\delta _{RR^{}}\delta _{LL^{}}}{\epsilon n_l\left(\kappa a_R\right)n_l^{}\left(\kappa a_R^{}\right)},$$
(5)
where $`\eta _{RL}^\alpha \left(\epsilon \right)`$ is the phase shift with respect to the hard-sphere medium,

$$\mathrm{tan}\eta _{RL}^\alpha \left(\epsilon \right)\mathrm{tan}\eta _{RL}\left(\epsilon \right)\mathrm{tan}\alpha _{RL}\left(\epsilon \right).$$

The rows and columns of the kink matrix run merely over active channels. In the fomalism above, we have for simplicity used the notation of scattering theory, which is analytical for $`\epsilon >0,`$ and for $`\epsilon <0.`$ Screened scattering theory with the normalization (5) is however analytical in a region of interest around $`\epsilon =0.`$

Finally, the Lagrange matrix which gives the MTO set (1) in terms of the KPW set (4), is given solely in terms of the values of the Green matrix, $`G\left(\epsilon \right)K\left(\epsilon \right)^1`$, on the energy mesh $`\epsilon =ϵ_0,ϵ_1,\mathrm{},ϵ_N.`$ The Hamiltonian and overlap matrices in the MTO representation are given in terms of the same values of the Green matrix, plus the values of its first energy derivative, $`\dot{G}\left(\epsilon \right).`$

## 4 Further downfoldings

In the valence and lowest conduction bands of Si, there are only $`s`$\- and $`p`$-, but no $`d`$-electrons. To describe these bands, we should therefore be able to use a basis with only Si $`s`$\- and $`p`$-MTOs, that is, with only 4 orbitals per atom. We thus let the Si $`s`$\- and $`p`$-partial waves remain active, while the Si $`d`$-waves are now included among the passive ones, i.e. those ’folded down’ into the tails of the screened-spherical waves in (4). The results for the bands and the $`p_{111}`$-QMTO are shown in figures 3 and 4. These bands are indistinguishable from those obtained with the Si $`spd`$-set, on the scale of the figure, although between the energies of the mesh, the former bands do lie slightly above the latter. However, by making the mesh denser (increasing $`N),`$ the accuracy can be increased arbitrarily. The KPW of the $`sp`$-set is seen to have $`d`$-character on the nearest Si neighbor, and the QMTO and the KPW, particularly the one at the highest energy, are seen to be somewhat less localized than those for the $`spd`$-set.

It is even possible to construct an arbitrarily accurate MTO-basis which spans merely the occupied orbitals, that is, which spans the valence band with a basis of merely 2 orbitals per atom. For tetrahedrally coordinated covalent semiconductors like Si, it is customary to take the valence-band orbitals as the bond-orbitals, which are the bonding linear combinations of directed $`sp^3`$-hybrids of orthonormal orbitals. It is, however, far simpler and more general, e.g. not limited to elemental semiconductors and tetrahedral structures, to take the valence-band orbitals as the $`s`$\- and $`p`$-MTOs on every second Si atom, all partial waves on the nearest neighbors being downfolded. This corresponds to a Si4+Si4- ionic picture. This QMTO-set turns out to describe merely the valence band, and to do so surprisingly well considering the fact that the two silicons are treated differently, so that the degeneracy along the XW-line is, in fact, slightly broken. The error between the energy points is proportional to $`\left[\epsilon _i\left(𝐤\right)ϵ_0\right]\left[\epsilon _i\left(𝐤\right)ϵ_1\right]\left[\epsilon _i\left(𝐤\right)ϵ_2\right],`$ exactly as for the basis with 4 orbitals per atom shown in figure 3, because we use QMTOs in both cases, but the prefactors are larger for the smaller basis: As the number of active channels decreases, the KPWs attain longer range and stronger energy dependence. However, by making the energy mesh finer, the errors of the MTO set can be made arbitrarily small. In the bottom line of figure 5, we show the result of such a valence-band-only calculation with $`N=3`$ for Ge, together with the $`p_{111}`$ cubic MTO (CMTO) centered on the Ge atom to the right. The accuracy of the valence band is superb, and the MTO is seen to spill over onto the nearest-neighbor atom(s) which were chosen not to have orbitals associated with them.

Since this basis is complete for the occupied states, we may compute the density-functional ground-state properties in real space by taking traces, provided that we first Löwdin orthonormalize the basis in real space. The sum of the one-electron energies is then computed as the trace of the Hamiltonian, i.e. as the sum of the energies of the orthonormal orbitals, and the charge density is computed as the sum of the squares of these orbitals. This is a method where the amount of computation increases merely linearly with the size of the system, a so-called order-$`𝒩`$ method. Here, $`𝒩`$ refers to the number of atoms in the system and not to the order $`N`$ of the MTOs. This NMTO method, which generates the complete basis for the occupied states a priori, should be superior to current order-$`𝒩`$ methods, which either use inaccurate empirical tight-binding models or project onto the occupied states during the course of a large-basis-set calculation.

In order to demonstrate in further detail that our method works, we consider Si in the diamond structure for which the valence-band Wannier functions can be taken as bond orbitals. First, we orthonormalize our symmetry-breaking Si4+ Si $`s,p_x,p_y,p_z`$ QMTO set. The resulting Si $`s`$ and Si $`p_x`$ orbitals are shown in the (110)-plane in figures 6 and 7. These orthogonalized orbitals are seen to remain fairly localized. Then, we transform to the four congruent $`sp^3`$-hybrids centered nominally on every second Si atom. As figure 8 shows, such an $`sp^3`$-hybrid is, in fact, the bond orbital. Hence, folding all partial waves of the atoms chosen not to carry orbitals, into the tail of the $`sp^3`$-directed orbital on one of the other Si atoms, has made that orbital look like a bond orbital. The reason why the figure does not show exact symmetry between the two sites is caused by our use of an energy mesh with only 3 points in the valence band. Making the energy mesh finer will generate the exact symmetry.

Now, Sn is metal because the bonding and antibonding bands overlap. This should, however, not prevent our method from working for the occupied states only, because NMTOs are energy selective. Since the orbitals will be shaped in such a way that the basis set solves Schrödinger’s equation exactly for the energies on the mesh, we may merely have to choose several energy points in the region of band overlap below the Fermi level –and, of course, no energy points above. Remember that our ionic prescription does not make use of the fact that the Wannier functions for the valence and conductions bands are respectively bonding and antibonding.

Since the ionic Si Si$`\left(sp\right)`$ set gives the occupied states in diamond-structured Si with arbitrary accuracy, the same procedure with the $`sp`$-orbitals placed exclusively on the anion, and the $`d`$-orbitals on the cation, will of course work for any IV-IV, III-V, II-VI, and I-VII semiconductor and insulator. CuBr, for instance, would be thought of as an ionic compound Cu+Br- with the closed-shell configuration Cu$`3d^{10}`$ Br$`4s^24p^6,`$ and the basis should therefore have the $`d`$-MTOs on the Cu atoms and the $`s`$ and $`p`$-MTOs on the Br atoms. This is illustrated in the upper line of figure 5. Being a single-site prescription, this works for CuBr in any structure. What we have seen is thus, that the MTO basis can be designed a priori to span the Hilbert space of the occupied states only. That is, there is one, and only one, such MTO per electron. To specify such a set, one would, in order to get the electron-count right, put the orbitals where the electrons are though to be, and leave it to the method to shape the tails of these orbitals in such a way that the basis solves Schrödinger’s equation exactly for occupied states of the given static mean-field (e.g. LDA) potential. Such ionic MTO basis sets, which ’automatically’ span the occupied –and no further– states of any band insulator, could make density-functional molecular-dynamics calculations highly efficient for such systems.

How could one imagine to treat a chemical reaction like: 2$`H_2O2H_2+O_2`$ with MTO bases of occupied states only? For water, it is natural to use the ionic description H$`{}_{}{}^{+}{}_{2}{}^{}`$O$`^{},`$ according to which the orbital configuration is O $`p^6,`$ i.e. one would put the $`p`$ orbitals on oxygen and fold down all partial waves centered on the hydrogens. In principle, one might stick to this configuration throughout the reaction, because it keeps the electron-count right. However, the oxygen-centered orbitals would eventually look strange and have long range, because they would have to separate off pieces of wave functions sitting on the hydrogens. Such oxygen orbitals might be more time-consuming to generate. At some stage in the reaction, it might therefore be appropriate to switch to configurations such as H$`s`$ H$`s,`$ or H H$`s^2`$ for the hydrogen molecule; this is analogous to our treatment of the occupied states in tetrahedrally coordinated Si. For O2 with the open-shell molecular configuration $`pp\sigma ^2`$ $`pp\pi ^4`$ $`pp\pi ,`$ we might use an ’ionic’ configuration like: O $`(zxy)`$ O $`\left(zxy\right)`$ with $`z`$ referring to the local $`z`$-direction of the molecule.

This example immediately leads to a treatment of open-shell systems by means of spin and possibly orbital polarizations. We have seen that the 3rd-generation MTO method offers the possibility of designing single-electron bases of atom-centered localized orbitals, which span the wave functions in a given energy region of a given mean-field potential. These orbitals can even be symmetry-breaking, as in the case of diamond, Si, Ge, and Sn, without the generating mean field having to be so. These orbitals thus seem to have great potential in the design of many-electron wave functions which describe correlated electron systems in a realistic way.

## 5 Conclusion

We have solved the long-standing problem of deriving energy-independent, short-ranged orbitals from scattering theory (Hubbard 1967). The present formalism contains exactly the right ’physics and chemistry,’ we feel. This should give the computational method great speed and accuracy, and make it a vehicle for discovery and understanding.

## 6 References

Andersen O K and Jepsen O 1984 Phys. Rev. Lett. 53 2571

Andersen O K, Postnikov A V, and Savrasov S Yu 1992 Applications of Multiple Scattering Theory to Materials Science eds. W H Butler, P H Dederichs, A Gonis and R L Weaver, Mat. Res. Soc. Symp. Proc. Vol. 253 (Pittsburgh: Materials Research Society) 37-70

Andersen O K, Jepsen O, and Krier G 1994 Lectures on Methods of Electronic Structure Calculations eds V. Kumar, O.K. Andersen, and A. Mookerjee (Singapore: World Scientific Publishing Co) 63-124

Andersen O K, Arcangeli C, Tank R W, Saha-Dasgupta T, Krier G, Jepsen O, and Dasgupta I 1998 Tight-Binding Approach to Computational Materials Science eds L Colombo, A Gonis, and P Turchi, Mat. Res. Soc. Symp. Proc. Vol. 491 (Pittsburgh: Materials Research Society) 3-34

Andersen O K, Saha-Dasgupta T, Tank R W, Arcangeli C, Jepsen O, and Krier G 2000 Electronic Structure and Physical Properties of Solids. The Uses of the LMTO Method ed H Dreysse (New York: Springer Lecture Notes in Physics) 3-84

Andersen O K and Saha-Dasgupta T 2000 Phys. Rev. B 62 R16 219

Arcangeli C and Andersen O K to be published

Dasgupta I et al 2002 Bull Mater Sci this issue

Dasgupta I, Saha-Dasgupta T, Andersen O K, Pavarini E, and Jepsen O to be published

Hubbard J 1967 Proc. Phys. Soc. London 92 921, and references therein

Kohn W and Rostoker J 1954 Phys. Rev. 94 111

Korotin M A, Anisimov V I, Saha-Dasgupta T, Dasgupta I 2000 Journal of Physics: Condensed Matter 12 113-124

Korringa J 1947 Physica 13 392

Müller T F A, Anisimov V, Rice T M, Dasgupta I, Saha-Dasgupta T 1998 Phys. Rev. B 57 R12655

Pavarini E, Dasgupta I, Saha-Dasgupta T, Jepsen O, and Andersen O K 2001 Phys. Rev. Lett. 87 047003

Sarma D D, Mahadevan P, Saha-Dasgupta T, Ray S, and Kumar A 2000 Phys. Rev. Lett. 85 2549

Savrasov D and Andersen O K to be published

Tank R W and Arcangeli C 2000 phys. stat. sol. (b) 217 89-130

Valenti R, Saha-Dasgupta T, Alvarez J V, Pozgajcic K, and Gros C 2001 Phys. Rev. Lett. 86, 5381 (2001)

Zeller R, Dederichs P H, Ujfalussy B, Szunyogh L, and Weinberger P 1995 Phys. Rev. B 52 8807


format: text

text:
# Theory and simulation of the nematic anchoring coefficient

## I Introduction

The anchoring phenomenon is the tendency of a liquid crystal to orient in a particular direction when in contact with the container walls. The equilibrium director orientation set by the interaction of the liquid crystal with the aligning surface is called an *easy orientation axis*. The simplest, strong anchoring, assumption is that the director has a fixed orientation at the boundaries along the easy orientation axis. However, it has been discovered that the coupling of the director with the orienting surface can be rather weak. This results in deviation of the surface director from the easy axis in response to small perturbations.

On a phenomenological level, weak anchoring can be described by adding an appropriate *surface potential* to the free energy of the system. Then minimization of the free energy functional gives both the equations for the director in the bulk and the appropriate boundary conditions de Gennes and Prost (1995). The simplest form of the surface potential has been proposed by Rapini and Papoular Rapini and Papoular (1969)

$$f_\mathrm{s}=\frac{1}{2}W(𝒏𝒆)^2,$$
(1)
where the parameter, $`W`$, is termed an anchoring energy.

Since then, numerous experimental methods have been used to measure the surface anchoring coefficient $`W`$ Subacius et al. (1995); Yokoyama and van Sprang (1985); Nastishin et al. (1999); its value has turned out to be extremely important for liquid crystal devices, i.e. displays, optical switches. However, comparatively little work on systematic experimental investigation of the anchoring phenomenon has been presented up till now. From the available experimental data, one can say that the extrapolation length $`\lambda =K_{33}/W`$ is inversely proportional to the squared value of the bulk order parameter $`\lambda Q^2`$. Taking into account that the elastic constant $`K_{33}`$ is typically proportional to $`Q^2`$, this gives for the anchoring parameter $`WQ^4`$ Mertelj and Čopič (2000).

There have been several attempts to estimate the anchoring coefficient theoretically Tjipto-Margo and Sullivan (1988) and by combining molecular simulation with a local density functional approach Stelzer et al. (1997); Teixeira et al. (2001). The main difficulty here is that one needs to know the direct pair correlation function of the nematic state, which is usually unknown, and must hence be estimated with some uncontrolled approximations or assumptions (see, however, Ref. Phuong et al. (2001)). Moreover, it is often assumed that the director in the interface region is fixed Stelzer et al. (1997). However, it has already been noticed that these approximations may give incorrect results, by at least an order of magnitude, for example in the calculation of bulk elastic constants Phuong et al. (2001).

Similar studies have also been done using lattice models. The existence of subsurface deformations and effective intrinsic anchoring arising from the incomplete molecular interaction close to the surface has been studied using a hexagonal lattice approach Skačej et al. (1997). Monte-Carlo simulations of the Lebwohl-Lasher model have shown that the extrapolation length is not in general equal to the ratio of the bulk to surface couplings Priezjev and Pelcovits (2000). However, the results were not sufficiently robust to determine the dependence of the anchoring energy on the order parameter.

The present work attempts to remedy the situation by combining the elastic description with Onsager theory and Monte Carlo simulation results. We study the dependence of the elastic ($`K_{33}`$) and surface anchoring ($`W`$) coefficients on the liquid crystal state point, which is defined by the bulk value of the order parameter $`Q`$.

The paper is organized as follows. We define the geometry and derive director profiles using elastic theory in section II. In section III we discuss the Onsager approach which allows us to calculate the single particle density of the bulk and confined systems, while section IV outlines the fluctuation approach to calculating elastic coefficients. Section V gives details of the technique used in Monte Carlo simulations and the results are summarized in section VI.

## II Elastic description

The easiest way to obtain the surface anchoring coefficient $`W`$ is to create a director deformation far from the surface. Then, measuring the response of the director near the cell surface, or fitting the director profile with a theoretically predicted profile, yields an estimate of the anchoring extrapolation length and ratios of the elastic constants.

Indeed, consider one of the possible geometries suitable for the measurement of the zenithal anchoring strength. Let the director have fixed orientation at the boundary $`z=L`$. The surface at $`z=0`$ is assumed to provide homeotropic anchoring of strength $`W`$. In the elastic description, deformations of the director field $`𝒏`$ are described by the total free energy de Gennes and Prost (1995)

$$_{\text{el}}=_Vf_\text{b}𝑑V+_Sf_\text{s}𝑑S.$$
(2)
Here $`f_\text{b}`$ is the Frank-Oseen elastic free energy density

$`f_\text{b}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{K_{11}(𝒏)^2+K_{22}(𝒏[\times 𝒏])^2+`$
$`+K_{33}[𝒏\times [\times 𝒏]]^2\},`$
where $`K_{11}`$, $`K_{22}`$, and $`K_{33}`$ are elastic constants, and the integration extends over the sample volume $`V`$. $`f_\text{s}`$ is the surface anchoring energy density, which we assume to be of the Rapini-Papoular form (1), and it is integrated over the boundary surface $`S`$.

In slab geometry, the director $`𝒏`$ is assumed to lie in the $`xz`$ plane and depends only on the $`z`$ coordinate. Then it can be parametrized as $`𝒏=(\mathrm{sin}\theta ,0,\mathrm{cos}\theta )`$ which transforms eqn (2) into a free energy per unit area

$$_{\text{el}}/S=\frac{1}{2}K_{33}_0^Lf_{13}(\theta )\left(\frac{\theta }{z}\right)^2𝑑z\frac{1}{2}W\mathrm{cos}^2\theta _0,$$
(4)
where $`f_{13}(\theta )=1\delta \mathrm{sin}^2\theta `$, $`\delta =(K_{33}K_{11})/K_{33}`$, $`\theta _0=\theta (z=0)`$.

The absence of explicit $`z`$-dependence in the free energy (4) implies the first integral

$$f_{13}(\theta )\left(\frac{\theta }{z}\right)^2=\text{const}.$$
(5)
Boundary conditions read

$`K_{33}f_{13}(\theta _0){\displaystyle \frac{\theta }{z}}|_{z=0}`$ $`=`$ $`{\displaystyle \frac{1}{2}}W\mathrm{sin}2\theta _0,`$ (6)
$`\theta (z=L)`$ $`=`$ $`\theta _L.`$
Here we assumed that the director angle at the boundary $`z=L`$ is fixed.

Integrating eqn (5), together with the boundary conditions, yields

$`E(\theta ,\delta )=E(\theta _0,\delta )+[E(\theta _L,\delta )E(\theta _0,\delta )]z/L,`$ (7)
$`[E(\theta _L,\delta )E(\theta _0,\delta )]\sqrt{f_{13}(\theta _0)}=(L/2\lambda )\mathrm{sin}2\theta _0,`$
where $`E(\theta ,\delta )=_0^\theta \sqrt{f_{13}(x)}𝑑x`$ is the incomplete elliptic integral of the second kind, and $`\lambda =K_{33}/W`$ is the anchoring extrapolation length. For small angles $`\theta _0`$ eqn (7) can be simplified and has the form

$$E(\theta ,\delta )=E(\theta _L,\delta )\frac{z+\lambda }{L+\lambda }.$$
(8)
Note, that for small angles $`\theta _L`$ and, correspondingly, for small $`\theta `$, $`E(\theta ,\delta )=\theta `$ and we have linear dependence of the director angle on the $`z`$ coordinate

$$\delta n_x(z)\theta (z)=\theta _L\frac{z+\lambda }{L+\lambda }.$$
(9)
Using eqn (8), one can fit the director profiles $`\theta (z)`$ with $`\lambda =K_{33}/W`$ and $`\delta =(K_{33}K_{11})/K_{33}`$ as adjustable parameters. To simplify the procedure, it is more appropriate to fit $`z(\theta )`$ rather than $`\theta (z)`$.

## III Onsager approach

The Helmholtz free energy in the Onsager approach is expressed in terms of the single-particle density, $`\rho (\mathrm{𝟏})`$, where $`(\mathrm{𝟏})=(𝒓_1,𝛀_1)`$ represents both position $`𝒓_1`$ and orientation $`𝛀_1`$. It has the following form Onsager (1949); Evans (1992)

$`\beta [\rho ]=`$ $`{\displaystyle \rho (\mathrm{𝟏})\left\{\mathrm{ln}\rho (\mathrm{𝟏})\mathrm{\Lambda }^31\beta \mu +\beta U(\mathrm{𝟏})\right\}d(\mathrm{𝟏})}`$ (10)
$`{\displaystyle \frac{1}{2}}{\displaystyle f(\mathrm{𝟏},\mathrm{𝟐})\rho (\mathrm{𝟏})\rho (\mathrm{𝟐})d(\mathrm{𝟏})d(\mathrm{𝟐})}.`$
Here $`\beta =1/k_\text{B}T`$, $`\mathrm{\Lambda }`$ is the de Broglie wavelength, $`\mu `$ is the chemical potential, $`U`$ is the external potential energy (including the surface potential), and $`f(\mathrm{𝟏},\mathrm{𝟐})`$ is the Mayer $`f`$-function

$$f(\mathrm{𝟏},\mathrm{𝟐})=\mathrm{exp}\left[𝒱(\mathrm{𝟏},\mathrm{𝟐})/k_\mathrm{B}T\right]1,$$
(11)
where elongated particles interact pairwise through the potential $`𝒱(\mathrm{𝟏},\mathrm{𝟐})`$.

The equilibrium single-particle density that minimizes the free energy (10) is a solution of the following Euler-Lagrange equation

$$\mathrm{ln}\rho (\mathrm{𝟏})\mathrm{\Lambda }^3\beta \mu +\beta U(\mathrm{𝟏})f(\mathrm{𝟏},\mathrm{𝟐})\rho (\mathrm{𝟐})d(\mathrm{𝟐})=0,$$
(12)
which can be obtained from the variation of the functional (10). In practice, we find it more convenient to minimize the functional (10) instead of solving the integral equation (12).

### III.1 Bulk problem

In the bulk problem, the single particle density is independent of position, $`\rho (\mathrm{𝟏})=\rho (𝛀_1)`$. Then the integrals over position may be performed directly. To perform the integration over the angles, we expand the Mayer $`f`$-function and the single-particle density in spherical harmonics

$$f(\mathrm{𝟏},\mathrm{𝟐})=\underset{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}{}f^{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}(r_{12})\mathrm{\Phi }^{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}(𝛀_1,𝛀_2,𝛀_r),$$
(13)
$$\rho (𝛀)=\underset{\mathrm{}}{\overset{\text{even}}{}}\rho _{\mathrm{}}Y_\mathrm{}0(𝛀),$$
(14)
where $`\mathrm{\Phi }^{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}(𝛀_1,𝛀_2,𝛀_r)`$ is a rotational invariant Gray and Gubbins (1984)

$`\mathrm{\Phi }^{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}=4\pi {\displaystyle \underset{m_1,m_2,m_r}{}}\left(\begin{array}{ccc}\mathrm{}_1& \mathrm{}_2& \mathrm{}_r\\ m_1& m_2& m_r\end{array}\right)\times `$ (17)
$`Y_{\mathrm{}_1m_1}(𝛀_1)Y_{\mathrm{}_2m_2}(𝛀_2)C_{\mathrm{}_rm_r}(𝛀_r).`$ (18)
Here, $`\left(\begin{array}{ccc}\mathrm{}& \mathrm{}^{}& \mathrm{}^{\prime \prime }\\ m& m^{}& m^{\prime \prime }\end{array}\right)`$ is a $`3j`$ symbol, $`𝛀_r`$ is the direction of the unit vector $`\widehat{𝒓}_{12}=𝒓_{12}/r_{12}`$ where $`𝒓_{12}=𝒓_1𝒓_2`$, $`C_\mathrm{}m`$ is a reduced spherical harmonic. In writing eqn (14), we assumed that the director $`𝒏_0`$ is pointing along the $`z`$ axis.

To minimize the grand potential it is convenient to expand the logarithm of the density in spherical harmonics

$$\mathrm{ln}\rho (𝛀)=\underset{\mathrm{}}{\overset{\text{even}}{}}c_{\mathrm{}}Y_\mathrm{}0(𝛀).$$
(19)
The grand potential then can be rewritten in the form

$`{\displaystyle \frac{\beta }{V}}=\sqrt{4\pi }(1+\beta \mu )\rho _0+{\displaystyle \underset{\mathrm{}=0}{\overset{\text{even}}{}}}c_{\mathrm{}}\rho _{\mathrm{}}+{\displaystyle \frac{2\pi V_{\mathrm{}\mathrm{}}}{\sqrt{2\mathrm{}+1}}}\rho _{\mathrm{}}\rho _{\mathrm{}},`$ (20)
where the coefficients $`\rho _{\mathrm{}}`$ may be expressed in terms of the parameters $`c_{\mathrm{}}`$, and where the pair-excluded volume is expanded in coefficients

$$V_{\mathrm{}\mathrm{}}=4\pi _0^{\mathrm{}}r^2𝑑rf^\mathrm{}\mathrm{}0(r).$$
(21)
The grand potential was numerically minimized with respect to variation of the parameters $`c_{\mathrm{}}`$, by the conjugate gradient method Press et al. (1992). The resulting single-particle density was used to calculate elastic constants, for different values of the order parameter.

### III.2 Elastic constants

To evaluate the elastic constants we used the expressions derived by Poniewierski and Stecki Poniewierski and Stecki (1979); Lipkin et al. (1985); Yokoyama (1997)

$`K_{11}=M_{xxxx}=M_{yyyy},`$ (22)
$`K_{22}=M_{xxyy}=M_{yyxx},`$
$`K_{33}=M_{zzxx}=M_{zzyy},`$
where

$`M_{\alpha \beta \gamma \delta }`$ $`=`$ $`{\displaystyle \frac{1}{2}}k_BT{\displaystyle }d𝒓_{12}d𝛀_1d𝛀_2f(𝒓_{12},𝛀_1,𝛀_2)\times `$ (23)
$`r_{12}^\alpha r_{12}^\beta \rho ^{}(\mathrm{cos}\theta _1)\rho ^{}(\mathrm{cos}\theta _2)\mathrm{\Omega }_{1\gamma }\mathrm{\Omega }_{2\delta }.`$
All integrals are evaluated in a local coordinate frame with the $`z`$ axis parallel to the director at point $`𝒓`$.

As discussed in appendix A, performing the integrations over the angles and making use of the properties of $`3j`$ symbols and spherical harmonics, one obtains a simplified expression previously given by Stelzer et al Stelzer et al. (1995):

$`\beta K_{11}`$ $`=`$ $`K(1,1,3),`$ (24)
$`\beta K_{22}`$ $`=`$ $`K(1,5,1),`$
$`\beta K_{33}`$ $`=`$ $`K(2,4,4),`$
where

$`K(a_1,a_2,a_3)`$ $`={\displaystyle \frac{4\pi ^2}{3}}{\displaystyle \underset{\mathrm{}}{\overset{\text{even}}{}}}\mathrm{}(\mathrm{}+1)\rho _{\mathrm{}}^2\left[{\displaystyle \frac{I^{\mathrm{},\mathrm{},0}}{\sqrt{2\mathrm{}+1}}}{\displaystyle \frac{6a_1+a_2\mathrm{}(\mathrm{}+1)}{5}}\sqrt{{\displaystyle \frac{(2\mathrm{}2)!}{(2\mathrm{}+3)!}}}I^{\mathrm{},\mathrm{},2}\right]+`$
$`{\displaystyle \frac{a_3}{5}}\sqrt{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{(\mathrm{}+3)!}{(\mathrm{}1)!}}\sqrt{{\displaystyle \frac{(2\mathrm{})!}{(2\mathrm{}+5)!}}}\rho _{\mathrm{}}\rho _{\mathrm{}+2}I^{\mathrm{},\mathrm{}+2,2},`$
and where $`I^{\mathrm{}_1,\mathrm{}_2,\mathrm{}}`$ are radial integrals over the expansion coefficients of the Mayer $`f`$-function.

### III.3 Slab geometry

In slab geometry, with the $`z`$ direction normal to the surfaces, the single-particle density and the external potential depend on the $`z`$ coordinate only. To perform the integrations over the angles, we expand the Mayer $`f`$-function into rotational invariants, similar to the bulk system (13). The single-particle density and its logarithm can be also expanded in spherical harmonics

$`\rho (z,𝛀)={\displaystyle \underset{\mathrm{}m}{}}\rho _\mathrm{}m(z)Y_\mathrm{}m^{}(𝛀),`$
$`\mathrm{ln}\rho (z,𝛀)={\displaystyle \underset{\mathrm{}m}{}}c_\mathrm{}m(z)Y_\mathrm{}m(𝛀).`$ (25)
The difference from the bulk case, eqns (14,19), is that the director is allowed to vary in the $`xz`$ plane, so an expansion in Legendre polynomials ($`m=0`$) is not sufficient. Conducting angular and $`x`$\- and $`y`$-integrations gives

$`{\displaystyle \frac{\beta }{S}}={\displaystyle 𝑑z\sqrt{4\pi }(1+\beta \mu )\rho _{00}(z)}+`$ (26)
$`{\displaystyle \underset{\mathrm{},m}{}}\left(c_\mathrm{}m(z)+\beta U_\mathrm{}m(z)\right)\rho _\mathrm{}m+`$
$`2\pi {\displaystyle 𝑑z_1𝑑z_2A_{\mathrm{}_1,\mathrm{}_2,m}(z_{12})\rho _{\mathrm{}_1m}(z_1)\rho _{\mathrm{}_2m}(z_2)},`$
where the pair-excluded area at given $`z`$-separation $`z_{12}`$ is expanded in coefficients

$`A_{\mathrm{}_1,\mathrm{}_2,m}(z_{12})=2\pi {\displaystyle \underset{\mathrm{}}{\overset{\text{even}}{}}}\left(\begin{array}{ccc}\mathrm{}_1& \mathrm{}_2& \mathrm{}\\ m& \overline{m}& 0\end{array}\right)\times `$ (29)
$`{\displaystyle _0^{\mathrm{}}}s𝑑sf_{\mathrm{}_1,\mathrm{}_2,\mathrm{}}(r_{12})P_{\mathrm{}}(\mathrm{cos}\theta _r).`$ (30)
Here $`s^2=(x_1x_2)^2+(y_1y_2)^2`$, $`\mathrm{tan}\theta _r=s/(z_1z_2)`$, $`\overline{m}=m`$. Integrating, we took into account that the $`3j`$ symbol vanishes unless $`m_1+m_2+m_r=0`$.

To obtain equilibrium single-particle density profiles, the grand potential then was tabulated on a regular grid of points $`z_i`$ and numerically minimized with respect to variation of the parameters $`c_\mathrm{}m(z_i)`$, by the conjugate gradient method.

### III.4 Anchoring energy

To obtain the microscopic expression for the extrapolation length $`\lambda `$, we start from the equation for the single-particle density (12). Assume that the solution for a ground state, i.e. for homeotropically aligned liquid crystal in slab geometry, is given by the single-particle density $`\rho _0`$. Consider a small perturbation around the ground state, $`\rho =\rho _0+\delta \rho `$. To first order in $`\delta \rho `$, eqn (12) can be written as

$$\frac{\delta \rho (\mathrm{𝟏})}{\rho _0(\mathrm{𝟏})}=f(\mathrm{𝟏},\mathrm{𝟐})\delta \rho (\mathrm{𝟐})d(\mathrm{𝟐}).$$
(31)
In the case of slowly-varying director fields we assume that the free energy functional is locally in equilibrium. This is equivalent to the mathematical simplification Stelzer et al. (1997); Somoza and Tarazona (1991)

$$\rho (𝒓,𝛀)=\rho _0(𝒓,𝒏(𝒓)𝛀).$$
(32)
Then the density perturbation can be written in terms of the perturbation of the director

$$\delta \rho (𝒓,𝛀)=\rho _0^{}(𝒓,𝒏_0𝛀)\delta 𝒏(𝒓)𝛀,$$
(33)
where the prime denotes a partial derivative with respect to $`(𝒏_0𝛀)`$.

In slab geometry, with the $`z`$ axis normal to the surfaces, the single particle density $`\rho _0`$ depends on the $`z`$ coordinate only and can be expanded in spherical harmonics similar to eqn (25). Note, that $`\rho _0`$ does not depend on $`\varphi `$ which implies $`m=0`$ in expressions (25). We also assume that the director is parallel to the $`xz`$ plane, $`\delta 𝒏=(\delta n_x,0,0)`$. Conducting angular integrations in eqn (31) and making use of the properties of $`3j`$ symbols we obtain

$`c_\mathrm{}_1(z_1)\delta n_x(z_1)=4\pi {\displaystyle \underset{\mathrm{}_2=2}{\overset{\text{even}}{}}}\sqrt{{\displaystyle \frac{\mathrm{}_2(\mathrm{}_2+1)}{\mathrm{}_1(\mathrm{}_1+1)}}}\times `$
$`{\displaystyle _0^{\mathrm{}}}A_{\mathrm{}_1,\mathrm{}_2,1}(z_2z_1)\rho _\mathrm{}_2(z_2)\delta n_x(z_2)𝑑z_2.`$ (34)
Equation (III.4) is a homogeneous Fredholm equation of the second kind. It allows one to calculate the director profile (for small director deviations from the ground state, $`𝒏_0=𝒆_z`$) once the single-particle density $`\rho _0`$ of the ground state is known. Eqn (III.4) is not valid for every $`\mathrm{}_1`$, in spite of the derivation. This is because we are trying to map the single-particle density variation onto the director variation (eqn 32). However, this equation should be valid for the leading term in the density expansion, $`\mathrm{}_1=2`$, which we consider below.

First, we construct the solution to the eqn (III.4) in the cell bulk. The kernel $`A_{\mathrm{}_1,\mathrm{}_2,1}(z_2z_1)`$ is a short-ranged function. It equals zero for $`|z_2z_1|>a`$, where $`a`$ is the molecular length. The bulk director is a slowly-varying function on this length scale. Therefore we can expand it in a Taylor series

$$\delta n_x(z_2)=\delta n_x(z_1)+\delta n_x^{}(z_1)z_{12}+\delta n_x^{\prime \prime }(z_1)z_{12}^2/2+\mathrm{},$$
(35)
and restrict ourselves to second order in the expansion. Then the equation for the director (III.4) can be rewritten as a second-order linear differential equation

$$a_2(z)\delta n_x^{\prime \prime }(z)+a_1(z)\delta n_x^{}(z)+a_0(z)\delta n_x(z)=0,$$
(36)
where

$`a_n(z)`$ $`=`$ $`{\displaystyle \frac{4\pi }{\sqrt{6}}}{\displaystyle \underset{\mathrm{}_2=2}{\overset{\text{even}}{}}}\sqrt{\mathrm{}_2(\mathrm{}_2+1)}{\displaystyle _0^L}A_{2,\mathrm{}_2,1}(z_2z)\rho _\mathrm{}_2(z_2){\displaystyle \frac{z_{12}^n}{n!}}𝑑z_2c_2(z)\delta _{n,0}.`$ (37)
In the bulk, the expansion coefficients $`c_2`$ and $`\rho _\mathrm{}_2`$ do not depend on the $`z`$ coordinate. Expansion coefficients of the excluded area, $`A_{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}(z)`$, are even functions of $`z`$, which immediately implies that $`a_1(z)=0`$ in the bulk, since the integrated function $`zA_{\mathrm{}_1,\mathrm{}_2,\mathrm{}_r}(z)`$ is odd (because $`a_1(z)`$ is zero in the bulk, we expanded $`\delta n_x(z)`$ up to second order).

To prove that $`a_0(z)=0`$ in the bulk, we consider again the equation for the single particle density (12). Performing the usual expansion of density and $`f`$-function (13,14) and integrating over the angles and $`x,y`$ coordinates we obtain

$`c_0(z_1)`$ $`=`$ $`\sqrt{4\pi }\beta \mu 4\pi {\displaystyle _0^L}A_{0,0,0}(z_2z_1)\rho _0(z_2)𝑑z_2,`$
$`c_\mathrm{}_1(z_1)`$ $`=`$ $`4\pi {\displaystyle \underset{\mathrm{}_2=0}{\overset{\text{even}}{}}}{\displaystyle _0^L}A_{\mathrm{}_1,\mathrm{}_2,0}(z_2z_1)\rho _\mathrm{}_2(z_2)𝑑z_2,(\mathrm{}_10)`$ (38)
where $`A_{\mathrm{}_1,\mathrm{}_2,x}`$ are the expansion coefficients of the excluded area in a series of spherical harmonics. It is easy to show that the left hand side of the second equation in (38) equals $`a_0(z)`$ in the cell bulk. Indeed, using the symmetry of the excluded volume expansion coefficients we can write $`4\pi _0^LA_{2,\mathrm{}_2,1}(z_2z_\mathrm{b})𝑑z_2=4\pi _{\mathrm{}}^{\mathrm{}}A_{2,2,0}(z)𝑑z\delta _{2,\mathrm{}_2}=V_{22}\delta _{2,\mathrm{}_2}`$, which converts the expression for $`a_0(z)`$ to the left hand side of the second equation in (38). In fact, the conclusion $`a_0(z)=0`$ in the cell bulk is a consequence of the invariance of the grand potential with respect to director rotations.

Therefore, eqn (36) in the cell bulk reduces to $`\delta n_x^{\prime \prime }(z)=0`$ which means we have a linear dependence of the director on $`z`$ coordinate, in agreement with the result of elastic theory (9). This also means that, in the cell bulk, $`c_0+c_1z`$ is an eigenvector of the Fredholm equation (III.4).

To solve eqn (III.4) near the cell boundary is a much more challenging task. Numerically it can be done using, for example, an iterative method Press et al. (1992). Here we try to construct a crude analytical solution which will give us some qualitative understanding of what is happening in the interface region.

To begin with, we simplify eqn (III.4). The sum over $`\mathrm{}_2`$ in the right-hand side of eqn (III.4) is converging very fast (for ellipsoids with the elongation $`e=15`$, every term is about ten times smaller than the previous one). Therefore, we truncate the sum leaving only the $`\mathrm{}_2=2`$ term. Then, the kernel of the Fredholm equation $`A_{2,2,1}(z_1z_2)`$ can be expanded in a Taylor series

$$A_{2,2,1}(z_1z_2)=\underset{n=0}{\overset{\mathrm{}}{}}z_1^n\mathrm{\Phi }_n(z_2),$$
(39)
where $`\mathrm{\Phi }_n(z)`$ are some functions. A practical example of such an expansion is given in appendix B. The kernel is then separable, and the problem is reduced to the solution of a set of linear algebraic equations. Indeed, substituting eqn (39) into eqn (III.4) we obtain

$$\delta n_x(z)=c_2^1(z)\underset{n=0}{\overset{\mathrm{}}{}}b_nz^n,$$
(40)
where the $`b_n`$ are some constants. Substituting (40) back into eqn (III.4) we obtain an infinite set of linear equations for the coefficients $`b_n`$. To obtain an analytical expression, we perform further simplifications. First we note that, in the cell bulk, the director is a linear function of the $`z`$ coordinate. Therefore, to a good approximation, we may retain only the first two terms in eqn (III.4) since in the cell bulk $`c_2(z)=\text{const}`$. Using again the director profile given by elastic theory, eqn (9), we obtain an expression for the extrapolation length

$$\lambda =\frac{b_0}{b_1}=\frac{4\pi _0^{\mathrm{}}zA_{2,2,1}(z)\rho _2(z)c_2^1(z)𝑑z}{14\pi _0^{\mathrm{}}A_{2,2,1}(z)\rho _2(z)c_2^1(z)𝑑z}.$$
(41)
This expression is able to give qualitative explanations of the anchoring phenomenon in our system. In the ideal case, often considered in phenomenological approaches Yokoyama (1997), the density and order parameter are assumed to be constant in the cell, i.e. $`\rho _2(z)c_2^1(z)=\mathrm{const}`$. The anchoring appears only due to the presence of the interface, which breaks translational symmetry. According to eqn (38), $`c_2=V_{22}\rho _2`$, therefore $`4\pi _0^{\mathrm{}}A_{2,2,1}(z)\rho _2(z)c_2^1𝑑z=1/2`$. Then the anchoring coefficient is proportional to the first moment of the excluded area coefficient $`A_{2,2,1}(z)`$

$$\lambda =\frac{8\pi }{V_{22}}_0^{\mathrm{}}zA_{2,2,1}(z)𝑑z,$$
(42)
and does not depend on the value of density or order parameter, since excluded area, as well as excluded volume, are completely defined by the geometry of the overlapping molecules. This means that in the ideal case of a uniform nematic, the anchoring coefficient $`W=K_{33}/\lambda `$ has the same dependence on the order parameter as the elastic constant $`K_{33}`$.

In reality, one observes rather strong subsurface variations of the density and order parameter. The ratio $`\rho _2(z)/c_2(z)`$ then comes into play and contributes to the overall dependence of the extrapolation length on the order parameter. The numerical results which show this dependence are presented in sec. VI.2. A simple physical explanation is also possible. Presmectic ordering and higher density of the nematic phase in the interface region indicate that the nematic liquid crystal is at different state point. The mesophase is more ordered at this state point and this ordering is less sensitive to the density variation in the cell bulk. Since this ordering defines the director profile at the interface, it also affects the dependence of the extrapolation length on the state point.

## IV Thermal fluctuations

Another simulation method to measure bulk elastic constants $`K_{ii}`$ and the zenithal anchoring energy $`W`$ is based on the measurement of the director fluctuation amplitudes $`\delta 𝒏`$ in the liquid crystal cell Allen and Frenkel (1988, 1990); Allen (1999); Andrienko et al. (2000). Consider again slab geometry with homeotropic anchoring of the director at both cell surfaces. Consider a small perturbation of the director around the equilibrium distribution

$$𝒏(𝒓)=𝒆_z+\delta 𝒏(𝒓).$$
(43)
Minimizing the free energy (2) and linearizing the equations for the director and boundary conditions with respect to $`\delta 𝒏`$, we obtain

$`K_{11}(\delta n_{x,xx}+\delta n_{y,yx})+K_{22}(\delta n_{x,yy}\delta n_{y,yx})+K_{33}\delta n_{x,zz}=0,`$
$`K_{11}(\delta n_{x,xy}+\delta n_{y,yy})+K_{22}(\delta n_{x,xx}\delta n_{y,yx})+K_{33}\delta n_{x,zz}=0,`$
$$W\delta 𝒏+K_{33}\frac{}{z}\delta 𝒏|_{z=L}=\mathrm{𝟎},W\delta 𝒏K_{33}\frac{}{z}\delta 𝒏|_{z=0}=\mathrm{𝟎}.$$
(44)
Here we adopt the notation $`\delta n_{\alpha ,\beta \gamma }=_\beta _\gamma (\delta n_\alpha )`$.

We now expand $`\delta 𝒏(𝒓)`$ in a series of orthogonal functions

$`\delta 𝒏(𝒓)={\displaystyle \frac{1}{V}}{\displaystyle \underset{𝒒_{},q_z}{}}\mathrm{e}^{\mathrm{i}𝒒_{}𝒓_{}}\times `$
$`\left[\delta 𝒏^{(+)}(𝒒_{},q_z)\mathrm{e}^{\mathrm{i}q_zr_z}+\delta 𝒏^{()}(𝒒_{},q_z)\mathrm{e}^{\mathrm{i}q_zr_z}\right],`$ (45)
where $`𝒒_{}=(q_x,q_y)`$, and $`\delta 𝒏^{()}=(\mathrm{i}\chi \xi )/(\mathrm{i}\chi +\xi )\delta 𝒏^{(+)}`$ to satisfy the boundary conditions. Here we have introduced the dimensionless wave vector $`\chi =q_zL`$ and anchoring parameter

$$\xi =\frac{WL}{K_{33}}=\frac{L}{\lambda },$$
(46)
where $`\lambda `$ is the extrapolation length de Gennes and Prost (1995). The wave vectors $`q_z`$ form a discrete spectrum because of confinement in the $`z`$ direction which depends on the anchoring energy $`W`$. The explicit form of the $`q_z`$ spectrum is given by the secular equation:

$$\left(\xi ^2\chi ^2\right)\mathrm{sin}\chi +2\xi \chi \mathrm{cos}\chi =0.$$
(47)
In molecular simulations, rather than measuring director fluctuations, it is more convenient to measure fluctuations of the second-rank order tensor components (following Ref. Forster (1975)). We define the real-space order tensor density

$`Q_{\alpha \beta }(𝒓)`$ $`=`$ $`{\displaystyle \frac{V}{N}}{\displaystyle \underset{i}{}}\delta (𝒓𝒓_i)Q_{\alpha \beta }^i,`$
$`Q_{\alpha \beta }^i`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left(u_{i\alpha }u_{i\beta }{\displaystyle \frac{1}{3}}\delta _{\alpha \beta }\right),`$
where $`\alpha ,\beta =x,y,z`$, in terms of the orientation vectors $`𝒖_i`$ of each molecule $`i`$ (we consider only uniaxial molecules). If we assume that there is no variation in the degree of ordering, and neglect biaxiality of the order tensor, we may write

$$Q_{\alpha \beta }(𝒓)=\frac{3}{2}Qn_\alpha (𝒓)n_\beta (𝒓)\frac{1}{2}Q\delta _{\alpha \beta },$$

where $`Q`$ is the order parameter, i.e. the largest eigenvalue of $`Q_{\alpha \beta }(𝒓)`$.

Measurements are performed directly in reciprocal space. The Fourier transform of the real-space order tensor is

$$Q_{\alpha \beta }(𝒌)=_VQ_{\alpha \beta }(𝒓)\mathrm{e}^{\mathrm{i}𝒌𝒓}𝑑𝒓=\frac{V}{N}\underset{i}{}Q_{\alpha \beta }^i\mathrm{e}^{\mathrm{i}𝒌𝒓_i}.$$

Then the fluctuations $`\left|Q_{\alpha \beta }(𝒌)\right|^2`$ can be easily measured from simulations

$`|Q_{\alpha \beta }(k_z)|^2={\displaystyle \frac{V^2}{N^2}}[\left({\displaystyle \underset{i}{}}Q_{\alpha \beta }^i\mathrm{cos}(k_zz_i)\right)^2+`$
$`\left({\displaystyle \underset{i}{}}Q_{\alpha \beta }^i\mathrm{sin}(k_zz_i)\right)^2],`$ (48)
and compared with the prediction of elastic theory

$`|Q_{\alpha z}(k_z)|^2={\displaystyle \frac{9}{8}}k_\text{B}T{\displaystyle \frac{Q^2V}{K_{33}}}{\displaystyle \underset{q_z}{}}{\displaystyle \frac{\chi ^2+\xi ^2}{q_z^2\left(2\xi +\chi ^2+\xi ^2\right)}}\times `$
$`\left|{\displaystyle \frac{\mathrm{e}^{\mathrm{i}(\kappa +\chi )}1}{\kappa +\chi }}+\left({\displaystyle \frac{\mathrm{i}\chi \xi }{\mathrm{i}\chi +\xi }}\right){\displaystyle \frac{\mathrm{e}^{\mathrm{i}(\kappa \chi )}1}{\kappa \chi }}\right|^2`$ (49)
where $`\mathrm{}`$ denotes an ensemble average, $`\kappa =k_zL`$.

We measure $`Q`$ and $`\left|Q_{\alpha z}(k_z)\right|^2`$ from simulations, eqn (48), and then compare with the theoretical prediction, eqn (49), which is parametrized by $`L`$, $`\lambda `$ and $`K_{33}`$. Both the permitted $`q_z`$ spectrum, and the variation of $`\left|Q_{\alpha z}(k_z)\right|^2`$ with $`k_z`$, are sensitive to the anchoring strength parameter $`\xi =L/\lambda `$ .

## V Molecular model and simulation methods

We performed Monte Carlo (MC) simulation of a liquid crystal confined between parallel walls, with finite homeotropic anchoring at the walls. We used a molecular model which has been studied earlier in this geometry Allen (1999). The molecules in this study were modeled as hard ellipsoids of revolution of elongation $`e=a/b=15`$, where $`a`$ is the length of the major axis and $`b`$ the length of the minor axis. Such systems cannot form smectic or columnar phases Frenkel et al. (1984), and the phase transitions are not thermally driven as they are for most mesogens. Therefore, it should be borne in mind that simulation results are model specific, however the advantage of using systems of hard ellipsoids is that they exhibit a simple phase behaviour with some resemblance of that of real nematogens. In addition, as the elongation increases, the nematic-isotropic phase boundary approaches the Onsager limit. This eases the comparison between simulation and density-functional Onsager-type theories.

The phase diagram and properties of this family of models are well studied Frenkel et al. (1984); Frenkel and Mulder (1985); Tjipto-Margo et al. (1992); Allen (1993); Allen et al. (1993). It is useful to express the density as a fraction of the close-packed density $`\rho _{\text{c}p}`$ of perfectly aligned hard ellipsoids, assuming an affinely-transformed face-centered cubic or hexagonal close packed lattice. For this model, temperature is not a significant thermodynamic quantity, so it is possible to choose $`k_\text{B}T=1`$ throughout.

The slab geometry is defined by two hard parallel confining walls, which cannot be penetrated by the *centers* of the ellipsoidal molecules. Packing considerations generate homeotropic ordering at the surface, since the molecules have more free volume if their axes are normal to the wall. Surface anchoring in this system has been studied recently by applying an orienting perturbation at one of the walls and observing the response at the other Allen (1999) and by measuring amplitudes of the director fluctuations Andrienko et al. (2000). This yielded an estimate of the extrapolation length $`\lambda 17.5b1.16a`$ at one particular state point corresponding to the value of the order parameter $`Q0.85`$.

A sequence of runs was carried out for systems of $`N=2000`$ particles using the constant-$`NVT`$ ensemble, allowing typically $`10^5`$ MC sweeps for equilibration and $`4\times 10^7`$ sweeps for accumulation of averages (one sweep is one attempted move per particle). The wall separation was fixed and equal to $`L_z4.93a`$. Note, that the size of the simulation box $`L_z`$ is not equal to the liquid crystal cell thickness $`L`$ appearing in the elastic theory. The difference can be ascribed to partial penetration of the walls by the liquid crystal molecules, and/or formation of a solid layer near the walls. It is possible to write $`L=L_z+2L_w`$, where $`L_w`$ characterizes the wall. In our previous publication we found $`L_w0.3a`$ Andrienko et al. (2000).

The same molecular model and interaction of the molecules with the walls was adopted for the Onsager calculations. It was found sufficient to include terms with $`0l8`$ in the expansion of $`\mathrm{ln}\rho `$, while terms with $`0l10`$ were included in the expansion of the pair excluded volume coefficients $`V_{\mathrm{}\mathrm{}}`$ and the excluded area coefficients $`A_{\mathrm{}_1,\mathrm{}_2,m}`$.

## VI Results and discussion

### VI.1 Onsager theory, bulk

Minimization of the free energy functional (10) was carried out for several values of chemical potential $`\mu `$. From the single-particle density we evaluated Frank elastic constants, which, together with the values of the bulk density as a fraction of the close-packed density and the value of the order parameter are given in table 1.

The results are typical for the elastic constants of prolate bodies. They increase with fluid density (order parameter) and $`K_3>K_1>K_2`$ Tjipto-Margo et al. (1992).

It is essential to carry out bulk calculations if we want to know both anchoring extrapolation length $`\lambda =K_{33}/W`$ and anchoring strength $`W`$. The bulk problem provides us with the absolute values of the elastic constants; the elastic theory includes only ratios of elastic constants in the expression for the director profile.

### VI.2 Onsager theory, slab geometry

Minimization of the grand potential was carried out in slab geometry for the same values of chemical potential $`\mu `$ as considered in bulk. From the density and order parameter profiles we were able to extract values of these quantities in the central part of the cell which agreed with the bulk Onsager calculations.

Together with the simulation results, we plot dependences of the bulk order parameter versus density in Fig. 1.

The Onsager theory does not describe the bulk equation for the state perfectly: the predicted bulk density for a given order parameter value is larger than the value obtained in the simulation.

With this choice, the same values of the chemical potential and the same slab dimensions, minimization of the grand potential was carried out for the system with an external field applied near the right wall. Profiles of director angle $`\theta (z)`$ are compared with elastic theory (eqn (8)) in Fig. 2.

The elastic theory has been fitted to the director angle profiles predicted by the Onsager theory using two adjustable parameters, the extrapolation length $`\lambda `$ and the elastic constant ratio $`\delta `$. Note that only part of the bulk region has been used for fitting, $`20<z/b<80`$ where the elastic theory is applicable.

The anisotropy of elastic constants $`\delta `$ obtained from fitting is shown in Fig. 3, together with the results of calculations using the Poniewierski and Stecki expressions.

The value of $`\delta `$, and hence the splay constant $`K_1`$, comes into play only as the deformation angle $`\theta (z)`$ increases. Therefore, for the external field with easy axis $`\theta _L=\pi /4`$, the error in the determination of $`\delta `$ from simulation data is quite large.

The dependence of the extrapolation length, $`\lambda `$, on the order parameter $`Q`$ is shown in Fig 4. It is seen that Onsager theory predicts the extrapolation length to grow linearly with the order parameter, which is completely different from the experimental results in thermotropic liquid crystals, where $`\lambda `$ decreases with the increase of the order parameter as $`Q^2`$.

We have also carried out minimization of the grand potential for the system without an external field. As a result, we obtained the single-particle density with a homeotropic distribution of the director in the cell. The dependence of the extrapolation length on the order parameter was then calculated using eqn (41) and is also shown in Fig 4. The results qualitatively agree with the results of fitting obtained by combining Onsager theory and elastic theory. The extrapolation length tends to grow with increase in the order parameter and has the same order of magnitude. Plugging this single-particle density into all equations preceding eqn (41) we were able to check the approximations we did deriving this equation. We found that the most crude approximation is the truncation of the sum in eqn (40). This is not a fundamental problem and can be easily corrected by taking into account a sufficient number of expansion coefficients. However, this points out that the dependence of the director on the $`z`$ coordinate near the cell surface is different from the linear dependence in the cell bulk. Moreover, to obtain correct quantitative values of the anchoring coefficient, or extrapolation length, one needs to know the director distribution at the surface. The assumption $`\delta 𝒏=\text{const}`$ in the interface region, which has been made in Stelzer et al. (1997), may lead to absolutely incorrect estimates of the anchoring coefficient.

### VI.3 Simulation

Simulations were carried out in slab geometry for several values of the number density. The density variation $`0.25\rho /\rho _{\text{cp}}0.34`$ provides a sufficient range of order parameter variation in the nematic phase, $`0.68Q0.88`$ for us to study the effect of ordering on elastic behaviour.

The order tensor fluctuations in reciprocal space were calculated using expression (48). To fit the simulation results we used expression (49). Recall that the size of the simulation box $`L_z`$ is not necessarily equal to the liquid crystal cell thickness $`L`$ appearing in the elastic theory. We found that $`L=L_z+2L_\text{w}`$, with $`L_\text{w}4.5b=0.3a`$, almost independent of the density.

Using this value of $`L_\text{w}`$ we adjusted the extrapolation length $`\lambda `$ and the elastic constant $`K_{33}`$ to obtain the best fit to the fluctuation amplitudes. Typical order fluctuation amplitudes together with the corresponding fitting curves are plotted in Fig 5.

From this plot one can see that the fitting curves agree quite well with the simulation results for small values of the wave-vector $`k_z`$. At higher $`k_z`$, the structure is not perfectly reproduced, as one would expect for a theory valid for long-wavelength fluctuations, but the agreement is satisfactory.

The elastic constant, $`K_{33}b/k_\text{B}T`$ versus order parameter $`Q`$ is plotted in Fig 6. The agreement with the Onsager theory is satisfactory, especially if we take into account that the equation of state is not perfectly reproduced by the Onsager theory.

The dependence of the extrapolation length, $`\lambda `$, on the order parameter $`Q`$ is shown in Fig 4, together with the results from the Onsager theory. It should be noted that combining the elastic approach with the Onsager calculations does not allow one to determine, separately, $`L_\text{w}`$ and $`\lambda `$. Therefore, the results of the Onsager theory in Fig 4 really represents $`\lambda +L_\text{w}`$, which is one of the possible origins of the systematic difference between the two approaches.

To double check the results obtained by examining the director fluctuation amplitudes, we performed the same type of experiments as in the Onsager slab system. Within a range $`7.5b`$ of the right-hand wall, a strong coupling field was applied to molecular orientations, $`V^{\text{ext}}(𝒖_i𝒆_z)^2`$, aligning the molecules near the right wall parallel to it. After the system was equilibrated, the director profile was fitted to the result of the elastic theory, eqn (8). The dependence of the extrapolation length $`\lambda `$ on the order parameter $`Q`$ is also shown in Fig 4. One can see that the agreement between the two methods is quite good.

Finally, we plot the dependence of the anchoring energy coefficient, $`W=K_{33}/\lambda `$, which is an intrinsic charachteristic of the interface region, in Fig. 7. For Onsager theory, $`\lambda (Q)`$ is given by the results in slab geometry fitted with the elastic theory (Fig. 4), $`K_{33}(Q)`$ \- by the Poniewierski - Stecki expressions (Fig. 6). For MC simulations, we used the elastic constant obtained from the analysis of the director fluctuation amplitudes. All methods predict that the anchoring coefficient is a non-monotonic function of the order parameter, even though the actual variation is small.

## VII Conclusions

We have studied the dependence of the zenithal surface anchoring coefficient on the order parameter of a lyotropic nematic liquid crystal modeled by hard ellipsoids. Several techniques have been used: Onsager theory combined with elastic theory; Monte Carlo simulations fitted to elastic theory; analysis of the director fluctuation amplitudes obtained from Monte Carlo simulations. The results of these methods agreed qualitatively with each other.

Perhaps the most interesting aspect of this study is the increase of the anchoring extrapolation length with the increase of the nematic order parameter. This implies that the bulk elastic moduli in our system grow faster than the surface anchoring strength. This is opposite to the experimentally observed behaviour in thermotropic nematics, where the extrapolation length goes down with increase of the order parameter.

A microscopic semi-qualitative expression for the extrapolation length allowed us to conclude that subsurface variations of the single-particle density, mainly defined by the nematic order parameter and density variation, contribute substantially to the anchoring phenomenon. We showed that for an ideal system, in which the single-particle density in the cell is assumed to be uniform, the extrapolation length does not depend on the nematic order parameter. This dependence is therefore associated with the subsurface variations of the single-particle density.

We now turn to a brief discussion about possible generalizations of this work. First, it would be interesting to go beyond the limits of Onsager theory and use the direct pair correlation function in the nematic liquid instead of the Mayer $`f`$-function. Second, as was shown in sec. III.4, the anchoring coefficient can be calculated if we know the single-particle density and the direct pair correlation function (or Mayer $`f`$-function in case of Onsager theory). This has also been done by Stelzer et al Stelzer et al. (1997) with numerous approximations. Combining elastic theory and local density functional theory one can avoid these approximations, or at least check their validity. This work is in progress.

###### Acknowledgements.

This research was supported by EPSRC grants GR/L89990 and GR/M16023. D.A. acknowledges the support of the Overseas Research Students Grant; M.P.A. is grateful to the Alexander von Humboldt foundation, the British Council ARC programme, and the Leverhulme Trust.

## Appendix A Elastic constants

Performing the integrations over the angles in expression (23) and making use of the properties of $`3j`$ symbols and spherical harmonics one obtains

$`\beta K_{11}`$ $`=`$ $`\kappa (1,1,1,1,1,1),`$ (50)
$`\beta K_{22}`$ $`=`$ $`\kappa (1,1,1,1,1,1),`$
$`\beta K_{33}`$ $`=`$ $`\kappa (2,0,2,0,2,0),`$
where

$`\kappa (b_1,b_2,b_3,b_4,b_5,b_6)=`$ (67)
$`{\displaystyle \frac{4\pi ^2}{3}}{\displaystyle \underset{\mathrm{}}{\overset{\text{even}}{}}}\mathrm{}(\mathrm{}+1)\rho _{\mathrm{}}^2\left[\left(\begin{array}{ccc}\mathrm{}& \mathrm{}& 0\\ 1& 1& 0\end{array}\right)I^{\mathrm{},\mathrm{},0}+{\displaystyle \frac{b_1}{5}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}& 2\\ 1& 1& 0\end{array}\right)I^{\mathrm{},\mathrm{},2}+{\displaystyle \frac{b_2}{5}}\sqrt{{\displaystyle \frac{3}{2}}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}& 2\\ 1& 1& 2\end{array}\right)I^{\mathrm{},\mathrm{},2}\right]+`$
$`\sqrt{{\displaystyle \frac{(\mathrm{}+1)!}{(\mathrm{}3)!}}}\rho _{\mathrm{}}\rho _\mathrm{}2\left[{\displaystyle \frac{b_3}{5}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}2& 2\\ 1& 1& 0\end{array}\right)+{\displaystyle \frac{b_4}{5}}\sqrt{{\displaystyle \frac{3}{2}}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}2& 2\\ 1& 1& 2\end{array}\right)\right]I^{\mathrm{},\mathrm{}2,2}+`$
$`\sqrt{{\displaystyle \frac{(\mathrm{}+3)!}{(\mathrm{}1)!}}}\rho _{\mathrm{}}\rho _{\mathrm{}+2}\left[{\displaystyle \frac{b_5}{5}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}+2& 2\\ 1& 1& 0\end{array}\right)+{\displaystyle \frac{b_6}{5}}\sqrt{{\displaystyle \frac{3}{2}}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}+2& 2\\ 1& 1& 2\end{array}\right)\right]I^{\mathrm{},\mathrm{}+2,2},`$
and $`I^{\mathrm{}_1,\mathrm{}_2,\mathrm{}}`$ are radial integrals over the expansion coefficients of the Mayer $`f`$-function

$$I^{\mathrm{}_1,\mathrm{}_2,\mathrm{}}=𝑑rr^4f^{\mathrm{}_1,\mathrm{}_2,\mathrm{}}(r).$$
(68)
Expression (50) is the same as obtained before by Stelzer et al Stelzer et al. (1995), except that we used $`3j`$ symbols instead of Clebsch-Gordan coefficients and a different normalization of the single-particle density. To simplify (67) we use the relation between the $`3j`$ symbols

$$\sqrt{\frac{3}{2}}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}+2& 2\\ 1& 1& 2\end{array}\right)=\frac{1}{2}\left(\begin{array}{ccc}\mathrm{}& \mathrm{}+2& 2\\ 1& 1& 0\end{array}\right),$$
(69)
and take into account that sums with $`b_3`$, $`b_5`$ and $`b_4`$, $`b_6`$ are identical, since

$$\left(\begin{array}{ccc}\mathrm{}_1& \mathrm{}_2& 2\\ 1& 1& 0\end{array}\right)=(1)^{\mathrm{}_1+\mathrm{}_2}\left(\begin{array}{ccc}\mathrm{}_2& \mathrm{}_1& 2\\ 1& 1& 0\end{array}\right)=\left(\begin{array}{ccc}\mathrm{}_2& \mathrm{}_1& 2\\ 1& 1& 0\end{array}\right),$$

and $`I^{\mathrm{}_1,\mathrm{}_2,2}=I^{\mathrm{}_2,\mathrm{}_1,2}`$. Performing simplifications and taking into account explicit expressions for the $`3j`$ coefficients Edmonds (1974) we obtain expression (24).

## Appendix B Expansion of the kernel

From the numerical data we found that the kernel of the integral equation $`A_{221}(z)`$ can be accurately approximated with a Gaussian function

$$A_{2,2,1}(z)=\frac{V_{22}}{4\pi \sqrt{\pi }l}\mathrm{exp}\left[(z/l)^2\right],$$
(70)
where we took into account that $`4\pi _{\mathrm{}}^{\mathrm{}}A_{2,2,1}(z)𝑑z=V_{22}`$. Here $`l`$ is a geometrical parameter which depends only on the elongation of the molecules.

Using the generating function for the Hermite polynomials Gradshtein and Ryzhik (1994)

$$\mathrm{exp}(t^2+2tx)=\underset{n=0}{\overset{\mathrm{}}{}}H_n(x)\frac{t^n}{n!}$$
(71)
one can write the kernel in a separable form

$`A_{2,2,1}(z_2z_1)={\displaystyle \frac{V_{22}}{4\pi \sqrt{\pi }l}}\mathrm{exp}[(z_2/l)^2]\times `$
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}\left({\displaystyle \frac{z_1}{l}}\right)^nH_n\left(z_2/l\right),`$ (72)
similar to the general case we used in sec. III.4.

Using this approximation to the kernel one can show that the expression for the anchoring coefficient (42) reads

$$\lambda =l/\sqrt{\pi }.$$
(73)
For ellipsoids with elongation $`e=15`$, we found $`l5.46b`$. This results in an anchoring coefficient $`\lambda 3.1b`$ which is much lower than the actual value of the anchoring in the system with the density and order parameter variations at the surfaces. Hence, changes in the density, order parameter, and, as a result, in the director profile in the interface region cannot be neglected (see also sec. VI.2).