The role of counterions in ionic liquid crystals

Previous theoretical studies of calamitic (i.e., rod-like) ionic liquid crystals (ILCs) based on an effective one-species model led to indications of a novel smectic-A phase with a layer spacing being much larger than the length of the mesogenic (i.e., liquid-crystal forming) ions. In order to rule out the possibility that this wide smectic-A phase is merely an artifact caused by the one-species approximation, we investigate an extension which accounts explicitly for cations and anions in ILCs. Our present findings, obtained by grand canonical Monte Carlo simulations, show that the phase transitions between the isotropic and the smectic-A phases of the cation-anion system are in qualitative agreement with the effective one-species model used in the preceding studies. In particular, for ILCs with mesogenes (i.e., liquid-crystal forming species) carrying charged sites at their tips, the wide smectic-A phase forms, at low temperatures and within an intermediate density range, in between the isotropic and a hexagonal crystal phase. We find that in the ordinary smectic-A phase the spatial distribution of the counterions of the mesogens is approximately uniform, whereas in the wide smectic-A phase the small counterions accumulate in between the smectic layers. Due to this phenomenology the wide smectic-A phase could be interesting for applications which hinge on the presence of conductivity channels for mobile ions.


I. INTRODUCTION
Ionic liquid crystals (ILCs) are versatile materials which exhibit, on the one hand, properties of ionic systems, such as the capability of charge transport, and, on the other hand, they are able to form mesophases which is the distinctive feature of liquid crystals [1,2].
Recently, several numerical studies [3][4][5][6] aimed at gaining insight into the link between the underlying molecular features of ILCs and their resulting properties, like the phase behavior or the formation of nanostructures. The complicated interplay of anisotropy and ionic molecular properties renders ILCs to be challenging objects for theoretical studies. Molecular dynamics (MD) computer simulations can be performed with a rather detailed description of the underlying molecules [7]. Such simulations can provide a rich and detailed picture of the structures of the formed mesophases. But the complexity of the underlying models makes it still difficult to pinpoint basic characteristics of these systems which give rise to the observed properties (or at least being essential for their appearance).
In this regard, studying simplified models allows one to elucidate the generic interplay of the key features of ILCs. This might provide insight into elementary mechanisms on the molecular level, which are necessary ingredients to the rich phenomenology of ILC systems. The natural drawback of simplistic models is that not all of the observable features are captured. At the same time, one is able to rule out minimal models, which are insufficient to describe (or explain) particular phenomena * bier@is.mpg.de of the object of interest. Thereby, it is possible to successively approach and isolate the minimal physical requirements for the diverse phenomenology occurring for complex systems such as ILCs. Moreover, this procedure yields a generic understanding of the individual microscopic mechanisms involved in the complex, mesoscopic and macroscropic materials properties.
Typically the mesogenes of ILC systems are composed of long alkyl-chains in combination with charged groups, like imidazolium rings [1]. Due to the alkyl chains, ILC molecules exhibit a large aspect ratio, rendering the molecules highly anisotropic in shape. However, due to the flexibility of the molecules the formation of liquid crystalline phases (or mesophases) is not self-evident. Notwithstanding, in the spirit of the aforementioned simplified model descriptions, within the present study we refrain from focussing on the formation of orientationally ordered nanostructures, i.e., mesophases, by flexible molecules, which is typically driven by complex mechanisms like microphase segregation [1,2]. Instead, we consider a coarse-grained model of rigid particles, which gives generically rise to the formation of mesophases (in particular smectic phases) due to the underlying shape of the molecules and due to the presence of van der Waals interactions. We are aiming at elucidating the influence of charges on the formation of liquid crystalline structures. More specifically, we are concerned with the role of counterions on the structure formation. To this end, we apply a previously introduced model [3] of ILCs, which incorporates only one of the ion species explicitly, while the counterions are only considered as a homogeneous screening background. Within this one-species ILC model an interesting phase behavior could be observed which is very sensitive to the charge distribution within the mesogenic ions [5].
Inspired by ILCs which are based on positively charged imidazolium rings attached to long alkyl chains [2], we adopt this notion and refer to the mesogenic species as cations, while the spherical counterions are referred to as anions.
This study is structured as follows: In Sec. II the model used to describe the ILC molecules is presented, and additional information concerning the methods to analyze the simulation data is given. In Sec. III our present results, obtained by grand canonical Monte Carlo (MC) simulations, are shown and discussed, followed by our conclusions in Sec. IV.

A. Pair interactions
For common examples of ionic liquid crystals (ILCs), the cations (⊕) are characterized by a molecular structure exhibiting charged groups, e.g., in the form of imidazolium rings, attached to rather long alkyl-chains, whereas the anions ( ) are typically much smaller in size, e.g., in the form of iodide (I − ) [1,2]. While the charged compounds introduce ionic properties, it is the presence of the alkyl-chains as part of the cations which leads to the formation of liquid-crystalline phases, socalled mesophases, within ILCs. Nonetheless, the occurrence of mesophases within ILCs, such as smectic phases, is a highly non-trivial matter, because the internal flexibility of the alkyl-chains hinders the formation of layered structures. However, the interplay of the hydrophilic charged groups and the lipophilic alkyl-chains stabilizes smectic structures via microphase segregation of the two molecular compounds [1].
In line with the scope of the present study, we reduce the degree of intra-molecular complexity by considering a coarse-grained description of ILCs, in which the cations (⊕) are represented by rigid ellipsoids of length L ⊕ and width R ⊕ , while the anions ( ) are spherical particles of diameter R < R ⊕ . Although such an approach does not allow one to study the underlying mechanisms leading to the formation of smectic layers within the above described actual ILCs, here we are interested in the dependence of the smectic structures on the location of the charged groups. Thus, the model parameters (see below) are tuned such that, within the simplistic model, smectic phases are formed. We shall analyze how these layered structures depend on the intra-molecular location of the charges.
To this end, the pair interactions used here consist of (i) a hard-core contribution φ hc ij , (ii) an attractive energy contribution φ vdW ij , accounting for van der Waals forces , and (iii) the electrostatic interaction φ es ij due to the presence of charges.
Thus, the total pair potential φ ij between two particles i and j, where i, j ∈ {⊕, }, reads As mentioned above, the cations are rigid prolate ellipsoids of length-to-breadth ratio L ⊕ /R ⊕ . Thus, the orientation of a cation is fully described by the direction ω ω ω of its long axis. The total interaction potential between a pair of cations, i.e., Eq. (1) with i = j = ⊕, is of the following form: where r r r denotes the center-to-center distance vector between two cations with orientations ω 1 ω 1 ω 1 and ω 2 ω 2 ω 2 . The contact distance between two cations depends on the orientations of both cations and on the di-rection of the center-to-center distance vector, expressed by the unit vectorr r r := r r r/|r r r|. In Eq. (2), the contributions beyond the hard-core repulsion at contact, i.e., for |r r r| ≥ R ⊕ σ, are subdivided into two parts. The attractive interactions φ att ⊕⊕ , due to short-ranged van der Waals forces between the cations, are modeled by the Gay-Berne potential φ GB (r r r, ω 1 ω 1 ω 1 , ω 2 ω 2 ω 2 ) = φ att ⊕⊕ [8,9]. φ GB is a modification of the Lennard-Jones pair potential designed for ellipsoidal particles: with the anisotropic interaction strength Both the contact distance R ⊕ σ, i.e., Eq. (3), and the direction-and orientation-dependent interaction strength ε(r r r, ω 1 ω 1 ω 1 , ω 2 ω 2 ω 2 ), i.e., Eq. (5), depend on the cation length-to-breadth ratio where ε R⊕ /ε L⊕ is called the anisotropy parameter, defined as the ratio of the depth ε R⊕ of the Gay-Berne potential minimum for parallel particles positioned side by side, i.e., withr r r ·ω 1 ω 1 ω 1 =r r r ·ω 2 ω 2 ω 2 = 0, and of the depth ε L⊕ of the Gay-Berne potential minimum for parallel particles positioned end to end, i.e., withr r r · ω 1 ω 1 ω 1 =r r r · ω 2 ω 2 ω 2 = 1. The length-to-breadth ratio L ⊕ /R ⊕ and the anisotropy parameter ε R⊕ /ε L⊕ specify the molecular properties like the shape and the chemical structure of the underlying cation molecules within the present coarse-grained model. As we aim for comparing the present findings with those of our previous studies using a one-species model -comprising only the ellipsoidal cations while the anions were incorporated implicitly as a homogeneous screening background -we choose the following set of parameters for the Gay-Berne potential, which allows for such a comparison with Ref. [5]: It is worth mentioning, that in the case of spherical cations, i.e., for L ⊕ = R ⊕ , Eq. (4) reduces to the (isotropic) Lennard-Jones potential iff ε R⊕ /ε L⊕ = 1, because in that case σ(r r r, ω 1 ω 1 ω 1 , ω 2 ω 2 ω 2 ) = 1 and ε(r r r, ω 1 ω 1 ω 1 , ω 2 ω 2 ω 2 ) = ε ⊕⊕ . The relations L ⊕ = R ⊕ and ε R⊕ = ε L⊕ describe molecules of rather spherical shape, which, however, due to their internal chemical structure exhibit non-spherical van der Waals interactions.
The remaining contribution φ es ⊕⊕ in Eq. (2) is the electrostatic repulsion between the cations. Since for the present study the electrostatic interactions are of particular interest, their implementation is discussed in detail in the next section. At this point, we only point out that the charge sites are located symmetrically at a distance D from the cation center along the long axis. Thus, φ es ⊕⊕ (r r r, ω 1 ω 1 ω 1 , ω 2 ω 2 ω 2 ) depends not only on the distance r = |r r r| between the centers of two cations, but also on the orientations ω 1 ω 1 ω 1 and ω 2 ω 2 ω 2 , as well as on the relative direction of the centers of the cations.
Anions are modeled as hard spheres of diameter R with a negative charge site in their center. To account for the omnipresent attractive van der Waals forces, typically the Lennard-Jones potential is used to mimic these dispersion-induced interactions between spherical particles. However, here, we are neglecting any contributions arising from dispersion forces between anions, i.e., φ vdW = 0, because we focus on ionic effects originating from the much stronger electrostatic interaction. In particular we focus on the influence of the anion distribution on the liquid-crystalline structure of the cations. Thus, the pair potential φ (r) between two anions separated by distance r reads φ es is the residual anion-anion electrostatic interaction. Due to the spherical shape of the anions it is an isotropic function and depends only on the distance r. As mentioned above, further details of the electrostatic interactions are provided in the next section.
We do not consider contributions due to van der Waals forces between cations and anions, i.e., φ vdW ⊕ = 0. These contributions might be necessary in order to describe quantitatively reliably a specific type of ionic liquid crystal system. Here, however, we are not interested in such a quantitative analysis, but we are rather aiming at a general understanding of the mechanisms leading to structures and distinct phases in ILCs. In particular, we want to understand how a possibly non-uniform counterion distribution affects the phase behavior which is observed within the effective one-species model used in Ref. [5]. In practice, this means that the effectively treated electrostatic interaction is altered such that the valency dependence of the Coulomb potential is now explicitly incorporated. Since this is a key issue of the present study, the implementation of the electrostatic interactions among all particles is explicitly given in the following section.
Finally, we point out that the remaining independent parameters R ⊕ (which denotes the cation width, see Eq. (3)) and the cation-cation interaction constant ε ⊕⊕ (see Eq. (5)), are chosen as the length and the energy scale of the system.

B. Electrostatic energy contributions
Within the present model, both cations (⊕) and anions ( ) carry point-like charge sites. While each anion carries a single charge site in its center, cations exhibit two distinct charge sites, located at a distance D from their geometrical center. Thus, the electrostatic interactions among all types of particles are given by and φ es ⊕ = −2 γ φ (|r r r + D ω ω ω|) +φ(|r r r − D ω ω ω|) with the electrostatic interaction strength γ = q 2 /(4πε 0 ), where ε 0 denotes the vacuum permittivity, and q > 0 is the charge of a single site of the cations. The factors of 4 in Eq. (9) and of 2 in Eq. (11) occur, because the negative charge site in the center of an anion has to be twice as strong as a single cation charge site, such that the valency is the same for cations and anions. Thus, in order to guarantee global charge neutrality, the system contains the same number of cations and anions (this issue will be discussed in detail in the next section). We note, that the factor of 4 is also recovered in Eqs. (10) and (11) for D = 0. In Eqs. (9) and (10) the electrostatic energy contributions φ es ≥ 0 and φ es ⊕⊕ ≥ 0 are repulsive, while in Eq. (11) the negative sign indicates the electrostatic attraction of cations and anions, i.e., φ es ⊕ ≤ 0. For point-like charges q in d = 3 dimensions, the electrostatic interaction potential decays as the inverse of the distance, i.e.,φ(r) = 1/r (see Eq. (9)). Thus, it is longranged and this property is well-known to lead to a wide range of peculiarities of ionic systems [10]. In the present context it is important to note, that, in order to accurately account for the long-ranged character of the interactions, in computer simulations (MC or MD) of bulk systems, based on periodic boundary conditions, one has to resort to sophisticated methods like the Ewald summation [11]. The Ewald summation splits the full electrostatic contribution to the total energy of a given configuration into a short-ranged and a long-ranged part by expanding the actual charge density by a set of Gaussian screening charge clouds. While the first contribution is a sum over short-ranged interaction potentials and can be calculated in real-space, the second contribution contains the long-ranged part which can be calculated by Fourier transformation by expoiting the periodic boundary conditions. A different perspective on this method is, that the Ewald summation separates the electrostatic energy into two contributions, such that the first one expresses the valency dependence of the Coulomb interaction, while the second one is determined by the long-ranged part.
The motivation for the present study is to analyze the effects incorporated due to accounting for both ion species and to compare them with those occurring within the effective one-species model which has been used previously for studying the bulk phase behavior of ILCs [3,5]. First, it is interesting to study a system of cations and anions interacting via short-ranged potentials and analyze how the valency dependence affects the previous results. In a second step the full Ewald summation allows one to investigate the relevance of the long-ranged character. In this way one can gain insight into the influence of both these two fundamental properties of the Coulomb interaction on the phase behavior of ILCs.
The interaction potential, resembling the short-ranged contributions to the total electrostatic energy contribution, is described by a Yukawa potential (see Eq. (9)) with decay length λ = 5R ⊕ (Eqs. (9)-(11)). While Eq. (12) exhibits the same functional form as the one used in Refs. [3,5], it is important to emphasize that the decisive new aspect of the present study is the actual presence of counterions, i.e., positive (repulsive) and negative (attractive) electrostatic energy contributions. The decay length λ = 5R ⊕ has been chosen such as to match the parameters of the interaction potential in Ref. [5]. In addition, for the given cation length L ⊕ = 4R ⊕ and the anion diameter R < R ⊕ , the chosen decay length λ = 5R ⊕ is larger than the particle sizes, such that Eq. (12) corresponds to a weak artificial screening of the pure Coulomb potential 1/r. A previous study [12] suggests that the valency dependence rather than the longranged character of the electrostatic interaction is decisive for the phase behavior of ionic fluids, i.e., the weak artificial screening in Eq. (12) is expected to give rise to at most some quantitative consequences. Moreover, in Ref. [13], it has been shown recently that there are plenty of alternatives to the functional form of Eq. (12), which serve to describe the structure of actual ionic systems remarkably well. Thus we expect that incorporating the full Coulomb interaction via the Ewald method only leads to a quantitative change of the phase behavior, such as a shift of the phase transitions in the temperature-density plane, which, however, does not significantly alter the occurrence of the phases or their structural properties on a qualitative level. In Fig. 1 the full potentials for the cation-anion and the cation-cation interactions are presented. While panels (a) and (b) show the full cation-anion interaction for D = 0 and D/R ⊕ = 1.8, respectively, in panels (c) and (d) the cation-cation interaction potential is illustrated for the same types of charge distributions of the cations as in (a) and (b), respectively. For the considered particle sizes, i.e., L ⊕ /R ⊕ = 4 and R /R ⊕ = 1/4, the excluded volume (illustrated by the beige area with a solid black rim) between cations and anions is much smaller compared to the excluded volume between two cations. If D = 0, the electrostatic interaction is strongest for a side-by-side position of the two considered particles, while for D/R ⊕ = 1.8 it is strongest at the tips. Moreover, the additional electrostatic repulsion among cations leads to a narrowing of the most attractive region (which stems from the attractive Gay-Berne interaction) at close distances for a side-by-side configuration of the cations, if the cation charges are located at the center of the molecule.

C. Pair distribution functions
In the course of the simulations the structure of the fluid is analyzed via pair distribution functions. For simulations of bulk systems, the pair distribution functions g ij (r r r, r r r ) := ij (r r r|r r r )/ j can be defined as the ratio of the conditional density ij (r r r|r r r ) of particles of species i at position r r r, provided a particle of species j is located at r r r , and of the constant mean density j in the simulation box [10,14].
As we are mainly interested in observing (smectic) layer structures, we monitor the pair distribution in the direction of the layer normaln n n, as well as the particle distribution within the layers, i.e., in directions r r r ⊥ lateral to the layer normal. Parallel to the layer normaln n n, the statistics along the simulation trajectories invokes all pairs of particles at distances n := |(r r r − r r r ) ·n n n|: Additionally, for the planes perpendicular ton n n -associated with the vector r r r ⊥ := r r r − (r r r ·n n n)n n n -we monitor the radial distribution of cation pairs via where n ∈ {0, d} refers to the plane for which n = 0 and n = d, respectively. Thus g (0) ⊕⊕ (r ⊥ ) monitors the (cation) pair distribution in the 0-th plane, i.e., the plane which contains the reference cation at r r r , while g (d) ⊕⊕ (r ⊥ ) monitors the distribution of particles in the two neighboring layers, with respect to the cation at r r r .
We add the following remarks concerning the computation of Eqs. (13) and (14): • The direction of the layer normaln n n is determined by calculating the director [15] of each configuration along the simulated trajectories. For the relevant cases, which are analyzed within the scope of the present study, the director and the layer normal point into the same direction, i.e., they are (almost) parallel.
• The evaluation of the conditional densities || ij (n) and (n) ⊕⊕ (r ⊥ ) requires to count the number of particles which are a distance n and a distance r ⊥ , respectively, apart from the central reference particle. To this end, one considers small but nonzero volumina at n and r ⊥ , which are given by straight slices of width ∆n and by annuli of width ∆r ⊥ , respectively. The straight slices for calculating || ij (n) extend in lateral direction up to the boundaries of the simulation box (see the illustration of such a slice in Fig. 2(a)), while the annuli for calculating (n) ⊕⊕ (r ⊥ ) have an extent in the direction of the layer normal from ∆ || = 1 × R ⊕ to ∆ || = 2 × R ⊕ .
• While the volume of an annulus is given by ∆V ⊥ = π∆ || ∆r ⊥ (∆r ⊥ + 2r ⊥ ), the volumina ∆V || of the straight slices cannot be calculated straightforwardly, as they are cut off at the boundaries of the simulation box. However, by considering a reference configuration with a homogeneous and isotropic distribution of particles in the simulation box, the volumina of the straight slices can be approximated by ∆V || ≈N / ref , whereN denotes the number of particles counted within the considered slices for the isotropic and homogeneous reference configuration of mean density ref .
Here, βH(ζ) denotes the total (potential) energy of a given configuration ζ (see Sec. II A and Eq. (1)) in units of the thermal energy k B T = β −1 , βµ is the chemical potential, and N (ζ) is the total number of particles. Since we are considering 1:1-ionic mixtures, there is an equal number of cations and anions, i.e., N ⊕ (ζ) = N (ζ) = N (ζ)/2.
Since we are mainly concerned with (smectic) structures formed by the mesogenic cations, the number density is given in terms of the cation packing fraction η := π 6 L ⊕ R 2 ⊕ N ⊕ /L 3 , where π 6 L ⊕ R 2 ⊕ denotes the cation volume and N ⊕ refers to the thermally averaged total number of cations. Temperature is measured in terms of the ratio of the thermal energy k B T and the interaction strength ε ⊕⊕ of the Gay-Berne potential (Eq. (5)), i.e., If not stated otherwise, for the simulations we used the following model parameters: L ⊕ /R ⊕ = 4, ε R⊕ /ε L⊕ = 3, R /R ⊕ = 1/4, γ/(ε ⊕⊕ R ⊕ ) = 0.045, and λ/R ⊕ = 5. Furthermore, for calculating the total energy, all pair interactions (Eq. (1)) have been truncated beyond the range R cut /R ⊕ = 6.
The phase diagrams displayed in Fig. 4 are obtained by performing simulations for numerous state points (T * , βµ), whereby each run is initialized with an isotropic configuration. By performing short additional simulation runs initialized with smectic-A and crystalline configrations, it has been checked for all considered state points that the simulation results do not depend on the initialization. The phase transitions in Fig. 4 are resolved with an accuracy of ∆(βµ) = 0.1 in terms of the chemical potential βµ. Taking into account the values of ∂η ∂(βµ) obtained from the simulations, this accuracy is sufficient to resolve the white two-phase regions in Fig. 4 with an accuracy of ∆η ≈ ∂η ∂(βµ) ∆(βµ) ≤ 0.01 in terms of the packing fraction η.
A. Dependence of smectic structures on the charge distribution within the cations

Formation of the phase SAW
The present model (see Sec. II) can be understood as an extension of the effective one-species model of ILCs used in previous theoretical studies [3,5,12], such that the present, extended model accounts for the explicit presence of both ionic species. The present approach allows us to study explicitly the effect of incorporating the valency dependence of the Coulomb interaction on top of the previously studied one-species model. We note, that yet it is necessary to introduce hard-core interactions between the cations and anions, as well as among the anions, in order to avoid divergences in the electrostatic interactions (caused by mutual penetration).
One of the most striking findings within the effective one-species model in Ref. [5] is the sensitive dependence of the occurring smectic structures on the location of charges within the ellipsoidal cations. Therefore, we first analyze how the presence of counterions affects this dependence.
We start our analysis at fixed temperature T * = 0.55 and discuss the structures which are formed at sufficiently high densities, such that the isotropic phase becomes thermodynamically unstable (or at least metastable) with respect to smectic phases. We note that, for the considered parameters, no nematic phase is observed in the density regime between the isotropic and the smectic phase. (For larger values of the cation length-to-breadth ratio L ⊕ /R ⊕ this might, however, be the case.) Performing Monte Carlo simulations for D = 0 and D/R ⊕ = 1.8 at η > 0.35, two distinct structures, which are shown in Figs. 2(a) and (b), can be observed. The snapshots show an ordinary smectic-A phase for D = 0 (a) and the phase S AW for D/R ⊕ = 1.8 (b). While in panel (a) one recognizes a typical smectic layer structure with a layer spacing comparable to the cation length L ⊕ , in panel (b) alternating layers are observed, in which the cations (illustrated in red) are well-aligned with the layer normal n n n, as well as intermediate layers of cations (depicted in blue) which are oriented almost perpendicularly to the layer normal. Interestingly, the intra-layer structure is also different for the two phases. For example for the two layers of the common phase S A and the phase S AW in FIG. 2. In (a) and (b) we show snapshots of two configurations, belonging to the phases SA and SAW , respectively. While for the phase SA a layer spacing of the size of the cation length L⊕ can be observed, and all cations tend to be aligned with the smectic layer normaln n n, in (b) the alternating layer structure of the phase SAW is clearly visible. In between the layers of cations, being well-aligned with the layer normaln n n (red-colored ellipsoids), secondary layers are observed, in which the cations are preferentially perpendicular to the layer normal (blue ellipsoids). Due to the alternating layer structure the layer spacing is significantly increased. In panels (a) and (b) the anions are depicted as small black dots, which, however, in order to increase visibility, are three times larger than the actual anions (R /R⊕ = 1/4). In (c) and (d) the intra-layer structures of the phases SA and SAW , respectively, are shown. While in (c) for the phase SA a fluid-like structure is observed, the snapshot of a main layer of the phase SAW resembles a hexagonal structure (highlighted by the green hexagon and the thick violet lines in panel (d)). This observation is accompanied by a higher cation density within the main layers of the phase SAW as compared with the SA layers. Note, that the green slab in the upper left corner of panel (a) depicts a slice which is used to evaluate the pair distribution functions g || ij (n) in the direction of the layer normaln n n (green arrow). Similarly, in panel (c) the green concentric circles indicate the annulus for calculating the lateral pair distribution function g Figs. 2 (c) and (d), respectively, one observes a (typical) fluid-like structure for the phase S A , while a dense and fairly ordered structure is observed for the main layers of the phase S AW .
In order to discuss the structure of the two phases in more detail, the pair distribution functions in the direction of the layer normaln n n, i.e., g || ij (n), and in lateral directions perpendicular ton n n, i.e., g (n) ⊕⊕ (r ⊥ ), (compare Eqs. (13) and (14)) are analyzed in Fig. 3. In Figs. 3(a) and (i) the pair distribution functions are shown for the phase S A , as depicted in Fig. 2(a), for D = 0. According to the red curve in Fig. 3(a), the cation-cation correlations, i.e., g || ⊕⊕ (n), clearly show that a layer structure with layer spacing d/R ⊕ ≈ 3.5, comparable to the particle length L ⊕ = 4R ⊕ , is formed. Panel (i) shows, for this smectic-A structure, the lateral correlations among the cations, i.e., g (n) ⊕⊕ (r ⊥ ). Within the layer in which the reference cation is located (black curve, corresponding to n = 0), clearly a fluid-like structure with rapidly decaying correlations is observed. For neighboring layers (magenta curve, i.e., n = d) one finds no correlations at all. Thus, these findings confirm that the structure shown in Fig. 2(a) is an ordinary smectic-A phase (S A ).
The phase S AW , formed for D/R ⊕ = 1.8, is shown in Figs. 3 (c) and (iii). The cation-cation correlations indicate the alternating layer structure consisting of main layers of high cation density (the peaks of the red curve at n/R ⊕ = 0 and n/R ⊕ ≈ 6) and secondary layers (at n/R ⊕ ≈ 3 and n/R ⊕ ≈ 9). Interestingly, analyzing the lateral structure of the phase S AW in panel (iii), we find a pronounced structure which is distinct from the (fluidlike) pair distribution function obtained for the ordinary phase S A . The peak positions yield that this resembles a hexagonal structure, which is also confirmed by the snapshot of a S AW layer, shown in Fig. 2 (d). However, the lateral correlations within the neighboring layer with respect to the reference cation (magenta curve in Fig. 3(iii)) are almost vanishing and therefore the phase S AW is indeed a smectic phase and not a crystal-like structure. Thus, neighboring smectic layers can be sheared without any cost of free energy. We note, that the weak oscillations, which are visible in the magenta curve in Fig. 3(iii), are artifacts of the periodic boundary conditions: If the layer normal is not parallel to one of the main axis of the cubic simulation box, i.e.,n n n / ∈ {x x x,ŷ ŷ y,ẑ ẑ z}, the smectic layers are not correctly continued by periodic images of the simulation box. For example, the smectic layer in the lower right corner of Fig. 2(b) is continued to below by the periodic image of the third smectic layer counted from the lower right corner. These artificial correlations (in lateral directions) between neighboring smectic layers occur in principle also for the ordinary phase S A . However, due to the short-ranged lateral correlations, they are not visible in Fig. 3(i).
Presumably the different structures of the phases S A and S AW are directly related to the slightly higher density within the main layers of the phase S AW as compared to the layers of the phase S A (see the values of the cation-cation pair distribution function g ⊕⊕ (n), i.e., the red curves in Figs. 3(a) and (c), at n = 0). The higher local cation densities within the S AW main layers are stabilized by the Gay-Berne attraction for parallel oriented cations. Yet, the charges at the tips are indispensable for the formation of the phase S AW as they provide a net repulsion of neighboring smectic layers and therefore make it energetically favorable to maintain a larger distance between the dense main layers separated by the    The alternating layer structure of cations with significantly larger layer spacing d/R⊕ ≈ 6 > L⊕/R⊕ is apparent from g || ⊕⊕ (n). Due to the enhanced density within the main cation layers (see the maxima of g || ⊕⊕ (n) at n = 0 and n/R⊕ ≈ 6) in lateral directions a hexagonal structure can be observed, unlike the fluid-like lateral structure of the ordinary phase SA. However, correlations between cations in neighboring layers are almost absent and thus the phase SAW is a genuine smectic phase and not a crystal. We note that the small variations in g  ⊕⊕ (r ⊥ ) at large distances r ⊥ > 5 is another artifact, due to insufficient statistics. Interestingly, for the phase SAW a considerable inhomogeneous distribution of anions (in normal direction) can be observed in panel (c). The anions prefer to be close to the locations of the cation charges in the main layers, e.g., at n/R⊕ ≈ 1.8, as can be inferred from g || ⊕ (n). Also for D/R⊕ = 1.8 the ordinary phase SA (at higher temperature T * = 0.58 and η ≈ 0.42) and the hexagonal crystal C (at T * = 0.55 and η ≈ 0.46) can be observed (see panels (d) and (iv), respectively (e) and (v)). However, in contrast to the findings for the corresponding phases for D = 0, for D/R⊕ = 1.8 the anions exhibit a considerably inhomogeneous spatial distribution. They prefer to be located in between the cation layers, similar to the phase SAW .
intermediate secondary layers. Given the relatively weak electrostatic interaction as compared to the Gay-Berne attraction (see Fig. 1), already small density differences within the (main) layers of the smectic-A phases decide on the stability of S A or S AW . The large layer spacing in combination with the slightly increased density in the main layers and the substantially lower density in the secondary layers rationalizes the intermediate mean density range in which the phase S AW is observed.
For D = 0, however, the repulsion between cations in neighboring layers is weaker, as compared to the case D/R ⊕ = 1.8 (see the interaction landscape for the two cases in Figs. 1(c) and (d)). Thus, the energetic benefit of a larger layer spacing is insufficient for stabilizing the phase S AW in favor of the phase S A with the layer spacing being comparable to the cation length. The comparison of the two cases, in which the cation charges are either localized in the center, i.e., D = 0, or the charges are positioned close to the tips, i.e., D/R ⊕ = 1.8, underscores the importance of the charge distribution within the cations for the formation of the phase S AW : In agreement with previous results [5], obtained within the effective one-species model, one can conclude, that, due to the cation charges at the tips, a considerable net repulsion of adjacent (main) layers occurs at small distances (see Fig. 1(d)), which drives the main layers apart to distances larger than the cation length L ⊕ .
In the case D/R ⊕ = 1.8, apparently the incorporation of explicit anions does not affect the formation of the phase S AW . However, as we shall present in the next subsection, the distribution of the anions does sensitively depend on the charge distribution within the cations.

Anion distribution
The different cation charge distributions for D = 0 and D/R ⊕ = 1.8 not only lead to remarkably different smectic structures, formed by the ellipsoidal cations, but moreover, for the two cases the distribution of anions is also distinct.
Revisiting Fig. 3(a), which depicts the pair distribution functions g || ij (n) along the layer normal for the ordinary smectic-A phase S A for D = 0, one finds that the anions are rather homogeneously distributed around the layers of cations. (See the green and blue curves which show only minor spatial variations.) The pair distribution function g || ⊕ (n) (blue curve) shows a sparse tendency of the anions to be located in between the cation layers. In contrast, for the phase S AW at the same temperature T * = 0.55 the anions are strongly pushed out of the main layers formed by the cations. The highest probability to find the anions is close to the location of the charges of the cations in the main layers, i.e., the maxima of g || ⊕ (n) are found at a distance n/R ⊕ ≈ D/R ⊕ = 1.8 away from the centers of the smectic layers.
A similarly strong inhomogeneous distribution of anions is found for the ordinary phase S A which is formed by cations with D/R ⊕ = 1.8 and at higher temperatures. In Fig. 3(d) the phase S A , for D/R ⊕ = 1.8 at the slightly higher temperature T * = 0.58, is analyzed. While the structure (in normal direction) of the cations (red curve in panel (d)) is very similar to the phase S A formed for D = 0 (see Figs. 3(a) and (i)), analogously to the case of the phase S AW , the anions are preferentially located in between the smectic layers formed by the cations (see the blue curve in panel (d)).
The same findings are obtained for the respective crystal-like phases which are observed in both cases, i.e., for D = 0 and D/R ⊕ = 1.8, at low temperatures and for large densities (see the detailed discussion of the phase behavior for the two cases in the following Subsec. III B). Thus we conclude that, for charges at the tips of the cations, not only layer structures of the cations are observed, but also for the anions, although the inhomogeneity in the distribution of the anions is not as pronounced .52 a first-order phase transition occurs from the isotropic fluid phase I (violet-colored area) to the ordinary smectic-A phase SA (blue), while at lower temperatures a direct first-order phase transition takes place to the crystalline phase C (yellow) with a hexagonal lattice structure. In contrast, for D/R⊕ = 1.8 (b) at low temperatures the phase SAW (green) is stable within an intermediate density range, such that, upon increasing η, first one observes a discontinuous phase transition from the isotropic fluid I to the phase SAW , followed at large densities by a first-order phase transition from the phase SAW to the hexagonal crystal C.
as for the cations. In line with these findings, for D = 0, i.e., if the charges of the cations are localized in the centers, layering is only observed for the cations. Supposedly, this structural behavior is driven by the interplay of, on the one hand, the electrostatic attraction of the anions towards the cation centers and, on the other hand, with the steric hard-core repulsion, which hinders the anions to penetrate the cation layers.
This striking difference in the anion distributions for the two cases is also relevant for potential applications of ILCs, because some technologies incorporating ILCs like dye-sensitized solar cells (DSSCs) benefit from a higher density of counterions, as they are needed as charge carriers in these applications [17]. Thus, ILCs with charges at the tips, seem to be the best candidates for such applications, whereas cations with D = 0 seem to be good candidates if a smectic phase of cations in combination with a homogeneous distribution of anions is required.

B. Phase behavior
After having analyzed the differences of the smectic structures associated with the charge distributions D = 0 and D/R ⊕ = 1.8, we now focus on the phase behavior, i.e., identifying those regions in the (T * , η) plane, which correspond to the aforementioned distinct structures.
First, we analyze the phase behavior of ILCs consisting of cations with the charges localized in their center, i.e., for D = 0. Figure 4(a) shows the T * -η phase diagram for D = 0. Within the investigated temperature range T * ∈ [0.45, 0.6], the isotropic fluid phase I (violet) is thermodynamically stable up to (cation) packing fractions η ≈ 0.25 . . . 0.35. At sufficiently high temperatures T * ≥ 0.52 a first-order phase transition, indicated by a density gap, is observed to the ordinary smectic-A phase S A (blue). However, at lower temperatures, a direct phase transition to a crystal-like structure C (yellow) occurs. While the phase S A and the crystalline phase C show a rather similar layer spacing d/R ⊕ ≈ 3.5 in normal direction (compare the pair distribution functions g || ⊕⊕ (n) shown in Figs. 3(a) and (b)), their lateral structures are quite distinct. For the phase C a hexagonal ordering (see Fig. 3(ii)) can be observed and, moreover, there are strong correlations between neighboring layers. Indeed, comparing the observed peaks in the distribution functions with an ideal three-dimensional hexagonal lattice, one finds very good agreement concerning the peak positions (in the 0-th layer, as well as in the neighboring layers). Interestingly, although in general the anions have a tendency to be localized in between the cation layers in the phase C, here the distribution of anions is fairly homogeneous (relative to the pronounced peaks in the cation distribution). Thus, while at the considered high densities the cations form a well-marked hexagonal lattice, the anions are still rather motile. Now we turn to the ILCs with the cation charges at the tips, i.e., D/R ⊕ = 1.8. The corresponding phase diagram is shown in Fig. 4(b). At high temperatures T * 0.58 one also finds a first-order phase transition from the isotropic fluid I to the ordinary smectic phase S A . Besides the additional inhomogeneous distribution of the anions, which has not been observed for the phase S A for D = 0 (see Sec. III A 2), here the smectic-A phase is very similar to the phase S A forming for D = 0. The layer spacing d/R ⊕ ≈ 3.5 is comparable to the cation length L ⊕ and the cations are well-aligned with the layer normal. Furthermore, the lateral correlations (see Fig. 3(iv)) clearly exhibit a fluid-like structure within the smectic layers. If the temperature is lowered to T * 0.55, a different structure appears. From Figs. 3(c) and (iii) one infers that this distinct structure corresponds to the phase S AW (green). Since the phase S AW emerges in an intermediate density region, for D/R ⊕ = 1.8 the isotropic fluid I is stable only at small densities at low temperatures. An alternating structure of main and secondary layers occurs which leads to a significantly increased layer spacing d/R ⊕ ≈ 6.0. Due to this increased layer spacing, driven by the electrostatic repulsion of neighboring main layers, the bulk densities η of the phase S AW are lower than the bulk densities η at which the ordinary phase S A is observed. However, by further increasing the num-ber of particles in the system (via raising the chemical potential), at sufficiently high densities the phase S AW becomes metastable with respect to the hexagonal crystal C (see Fig. 4(b) for T * 0.55 and η 0.4).
Comparing the present findings with the theoretical predictions of Ref. [5] (see Fig. 5 therein), it is remarkable that the DFT results for the effective one-species model, on a qualitative level, predict the stability of the phase S AW in the same thermodynamic region as the present Monte Carlo simulations for the enhanced ILC model, i.e., at low temperatures and at intermediate densities.
Furthermore, in both approaches the charge distribution of the (ellipsoidal) cations turns out to be crucial for the formation of the phase S AW , i.e., in order to have a stable phase S AW it is indispensable that the charges are close to the tips of the cations.
Finally, it is worth mentioning, that although no stable state points of the phase C have been found at high temperatures, i.e., T * 0.52 for D = 0 and T * 0.58 for D/R ⊕ = 1.8, respectively, at sufficiently large densities a crystal-like phase is expected to be the stable configuration, even at high temperatures. Supposedly, with our MC simulations we have been unable to reach the very high density region, which requires a large number of particles in the simulation box and thus slows down the simulations.

IV. SUMMARY AND CONCLUSIONS
The objective of the present study is to shed light on the role of counterions in forming smectic structures in ionic liquid crystals (ILCs). In particular, our analysis aims at investigating how the phase behavior and the structural properties of an effective one-species model, which has been employed in previous theoretical studies of ILCs [3,5,6], are affected by taking the counterions explicitly into account. These previous models represent a simplified description of ILCs, which are composed of anisotropic mesogenic ions (for typical examples, these are cations [1,2]), which are embedded in a homogeneous screening background consisting of much smaller anions.
The present model, which incorporates both cations and anions on equal footing, can be understood as an improvement of the previous one, because it does not take at the outset any specific distribution of the anions. Thereby it allows one to investigate not only the influence of the anions on the previously observed liquid-crystalline structures [5], but, moreover, the anion distribution by itself can be analyzed, which is of particular interest for potential technological applications of ILCs, e.g., as electrolyte materials in solar cells [17].
The current coarse-grained model (see Sec. II and Fig. 1) exhibits rigid ellipsoidal cations with two charge sites, symmetrically located at a distance D from the molecular center and spherical anions with one central charge site. All results of this work have been obtained using grand canonical Monte Carlo (MC) simulations.
Depending on the intra-molecular cation charge distribution, distinct smectic structures are observed. They are formed upon increasing the density (expressed in terms of the cation packing fraction η) such that the isotropic fluid I, which is the stable phase at low densities, becomes metastable (and ultimately unstable) with respect to forming smectic layer structures. For T * = 0.55 a first-order phase transition from the phase I to an ordinary smectic-A phase S A is observed for η > 0.35, if the cations carry a single charge site in their center, i.e., for D = 0 (see Fig. 2(a)). The designation of the emerged structure as 'smectic-A' is based on the following observations. The ellipsoidal cations form layers in which they mostly orient parallel to the layer normaln n n and the layer spacing d ≈ 3.5R ⊕ (Fig. 3(a)) is of the size of the cation length L ⊕ . Moreover, in the directions which are lateral with respect ton n n, i.e., within the smectic layers, a fluid-like structure is observed (see Figs. 2(c) and 3(i)).
In contrast for D/R ⊕ = 1.8, i.e., for cations with charges at their tips, at the same temperature T * = 0.55 (and at sufficiently high densities, i.e., η 0.3) a layer structure, which is distinct from the ordinary phase S A , is found: Alternating layers of cations are observed, which are mostly parallel to the layer normaln n n, and of cations, which are oriented mostly perpendicular ton n n ( Fig. 2(b)). This structure can be identified as the wide smectic-A phase S AW , which has been found previously [5]. In agreement with the previous findings, the (main) layers, in which the cations are well aligned, show a much larger (local) density as compared with the (secondary) layers, in which the cations are mostly perpendicular tô n n n (see g || ⊕⊕ (n) in Fig. 3(c)). Moreover, the layer spacing d ≈ 6R ⊕ of the phase S AW is significantly larger than for the ordinary smectic-A phase S A . Due to the high density in the main layers, a lateral hexagonal structure is observed (see Figs. 2(d) and 3(ii)). However, there are no visible correlations among neighboring layers and therefore the phase S AW is a genuine smectic phase and not a crystal.
Interestingly, from comparing the distribution of anions in normal directionn n n, we infer that for D = 0 the anions are rather homogeneously distributed around the cation layers (Figs. 3(a) and (b)), while for D/R ⊕ = 1.8 a pronounced localization of anions in between the cation layers is observed. While for D = 0 the competing electrostatic attraction and the steric repulsion of cations and anions for small center-to-center distances presumably lead to the fairly homogeneous distribution of anions, for D/R ⊕ = 1.8 the anions are not strongly inhibited by steric repulsion to accumulate at the tips of the cations.
Concerning the phase behavior (see the phase diagrams in Fig. 4) for the currently studied model, we find a remarkable (qualitative) agreement with the previously studied one-species model description of ILC systems. The phase S AW is formed only if the cation charges are positioned at their tips. Furthermore, its stable region is found at lower temperatures, as compared to the ordinary smectic-A phase S A , and at intermediate densities, i.e., in particular at lower densities than the ones for the stable phase S A and at higher densities than the ones of the stable isotropic fluid I, which is in agreement with the phase behavior predicted by DFT [5]. In both cases, i.e., for D = 0 and D/R ⊕ = 1.8, at very high densities a three-dimensional hexagonal crystal C is formed (see These results are not only consistent with the previous findings for the effective one-species description of ILCs, but, moreover, they pinpoint the significance of the (intra-molecular) charge distribution for the phase behavior as well as for the structural properties of ILC systems.
Future studies might focus on the dependence of the thermal behavior and structural properties of ILCs on the anion size and shape, but also on the strength of the Gay-Berne potential as compared to the electrostatic interactions. This, in particular, is a subtle issue, because ILCs typically exhibit an effective charge which is even less than one elementary charge, due to effects like charge delocalization [4].