Theory of inhomogeneous calamitic Coulomb fluids

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 {\sl 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 calamitic counterion-only system in proximity of a charged planar wall. The nematic order parameter profile, the counterion density profile and the electrostatic potential profile are interpreted within the framework of a nematic-isotropic wetting phase 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 pH 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 calamitic shape 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 abundand. 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.
FIG. 1. A schematic presentation of the system with calamitic cations and simple anions in the bulk (left). An inhomogeneous system of calamitic cationsclose 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 calamitic 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 [13,[42][43][44] 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 [34,36,37], 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 calamitic 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 calamitic cations and point-like anions with microscopic interactions of 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 calamitic 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. Our theory presents a natural generalization of the Gouy-Chapman theory to the case of nematogenic calamitic counterions.

A. Order parameters
Let us consider a system composed of monovalent anions (charge −e) and polyvalent rod-like cations (charge +qe) with chemical potential µ, 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. [45]) with n(x n ) andq ij (x n ) the director and the orientational order parameter of the n − th rod-like cation, together with the microscopic charge densityρ so that the total charge density isρ(x) =ρ (+) (x) +ρ (−) (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 and where the corresponding thermodynamic densities of the rod cations and anions are c + (x), c − (x) and the nematic scalar order parameter is S(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 [46]. This constraint was relaxed in a recent analytical formulation [47]. By assumption then where n, m run over all the particles in the system and n n , n m are the unit directors of the n-th and m-th particles.
The interaction strengths are all expressed in thermal units. The angular part of the Onsager's nematic interaction between two spherocylinders of length L and diameter 2D is 2L 2 D sin(n, n ′ ) [22], where n, n ′ are the unit direction vectors of the two particles is also related to the above Ansatz.
In fact it was shown by Doi [48] that sin(n, n ′ ) = π 4 1 − 15 where the next order term corresponds to 4th rank irreducible tensors. The Onsager interactions are thus -at least to the lowest order -equivalent to the first term in Eq. 5, where the nematic interaction has a delta function range. 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 where we further assume that the orientational interaction is short-ranged and that the scalar part is due to the longrange electrostatic 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.

B. Grand canonical partition function
The grand canonical partition function for this two component system is defined standardly as Ξ[r n , n n ] = where the integral over the orientational degrees of freedom, [n m ], is only over the cationic species. Introducing the "decomposition of unit" in the form together with the functional integral representation with auxiliary potentials φ(x) and Φ ij (x) for the functional delta functions, one then remains with the following form of the partition function where we introduced the field coupling part, H[Q ij (x), ρ(x)], that stems from the two decompositions of unit as well as the interaction energy written with the collective coordinates, in the form while the configurational part in the internal space of the particles has the form where we explicitly used the definitions of the microscopic scalar charge density and the microscopic tensor nematic order parameter density, Eqs. 2 and 1.
H * [r n , n n ] is diagonal in the internal space of the particles so that after summing separately and explicitly over both N + and N − , with Here we introduced the orientational partition function of a single particle The orientational trace is defined in such a way that < 1 > Ω = 1. The log of this expression is actually the orientational entropy of the rods. Were the field Φ ij (x) pure imaginary, the above distribution would correspond to the Bingham distribution, and the Φ ij (x) field would be playing the role of the Bingham field. Since the Q ij (x) and ρ(x) functional integrals in Eq. 11 are obviously Gaussian, these variables can be integrated out explicitly, yielding a pure field theoretical representation of the original configurational partition function Eq. 8 in terms of the tensor, Φ ij (x), and scalar, φ(x), fields with where the effective field action is finally obtained in a highly non-linear form The formal identity of the configurational, Eq. 8, and auxiliary field representations, Eq. 17, can be recapitulated as The steps leading to this identity are analogous to the case of Edwards transform for the Coulomb fluid partition function [49], except for the orientational Lebwohl-Lescher part that leads to a tensor order parameter and a tensor auxiliary field. The equivalence is formal and further steps are needed to actually calculate the functional integrals of the very non-linear field action. We will eventually take recourse to the saddle-point approximation which is, in the case of the Coulomb fluid, equivalent to the mean-field Poisson-Boltzmann approximation.
Notably, the two fluctuational T rLog expressions at the end of the above equation pertain to the Casimir-type fluctuation terms [49] 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 and will omit them in what follows.

A. Order parameters and auxiliary fields
We now investigate the consequences of the field-theoretical formulation of the partition function for the microscopic model under consideration. First we need to find a physical meaning of the two formal auxiliary fields that were introduced above. This proceeds similarly to the field-theoretical formulation of the Coulomb fluid partition function.
One notes that in the case of an added ghost source charge density, , and/or an added ghost quadrupolar orientational order density, Q 0 ij (x), Eq. 7 acquires two additional terms Repeating now the steps leading to Eq. 18 one then obtains the effective field action for the two-component field theory with ghost sources as It then follows straighforwardly that and similarly Rearranging these equations obtained from Eq. 21, one can derive a set of two other equations stemming from the derivatives of Eq. 20, thus ending up with two general identities and analogously that we will use later in their saddle-point variant. The latter identity establishes the thermodynamic average i φ(x) as the average of the electrostatic potential, and the former establishes the average i Φ ij (x) as the average of the nematic potential. One should note that these are general expressions valid exactly without any approximations. Clearly, the auxiliary fields and the order parameters are dual to one another.

B. Thermodynamic relations
The fields and order parameters need to satisfy several thermodynamic relations. There is first the thermodynamic identity involving the fugacity and the average number of molecules in the system and Another general relation is derived from the invariance of the functional integral with respect to the linear transformation of the fluctuating fields. This immediately yields two relations and where the averages stand for Distinguishing between the normalized .. .
and unnormalized orientational averages with respect to the distribution Eq. 16, one derives from Eqs. 28 and 29 that where we took into account the definition of the double brackets Eq. 31. Furthermore by analogy Note the difference between the thermodynamic averages . . . and the orientational averages . . .

Ω
. 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.

A. Maier-Saupe and Poisson-Boltzmann equations
Since the field action Eq. 18 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 fields.
At the saddle-point the fields are pure imaginary, so that one can transform The thermodynamic averages . . . are given at the value of the mean-field, and the two self-consistent field equations, Eqs. 24, 25, are then reduced to as well as with the orientational partition function Eq. 16 as In what follows we will assume a short range attractive orientational potential and a Coulombic positional potential. This implies and where ε ≡ ǫ 0 ǫ, with ǫ the dielectric permittivity of the solvent and ǫ 0 the permittivity of space.
With this in mind we then derive the tensorial part in the form of a modified Meier-Saupe equation 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. To describe orientational relaxation effects one needs also the elastic deformation energy which would stem from the expansion of the orientation interaction potential w.r.t. the gradient of the tensorial order. The scalar part of the mean-field equations can be written in the form of the Poisson-Boltzmann equation as where we introduced the mean-field cationic and anionic densities as The two mean-field equations Eqs. 40 and 41 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 [50,51]. Consequently any other degree of freedom would introduce its own mean-field equation.
Inserting the mean-field Ansatz into the free energy Eq. 18 we obtain the general form of the mean-field free energy F [Φ * ij (x); φ * (x)] as a functional of the mean-field nematic potential, Φ * ij (x), and the mean-field electrostatic potential, φ * (x) as The mean field free energy can be derived also from the Onsager theory as was done by Doi [48]. We recapitulate this alternative derivation as it applies to our model in the Appendix. The mean-field theory for calamitic cations thus decouples into the MS equation for simple liquid crystals and the PB equation for simple Coulomb fluids, except that formally the electrostatic potential of the calamitic cations is transformed to qφ * (x) −→ qφ * (x) − ln P(−iΦ * ij (x)). 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.

B. Solution of the mean-field equations
There are of course countless cases that one can dwell on in order to apply the general theory, and therefore 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 so that the only spatial dependence is in the direction of the surface normal, which we choose to coincide with the z axis. This tallies with the paradigmatic case in the solution of the ordinary PB equation for a counterion-only system.
The scalar mean-field equation in this case is the modified Poisson-Boltzmann equation that is obtained as is the second derivative of the mean potential with respect to z and c (+) (φ * (z)) is given by Eq. 42. Comparing the above expression with the standard PB equation for a point-like counterion-only case, the difference is manifest in the sense that the mean-field single ion free energy includes also the orientational entropy of the rods, equal to ln P(−iΦ * ij (x)). We will solve this equation later. First we analyse the tensor mean-field equation, i.e., the modified Maier-Saupe equation. Since the tensor auxiliary field and the tensor order parameter are dual, we can assume a non-vanishing orientational order characterized with the directorn so that and s(z), proportional but not equal to the nematic order parameter, then simply represents the absolute value of Φ * ij (z). In infinite space without any anchoring fields,n is arbitrary but would be set by e.g. Rapini-Papoular surface free energy terms, which we did not invoke explicitly since the Maier-Saupe equation is local and thus valid for any local orientation.
One can recall, Eq. 3 and Eq. 45, which means that wherefrom we derive the exact connection between the nematic order parameter S and the parameter s introduced in Eq. 45 as where c + is given by Eq. 42. In the above formula s is a coefficient of Φ * ij (z) in Eq. 45, while S is an orientational order parameter, so that it is strictly limited to the interval 0 ≤ S ≤ 1.
Then we see that the exponent in Eq. 37 equals Multiplying both sides of Eq.40 by n inj − 1 3 δ ij this turns out to yield 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. Developing further, this eventually yields where γ = 9 4 s and from the definition Eq. 31 it follows that which can be indeed written as the logarithmic derivative in Eq. 51. Above the J(γ) function has been defined as with D(x) as the Dawson's integral [52]. The Maier-Saupe expression can then be cast definitively as where γ = 9 4 s and the nematic order parameter is extracted from Eq. 48 as S = γ/( 9 4 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. 38 we only considered local nematic interactions. A similar path can be pursued also in 2D except that instead of the Dawson integral one ends up with a modified Bessel function integral [53].
Eq. 54 corresponds exactly to the Kleinert formulation [52] of the Maier-Saupe theory, if one takes the coupling strength of nematic interactions equal to A 0 = 3 2 u QQ (0)c + , i.e., to the local strength of the nematic potential multiplied by the local density. Obviously, the larger the local density of the rods, the smaller needs to be the orientational interaction driving the isotropic nematic transition.
The solution of the Maier-Saupe equation for 0 ≤ S ≤ 1 exhibits a first order transition at a critical value of c (+) , where the order parameter makes a discontinuous jump from an isotropic phase to a nematic phase. More about this solution below.

C. First integral and the phase portrait analysis
Introducing now the generalized van't Hoff osmotic pressure as it is possible to write the mean-field equations in a "Lagrangian" form (see Ref. [54]) Proceeding now as in the case of the first integral of the Gouy-Chapman theory [55], by multiplying the first equation by φ * ′ and the second one by s ′ (z) and then summing them up, we obviously remain with which can be cast into the form of the first integral that generalizes the standard Poisson-Boltzmann result 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. 56 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 [56]. The phase portrait method allows for an easy and physical visualization of the types of solutions of the mean-field theories and has been successfully applied also to the case of the Poisson-Boltzmann-type theories [57].
One can now obtain the full implicit form of the solution of the mean-field equations by introducing two new variables where c ≡ c + is the counterion density, Eq. 55. Q can thus be interpreted as the dimensionless density of the calamitic counterions. It follows from the first integral Eq. 58 that 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 γ(Q) is a solution of Numerical solutions of the above equations are presented in Fig. 2 in the form of the dependence of the nematic order parameter S(Q) or equivalently the parameter s on the dimensionless calamitic counterion density: s = s(Q), with the jump from zero to 0.604 at the critical value Q = 11.15, obviously corresponding to a first order isotropic-nematic transition. Fig. 2 also shows the phase portrait of the system P = P (Q), Eq. 59. The dotted line in the inset corresponds to the pure PB branch (isotropic branch) while the solid line MS branch (nematic branch) separates from the PB branch at the nematic transition point.
The surface charge density at the bounding surface sets the boundary value to which follows directly from the Gauss boundary condition for the electrostatic field at a surface with surface charge density σ, i.e. εφ * ′ (z = 0) = σ. The solution then amounts to the P (Q) curve starting at P 0 and then moving along the solution line to P = 0. For P 0 large enough (see the P 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 an ordered wetting phase close to the surface, with a phase boundary at the transition point corresponding to a finite thickness of the wetting layer. After that the solution migrates to the PB branch of a disordered system that it follows 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 P B intercept in Fig. 2) there is no 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.
The final dependence of the mean-field calamitic counterion density as well as the electrostatic potential on the separation from the bounding surface, z, is obtained by integrating the first integral, Eq. 58, as on the MS branches (from left to right, full curves). As c0 increases (not shown) the curves are systematically displaced towards larger values of the z coordinate, viz., the I-N transition is displaced away from the surface. At the I-N transition, indicated by the gray vertical lines, there is a discontinuous jump corresponding to the Donnan potential difference, qϕ * D ≃ 0.99, corresponding to the phase boundary in the system. Note that the transition lines seem to be skewed and not vertical, but this is purely an artifact of the discrete sampling of the curve.
where φ * 0 is the value at the charged surface. This expression can be modified by taking into account the connection between Q and φ * (z) from Eq. 59, leading to The value of Q 0 is obtained from the boundary condition and Eq. 62 in the form which yields Q 0 = Q 0 (σ, u QQ (0)). 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 z will depend on the Q 0 , too. In the case of an isotropic solution, the equation Eq. 64 yields correctly the standard counterion-only Gouy-Chapman result, as it should. Evaluating now the integral in Eq. 64 numerically we obtain the general solution to the problem, including the migration from the PB to the MS branch. First we make the substitution Q = u QQ (0)c, where c is the calamitic counterion density. Then with 2q 2 /ε = 8πℓ B λ + and ℓ B the Bjerrum length. The dimensionless distance is thus 8πℓ B λ + z. c 0 is the density of the counterions at the bounding charged wall at z = 0. c 0 and u QQ (0), implicit in γ(c), are the two parameters of the density profile. Inverting the dependence z = z(c) one finally obtains c(z), the dependence of the calamitic counterion density on the separation from the charged wall. One notices, see Fig. 3, that the c(z) dependence is continuous but displays a discontinuity in the derivative at the isotropic -nematic transition. Nevertheless we will see that the underlying electrostatic potential does display a discontinuity at the wetting phase boundary.
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 calamitic counterions also contain the orientational entropy as part of the mean field energy, see Eq. 59, and thus the electrostatic potential is given rather by −qφ * (z) = log c(z) λ (+) J γ(z) .
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. The Donnan potential difference, φ * D , can be obtained straightforwardly from Eq. 67 as where s I−N is the jump of the orientational order parameter at the I − N transition, P A case in Fig. 2, and is therefore universal for all the electrostatic potential curves. The value of the Donnan potential difference across the phase boundary can be read off the graph, Fig. 4, as 0.99, which equals exactly log J 9 4 × 0.604 according to Eq. 68. On the other hand, should the Q 0 that corresponds to the boundary condition Eq. 62 fall below the critical value of 11.1, P B case in Fig. 2, the solution for the electrostatic potential remains on the PB branch for all values of z, displaying no discontinuities.

V. DISCUSSION AND CONCLUSIONS
While one can formulate the theory of inhomogeneous charged calamitic 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.
The main feature of the theory presented is the two mean-field equations which present generalizations of the standard Maier-Saupe and Poisson-Boltzmann equations. Their solution for a single charged surface with calamitic 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 a Donnan potential difference corresponding to the nematic order parameter jump.
Among the possible and obvious generalizations of the present theory we should mention two explicitly. The first one is the non-locality of the nematic interactions which implies the following form for the u −1 where ξ 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. 40, and one ends up with a system of two coupled non-linear differential equations. The no doubt complicated solutions would reduce to those studied above when the nematic correlation length ξ 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. 16 of the form where t is the strength of the quadrupolar moment of the charged rod. This generalization would treat the rodrod 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 of the type that was derived in Refs. [12,31,40]. In addition, one can also use the present theory to calculate the effective macroscopic dielectric function in the same way as was already done by Andelman et al. [38,58,59] in a different context. This can be pursued either from writing the Poisson-Boltzmann equation in the standard form, but with a modified dielectric constant, or by analyzing the second variational derivative of the free energy with respect to the electrostatic fields. Both methodologies yield an explicit expression for the dielectric constant, either local or a non-local one [60].
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.