# Theory of inhomogeneous rod-like Coulomb fluids

Rudolf Podgornik  School of Physical Sciences and Kavli Institute for
Theoretical Sciences, University of Chinese Academy of Sciences, Beijing
100049, China  Wenzhou Institute of the University of Chinese Academy of
Sciences, Wenzhou, Zhejiang 325000, China  Institute of Physics, Chinese
Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>Also
affiliated with Department of Physics, Faculty of Mathematics and Physics,
University of Ljubljana, SI-1000 Ljubljana, Slovenia and Department of
Theoretical Physics, Jožef Stefan Institute, SI-1000 Ljubljana, Slovenia.

###### Abstract

A field theoretic representation of the classical partition function is
derived for a system composed of a mixture of anisotropic and isotropic mobile
charges that interact via long range Coulomb and short range nematic
interactions. The field theory is then solved on a saddle-point approximation
level, leading to a coupled system of Poisson-Boltzmann and Maier-Saupe
equations. Explicit solutions are finally obtained for a rod-like counterion-
only system in proximity of a charged planar wall, generalizing the standard
Gouy-Chapman results. The nematic order parameter profile, the counterion
density profile and the electrostatic potential profile are interpreted within
the framework of a nematic wetting layer with a (Donnan) potential difference.

##  I Introduction

Coulomb fluids composed of anisotropic charge carriers are ubiquitous in many
contexts. In fact it is worth noting that strong electrostatic interactions
between rod-like charges were already invoked in the case of nematic ordering
of tobacco mosaic virus (TMV) in the seminal work of Bernal and Fankuchen [1],
which is also one of the first cases of the application of the Poisson-
Boltzmann theory to biological systems [2]. Apart from the TMV, other charged
rod-like viruses and virus-like nanoparticles have been used in functional
materials assembly [3] whose formation is controlled by the electrostatic
interaction with the substrate [4].

A strong electrostatic attraction between the substrate and the deposited
filamentous viruses enhances a stable film formation, as is clear from the
effect of the ionic strength and the $p$H of the solution [5, 6]. Different
types of filamentous viruses [7], as well as short fragments of DNA [8],
F-actin [9], and cellular scaffold microtubules [10] all exhibit also
properties of polyvalent rod-like ions as do also many other multivalent
strongly anisotropic biological polyvalent ions [11] that can be either
modeled as spatially distributed point charges or as higher order point
multipoles [12, 13, 14]. Last but not least, ionic liquid crystals (ILCs) [15]
are solvent-less ionic systems with a dual ionic and organic nature [16], that
are composed of cations and anions with at least one ionic species and
characterized by highly anisotropic molecular shape [17].

The rod-like shape of ions leads to ordered structures whose formation
exhibits features of liquid crystal ordering as well as long-range
electrostatic interactions. It is this latter example that has recently
witnessed a real proliferation of theoretical works set to illuminate the
basic principles of order formation in these complex Coulomb systems.

The theoretical approaches to charged anisotropic systems in the bulk are
varied and abundant. Here we delimit to a few that are directly relevant for
subsequent developments. The work of Deutsch and Goldenfeld [18, 19] for thin
charged rods relies on a collective coordinate transformation method applied
to the ordering of charged rods in 3D. In addition, for rod-like charged
cylinders a generalized Onsager theory could also be used to describe the
ordering transition with electrostatic interactions strongly modifying the
hard core diameter of the rods as well as providing a mechanism for twisting
interaction as first described by Odijk [20, 21, 22]. This approach has seen
many further developments with different level modifications and extensions
[23, 24, 25, 26, 27]. A generalized variational field theory of particles with
rigid charge distributions [28] and an order parameter based mean-field
approximation of rod-like polyelectrolytes [29] both lead to an ordering
transition in 3D.

The properties of bulk systems composed of a mixture of multipolar charges
have also been analyzed in detail based on a field theoretic approach that
naturally incorporates also non-local dielectric response [30, 31, 32], while
the nematic order and electrostatics in the case of ion-doped nematic
electrolytes, with an anisotropic dielectric response, have been considered on
a phenomenological level [33]. Bulk properties as well as electrical double
layers in ionic liquid crystals have been analyzed in the work of Bier within
the density functional approach [34, 35] that was formulated for homogeneous
as well as inhomogeneous systems with interfaces [36, 37]. It is the latter
case that is particularly interesting as it should exhibit features of both,
the nematic ordering as well as the Gouy-Chapman-type electrostatic double
layers.

Figure 1: A schematic presentation of the system with rod-like cations and
simple anions in the bulk (left). An inhomogeneous system of rod-like cations
close to a charge surface (right). Both the cations as well as the anions are
charged, though in general differently. Apart from long range Coulomb
interactions the rod-like cations interact also via short range nematic
interactions described by a Lebwohl-Lasher interaction potential.

The problem of an inhomogeneous system with multipolar anisotropic charges has
also been addressed on various levels of modelling and approximations. The
case of a system bounded by charged wall(s) with mobile dipolar charges has
been analyzed with field theoretical approach [38, 39, 40], while quadrupolar
charges [12, 41] and finite size dumbbell charges [42, 13] have been
investigated on different levels of approximations, including field
theoretical approach. Density functional theory was formulated for the case of
inhomogeneous systems of charged anisotropic particles with interfaces [36,
37, 34], specifically a semiinfinite isotropic or nematic bulk system in
contact with a charged hard wall exhibiting nematic wetting of the substrate,
which is close to our point of departure.

Our focus here will be on how to formulate a theory of an inhomogeneous system
of anisotropic charges with Coulomb and nematic interactions \- as in the case
of rod-like cations close to an oppositely charged wall - based on a formalism
that could be seen as a straightforward generalization of the Gouy-Chapman
theory for simple ions. Thus, in what follows we will derive a statistical
theory of a Coulomb system composed of anisotropic rod-like cations and point-
like anions with microscopic interactions of the simple Lebwohl-Lasher nematic
as well as Coulomb type and derive its mean-field approximation form. We will
then apply the general theory to the case of an inhomogeneous system, composed
of a negatively charged planar interface in the presence of rod-like
counterions, thus generalizing the standard counterion-only Gouy-Chapman
problem of colloid physics. We will derive the equations governing the density
distribution and electrostatic potential on the mean-field level, in the case
of a one dimensional system with an electrified interface. By solving the
mean-field level equations that emerge as a coupled system of the Maier-Saupe
and the Poisson-Boltzmann equation, we are able to derive some salient
properties of inhomogeneous nematic ordering induced by the charged interface
as well as the modifications in the electric double layer distribution wrought
by the presence of nematic order.

##  II Collective description and field theoretic representation

###  II.1 Order parameters

Let us consider a system composed of monovalent anions (charge $-e$) and
polyvalent rod-like cations (charge $+qe$) with chemical potential $\mu$, a
situation often encountered in many ionic liquid system with cations being
extended stiff rods, while anions are considered to be point-like particles,
see Fig. 1.

The respective microscopic order parameters are the nematic order parameter
density defined as (see e.g. Ref. [43])

$\displaystyle\hat{\cal Q}_{ij}({\bf x})$ $\displaystyle\equiv$
$\displaystyle\sum_{(+)}{\textstyle\frac{3}{2}}\left({\bf n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\delta\left({\bf x}-{\bf
x}_{n}\right),$ (1)

with ${\bf n}({\bf x}_{n})$ the director of the $n-th$ rod-like cation,
together with the microscopic charge density

$\displaystyle\hat{\rho}_{(\pm)}({\bf
x})\equiv(+q,-1)\sum_{(\pm)}\delta\left({\bf x}-{\bf x}_{n}\right)$ (2)

so that the total charge density is $\hat{\rho}({\bf x})=\hat{\rho}_{(+)}({\bf
x})+\hat{\rho}_{(-)}({\bf x})$. The cations are polyvalent, with valency $q$,
and the anions are univalent $q=1$. The ensemble averaged forms of the above
order parameters are defined as

$\displaystyle\Big{<}\hat{\cal Q}_{ij}({\bf x})\Big{>}=c_{+}({\bf x})~{}{\cal
S}({\bf x})~{}{\textstyle\frac{3}{2}}\left({\bf n}_{i}({\bf x}){\bf
n}_{j}({\bf x})-{\textstyle\frac{1}{3}}\delta_{ij}\right)$ (3)

and

$\displaystyle\Big{<}\hat{\rho}_{+}({\bf x})\Big{>}=qc_{+}({\bf
x}),\qquad\Big{<}\hat{\rho}_{-}({\bf x})\Big{>}=-c_{-}({\bf x}),$ (4)

where the corresponding thermodynamic densities of the rod cations and anions
are $c_{+}({\bf x}),c_{-}({\bf x})$ and the nematic scalar order parameter is
${\cal S}({\bf x})$.

The microscopic Hamiltonian is assumed to be of the general Lebwohl-Lasher
type, though in the original formulation it was assumed to act only between
nearest neighbors [44]. This constraint was relaxed in a recent analytical
formulation [45]. By assumption then

$\displaystyle{\cal H}[{\bf r}_{n},{\bf n}_{n}]$ $\displaystyle=$
$\displaystyle{\textstyle\frac{1}{2}}\sum_{n,m}~{}u_{QQ}({\bf x}_{n}-{\bf
x}_{m})\left(\Big{(}{\bf n}_{n}{\bf
n}_{m}\Big{)}^{2}-{\textstyle\frac{1}{3}}\right)+$ (5)
$\displaystyle+{\textstyle\frac{1}{2}}\sum_{n,m}~{}u_{\rho\rho}({\bf
x}_{n}-{\bf x}_{m}),$

where $n,m$ run over all the particles in the system and ${\bf n}_{n},~{}{\bf
n}_{m}$ are the unit directors of the $n$-th and $m$-th particles. The
interaction strengths are all expressed in thermal units. The Lebwohl-Lasher
interaction type is chosen in order to simplify the derivation but is not
essential and can be easily generalized. The Onsager interactions are at least
to the lowest order equivalent to the first term in Eq. 5, where the nematic
interaction has a delta function range [46]. The theory presented here
therefore incorporates to the lowest order also the standard Onsager results.

The interaction energy Eq. 5 can be clearly recast as

$\displaystyle{\cal H}[{\bf r}_{n},{\bf n}_{n}]$ $\displaystyle=$
$\displaystyle{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\hat{\cal Q}_{ij}({\bf x})u_{QQ}({\bf x}-{\bf
x}^{\prime})\hat{\cal Q}_{ij}({\bf x}^{\prime})+$ (6)
$\displaystyle+{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\hat{\rho}({\bf x})u_{\rho\rho}({\bf x}-{\bf
x}^{\prime})\hat{\rho}({\bf x}^{\prime}),$

where we further assume that the scalar part $u_{\rho\rho}({\bf x}-{\bf
x}^{\prime})$ is due to the long-range Coulomb interactions. The coupling
between the two is omitted to the lowest order but can be considered for e.g.
dipolar or quadrupolar rods. The self-energies which were also omitted from
the formula above can be absorbed into the chemical potential when defining
the grand canonical partition function.

###  II.2 Field theoretical representation and thermodynamic relations

In the Appendix we derive the field theoretical representation of the grand
canonical partition function, $\Xi$, in terms of the tensor, $\Phi_{ij}({\bf
x})$, and scalar, $\phi({\bf x})$, auxiliary fields with

$\displaystyle{\Xi}[\Phi_{ij}({\bf x});\phi({\bf x})]$ $\displaystyle=$
$\displaystyle\int\\!\\!{\cal D}[\Phi_{ij}({\bf x})]\int\\!\\!{\cal
D}[\phi({\bf x})]~{}e^{-\beta{\underline{\cal H}}[\Phi_{ij}({\bf x});\phi({\bf
x})]},$ (7)

with the effective field action as

$\displaystyle-\beta{\underline{\cal H}}[\Phi_{ij}({\bf x});~{}\phi({\bf x})]$
$\displaystyle=$
$\displaystyle~{}-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}{\Phi}_{ij}({\bf x})u_{QQ}^{-1}({\bf x}-{\bf
x}^{\prime}){\Phi}_{ij}({\bf
x}^{\prime})-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}{\phi}({\bf x})u_{\rho\rho}^{-1}({\bf x}-{\bf
x}^{\prime}){\phi}({\bf x}^{\prime})+$ (8)
$\displaystyle+~{}\lambda_{(+)}\int_{V}d{\bf x}~{}e^{i\phi({\bf x})+\ln{{\cal
P}(\Phi_{ij}({\bf x}))}}+\lambda_{(-)}\int_{V}d{\bf x}~{}e^{-i\phi({\bf x})},$

where $\lambda_{(\pm)}$ are the fugacities of the two ionic species and ${\cal
P}(\Phi_{ij}({\bf x}))$ is the orientational partition function of a single
particle, see Eq. 55,

$\displaystyle{\cal P}(\Phi_{ij}({\bf
x}))\equiv\Big{<}e^{i{\textstyle\frac{3}{2}}\left({\bf n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Phi_{ij}({\bf
x})}\Big{>}_{\Omega}.$ (9)

The orientational trace is defined in such a way that $<1>_{\Omega}=1$. From
the expressions above it is clear that for the cations the interaction with
the fields has a scalar, electrostatic component, and a tensor nematic
component given by the orientational entropy. While the field theoretical
representation of the partition function cannot be solved analytically, it
suggests a number of developments that lead to meaningful approximations, such
as the saddle-point approximation that we will detail in what follows. In
addition other approximations that turned out to be useful in the context of
isotropic Coulomb fluids could be exploited as well [47].

There is a number of thermodynamic relations that need to be satisfied by the
derivatives of the partition function. In the case of a set chemical
potential, the fugacity and the average number of molecules in the system are
given by the standard relations. In addition, the invariance of the functional
integral with respect to the linear transformation of the fluctuating fields
yields two relations

$\displaystyle\Big{<}\frac{\delta{\underline{\cal H}}[\Phi_{ij}({\bf
x});~{}\phi({\bf x})]}{\delta\Phi_{ij}({\bf x})}\Big{>}=0\qquad{\rm
and}\qquad\Big{<}\frac{\delta{\underline{\cal H}}[\Phi_{ij}({\bf
x});~{}\phi({\bf x})]}{\delta\phi({\bf x})}\Big{>}=0,$ (10)

where the averages stand for

$\displaystyle\\!\\!\\!\\!\\!\Big{<}\dots\Big{>}=\frac{\int{\cal
D}[\Phi_{ij}({\bf x})]\int{\cal D}[\phi({\bf
x})]~{}\dots~{}e^{-\beta{\underline{\cal H}}[\Phi_{ij}({\bf x});\phi({\bf
x})]}}{\int{\cal D}[\Phi_{ij}({\bf x})]\int{\cal D}[\phi({\bf
x})]~{}e^{-\beta{\underline{\cal H}}[\Phi_{ij}({\bf x});\phi({\bf
x})]}}.~{}~{}~{}~{}$ (11)

These two relations, Eq. 10, furthermore yield the following exact connections
between the average values of the order parameters and the auxiliary fields

$\displaystyle\Big{<}\hat{\cal Q}_{ij}({\bf x})\Big{>}=-i\int_{V}d{\bf
x}^{\prime}~{}u^{-1}_{QQ}({\bf x}-{\bf x}^{\prime})\Big{<}\Phi_{ij}({\bf
x}^{\prime})\Big{>}=$
$\displaystyle=-i\lambda_{(+)}~{}\Big{<}\frac{\partial\ln{{\cal
P}(\Phi_{ij}({\bf x}))}}{\partial\Phi_{ij}({\bf x})}e^{i\phi({\bf
x})+\ln{{\cal P}(\Phi_{ij}({\bf x}))}}\Big{>}=$
$\displaystyle=\lambda_{(+)}~{}\Big{<}\Big{<}\\!\\!\\!\Big{<}{\textstyle\frac{3}{2}}\left({\bf
n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Big{>}\\!\\!\\!\Big{>}_{\Omega}e^{i\phi({\bf
x})+\ln{{\cal P}(\Phi_{ij}({\bf x}))}}\Big{>},$ (12)

where we took into account the definition of the double brackets Eq. 14.
Furthermore by analogy

$\displaystyle\Big{<}{\rho}({\bf x})\Big{>}=-i\int_{V}d{\bf
x}^{\prime}~{}u^{-1}_{\rho\rho}({\bf x}-{\bf x}^{\prime})\Big{<}\phi({\bf
x}^{\prime})\Big{>}=$
$\displaystyle~{}~{}~{}~{}+q\lambda_{(+)}\Big{<}~{}e^{i\phi({\bf x})+\ln{{\cal
P}(\Phi_{ij}({\bf x}))}}\Big{>}-\lambda_{(-)}~{}\Big{<}e^{-i\phi({\bf
x})}\Big{>}.$ (13)

Note the difference between the orientational averages
$\Big{<}\\!\\!\\!\Big{<}...\Big{>}\\!\\!\\!\Big{>}_{\Omega}$ ,
$\Big{<}\dots\Big{>}_{\Omega}$ defined as

$\displaystyle\Big{<}\\!\\!\\!\Big{<}...\Big{>}\\!\\!\\!\Big{>}_{\Omega}=\frac{\Big{<}(\dots)e^{i{\textstyle\frac{3}{2}}\left({\bf
n}_{i}{\bf n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Phi_{ij}({\bf
x})}\Big{>}_{\Omega}}{\Big{<}e^{i{\textstyle\frac{3}{2}}\left({\bf n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Phi_{ij}({\bf
x})}\Big{>}_{\Omega}}$ (14)

where the unnormalized orientational average, $\Big{<}\dots\Big{>}_{\Omega}$,
is with respect to the distribution Eq. 9.

On the saddle-point level we will soon see that the above two equations are
actually the modified Maier-Saupe (MS) self-consistent equation and the
modified Poisson-Boltzmann (PB) self-consistent equation for the tensor and
scalar fields, respectively.

###  II.3 Saddle-point approximation

Since the field action Eq. 54 is non-linear, no further exact developments are
feasible and one needs to resort to the saddle-point approximation that yields
two mean-field equations for the two auxiliary fields. At the saddle-point the
fields are pure imaginary, so that one can transform

$\Phi_{ij}({\bf x})\longrightarrow-i\Phi_{ij}^{*}({\bf x})\quad{\rm
and}\quad\phi({\bf x})\longrightarrow i\phi^{*}({\bf x}).$ (15)

The thermodynamic averages $\Big{<}\dots\Big{>}$ are given at the value of the
mean-field, and the two self-consistent field equations, Eqs. 12, 13, are then
reduced to

$\displaystyle\Big{<}\hat{\cal Q}_{ij}({\bf
x})\Big{>}=\lambda_{(+)}~{}\Big{<}\\!\\!\\!\Big{<}{\textstyle\frac{3}{2}}\left({\bf
n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Big{>}\\!\\!\\!\Big{>}_{\Omega}e^{-\phi^{*}({\bf
x})+\ln{{\cal P}(-i\Phi^{*}_{ij}({\bf x}))}}=c_{(+)}(\phi^{*}({\bf
x}))\Big{<}\\!\\!\\!\Big{<}{\textstyle\frac{3}{2}}\left({\bf n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Big{>}\\!\\!\\!\Big{>}_{\Omega},$
(16)

with the orientational partition function ${\cal P}(-i\Phi^{*}_{ij}({\bf x}))$
defined in Eq. 9 taken at the imaginary value of the tensor auxiliary field,
and

$\displaystyle\Big{<}{\rho}({\bf
x})\Big{>}=+q\lambda_{(+)}~{}e^{-\phi^{*}({\bf x})+\ln{{\cal
P}(-i\Phi^{*}_{ij}({\bf x}))}}-\lambda_{(-)}~{}e^{+\phi^{*}({\bf
x})}=+qc_{(+)}(\phi^{*}({\bf x}))-c_{(-)}(\phi^{*}({\bf x})),$ (17)

with obvious definitions of the cation and anion densities,
$c_{(+)}(\phi^{*}({\bf x}))$ and $c_{(-)}(\phi^{*}({\bf x}))$. In what follows
we will assume a short range attractive orientational potential and a
Coulombic positional potential. This implies

$\displaystyle u_{QQ}^{-1}({\bf x}-{\bf
x}^{\prime})=-u_{QQ}^{-1}(0)\delta({\bf x}-{\bf x}^{\prime})$ (18)

and

$\displaystyle u_{\rho\rho}^{-1}({\bf x}-{\bf
x}^{\prime})=-\varepsilon\nabla^{2}\delta({\bf x}-{\bf x}^{\prime})$ (19)

where $\varepsilon\equiv\epsilon_{0}\epsilon$, with $\epsilon$ the dielectric
permittivity of the solvent and $\epsilon_{0}$ the permittivity of space.
Other forms of the nematic interaction potential are of course possible but
would lead to more complicated tensorial saddle-point equations.

With this in mind we then derive the tensorial part in the form of a modified
Meier-Saupe (MS) equation

$\displaystyle\\!\\!\\!\\!\\!\\!{-u_{QQ}^{-1}(0)~{}\Phi^{*}_{ij}({\bf
x})+c_{(+)}~{}\Big{<}\\!\\!\\!\Big{<}{\textstyle\frac{3}{2}}\left({\bf
n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Big{>}\\!\\!\\!\Big{>}_{\Omega}=0}.$
(20)

We will see that the Maier-Saupe equation determines only the nematic order
parameter but not the orientation in the ordered phase, which is assumed to be
homogeneous. The orientation of the ordering would be determined by the
substrate-nematic interaction potential in line with the model described in
[48]. Without any loss of generality we can assume the ordering is
perpendicular to the bounding surface. To describe the orientational
relaxation effects one would need also the elastic deformation energy which
would stem from the expansion of the nematic interaction potential w.r.t. the
gradient of the tensorial order parameter as in the general inhomogeneous
Landau-de Gennes Ansatz [48]. We assume that the that the characteristic
length of the nematic order relaxation is much shorter then the electrostatic
relaxation length, either the Gouy-Chapman length or the Debye length, and can
thus be neglected. We discuss the possible generalizations in the Discussion
section.

The scalar part of the mean-field equations can be written in the form of the
Poisson-Boltzmann (PB) equation as

$\displaystyle{\varepsilon\nabla^{2}\phi^{*}({\bf x})+qc_{(+)}-c_{(-)}=0}.$
(21)

The two mean-field equations Eqs. 20 and 21 correspond to the nematic and
electrostatic degrees of freedom in a similar manner that the Edwards and the
polymer Poisson-Boltzmann equation correspond to polymer and electrostatic
degrees of freedom for charged flexible polymers [49, 50]. Consequently any
other degree of freedom would introduce its own mean-field equation. Clearly
the approximations entailed in the derivation of the mean-field equations
imply that the spatial relaxation of both fields is given solely by the
electrostatic component, so that in this respect the nematic field is
subservient to the electrostatic field.

Inserting the mean-field Ansatz into the free energy Eq. 54 we obtain the
general form of the mean-field free energy ${\cal F}[\Phi_{ij}^{*}({\bf
x});~{}\phi^{*}({\bf x})]$ as a functional of the mean-field nematic
potential, $\Phi_{ij}^{*}({\bf x})$, and the mean-field electrostatic
potential, $\phi^{*}({\bf x})$ as

$\displaystyle\beta{\cal F}[\Phi_{ij}^{*}({\bf x});~{}\phi^{*}({\bf x})]=$
$\displaystyle=$
$\displaystyle~{}-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\Phi_{ij}^{*}({\bf x})u_{QQ}^{-1}({\bf x}-{\bf
x}^{\prime})\Phi_{ij}^{*}({\bf
x}^{\prime})-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}\phi^{*}({\bf x})u_{\rho\rho}^{-1}({\bf x}-{\bf
x}^{\prime})\phi^{*}({\bf x}^{\prime})-$ (22)
$\displaystyle-~{}\lambda_{(+)}\int_{V}d{\bf x}~{}e^{-\phi^{*}({\bf
x})+\ln{{\cal P}(-i\Phi_{ij}({\bf x}))}}-\lambda_{(-)}\int_{V}d{\bf
x}~{}e^{+\phi^{*}({\bf x})}.$

The mean-field theory for rod-like cations thus couples the MS equation for
simple liquid crystals to the PB equation for simple Coulomb fluid, except
that formally the electrostatic potential of the rod-like cations is
transformed to $\phi^{*}({\bf x})\longrightarrow\phi^{*}({\bf x})-\ln{{\cal
P}(-i\Phi^{*}_{ij}({\bf x}))}$. Clearly one could expand the above free energy
in terms of the nematic order parameter obtaining a generalized Landau-de
Gennes free energy but then lose the direct connection with the Gouy-Chapman
theory. As already stated in the Introduction we aim to keep and develop this
connection.

The theory was formulated specifically for a rod-like cationic species and a
point-like anionic species, but based on the methodology other possibilities
are just as amenable to the same procedure of deriving the field theory as
well as the mean-field equations.

##  III Rod-like counterion-only system in one dimension

There are of course countless cases that one can dwell on in order to apply
the general theory, and thus some selectivity is in order. In what follows we
shall delimit ourselves to the case of a rod-like counterion-only system in
the presence of a single electrified surface, which is a direct generalization
of the standard Gouy-Chapman paradigm of the colloid electrostatics. The only
spatial dependence is in the direction of the surface normal, which we choose
to coincide with the $z$ axis. One also needs to remember that the counterion
only system has no reservoir and thus no chemical potential, so that the
number of counterions is set only by the neutralizing charge on the bounding
surface. Nevertheless this does not affect the above derivation.

###  III.1 Coupled system of Maier-Saupe and Poisson-Boltzmann equations

The scalar mean-field equation in this case is the modified Poisson-Boltzmann
equation that is obtained as

$\displaystyle\varepsilon{\phi^{*}}^{\prime\prime}(z)+qc_{(+)}(\phi^{*}(z))=0,$
(23)

where ${\phi^{*}}^{\prime\prime}(z)=\frac{d^{2}{\phi^{*}}(z)}{dz^{2}}$ is the
second derivative of the mean potential with respect to $z$. Comparing the
above expression with the standard PB equation for a point-like counterion-
only case, the difference is manifest in the presence of the orientational
entropy of the rods, equal to $\ln{{\cal P}(-i\Phi^{*}_{ij}({\bf x}))}$. We
will solve this equation later.

The tensor mean-field equation, i.e., the modified Maier-Saupe equation, can
be solved by assuming a non-vanishing orientational order characterized with
the director $\hat{\bf n}$

$\displaystyle\Phi^{*}_{ij}(z)$ $\displaystyle=$
$\displaystyle{s}(z)~{}{\textstyle\frac{3}{2}}\left(\hat{\bf n}_{i}\hat{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right),$ (24)

so that

$\displaystyle\Phi^{*}_{ij}(z)^{2}={\textstyle\frac{3}{2}}~{}{s}^{2}(z),$ (25)

and ${s}(z)$, proportional but not equal to the nematic order parameter, then
simply represents the absolute value of $\Phi^{*}_{ij}(z)$. In infinite space
without any anchoring fields, $\hat{\bf n}$ is arbitrary but would be set by
surface free energy terms [48], which we did not invoke explicitly since the
Maier-Saupe equation is local and thus valid for any orientation. One can
recall, Eq. 3 and Eq. 24, which means that

$\displaystyle{\cal
Q}^{*}_{ij}=u^{-1}_{QQ}(0)\Phi^{*}_{ij}=u^{-1}_{QQ}(0){s}~{}{\textstyle\frac{3}{2}}\left(\hat{\bf
n}_{i}\hat{\bf n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)=c_{+}~{}{\cal
S}~{}{\textstyle\frac{3}{2}}\left(\hat{\bf n}_{i}\hat{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)$ (26)

wherefrom we derive the exact connection between the nematic order parameter
${\cal S}$, strictly limited to the interval $0\leq{\cal S}\leq 1$, and the
parameter $s$ introduced in Eq. 24 as

${\cal S}=\frac{s}{u_{QQ}(0)~{}c_{+}}.$ (27)

By multiplying both sides of Eq.20 by $\left(\hat{\bf n}_{i}\hat{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)$ and tracing over the
indices, we finally obtain

$\displaystyle{s}(z)$ $\displaystyle=$
$\displaystyle{\textstyle\frac{3}{2}}u_{QQ}(0)c_{(+)}(\phi^{*}(z))~{}\Big{<}\\!\\!\\!\Big{<}({\bf
n}\cdot\hat{\bf
n})^{2}-{\textstyle\frac{1}{3}}\Big{>}\\!\\!\\!\Big{>}_{\Omega}.$ (28)

Since the above relation is a local one, it remains the same at any position
$z$, and the explicit dependence on the coordinate can be dropped in what
follows. This is a consequence of the assumption that the nematic order
relaxation is much shorter then the electrostatic relaxation length.
Developing further, this eventually yields

$\displaystyle\Big{<}\\!\\!\\!\Big{<}({\bf n}\cdot\hat{\bf
n})^{2}-{\textstyle\frac{1}{3}}\Big{>}\\!\\!\\!\Big{>}_{\Omega}=\Big{<}\\!\\!\\!\Big{<}\cos{\theta}^{2}-{\textstyle\frac{1}{3}}\Big{>}\\!\\!\\!\Big{>}_{\Omega}=\frac{\partial}{\partial\gamma}\ln{J(\gamma)},$
(29)

where obviously $\gamma={\textstyle\frac{9}{4}}{s}$ and from Eq. 14 it
furthermore follows that the orientational partition function is obtained as

$\displaystyle
J(\gamma)=\Big{<}e^{\gamma\left(\cos{\theta}^{2}-{\textstyle\frac{1}{3}}\right)}\Big{>}_{\Omega}=\int_{0}^{1}dz~{}e^{\gamma\left(z^{2}-{\textstyle\frac{1}{3}}\right)}=\frac{e^{{\textstyle\frac{2}{3}}\gamma}}{\sqrt{\gamma}}D(\sqrt{\gamma}),$
(30)

with $D(x)$ as the Dawson’s integral [51]. The Maier-Saupe expression can then
be cast definitively as

$\displaystyle{\gamma}={\textstyle\frac{27}{8}}u_{QQ}(0)c_{+}\frac{\partial}{\partial\gamma}\ln{J(\gamma)}$
(31)

where the nematic order parameter is then extracted from Eq. 27 as ${{\cal
S}={\gamma}/({\frac{9}{4}u_{QQ}(0)~{}c_{+}}}).$ In 2D the same type of
analysis would just replace the Dawson integral by a modified Bessel function
integral [52]. Eq. 31 corresponds exactly to the Kleinert formulation [51] of
the Maier-Saupe theory, if one takes his coupling strength of nematic
interactions equal to ${\cal A}_{0}={\textstyle\frac{3}{2}}u_{QQ}(0)c_{+}$.

Note again that the mean-field Maier-Saupe equation is a local equation that
pertains to every point in the domain, which is a consequence of the fact that
in the Ansatz Eq. 18 we only considered local nematic interactions.

The solution of the Maier-Saupe equation for $0\leq{\cal S}\leq 1$ exhibits a
first order isotropic-nematic transition at a critical value of $c_{(+)}$,
where the order parameter makes a discontinuous jump from an isotropic phase
to a nematic phase. It is important to realize that the solutions of the MS
theory for this counterion only system differ from the standard case, where
the two phases are kept under the same chemical potential that allows for the
density jump at the transition. The counterion only system is not regulated by
a chemical potential and exhibits no density jump.

###  III.2 First integral and the phase portrait analysis

Figure 2: Phase portrait solution of the Maier-Saupe equation. Phase portrait
$P(Q)$ obtained from the mean-field equation, Eq. 37. The solution starts with
$P_{0}$ obtained from the boundary condition Eq. 39 and then moves towards
$P=0$ far away from the surface. If $P_{0}$ is above the IN transition point,
e.g. ”A”, the solution first follows the MS branch (full line) and then at the
transition migrates to the PB branch that it follows until the electrostatic
field levels off to zero far away from the surface, just as in the case of
point-like counterions. On the other hand, if $P_{0}$ is below the transition
point, e.g. ”B”, the solution simply follows the PB line as if the counterions
are point-like. (Inset) The dependence of the nematic order parameter ${\cal
S}$ on the dimensionless density $Q$, defined in Eq. 36. The isotropic-nematic
(IN) transition takes place at a critical value $Q_{c}\simeq 3.$ and the jump
in the order parameter amounts to $0.313$.

Introducing now the generalized van’t Hoff osmotic pressure as

$\displaystyle p(\phi^{*}(z),s(z))$ $\displaystyle=$ $\displaystyle
c_{+}(\phi^{*}(z),s(z)),$ (32)

it is possible to write the mean-field equations in a ”Lagrangian” form (see
Ref. [53])

$\displaystyle\varepsilon\phi^{*}(z)^{\prime\prime}=\frac{\partial
p(\phi^{*}(z),s(z))}{\partial\phi^{*}(z)},$ $\displaystyle
s(z)=-{\textstyle\frac{2}{3}}u_{QQ}(0)\frac{\partial
p(\phi^{*}(z),s(z))}{\partial s(z)}.$ (33)

Proceeding now as in the case of the first integral of the Gouy-Chapman theory
[54], by multiplying the first equation by ${\phi^{*}}^{\prime}$ and the
second one by $s^{\prime}(z)$ and then summing them up, we obviously remain
with

$\frac{\partial p}{\partial\phi^{*}}{\phi^{*}}^{\prime}+\frac{\partial
p}{\partial
s(z)}s^{\prime}(z)=\Big{(}\varepsilon\phi^{*}(z)^{\prime\prime}\phi^{*}(z)^{\prime}+{\textstyle\frac{3}{2}}u_{QQ}(0)^{-1}s(z)s^{\prime}(z)\Big{)}=\frac{d}{dz}\Big{(}{\textstyle\frac{1}{2}}\varepsilon({\phi^{*}}^{\prime}(z))^{2}-{\textstyle\frac{3}{4}}~{}{u_{QQ}(0)}^{-1}s^{2}(z)\Big{)},$
(34)

which can be cast into the form of the first integral that generalizes the
standard Poisson-Boltzmann result

$\displaystyle\\!\\!\\!\\!\\!\\!\\!{\textstyle\frac{1}{2}}\varepsilon({\phi^{*}}^{\prime}(z))^{2}-{\textstyle\frac{3}{4}}\frac{s^{2}(z)}{{u_{QQ}(0)}}-p(\phi^{*}(z),s(z))=const.$
(35)

Because of the assumption of the short range nematic interactions, the field
$s(z)$ has no associated ”dynamics”, i.e., the first integral contains no
derivatives of $s(z)$.

The ”Lagrange equations” Eq. 33 can be solved and plotted in a phase portrait
mode that has been invoked previously by Pandit and Wortis to describe the
phase equilibria in lattice models with surfaces and interfaces [55]. The
phase portrait method allows for an easy and physical visualization of the
mean-field theories and has been successfully applied also to the case of the
Poisson-Boltzmann case [56].

Figure 3: The dimensionless counterion density profile $Q$ obtained from Eq.
40 as a function of the dimensionless separation from the bounding surface for
different values of $Q_{0}$ as indicated on the figure. For $Q_{0}\leq 3.$ the
dimensionless density dependence on $\tilde{z}$ coincides with the PB branch.
For $Q_{0}>3.$ it stays on the PB branch for $Q\leq Q_{0}$ but then exhibits a
jump at the transition point and continues on the MS branch afterwards. The
dimensionless density profile is continuous everywhere, but displays a
discontinuity of the derivative at the wetting layer. In addition for
$Q_{0}>3.$ the PB branch does not correspond to $Q_{PB}(\tilde{z}=0)=Q_{0}$ at
the surface but to a rescaled value such that $Q_{MS}(\tilde{z}=0)=Q_{0}$, the
difference stemming from the form of the boundary condition Eq. 39 on the MS
and PB branches.

One can now obtain the full implicit form of the solution of the mean-field
equations by introducing two new variables

$\displaystyle P$ $\displaystyle=$
$\displaystyle\sqrt{u_{QQ}(0)\varepsilon}~{}{\phi^{*}}^{\prime}(z)$
$\displaystyle Q$ $\displaystyle=$ $\displaystyle
u_{QQ}(0)\lambda_{(+)}~{}e^{-\phi^{*}(z)+\ln{J\big{(}\gamma(z)\big{)}}}=u_{QQ}(0)~{}c,$
(36)

where $c\equiv c_{+}$ is the counterion density, Eq. 32. $Q$ can thus be
interpreted as the dimensionless density of the rod-like counterions and $P$
is the dimensionless electrostatic field. It follows from the first integral
Eq. 35 that

$\displaystyle P=\pm\sqrt{2Q+{\textstyle\frac{8}{27}}{\gamma^{2}(Q)}},$ (37)

where for a single charged surface the constant in the first integral can be
obtained as vanishing, just as in the standard Gouy-Chapman case, meaning that
the osmotic pressure of the system is zero. The solution of the problem is
therefore completely specified by the dependence $P=P(Q)$, while $\gamma(Q)$,
Eq. 31, becomes a solution of

$\displaystyle\gamma(Q)={\textstyle\frac{27}{8}}Q~{}\frac{\partial}{\partial\gamma(Q)}{\ln
J\big{(}\gamma(Q)\big{)}}.$ (38)

Numerical solutions of the above equations are presented in Fig. 2 in the form
of the dependence of the nematic order parameter ${\cal S}(Q)$ on the
dimensionless rod-like counterion density, with the jump from zero to $0.313$
at the critical value $Q=3.$, obviously corresponding to a first order
isotropic-nematic transition. The phase portrait of the system, $P=P(Q)$,
presents the dependence of the dimensionless electrostatic field $P$ on the
dimensionless density showing the pure PB branchand the MS branch that
separates from the PB branch at the nematic transition point.

The surface charge density at the bounding surface sets the boundary value to

$\displaystyle
P_{0}^{2}=2Q_{0}+{\textstyle\frac{8}{27}}{\gamma^{2}(Q_{0})}={u_{QQ}(0)}\frac{\sigma^{2}}{\varepsilon},$
(39)

which follows directly from the Gauss boundary condition for the electrostatic
field at a surface with surface charge density $\sigma$, i.e.
$\varepsilon{\phi^{*}}^{\prime}(z=0)=\sigma$. The solution,
$Q_{0}=Q_{0}(\sigma)$, then implies that the $P(Q)$ curve starts at $P_{0}$
and then moves along the solution line to $P=0$, corresponding to the
vanishing of the electrostatic field far from the surface.

For $P_{0}$ large enough (see the $A$ intercept in Fig. 2) the solution first
follows the MS branch in the ordered phase until it reaches the transition
point. The system thus exhibits a finite thickness nematic wetting layer close
to the surface. After that the solution migrates to the PB branch, following
it until the electrostatic field levels off to zero far away from the surface.
On the other hand for $P_{0}$ smaller then a critical value (see the $B$
intercept in Fig. 2) there is no nematic wetting and the system remains
disordered along the whole $P(Q)$ solution, following the PB line as if the
counterions are point-like.

Clearly far away from the charged surface the system is disordered while in
the proximity to the surface, where electrostatic attraction between the
cations and the negatively charged surface increases their local
concentration, it orders up, creating a surface wetting layer of the nematic
phase.

A note of caution is in order here: in the nematic wetting context the wetting
layer and specifically its thickness is a result of the competition between
the order parameter inhomogeneities, described by the order parameter gradient
terms in the free energy, and the bulk free energy [57, 58, 48]. The wetting
layer and its properties in our case is related but does not coincide with the
the standard definition of nematic wetting. As already stated, we assumed that
the spatial relaxation of both the nematic as well as the electrostatic fields
is governed solely by the electrostatic component, so that in this respect the
nematic field is subservient to the electrostatic field. A more complete but
unfortunately also much more complicated theory would have to encompass both
relaxations separately.

Figure 4: Dimensionless mean-field electrostatic potential difference
$\phi^{*}(Q)-\phi^{*}(Q_{0})$ from Eq. 41 as a function of the dimensionless
distance from the surface $\tilde{z}$ for different values of $Q_{0}$, Eq. 39,
as indicated in the figure. For $Q_{0}\leq 3.$ the electrostatric potential
stays completely on the PB branch. For $Q_{0}>3.$ it stays on the PB branch
for $Q\leq 3.$ but then exhibits a jump at the transition point and continues
on the MS branch afterwards. The inset shows the details of the dependence for
$Q_{0}=5.$. We refer to this jump at the transition point as the ”Donnan
potential difference”. The value of $\tilde{z}$ at this jump, that can be read
off Fig. 3, also corresponds to the thickness of the nematic layer. In
addition for $Q_{0}>3.$, $\phi^{*}(Q)-\phi^{*}(Q_{0})=0$ at the surface only
for the MS branch while th ePB branch is suitably rescaled due to the form of
the boundary condition Eq. 39, on the MS and PB branches.

###  III.3 Dimensionless counterion density and electrostatic potential

The final dependence of the mean-field rod-like counterion density as well as
the electrostatic potential on the dimensionless distance from the bounding
surface, $\tilde{z}=q\sqrt{\frac{2}{\varepsilon}}z$, is obtained implicitly
from two equations. First by integrating the Poisson equation, Eq. 23, that
can be rewritten in the form

$\displaystyle\tilde{z}=\int^{Q_{0}}_{Q(z)}\frac{~{}dQ}{Q}\left(\frac{\partial
P(Q)}{\partial Q}\right),$ (40)

giving the dimensionless counterion density $Q$ as a function of $\tilde{z}$,
and then obtaining the dependence of the dimensionless potential on the
dimensionless density, $\phi^{*}(Q)$, that can be derived from Eq. 36, in the
form

$\displaystyle\phi^{*}(Q)-\phi^{*}(Q_{0})=\int^{Q_{0}}_{Q(z)}\frac{~{}dQ}{Q}\left(P(Q)\frac{\partial
P(Q)}{\partial Q}\right),$ (41)

so that $\phi^{*}(\tilde{z})$ is obtained form $\phi^{*}(Q(\tilde{z}))$. The
value of $Q_{0}$ is obtained from the boundary condition and Eq. 39 and yields
$Q_{0}=Q_{0}(\sigma)$. Clearly $Q_{0}$ can be either on the MS branch or the
PB branch, as is clear from Fig. 3, and consequently the functional dependence
on $\tilde{z}$ will depend on the $Q_{0}$, too. Note that the functional
dependence $Q_{0}=Q_{0}(\sigma)$ changes on the PB and the MS branch.

It is straightforward to obtain the limiting behavior for the dimensionless
density from Eqs. 40, 41 which leads to

$\displaystyle\lim_{Q\gg 1}Q(\tilde{z})=Q_{0}~{}e^{-\tilde{z}/a}\simeq
Q_{0}(1-\frac{\tilde{z}}{a})\qquad{\rm and}\qquad\lim_{Q\ll
1}Q(\tilde{z})=\frac{1}{(\tilde{z}+\tilde{z}_{0})^{2}},$ (42)

where $a=\lim_{Q\gg 1}\frac{\partial P(Q)}{\partial Q}=\sqrt{\textstyle 3/2}$
and $\tilde{z}_{0}^{2}=2/Q_{0}$. The former limit is valid for small
separations, while the latter is valid for large separations, coinciding with
the standard counterion-only Gouy-Chapman result, as it should.

Naively one would assume that the mean field electrostatic potential is
proportional to the log of the density, just as in the standard GC case.
However, the rod-like counterions also contain the orientational entropy as
part of the mean field energy, see Eq. 36, and thus the electrostatic
potential is given rather by

$\displaystyle-\phi^{*}(Q)=\log{Q/\gamma(Q)\big{)}}.$ (43)

While the spatial density profile is itself continuous with the derivative
being discontinuous at the I-N transition, the electrostatic potential is
discontinuous and displays a Donnan potential difference at the transition
point. This Donnan potential difference, $\phi^{*}_{D}$, can be obtained
straightforwardly from Eq. 43 as

$\displaystyle\phi^{*}_{D}=\log{J\big{(}\gamma_{I-N}\big{)}}=\log{J\big{(}{\textstyle\frac{9}{4}}s_{I-N}\big{)}},$
(44)

where $s_{I-N}$ is the jump of the orientational order parameter at the $I-N$
transition, and is therefore universal for all the electrostatic potential
curves. It can be viewed as a Lagrange multiplier that ensures the
electroneutrality of the system. The value of the Donnan potential difference
across the phase boundary can be read off the graph, Fig. 4, as $0.313$, which
equals exactly $\log{\big{(}J({\textstyle 2.53})\big{)}}$ according to Eq. 44.

##  IV Discussion and conclusions

While one can formulate the theory of inhomogeneous charged rod-like systems
on different levels of approximations, we were specifically motivated to
remain as close as possible to the Gouy-Chapman theory of point-like ions, the
reason being that the Poisson-Boltzmann mean-field framework presents the
foundation for the soft matter electrostatics and serves as a standard against
which the new developments are usually compared with. Of course this
necessarily implies also that the present theory shares at least the same
weaknesses as those well known for the Poisson-Boltzmann theory [47].

The main feature of the theory presented is the two coupled mean-field
equations which present generalizations of the standard Maier-Saupe and
Poisson-Boltzmann equations. Their solution for a single charged surface with
rod-like counterions leads to the existence of a nematic wetting layer, driven
by the interplay of nematic and electrostatic interactions between charged
rods and the bounding surface charges. The phase boundary at a finite distance
from the surface is in addition characterized by an electrostatic potential
jump corresponding to the nematic order parameter jump. While there are
obvious similarities with the standard nematic-isotropic transition and the
existence of nematic wetting, one needs to remember that in our case there is
no chemical potential driving the density changes as well as no order
parameter elastic energy driving the nematic wetting. The only driving force
is electrostatics.

Among the possible and obvious generalizations of the present theory we should
mention two explicitly. The first one is the tensorial nature of the
interaction potential

$\displaystyle u_{QQ}^{-1}({\bf x}-{\bf x}^{\prime})\longrightarrow
u_{iklm}^{-1}({\bf x}-{\bf x}^{\prime}),$ (45)

which would correspond to a more complicated liquid crystal elastic energy.
The second one is the non-locality of the nematic interactions which implies
the following form for the $u_{QQ}^{-1}({\bf x}-{\bf x}^{\prime})$ interaction

$\displaystyle u_{QQ}^{-1}({\bf x}-{\bf x}^{\prime})\longrightarrow
u_{QQ}^{-1}(0)\left(1+\xi^{2}\nabla^{2}\right)\delta({\bf x}-{\bf
x}^{\prime}),$ (46)

or the corresponding tensorial expression consistent with Eq. 46. Above $\xi$
is the nematic order correlation length. The non-local form of the interaction
potential in its turn leads to a non-local form of the Maier-Saupe equation,
Eq. 20, and one ends up with a system of two coupled non-linear differential
equations. This would in its turn require also the introduction of the surface
energy that would pin the surface value of the order parameter. The no doubt
complicated solutions would reduce to those studied above when the nematic
correlation length $\xi$ is much smaller then the electrostatic correlation
length, i.e., either the Gouy-Chapman length or the Debye length, depending on
the composition of the system.

The other possible and obvious generalization would be to include the higher
multipolar moments into the interaction energy, Eq. 5, such as the quadrupolar
electrostatic term. Formally this can be seen as leading to a modification in
Eq. 55 of the form

$\Phi_{ij}({\bf x})\longrightarrow\Phi_{ij}({\bf
x})+t\nabla_{i}\nabla_{j}\phi({\bf x}),$ (47)

where $t$ is the strength of the quadrupolar moment of the charged rod. This
generalization would treat the rod-rod electrostatic interactions more
accurately, allowing for the existence of the Odijk effect (preferred
perpendicular orientation of the rods) but would again imply a more
complicated form of the Poisson-Boltzmann equation.

A variation on the geometry of the model could be pursued for a system
confined between two charged surfaces with point-like co-ions. In this case
one can either expect a surface nematic wetting transition or indeed a
Fredericksz-type transition with a nematic phase between the surfaces and
isotropic layers vicinal to the surfaces. These variations in the geometry
setup would allow for interesting phenomena also in terms of the effective
electrostatic interactions between the bounding surfaces that would no doubt
deviate from the standard expectations.

##  V Acknowledgement

This work was funded by the Key project #12034019 of the National Science
Foundation of China. I would also like to acknowledge the support of the
School of physics, University of the Chinese Academy of Sciences, Beijing and
of the Institute of physics, Chinese Academy of Sciences, Beijing. Finally I
would like to thank D. Andelman, S. Buyukdagli and J. Everts for illuminating
and helpful discussions and in particular M. Bier for valuable comments on an
earlier version of this manuscript.

##  VI Appendix

Here we will give a short derivation of the field theoretic representation of
the partition function as a functional integral over the scalar and tensor
auxiliary fields based on the respective order parameters. The grand canonical
partition function for the system with Hamiltonian Eq. 6 is defined standardly
as

$\displaystyle{\Xi}[{\bf r}_{n},{\bf
n}_{n}]=\sum_{N^{+}}\sum_{N_{(-)}}\frac{\lambda_{(+)}^{N^{+}}\lambda_{(-)}^{N^{-}}}{N^{+}!N^{-}!}~{}\int\\!\\!\dots\\!\\!\int{\cal
D}[{\bf x}_{n}]{\cal D}[{\bf n}_{m}]~{}e^{-\beta{\cal H}[{\bf r}_{n},{\bf
n}_{n}]},$ (48)

where the integral over the orientational degrees of freedom, $[{\bf n}_{m}]$,
is only over the cationic species. Introducing the ”decomposition of unit” in
the form

$\displaystyle 1$ $\displaystyle\equiv$ $\displaystyle\int{\cal D}[{\cal
Q}_{ij}({\bf x})]~{}\Pi_{{\bf x}}\delta({\cal Q}_{ij}({\bf x})-\hat{\cal
Q}_{ij}({\bf x}))\times$ (49) $\displaystyle\int{\cal D}[\rho({\bf
x})]~{}\Pi_{{\bf x}}\delta(\rho({\bf x})-\hat{\rho}({\bf x})),$

together with the functional integral representation with auxiliary potentials
$\phi({\bf x})$ and $\Phi_{ij}({\bf x})$ for the functional delta functions,
one then remains with the following form of the partition function

$\displaystyle{\Xi}[{\bf r}_{n},{\bf n}_{n}]$ $\displaystyle=$
$\displaystyle\int{\cal D}[{\cal Q}_{ij}({\bf x})]{\cal D}[\Phi_{ij}({\bf
x})]\times$ (50) $\displaystyle~{}\int{\cal D}[\rho({\bf x})]{\cal
D}[\phi({\bf x})]~{}e^{-\beta\tilde{\cal H}[{\cal Q}_{ij}({\bf
x}),~{}\rho({\bf x})]}\times$
$\displaystyle\sum_{N^{+}}\sum_{N_{(-)}}\frac{\lambda_{(+)}^{N^{+}}\lambda_{(-)}^{N^{-}}}{N^{+}!N^{-}!}~{}\int{\cal
D}[{\bf x}_{n}]{\cal D}[{\bf n}_{m}]~{}e^{-\beta\tilde{\cal H}^{*}[{\bf
x}_{n},{\bf n}_{m}]},$

where we introduced the field coupling part, ${\cal H}[{\cal Q}_{ij}({\bf
x}),~{}\rho({\bf x})]$, that stems from the two decompositions of unit as well
as the interaction energy written with the collective coordinates, in the form

$\displaystyle-\beta\tilde{\cal H}[{\cal Q}_{ij}({\bf x}),~{}\rho({\bf x})]$
$\displaystyle=$ $\displaystyle-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf
x}d{\bf x}^{\prime}~{}{\cal Q}_{ij}({\bf x})u_{QQ}({\bf x}-{\bf
x}^{\prime}){\cal Q}_{ij}({\bf
x}^{\prime})-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}{\rho}({\bf x})u_{\rho\rho}({\bf x}-{\bf x}^{\prime}){\rho}({\bf
x}^{\prime})-$ (51) $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}-i\int_{V}d{\bf
x}~{}{\cal Q}_{ij}({\bf x})\Phi_{ij}({\bf x})-i\int_{V}d{\bf x}~{}\rho({\bf
x})\phi({\bf x}),$

while the configurational part in the internal space of the particles has the
form

$\displaystyle-\beta\tilde{\cal H}^{*}[{\bf r}_{n},{\bf n}_{n}]$
$\displaystyle=$ $\displaystyle i\int_{V}\\!\\!\\!d{\bf x}~{}\hat{\cal
Q}_{ij}({\bf x})\Phi_{ij}({\bf x})+i\int_{V}\\!\\!\\!d{\bf
x}~{}\hat{\rho}({\bf x})\phi({\bf
x})=i~{}\sum_{(+)}{\textstyle\frac{3}{2}}\left({\bf n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Phi_{ij}({\bf
x}_{i})+~{}iq\sum_{(+)}~{}\phi({\bf x}_{i})-i\sum_{(-)}~{}\phi({\bf x}_{i}),$
(52)

where we explicitly used the definitions of the microscopic scalar charge
density and the microscopic tensor nematic order parameter density, Eqs. 2 and
1. Since the ${\cal Q}_{ij}({\bf x})$ and $\rho({\bf x})$ functional integrals
in Eq. 50 are obviously Gaussian, these variables can be integrated out
explicitly, yielding a pure field theoretical representation of the original
configurational partition function Eq. 48 in terms of the tensor,
$\Phi_{ij}({\bf x})$, and scalar, $\phi({\bf x})$, fields with

$\displaystyle{\Xi}[\Phi_{ij}({\bf x});\phi({\bf x})]$ $\displaystyle=$
$\displaystyle\int\\!\\!{\cal D}[\Phi_{ij}({\bf x})]\int\\!\\!{\cal
D}[\phi({\bf x})]~{}e^{-\beta{\underline{\cal H}}[\Phi_{ij}({\bf x});\phi({\bf
x})]},$ (53)

where the effective field action is finally obtained in the form

$\displaystyle-\beta{\underline{\cal H}}[\Phi_{ij}({\bf x});~{}\phi({\bf x})]$
$\displaystyle=$
$\displaystyle~{}-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}{\Phi}_{ij}({\bf x})u_{QQ}^{-1}({\bf x}-{\bf
x}^{\prime}){\Phi}_{ij}({\bf
x}^{\prime})-{\textstyle\frac{1}{2}}\int\\!\\!\int_{V}d{\bf x}d{\bf
x}^{\prime}~{}{\phi}({\bf x})u_{\rho\rho}^{-1}({\bf x}-{\bf
x}^{\prime}){\phi}({\bf x}^{\prime})+$ (54)
$\displaystyle+~{}\lambda_{(+)}\int_{V}d{\bf x}~{}e^{i\phi({\bf x})+\ln{{\cal
P}(\Phi_{ij}({\bf x}))}}+\lambda_{(-)}\int_{V}d{\bf x}~{}e^{-i\phi({\bf x})}-$
$\displaystyle-{\textstyle\frac{1}{2}}{\ln}~{}{\rm Det}{\left(u_{QQ}({\bf
x}-{\bf x}^{\prime})\right)}-{\textstyle\frac{1}{2}}{\ln}~{}{\rm
Det}{\left(u_{\rho\rho}({\bf x}-{\bf x}^{\prime})\right)}.$

Here we introduced the orientational partition function of a single particle

$\displaystyle{\cal P}(\Phi_{ij}({\bf
x}))\equiv\Big{<}e^{i{\textstyle\frac{3}{2}}\left({\bf n}_{i}{\bf
n}_{j}-{\textstyle\frac{1}{3}}\delta_{ij}\right)\Phi_{ij}({\bf
x})}\Big{>}_{\Omega}.$ (55)

The orientational trace is defined in such a way that $<1>_{\Omega}=1$. The
log of this expression is actually the orientational entropy of the rods. Were
the field $\Phi_{ij}({\bf x})$ pure imaginary, the above distribution would
correspond to the Bingham distribution, and the $\Phi_{ij}({\bf x})$ field
would be playing the role of the Bingham field.

The formal identity of the configurational, Eq. 48, and auxiliary field
representations, Eq. 53, can be recapitulated as

$\displaystyle{\Xi}[{\bf r}_{n},{\bf n}_{n}]={\Xi}[\Phi_{ij}({\bf
x});\phi({\bf x})].$ (56)

The steps leading to this identity are analogous to the case of Edwards
transform for the Coulomb fluid partition function [47], except for the
orientational Lebwohl-Lescher part that leads to a tensor order parameter and
a tensor auxiliary field.

Notably, the two fluctuational ${\rm Tr}\log$ expressions at the end of the
above equation pertain to the Casimir-type fluctuation terms [47] and would be
combined with the fluctuational determinant of the second order expansion of
the above field theory. We will not delve into these details here.

## References

  * Bernal and Fankuchen [1941] J. Bernal and I. Fankuchen, The isotropic-nematic transition in charged liquid crystals, J. Gen. Physiol. 25, 147 (1941).
  * Šiber  _et al._ [2012] A. Šiber, A. L. Božič, and R. Podgornik, Energies and pressures in viruses: contribution of nonspecific electrostatic interactions, Phys. Chem. Chem. Phys. 14, 3746 (2012).
  * Busseron  _et al._ [2013] E. Busseron, Y. Ruff, E. Moulin, and N. Giuseppone, Supramolecular self-assemblies as functional nanomaterials, Nanoscale 5, 7098 (2013).
  * Culver  _et al._ [2015] J. N. Culver, A. D. Brown, F. Zang, M. Gnerlich, K. Gerasopoulos, and R. Ghodssi, Plant virus directed fabrication of nanoscale materials and devices, Virology 479-480, 200 (2015), 60th Anniversary Issue.
  * Rong  _et al._ [2011] J. Rong, Z. Niu, L. A. Lee, and Q. Wang, Self-assembly of viral particles, Current Opinion in Colloid & Interface Science 16, 441 (2011).
  * Yang  _et al._ [2013] S. H. Yang, W.-J. Chung, S. McFarland, and S.-W. Lee, Assembly of bacteriophage into functional materials, The Chemical Record 13, 43 (2013), https://onlinelibrary.wiley.com/doi/ pdf/ 10.1002/tcr.201200012 .
  * Tang  _et al._ [1996] J. X. Tang, S. Wong, P. T. Tran, and P. A. Janmey, Counterion induced bundle formation of rodlike polyelectrolytes, Ber. Bunsenges. Phys. Chem. 100, 796 (1996), https://onlinelibrary.wiley.com/ doi/abs/ 10.1002/bbpc.19961000620 .
  * Bellini  _et al._ [2012] T. Bellini, R. Cerbino, and G. Zanchetta, Dna-based soft phases, Top Curr Chem 318, 225 (2012).
  * Sanders  _et al._ [2007] L. K. Sanders, W. Xian, C. Guáqueta, M. J. Strohman, C. R. Vrasich, E. Luijten, and G. C. L. Wong, Control of electrostatic interactions between f-actin and genetically modified lysozyme in aqueous media, Proc. Natl. Acad. Sci. USA 104, 15994 (2007), https://www.pnas.org/content/104/41/15994.full.pdf .
  * Needleman  _et al._ [2004] D. J. Needleman, M. A. Ojeda-Lopez, U. Raviv, H. P. Miller, L. Wilson, and C. R. Safinya, Higher-order assembly of microtubules by counterions: From hexagonal bundles to living necklaces, Proceedings of the National Academy of Sciences 101, 16099 (2004), https://www.pnas.org/content/101/46/16099.full.pdf .
  * Teif and Bohinc [2011] V. B. Teif and K. Bohinc, Condensed dna: condensing the concepts, Progress in biophysics and molecular biology 105, 208 (2011).
  * Kanduč  _et al._ [2009] M. Kanduč, A. Naji, Y. S. Jho, P. A. Pincus, and R. Podgornik, The role of multipoles in counterion-mediated interactions between charged surfaces: strong and weak coupling, J. Phys.: Condens. Matter 21, 424103 (2009).
  * Bohinc and Lue [2015] K. Bohinc and L. Lue, The interaction between like-charged nanoparticles mediated by rod-like ions, Journal of Nanoscience and Nanotechnology 15, 3468 (2015).
  * Kanduč  _et al._ [2018] M. Kanduč, A. Naji, Y. S. Jho, P. A. Pincus, and R. Podgornik, Charged nanorods at heterogeneously charged surfaces, J. Phys.: Condens. Matter 149, 134702 (2018).
  * Alvarez Fernandez and Kouwer [2016] A. Alvarez Fernandez and P. Kouwer, Key developments in ionic liquid crystals, Int J Mol Sci. 17, 731 (2016).
  * Shi and Wang [2016] R. Shi and Y. Wang, Dual ionic and organic nature of ionic liquids., Scientific Reports 6, 19644 (2016).
  * Goossens  _et al._ [2016] K. Goossens, K. Lava, C. W. Bielawski, and K. Binnemans, Ionic liquid crystals: Versatile materials, Chem. Rev. 116, 4643 (2016).
  * Deutsch and Goldenfeld [1982a] J. M. Deutsch and N. D. Goldenfeld, Ordering in charged rod fluids, J. Phys. A: Math. Gen. 15, L71 (1982a).
  * Deutsch and Goldenfeld [1982b] J. Deutsch and N. Goldenfeld, The isotropic-nematic transition in charged liquid crystals, J. Phys 43, 651 (1982b).
  * Odijk [1986] T. Odijk, Theory of lyotropic polymer liquid crystals, Macromolecules 19, 2313 (1986).
  * der Schoot and Odijk [1990] P. V. der Schoot and T. Odijk, Structure factor of a semidilute solution of rodlike macromolecules, Macromolecules 23, 4181 (1990).
  * Vroege and Lekkerkerker [1992] G. J. Vroege and H. N. W. Lekkerkerker, Phase transitions in lyotropic colloidal and polymer liquid crystals, Macromolecules 55, 1241 (1992).
  * Wensink and Trizac [2014] H. H. Wensink and E. Trizac, Generalized onsager theory for strongly anisometric patchy colloids, J. Chem. Phys. 140, 024901 (2014).
  * Constantinescu and Jho [2014] A. Constantinescu and Y. Jho, Nematic phase emergence in solutions of similarly charged rodlike polyelectrolytes, Journal of the Physical Society of Japan 83, 014002 (2014).
  * Bruinsma [2001] R. Bruinsma, Liquid crystals of polyelectrolyte networks, Phys. Rev. E 63, 061705 (2001).
  * Eggen  _et al._ [2009] E. Eggen, M. Dijkstra, and R. van Roij, Effective shape and phase behavior of short charged rods, Phys. Rev. E 79, 041401 (2009).
  * Potemkin and Khokhlov [2004] I. I. Potemkin and A. R. Khokhlov, Nematic ordering in dilute solutions of rodlike polyelectrolytes, J. Chem. Phys. 120, 10848 (2004), https://doi.org/10.1063/1.173939 .
  * Lue [2006] L. Lue, A variational field theory for solutions of charged, rigid particles, Fluid Phase Equilibria 241, 236 (2006), a Festschrift in Honor of John M. Prausnitz.
  * Carri and Muthukumar [1999] G. Carri and M. Muthukumar, Attractive interactions and phase transitions in solutions of similarly charged rod-like polyelectrolytes, Journal of Chemical Physics 111, 1765 (1999).
  * Levy  _et al._ [2013] A. Levy, D. Andelman, and H. Orland, Dipolar poisson-boltzmann approach to ionic solutions: A mean field and loop expansion analysis, J. Chem. Phys. 139, 164909 (2013).
  * Buyukdagli and Ala-Nissila [2013] S. Buyukdagli and T. Ala-Nissila, Microscopic formulation of nonlocal electrostatics in polar liquids embedding polarizable ions, Phys. Rev. E 87, 063201 (2013).
  * Frydel [2016] D. Frydel,  _Advances in Chemical Physics_ , edited by S. A. Rice and A. R. Dinner (Wiley Pub., 2016) pp. 209–260.
  * Ravnik and Everts [2020] M. Ravnik and J. Everts, Topological-defect-induced surface charge heterogeneities in nematic electrolytes, Phys. Rev. Letts. 125, 037801 (2020).
  * Bier [2007] M. Bier,  _Interfaces in fluids of charged platelike colloids_ (Max-Planck-Institut fuer Metallforschung, Stuttgart (Disertation), 2007).
  * Kondrat  _et al._ [2010] S. Kondrat, M. Bier, and L. Harnau, Phase behavior of ionic liquid crystals, J. Chem. Phys. 132, 184901 (2010).
  * Bier  _et al._ [2005] M. Bier, L. Harnau, and S. Dietrich, Free isotropic-nematic interfaces in fluids of charged platelike colloids, J. Chem. Phys. 123, 114906 (2005).
  * Bier  _et al._ [2006] M. Bier, L. Harnau, and S. Dietrich, Surface properties of fluids of charged platelike colloids, J. Chem. Phys. 125, 184704 (2006).
  * Abrashkin  _et al._ [2007] A. Abrashkin, D. Andelman, and H. Orland, Dipolar poisson-boltzmann equation: Ions and dipoles close to charge interfaces, Phys. Rev. Lett. 99, 077801 (2007).
  * Azuara  _et al._ [2008] C. Azuara, H. Orland, M. Bon, P. Koehl, and M. Delarue, Incorporating dipolar solvents with variable density in poisson-boltzmann electrostatics, Biophysical Journal 95, 5587 (2008).
  * Buyukdagli and Blossey [2014] S. Buyukdagli and R. Blossey, Dipolar correlations in structured solvents under nanoconfinement, The Journal of Chemical Physics 140, 234903 (2014), https://doi.org/10.1063/1.4881604 .
  * Slavchov and Ivanov [2014] R. I. Slavchov and T. I. Ivanov, Quadrupole terms in the maxwell equations: Born energy, partial molar volume, and entropy of ions, J. Chem. Phys. 140, 074503 (2014).
  * Kim  _et al._ [2008] Y. W. Kim, J. Yi, and P. A. Pincus, Attractions between like-charged surfaces with dumbbell-shaped counterions, Phys. Rev. Lett. 101, 208305 (2008).
  * Frenkel [1991] D. Frenkel,  _Liquides, Cristallisation et Transition Vitreuse/ Liquids, Freezing and Glass Transition_ , edited by D. L. J. P. Hansen and J. Zinn-Justin (Les Houches, Session LI, 1991) pp. 693–759.
  * Lebwohl and Lasher [1972] P. A. Lebwohl and G. Lasher, Nematic-liquid-crystal order—a monte carlo calculation, Phys. Rev. A 6, 426 (1972).
  * Lu [2017] B.-S. Lu, Multiscale approach to nematic liquid crystals via statistical field theory, Phys. Rev. E 96, 022709 (2017).
  * Doi [1981] M. Doi, Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid crystalline phases, Journal of Polymer Science: Polymer Physics Edition 19, 229 (1981).
  * Naji  _et al._ [2013] A. Naji, M. Kanduč, J. Forsman, and R. Podgornik, Perspective: Coulomb fluids -weak coupling, strong coupling, in between and beyond, J. Chem. Phys. 139, 150901 (2013).
  * Braun  _et al._ [1996] F. N. Braun, T. J. Sluckin, and E. Velasco, Director distortion in a nematic wetting layer, Journal of Physics: Condensed Matter 8, 2741 (1996).
  * Podgornik [1992] R. Podgornik, Self-consistent-field theory for confined polyelectrolyte chains, J. Phys. Chem. 96, 884 (1992).
  * Borukhov  _et al._ [1995] I. Borukhov, D. Andelman, and H. Orland, Polyelectrolyte solutions between charged surfaces, Europhysics Letters (EPL) 32, 499 (1995).
  * Kleinert [2018] H. Kleinert,  _Collective Classical and Quantum Fields_ (World Scientific, Singapore, 2018).
  * Kumari  _et al._ [2020] S. Kumari, F. Ye, and R. Podgornik, Ordering of adsorbed rigid rods mediated by the boussinesq interaction on a soft substrate, The Journal of Chemical Physics 153, 144905 (2020), https://doi.org/10.1063/5.0022556 .
  * Maggs and Podgornik [2016] A. C. Maggs and R. Podgornik, General theory of asymmetric steric interactions in electrostatic double layers, Soft Matter 12, 1219 (2016).
  * Markovich  _et al._ [2020] T. Markovich, D. Andelman, and R. Podgornik,  _Handbook of Lipid Membranes_ , edited by C. Safinya and J. Raedler (Taylor and Francis, 2020) Chap. Charged Membranes: Poisson-Boltzmann theory, DLVO paradigm and beyond.
  * Pandit and Wortis [1982] R. Pandit and M. Wortis, Surfaces and interfaces of lattice models: Mean-field theory as an area-preserving map, Phys. Rev. B 25, 3226 (1982).
  * Gavish and Promislow [2016] N. Gavish and K. Promislow, On the structure of generalized poisson-boltzmann equations, Euro. Journal of Applied Mathematics 27, 667 (2016).
  * Sheng [1976] P. Sheng, Phase transition in surface-aligned nematic films, Phys. Rev. Lett. 37, 1059 (1976).
  * Sluckin and Poniewierski [1985] T. J. Sluckin and A. Poniewierski, Novel surface phase transition in nematic liquid crystals: Wetting and the kosterlitz-thouless transition, Phys. Rev. Lett. 55, 2907 (1985).