Theory of Inhomogeneous Rod-like Coulomb Fluids

: A ﬁeld 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 ﬁeld 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 ﬁnally obtained for a rod-like counterion-only system in proximity to a charged planar wall. The nematic order parameter proﬁle, the counterion density proﬁle and the electrostatic potential proﬁle are interpreted within the framework of a nematic–isotropic wetting phase with a Donnan potential difference.


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 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 the 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 also incorporates 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.
The problem of an inhomogeneous system with multipolar anisotropic charges has also been addressed on various levels of modeling and approximations. The case of a system bounded by charged wall(s) with mobile dipolar charges has been analyzed with the field theoretical approach [38][39][40], while quadrupolar charges [12,41] and finite size dumbbell charges [13,42] 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 semi-infinite isotropic or nematic bulk system in contact with a charged hard wall exhibiting the 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.

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 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 also interact via short range nematic interactions, described by Lebwohl-Lasher interaction potential.
The respective microscopic order parameters are the nematic order parameter density, defined as (see, e.g., [43] with n(x n ) as the director 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 [44]. This constraint was relaxed in a recent analytical formulation [45]. 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 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 Equation (5), where the nematic interaction has a delta function range [46]. The theory presented here, therefore, also incorporates the standard Onsager results to the lowest order. The interaction energy Equation (5) can be clearly recast as where we further assume that the scalar part u ρρ (x − x ) 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.

Field Theoretical Representation and Thermodynamic Relations
In the Appendix A we derive the field theoretical representation of the grand canonical partition function, Ξ, in terms of the tensor, Φ ij (x), and scalar, φ(x), auxiliary fields with with the effective field action as where λ (±) are the fugacities of the two ionic species and P (Φ ij (x)) is the orientational partition function of a single particle, see Equation (A8), The orientational trace is defined in such a way that < 1 > Ω = 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 where the averages stand for These two relations, Equation (10), furthermore yield the following exact connections between the average values of the order parameters and the auxiliary fields where we took into account the definition of the double brackets Equation (14). Furthermore by analogy Note the difference between the orientational averages ... .

Ω
, is with respect to the distribution Equation (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.

Saddle-Point Approximation
Since the field action Equation (A7) 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 purely 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, Equations (12) and (13), are then reduced to with the orientational partition function P (−iΦ * ij (x)) defined in Equation (9) taken at the imaginary value of the tensor auxiliary field, and with obvious definitions of the cation and anion densities, ). 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. 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- 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 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 The two mean-field equations, Equations (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 Equation (A7), we obtain the general form of the mean-field free energy F [Φ * ij (x); φ * (x)] as a functional of the meanfield nematic potential, Φ * ij (x), and the mean-field electrostatic potential, φ * (x) as 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 ). 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.

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.

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 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 counteriononly case, the difference is manifest in the presence of the orientational entropy of the rods, equal to ln P (−iΦ * ij (z)). 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 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 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 Equations (3) and (24), which means that whereby we derive the exact connection between the nematic order parameter S, strictly limited to the interval 0 ≤ S ≤ 1, and the parameter s, introduced in Equation (24) as By multiplying both sides of Equation (20) by n inj − 1 3 δ ij , and tracing over the indices, we finally obtain 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 where, obviously γ = 9 4 s, and from Equation (14) it furthermore follows that the orientational partition function is obtained as with D(x) as the Dawson's integral [51]. The Maier-Saupe expression can then be cast definitively as where the nematic order parameter is then extracted from Equation (27) as S = γ/( 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]. Equation (31) corresponds exactly to the Kleinert formulation [51] of the Maier-Saupe theory, if one takes his coupling strength of nematic interactions equal to A 0 = 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 Equation (18), we only considered local nematic interactions.
The solution of the Maier-Saupe equation for 0 ≤ S ≤ 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.

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 [53]) Proceeding now, as in the case of the first integral of the Gouy-Chapman theory [54], by multiplying the first equation by φ * and the second one by s (z) and then summing them up, we 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" Equation (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 also been successfully applied to the case of the generalized Poisson-Boltzmann theory [56].
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, Equation (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 Equation (35) 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), Equation (31), becomes a solution of Numerical solutions of the above equations are presented in Figure 2 in the form of dependence of the nematic order parameter S(Q) on the dimensionless rod-like counterion density, with the jump from zero to 0.313 at the critical value Q c = 3.0, obviously corresponding to a first order isotropic-nematic transition. The phase portrait of the system, P = P(Q), in Figure 2 presents the dependence of the dimensionless electrostatic field P on the dimensionless density, showing the pure PB branch, and 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 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, Q 0 = Q 0 (σ), 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 Figure 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 Figure 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 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 standard 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 [48,57,58]. The wetting layer and its properties in our case are related but do not coincide with 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.

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,z = q 2 ε u QQ (0) z, is obtained implicitly from two equations. The first one follows by integrating the Poisson equation, Equation (23), as giving the dimensionless counterion density Q as a function of dimensionless distancez. The second one is obtained from the dependence of the dimensionless potential on the dimensionless density, φ * (Q), that can be derived from Equation (36), in the form so that φ * (z) is obtained from φ * (Q(z)). AS already stated, the value of Q 0 is obtained from the boundary condition and Equation (39), yielding Q 0 = Q 0 (σ). Clearly Q 0 can be either on the MS branch or the PB branch, as is clear from Figure 3, and consequently the functional dependence onz will depend on the Q 0 , too. Note that the functional dependence Q 0 = Q 0 (σ) changes on the PB and the MS branch.   (40) as a function of the dimensionless separationz from the bounding surface for different values of Q 0 , as indicated on the figure. For Q 0 ≤ 3, the dimensionless density dependence onz coincides with the PB branch. For Q 0 > 3, it stays on the PB branch for Q ≤ 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 (z = 0) = Q 0 at the surface, but to a rescaled value, such that Q MS (z = 0) = Q 0 , the difference stemming from the form of the boundary condition Equation (39) on the MS and PB branches.
It is straightforward to obtain the limiting behavior for the dimensionless density from Equations (40) and (41), which leads to where a = lim Q 1 ∂P(Q) ∂Q = √ 3/2 andz 2 0 = 2/Q 0 . The former limit is valid for small distances, while the latter is valid for large distances, coinciding with the standard counteriononly 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 Equation (36)and thus the electrostatic potential is, rather, given by 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, φ * D , can be obtained straightforwardly from Equation (43) as where s I−N is the jump in the orientational order parameter at the I − N transition, 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, Figure 4, as 0.313, which equals exactly log J(2.53) , according to Equation (44). For Q 0 ≤ 3.0, the electrostatric potential stays completely on the PB branch. For Q 0 > 3.0, it stays on the PB branch for Q ≤ 3.0, 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 ofz at this jump, that can be read off in Figure 3, also corresponds to the thickness of the nematic layer. In addition, for Q 0 > 3.0, φ * (Q) − φ * (Q 0 ) = 0 at the surface, only for the MS branch, while the PB branch is suitably rescaled due to the form of the boundary condition of Equation (39) on the MS and PB branches.

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. Of course this necessarily also implies 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 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 or the corresponding tensorial expression consistent with Equation (46). Above ξ 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, Equation (20), and one ends up with a system of two coupled non-linear differential equations. This would, in its turn, also require the introduction of the surface energy that would pin the surface value of the order parameter. These, 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, Equation (5), such as the quadrupolar electrostatic term. Formally, this can be seen as leading to a modification in Equation (A8) of the form 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 in 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 also allow for interesting phenomena in terms of the effective electrostatic interactions between the bounding surfaces that would no doubt deviate from the standard expectations.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
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 Equation (6)  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 −βH * [r n , n n ] = i where we explicitly used the definitions of the microscopic scalar charge density and the microscopic tensor nematic order parameter density, Equations (1) and (2). Since the Q ij (x) and ρ(x) functional integrals in Equation (A4) are obviously Gaussian, these variables can be integrated out explicitly, yielding a pure field theoretical representation of the original configurational partition function, Equation (A1), in terms of the tensor, Φ ij (x), and scalar, φ(x), fields with where the effective field action is finally obtained in the form 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. The formal identity of the configurational, Equation (A1), and auxiliary field representations, Equation (A6), can be recapitulated as Ξ[r n , n n ] = Ξ[Φ ij (x); φ(x)]. (A9) 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 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.