Two-dimensional pattern formation in ionic liquids confined between graphene walls

TL;DR

This study uses molecular dynamics and Monte Carlo simulations to show that ionic liquids confined between graphene walls form pattern structures like stripes and hexagons, which are influenced by charge density and can impact electrochemical properties.

Contribution

It demonstrates that pattern formation in ionic liquids at graphene interfaces is a general feature driven by charge balance and can be modeled with a simple lattice approach, linking to Landau-Brazovskii theory.

Findings

Pattern formation is a universal feature of ionic liquids at interfaces.

Charge density perturbations trigger pattern transitions.

Patterns influence electrochemical interface properties.

Abstract

We perform molecular dynamics simulations of ionic liquids confined between graphene walls under a large variety of conditions (pure ionic liquids, mixtures with water and alcohols, mixtures with lithium salts and defective graphene walls). Our results show that the formation of striped and hexagonal patterns in the Stern layer can be considered as a general feature of ionic liquids at electrochemical interfaces, the transition between patterns being controlled by the net balance of charge in the innermost layer of adsorbed molecules. This explains previously reported experimental and computational results and, for the first time, why these pattern changes are triggered by any perturbation of the charge density at the innermost layer of the electric double layer (voltage and composition changes, and vacancies at the electrode walls, among others), which may help tuning electrode-ionic liquid interfaces. Using Monte Carlo simulations we show that such structures can be reproduced by a simple two-dimensional lattice model with only nearest-neighbour interactions, governed by highly screened ionic interactions and short-range and excluded volume interactions. We also show that the results of our simulations are consistent with those inferred from the Landau-Brazovskii theory of pattern formation in self-assembling systems. The presence of these patterns at the ionic liquid graphene-electrode interfaces may have a strong impact on the process of ionic transfer from the bulk mixtures to the electrodes, on the differential capacitance of the electrode-electrolyte double layer or on the rates of redox reactions at the electrodes, among other physicochemical properties, and is therefore an effect of great technological interest.

Tables

Table 1. Table 1: Number of molecules used in the reported MD simulations.

System	IL molecules	Dissolved molecules
Pure IL	900	-
5% Water-IL	950	50
5% methanol-IL	950	50
5% ethanol-IL	950	50
10% \ceLiBF4-IL +
5% vacancies	900	100
10% \ceLiBF4-IL-IL+
8% vacancies	900	100

Table 2. Table 2: Values of the average fraction of surface covered by anions ( η 𝜂 \eta ) in mixtures of \ce [BMIM] \ce [BF4]. The vacancy percentages correspond to those in the positively charged graphene wall.

Mixture	$\eta$	Phase
Pure	0.34	Stripes
5% water	0.36	Molten Lamella
5% methanol	0.29	Stripes
5% ethanol	0.16	Hexagons
5% vacancies	0.19	Hexagons
8% vacancies	0.17	Hexagons
10% Li salt +
5% vacancies	0.30	Stripes
10% Li salt +
8% vacancies	0.19	Hexagons

Equations

$$E/k_{B}T=N_{pp}{I}_{pp}+N_{nn}{I}_{nn}+N_{np}{I}_{np},$$

$$\frac{N_{np}^2}{\left(z{N}_n-N_{np}\right)\left(z{N}_p-N_{np}\right)}=\exp \left(-2\beta \omega \right),$$

$$N_n=\frac{z\left(I_{pp}-I_{np}\right)+\left({\mu }_n-{\mu }_p\right)}{z\left(I_{pp}-2I_{np}+I_{nn}\right)}.$$

$${\mathcal{L}}_B\left[\varphi \right]=\int d\mathbf{r}$$

$$\frac{\beta V_4}{4!}{\left({\nabla }^2\varphi \left(\mathbf{r}\right)\right)}^2],$$

$$f(\varphi )=(A_2/2+\beta V_0){\varphi }^2+A_3{\varphi }^3/3!+A_4{\varphi }^4/4!$$

$${\mathcal{L}}_B[\varphi ]=\beta \Omega [{\bar{\rho }}_q+\varphi ]-\beta \Omega [{\bar{\rho }}_q].$$

$${\mathcal{L}}_B$$

$$+\frac{1}{2}\int d\mathbf{r}\int d\Delta \mathbf{r}\varphi (\mathbf{r})\beta V_c(\Delta \mathbf{r})g(\Delta \mathbf{r})\varphi (\mathbf{r}+\Delta \mathbf{r}),$$

$$A_{3}=\frac{1}{2\rho_{tot}}\Big{[}\frac{1}{(\rho-\bar{\rho}_{q})^{2}}-\frac{1}{(\rho+\bar{\rho}_{q})^{2}}\Big{]},$$

Full text

\SetWatermarkText

Author’s Postprint \SetWatermarkScale0.5 \SetWatermarkAngle45 \SetWatermarkColor[gray]0.8

\makeFNbottom

[TABLE]

††footnotetext: a* Grupo de Nanomateriales, Fotónica y Materia Blanda. Departamentos de Física de Partículas y Física Aplicada, Facultade de Física, Universidade de Santiago de Compostela, Campus Vida s/n, E-15782 Santiago de Compostela, Spain.††footnotetext: b Grupo Mesturas, Universidade da Coruña, Departamento de Física e Ciencias da Terra, Facultade de Ciencias, Campus da Zapateira, 15071 A Coruña, Spain. ††footnotetext: c Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warszawa, Poland. ††footnotetext: ∗ Corresponding author, to which correspondence should be addressed. E-mail: luismiguel.varela@usc.es*

1 Introduction

Ionic liquids (ILs) have attracted, and continue to attract, considerable research efforts both for theoretical reasons (the desire to understand the peculiar properties of these complex systems) and for their many potential applications, including as solvents for synthesis and catalysis, thermal fluids for sensible heat storage in thermosolar power systems, lubricants and, prominently, as advanced electrolyte materials for fuel cells, batteries or supercapacitors, among others (see, e.g., refs. 1, 2, 3, 4, 5, 6, 7, 8, 9). Beyond their behaviour at bulk level, the properties of ILs at interfaces are specially relevant since it is at the interface where redox reactions usually take place. Until now, most studies in this area have been focused on the oscillations of charge in the direction perpendicular to the interface and their effect in the adsorption of cosolvents at the electrode (see the recent review of Fedorov and Kornyshev 10 and references therein). However, despite it probably exerts a profound influence on charging dynamics, 11 the lateral structure of ILs nanoconfined in the adlayer has received much less attention. Merlet et al.12 and Rotenberg and Salanne13 reported molecular dynamics (MD) studies of the interface between buthylmethylimidazolium hexafluorophosphate and graphite and Au(111) electrodes at constant potential, showing the occurrence of in-plane structures and fluctuations beyond mean-field in the first adsorbed layer of the IL. The authors analysed voltage driven changes in the lateral charge distribution that induce transitions form ordered states into disordered nonstoichiometric states, leading to peaks in the capacitance. These configurations were similar to those experimentally detected by Liu et al.14 and more recently by Ebeling et al.15 Using atomic force microscopy, the latter authors investigated the lateral structure of the innermost layer of the IL-solid interface and the associated three-dimensional (3D) structure induced by molecule-substrate interactions, revealing quasi (4$\times$4)R0*∘* overlayer. On the other hand, Dudka et al.11 have used a statistical theory based on bipartite lattices and a 3D off-lattice Monte Carlo (MC) simulation of the structure and phase behaviour of superionic liquids in nonpolarized nanoconfinement, showing the existence of crystal-like (ordered) and homogeneous (disordered) phases separated by first or second order phase transitions. Also on the theoretical side, Limmer16 reported an effective field theory predicting a fluctuation-induced, first-order interfacial transition associated with spontaneous ordering of the ions in the layer nearest the electrode, which explains the “anomalous” behaviour of the differential capacitance of the IL-based electric double layer. The author showed that two opposing interactions are behind this first-order interfacial transition: short-range interactions between like species leading to the phase separation, and long-range electrostatic interactions between ionized species, leading to the neighbourhood of opposite charges. The combined action of these two interactions induces long-ranged order (nonvanishing charge for vanishing potential) in a slab parallel to the electrode whose width equals the bare correlation length. However, neither the theoretical analysis nor the MC simulations provide any direct evidence of the lateral patterns in the near-electrode layer. On the other hand, by means of MD simulations Docampo et al.17 have shown that the lateral structure of the first layer of anions close to a graphene surface undergoes a transition from ordered stripes to an ordered phase with hexagonal symmetry upon the addition of water. Hence, the authors pioneeringly showed the existence of concentration driven structural transitions in the adlayer structure, which were recently complemented by Gomez-Gonzalez et al.,18 who reported the occurrence of these transitions also in mixtures of ILs with monovalent and divalent salts at the graphene interface.

These structural changes at IL-graphene interfaces seem to be a quite general emergent pattern of nanoconfined ILs that could have a strong impact on the physical and chemical properties, such as the capacitance and the mechanical response, which could be of great importance for electrochemical devices, voltage- and composition-dependent lubricants and actuators. So, although the interface structure may depend on the nature of the electrode and the IL-electrode interaction, it is important to know if there are some general features, for which we will consider the particularly relevant IL-graphene interface, and, if this is so, to give a proper theoretical interpretation of that behaviour. Moreover, the elucidation of the actual role played by the surface on these patterns is another critical issue still to be appropriately addressed.

In this work we use MD simulations to show that the formation of ordered patterns at IL interfaces is caused not only by the addition of water or by changes in the electrode potential, but also by changes in the structure of the electrodes or by the addition of other cosolvents such as alcohols or inorganic salts; indeed, by any other perturbation which alters the balance of charge in the adlayer. Additionally, by means of NVT MC simulations we show that the main features of such structural behaviour can be reproduced by a simple two-dimensional (2D) lattice model with only nearest-neighbour interactions and report a schematic phase diagram. Moreover, we propose a theoretical explanation of this behaviour using the Landau-Brazovskii phenomenological theory, 19 which has been previously used to explain pattern formation in self-assembling systems. 20, 21, 22, 23, 24

The essential technical details of the computational methods used are sketched in the next section, our results are presented and discussed in the section 3, and finally we summarize our main conclusions.

2 Computational Details

MD simulations were carried out using Gromacs 5.025 in a NVT ensemble at room temperature (298 $K$), which was kept constant using a V-rescale26 thermostat. In order to avoid the system being trapped in a local metastable configuration, it was annealed up to 600 $K$ for 1 ns before the production run. The simulation system was formed by \ce[BMIM][BF4] (1-butyl-3-methylimidazolium tetrafluoroborate) with different cosolvents (water, methanol, ethanol and lithium tetrafluoroborate) between two graphene electrodes. The number of IL pairs was set to 900 for the pure IL and its mixtures with \ceLiBF4. For the other cases, it was set to 950. These numbers were chosen so the number of dissolved molecules was at least 50 in order to have reliable enough statistical averages of their magnitudes. The exact number of molecules used for each system is shown in Table 1.

All interactions were represented by means of the OPLS-AA force field.27 The electrodes were modeled as a rigid graphene sheet with a charge of $\pm 1\leavevmode\nobreak\ e\leavevmode\nobreak\ \textrm{nm}^{-2}$. The pristine graphene sheets were generated using the Visual Molecular Dynamics (VMD) software.28 The vacancy defects on the electrodes were introduced by randomly removing carbon atoms from the pristine graphene sheet and distributing their charge among their first neighbours. Furthermore, an additional constraint in order to avoid the appearance of isolated atoms was used: the first and second neighbours of a removed carbon atom were excluded from subsequent removals.

The initial configurations were generated using Packmol29 with a distance between the graphene electrodes ca. 10 nm. Their lateral sizes were then chosen in order to reproduce the average bulk density of the corresponding mixture in the center of the box. After that, the systems were annealed up to 600 K during 1 ns, equilibrated at room temperature for 10 ns and, finally, production runs of 5 ns with steps of 2 fs were performed.

On the other hand, MC simulations were performed for a model Coulomb system composed of anions and cations distributed in a triangular lattice of 400 nodes in the NVT ensemble using the standard Metropolis algorithm. In each simulation step, two nearest-neighbour ions are allowed to swap positions, the move being accepted according to the normal Metropolis algorithm. The simulations consisted of 250000 MC stabilization loops (where loop means a number of steps equal to the number of nodes) and a production run of 250000 additional loops. The interaction energy was calculated using a lattice model similar to that reported by Dudka* et al.*, 11 as described in detail in Sec. 3.

3 Results and Discussion

In order to address the challenges mentioned in the Introduction, we first performed MD simulations to analyse the adlayer structure of different mixtures of \ce[BMIM][BF4] confined between both pristine and defective (with vacancies) graphene electrodes. Specifically, we investigated the behaviour of the pure IL and its mixtures with water, ethanol or methanol confined between two graphene sheets with pristine structures. Moreover, we investigated both the pure IL and its mixtures with a Li salt (\ceLiBF4) when one of the graphene sheets has different percentages of vacancy defects, in order to test the influence of the reduction of the surface symmetry. We note that the presence of defects in graphene-like materials are sometimes unavoidable, but can also be intentionally added for specific practical purposes (see, e.g., refs. 30 and 31). [BMIM]\ce[BF4] was chosen because its cation is among the most commonly used and its anion has tetrahedral symmetry, which makes it specially suitable for structural studies. The vacancies were randomly distributed in the surface and the partial charge of the removed atoms was homogeneously distributed among their closest neighbours.

In order to analyse the lateral structure of the mixtures in the adlayers, we computed the surface number density of the species in the mixtures in the 5 Å  slab closest to the electrode in all cases, which corresponds to the first minimum of the density for mixtures with water,17 and, in general, is a good approximation for most mixtures. Illustrative results for a mixture of \ce[BMIM][BF4] and \ce[Li][BF4] next to the graphene cathode with different percentages of vacancy defects are displayed in Fig. 1. The MD results in gray show that the systems have an amorphous structure, so the Fourier transform (which gives the in-plane structure factor) is not very informative, at least for the electrode sizes achievable with our fully atomistic MD simulations. However, a more clear description of the patterns can be obtained by calculating the Minkowski functionals 32 after treating the MD images with Otsu’s algorithm, 33 which uses Fisher’s discriminant analysis to separate the image into two populations. It is one of the most frequently used methods to reduce a gray level image to a binary one by means of automatic clustering-based image thresholding. The core of the algorithm is finding the optimal assignment in a binary image of the pixels of a gray scale one. For that, the image is assumed to contain only two different populations of pixels, foreground and background, with a bi-modal histogram, so a threshold must be defined in order to univocally assign pixels of the original image to one population or the other. Afterwards, the threshold that minimizes the intra-class variance, i.e., the variance of the distribution of pixels in each side of the threshold, is to be found. This calculation can be equivalently and more efficiently done maximizing the inter-class variance, since the sum of pairwise squared distances is constant. The results of these calculations clearly show the presence of striped and hexagonal phases, similar to those reported in ref. 13, depending of the vacancy defect percentage. Similar patterns were obtained in all the studied cases as shown in Table 2.

An appropriate magnitude to characterise the two phases mentioned above is the surface fraction occupied by the anions, $\eta$, which takes values close to 0.3 for the striped phase and to 0.2 for the hexagonal phase. The values of $\eta$ for the different investigated systems were obtained using the previously described Otsu’s algorithm to automatically calculate the surface area occupied by the anions, which was later divided by the corresponding surface area of the graphene wall; the results are shown in Table 2. It must be mentioned that the mixture with water is an intermediate state between both structures (molten lamella34, 35), i.e., it corresponds with a limiting state between the striped and hexagonal patterns.

In order to see if these structural transitions have any actual influence on other physically relevant parameters, we analyzed the presence of \ceLi in the adlayer for mixtures of \ce[BMIM][BF4] and \ceLiBF4 as a function of the number of vacancies, and the results are shown in Fig. 2. There we observe that surface inhomogeneities can change the amount of cosolvent adsorbed at the positively charged wall: a clear change in the number of adsorbed \ceLi+ species in the mixture \ce[BMIM][BF4] + \ceLiBF4 is registered at 7% of vacancies in the graphene cathode, which is the concentration for which the transition from the striped to the hexagonal phase takes place. It should also be pointed out that, in general, the effect of these patterns does not end in the first layer of the IL, but it should determine the 3D structure of the electric double layer of the electrode, as pointed out by Ebeling et al.15

For further exploring, under the simplest possible conditions, the mechanism of pattern formation at IL interfaces and to build the phase diagram of the system, we performed NVT MC simulations using a 2D lattice model similar to that reported by Dudka et al. 11, but using only two species (anions and cations) since we restricted ourselves to the simplest case of pure IL with nearest-neighbour interactions only. This latter assumption is realistic for ILs since the high charge density in the studied nanoconfined systems gives rise to highly screened Coulomb interactions and, hence, to very short screening lengths, so, in practice, only nearest-neighbour interactions need to be taken into account.15 This kind of analytically tractable simplified lattice model has been quite commonly used for qualitative (and even quantitative) describing dense ILs, one-dimensional Coulomb systems as single file pores and ultrathin slits (see ref. 15 and references cited therein). Similar lattice calculations were performed by Limmer16 using a 3D charge frustrated Ising model. Of course, repulsive anion-anion and cation-cation interactions ($I_{nn}>0$ and $I_{pp}>0$) and attractive anion-cation interactions ($I_{np}<0$) define the physical region, being the total energy of the system

[TABLE]

where $N_{nn}$, $N_{pp}$ and $N_{np}$ are the number of anion-anion, cation-cation and anion-cation pairs, respectively, and $k_{B}$ is the Boltzmann constant. The species are considered to be distributed at the nodes of a triangular lattice with periodic boundary conditions. This lattice was chosen since its symmetry is that of the adsorption sites (hole sites) of the graphene sheet. However, as mentioned above, the actual in-layer structures observed for pristine graphene and for graphene with vacancies suggest that, expectedly, the results are quite independent of the actual electrode’s lattice. In each simulation step, one node is chosen randomly, as it is one of its first neighbours. Then, an attempt to swap the nodes is accepted according to the conventional Metropolis algorithm.

In the construction of the phase diagram of the system the interaction $I_{pp}$ was set to 0 and used as reference. The number of anions, $N_{n}$, was set to 0.3 times the number of total nodes, $N$. The various phases can be distinguished by the number of first neighbour anions of anions per node ($n_{nn}$). An average value of $n_{nn}$ close to 6 corresponds to a homogeneous phase, and a value of about 0 corresponds to a hexagonal phase, while a value close to 2 corresponds to a striped phase or to a random distribution of ions. In the latter case, $\langle n_{nn}\rangle=zN_{n}/N$, where $z$ is the coordination number, so $\langle n_{nn}\rangle=1.8$ for $N_{n}/N=0.3$ on a triangular lattice. The resulting phase diagram is shown in Fig. 3. As can be seen, three different phases are found. However, if we restrict ourselves to physically feasible values in common ILs (with positive $I_{nn}$ and negative $I_{np}$), only disordered and hexagonal phases are possible for this fraction of anion-occupied sites. We have investigated the dependence of the phase on temperature, and the results are shown in Fig. 4. It can be seen there that, at high temperatures, the stable phase is the disordered one and at low temperatures the hexagonal or the homogeneous phase (this latter for negative values of the attractive parameter).

The main features of the MC phase diagram can be understood by resorting to a statistical lattice theory. Specifically, it is possible to describe the patterns observed in the mixture using the well-known quasi-chemical approximation,36, 37 given by

[TABLE]

where $N_{n}$ and $N_{p}$ are the numbers of anions and cations, respectively, and $2\omega=-(I_{nn}+I_{pp}-2I_{np})$, i.e., the energy required to change an $nn$ pair and a $pp$ pair into two $np$ pairs. Formally, eq. 2 is the usual mass-action law for chemical reactions, which is the reason for the name quasi-chemical approximation, and goes beyond a purely random distribution introducing independent pairs.36 We note that, in eq. 2, $zN_{n}-N_{np}=2N_{nn}$ and $zN_{p}-N_{np}=2N_{pp}$. At high temperatures, eq. 2 gives $N_{np}=zN_{n}N_{p}/\left(N_{n}+N_{p}\right)$, which is the value for random mixing.37 On the other hand, at low temperatures there are two possibilities: if $\omega>0$, the two species are fully segregated ($N_{np}=0$); if $\omega<0$, the equilibrium distribution is that which maximizes the first term in eq. 2, with the constraints imposed by the lattice. The result can be written in terms of the number of particles of the minority component, $N_{m}$. For the triangular lattice, a system with $N_{m}=N/2$ (where $N=N_{n}+N_{p}$) is organized in stripes; with $N_{m}=N/3$ or $N_{m}=N/4$ in hexagons, and with $N_{m}=0$ a homogeneous configuration is found. For intermediate values the system will adopt mixed configurations of two of them.

On the other hand, the number of anions in the ground state of the mixture in the model of eq. 1 can be estimated minimizing $f=E-\left(\mu_{n}-\mu_{p}\right)N_{n}$, where $E$ is given by eq. 1 and $\mu_{n}$ and $\mu_{p}$ are the electrochemical potentials of the two ionic species. Assuming, for simplicity, a random distribution of ions the equilibrium condition $df/dN_{n}=0$ yields the expression

[TABLE]

This result shows that it is possible to achieve any concentration of anions by tuning the difference between the surface electrochemical potentials, independently of the temperature. Accordingly, it should be possible to control the surface pattern by carefully adjusting the electrochemical potentials, e.g., modifying the electrode’s potential. A scheme of the various possible 2D patterns as a function of the electrochemical potentials ($h$) is shown in Fig. 5. It is visible there that the patterns coincide with our MD and MC predictions.

The variety of systems in which our simulations show the presence of striped and hexagonal patterns indicates that we are probably in front of a general feature of surface-confined ILs, a question which requires a sound theoretical foundation. To address this issue, we have used the phenomenological Landau-Brazovskii (LB) model.19 According to it, systems with mesoscopic-scale inhomogeneities can be described by the functional

[TABLE]

where $\phi({\bf r})$ is the order parameter (OP). Typically, it is assumed that

[TABLE]

with $A_{2},A_{4}>0$. $A_{n}$ and $V_{n}$ are phenomenological parameters. At the mean-field (MF) level, the equilibrium OP corresponds to the minimum of ${\cal L}_{B}$. When $V_{2}<0$, the inhomogeneous structure is favoured by the second term in eq. 4, because in the presence of inhomogeneities ($\nabla\phi({\bf r})\neq 0$) this term leads to a decrease of the functional. When $V_{4}>0$, the inhomogeneous structure is disfavoured by the third term in eq. 4. In the case of $V_{2}<0$ and $V_{4}>0$, competition between these two terms leads to a finite length-scale of inhomogeneities, $2\pi/k_{b}$, with $k_{b}^{2}=-6V_{2}/V_{4}$.19 The effective field theory given by eq. 5 is, in some sense, a generalization of Limmer’s theory proposed in ref. 16, where the OP represents the local excess charge density, and the Landau functional for a phase-separating mixture is supplemented with the Coulomb interactions between the ionized species. In our LB theory the specific and Coulombic interactions are treated at the same footing, and the neighbourhood of opposite charges is favoured when $V_{2}<0$ and $V_{4}>0$. Moreover, we explicitly consider possible asymmetric solutions by taking into account the term $A_{3}\phi^{3}/3!$, neglected in ref. 16. As we show later, this term plays a crucial role in the symmetry of the ordered structure.

At high temperature the Brazovskii functional takes the minimum for $\phi({\bf r})=0$, i.e., for the disordered phase. Ordered patterns with spatially periodic $\phi({\bf r})$ can correspond to the minimum of ${\cal L}_{B}$ for sufficiently low $T$, since for decreasing $T$ (increasing $\beta$) the role of the last two terms in eq. 4 increases. In 2D systems, the stable ordered pattern corresponds to parallel stripes for relatively small values of $|A_{3}|$, or to a hexagonal arrangement of circular domains with $\phi({\bf r})>0$ or $\phi({\bf r})<0$ for negative or positive values of $A_{3}$ with intermediate magnitude, respectively. In addition to the ordered phases, a disordered phase is present for very large $|A_{3}|$. The values of $|A_{3}|$ at the transitions parallel stripes - hexagonal structure, and hexagonal structure - disordered phase are functions of the remaining parameters in eq. 4, and both increase with decreasing $T$.38

In the presence of fluctuations, there are some changes on the phase diagram described above. The stability region of the disordered phase enlarges (fluctuations destroy the order), and a striped phase with orientational and without translational order appears.34, 35 In this phase, defects destroy the translational order, but a preferred direction of the stripes remains (anisotropic structure). At sufficiently low $T$ the general sequence of the structures for increasing $|A_{3}|$: oriented stripes - hexagonal - disordered - is preserved beyond MF. The disordered phase is isotropic, but the correlation function exhibits oscillatory decay at relatively low $T$, signaling short-range periodic order.

Predictions of the LB theory agree with the patterns that are formed in thin magnetic films with competing ferromagnetic and dipolar interactions,39 in amphiphilic systems 20, 21, 22, 23, 24 and in colloidal systems with competing short-range attraction and long-range repulsion (SALR). 40, 41, 38, 42, 43, 44 In magnetic films, the OP is identified with local magnetization. In the latter two examples, $\phi({\bf r})$ corresponds to the local deviation from the space-averaged value of the concentration difference between the polar and apolar components, and to the volume fraction of particles, respectively.

The patterns predicted by the LB theory agree very well with the charge inhomogeneities observed in the layer near the graphene electrode and with the predictions of our lattice model. The sequence of phases agrees with that in Fig. 5, and the short-range order in the disordered phase agrees with the upper right snapshot in Fig. 3. In order to describe these inhomogeneities by means of the functional in eq. 4, the OP $\phi({\bf r})$ must be identified, as has also been previously done by Ebeling et al.,15 with the deviations of the local charge density, $\rho_{q}({\bf r})=e(\rho_{p}({\bf r})-\rho_{n}({\bf r}))$ from its space-averaged value, $\bar{\rho}_{q}=e(\bar{\rho}_{p}-\bar{\rho}_{n})$, where $e$ is the elementary charge and $\rho_{p}({\bf r})$ and $\rho_{n}({\bf r})$ denote the local volume fraction of cations and anions respectively (controllable by different means such as voltage, concentration of additives, etc.). In equilibrium, $\bar{\rho}_{q}$ is determined by the external conditions.

To apply the continuous theory eq. 4 to the structure formed at the microscopic scale, we adopt the concept of the mesoscopic volume fraction.45 Let us assume, for clarity, that the ions and the neutral molecules are roughly spherical and have nearly the same diameter $D$. We choose this size as the scale of coarse-graining or smearing of the microscopic density. The mesoscopic volume fraction of cations at the point ${\bf r}$, $\rho_{p}({\bf r})$, is equal to the fraction of the volume of the sphere with the center at ${\bf r}$ and diameter $D$ that is covered by the cations. $\rho_{p}({\bf r})$ decreases continuously from $1$ when the center of the cation is shifted from ${\bf r}$. The mesoscopic volume fractions of the anions and the neutral component are defined in an analogous way. Note that with this definition each microscopic state is characterized by continuous fields, so differential operators can be applied.

In order to predict which pattern will be stable under given thermodynamic conditions, we need to know the dependence of $A_{3}$ on the measurable quantities. Here we follow the analysis of ref. 42, where a link between the phenomenological theory and the statistical thermodynamics was established by a systematic coarse-graining procedure for the SALR system. We assume that in our case ${\cal L}_{B}$ is equal to the excess grand potential associated with lateral inhomogeneities in the layer near the graphene electrode. Next we assume that the neutral molecules and/or the ionic vacancies (if present) in the near-surface layer are homogeneously distributed. With this assumption, the sum of the volume fractions of the anions and the cations in the considered layer is approximately constant, and ${\cal L}_{B}$ depends only on the lateral charge distribution $\phi({\bf r})=\rho_{q}({\bf r})-\bar{\rho}_{q}$,

[TABLE]

The excess grand potential consists of the excess internal energy and the excess entropy contributions. The excess entropy, $S(\bar{\rho}_{q}+\phi,\rho)-S(\bar{\rho}_{q},\rho)$, can be Taylor expanded in terms of $\phi$, and within the local-density approximation we obtain

[TABLE]

where the first and the second term on the right hand side represent the entropy and the internal energy contributions respectively, with $V_{c}$ and $g$ denoting the interaction potential and the pair distribution function respectively. The excess energy can be approximated by the LB-type expression, as shown in ref 42. The parameters $V_{n}$ can be expressed in terms of the Coulomb interactions. When $g$ is such that $g(r)=0$ for $r<D$, then $V_{2}<0,V_{4}>0$, as in the LB theory. We assume that the entropy can be approximated by the entropy of mixing of the ions and the neutral molecules, $S/V=-k_{B}T\sum_{i=p,n,s}x_{i}\ln x_{i}$, where $x_{i}$ is the mole fraction of the $i$-th component in the considered volume $V$ and the subscript $s$ is for the neutral component. We also assume that the above form for the entropy is valid locally (i.e., in each region with the center at ${\bf r}$ and linear size $D$, with $x_{i}$ taking the local value $x_{i}({\bf r})$). In the simplest case of equal sizes of all the components, $x_{i}=\rho_{i}/\rho_{tot}$, where $\rho_{tot}=\sum_{i=p,n,s}\rho_{i}$. By assumption, $\rho_{s}=\rho_{tot}-(\rho_{p}+\rho_{n})$ is independent of the lateral position.

The parameters $A_{n}(\rho,\bar{\rho}_{q})$ in eq. 7 can be easily obtained. We get $A_{2},A_{4}>0$, and

[TABLE]

which can be positive or negative depending on the sign of $\bar{\rho}_{q}$. Note that when $|\bar{\rho}_{q}|\ll\rho$, then $A_{3}$ is very small and stripes are expected by the LB functional for a weakly charged electrode. When $\bar{\rho}_{q}<0$ (positively polarized electrode), then $A_{3}<0$, and the hexagonal pattern of the positive charges in the negatively charged background occurs for a range of $|A_{3}|$ depending on $T$. For $\bar{\rho}_{q}>0$ (a layer near the cathode) $A_{3}>0$, and the “negative” of the above pattern, i.e., circular domains of negative charge forming a hexagonal structure appear for some range of $A_{3}$. The second case (i.e., $\bar{\rho}_{q}>0$) corresponds to our simulations performed in the presence of the graphene cathode.

From eq. 8, it follows that $|A_{3}|$ increases when $|\bar{\rho}_{q}|$ increases with fixed $\rho$, or $\rho$ decreases with fixed $\rho_{q}$. It means that the transition between oriented stripes and hexagons can be induced by increasing the charge imbalance in the adlayer (increasing $|\bar{\rho}_{q}|$), and/or by adding a neutral component (decreasing $\rho$). Indeed, any physical perturbation which couples to the OP $\phi$ (differences in the charge density) is candidate to trigger structural transitions in the 2D confined IL close to the interface. This behaviour agrees with all simulation results. The change of the pattern from stripes to hexagons was observed in simulations when $|\bar{\rho}_{q}|$ was increasing, and/or $\rho$ was decreasing. In particular, in Table 1 hexagons are observed for large charge imbalances ($\eta<0.2$), while the oriented stripes are associated to small ones ($\eta>0.2$).

To conclude, it is also noteworthy that these patterns could be also seen as the time-averaged outcome of the 2D dynamics of strongly adsorbed IL mixtures confined in the region close to the electrode, dynamics which has been recently investigated by means of NMR relaxometry.46 On the other hand, striped and hexagonal structures have also been reported as Turing patterns driven by a reaction-diffusion model.47, 48 However, as we have shown, it seems that an equilibrium model is enough to explain the observed patterns at the electrochemical interface.34 The relation of the hereby reported structural patterns with the dynamics of ILs under confinement, as well as with the Turing-like mechanism and its possible equivalence to the formalism hereby reported, deserve further attention and will be developed elsewhere.

4 Conclusions

Our MD simulations show that, in agreement to previously reported works,12, 13, 15, 11, 17 there are, at least, two different possible patterns (striped and hexagonal) for ILs in the Stern layer formed at IL-graphene interfaces, corresponding to the disorder (homogeneous) and ordered (crystal-like) phases reported in the literature. We have shown that the formation of these 2D structures can be theoretically explained on very general thermodynamic grounds using the Landau-Brazovskii phenomenological theory, thus showing that they are the result of an equilibrium of forces in the adlayer. The fact that these structures can be replicated, in the case of a pure 2D IL, using a simple lattice model with only nearest-neighbour interactions shows that the most relevant interactions in the pattern formation in such dense ionic environments are short-ranged core interactions. Additionally, our results show that the structural patterns are also triggered by changes in the structure of the electrode (vacancies). Hence, we proved that, independently of its origin, any perturbation affecting the net balance of electric charge at the IL-graphene interface (described here by the order parameter ($\phi$(r)) can trigger transitions between the possible 2D patterns. Thus, given the potential influence of the lateral structure of the adlayer on the charging dynamics, our findings may open new ways to induce structural transitions at those interfaces, which can be of technological interest.

Of course, molecular details may have a significant effect on the results at the quantitative level. In particular, a more accurate approximation for the entropy in specific systems will lead to some corrections to our expression for $A_{3}$ in eq. (8). Thus, the values of the average charge $\bar{\rho}_{q}$ and concentration $\rho$ of ions at the transitions between the different structures will be somewhat different. However, our results, as well as others previously reported, indicate that the sequence of structures should remain the same for a large class of systems, suggesting a universal behaviour in some way analogous to the universal sequence of gas, liquid and solid phases that occurs in majority of substances, independently of the molecular details. These universal features are predicted by the generic model that successfully describes many other systems exhibiting pattern formation, including systems composed of molecules as complex as lipids or block copolymers.22, 23

Finally, we must mention that although in this work we have focused entirely on 2D pattern formation in ILs confined between graphene walls, it is evident, as indicated above, that the formation of these 2D structures greatly influences (and is influenced by) the 3D structure of the electric double layer. This, indeed, has been reported by Ebeling et al.,15 who found a 3D zincblende-type ionic crystal-like organization in the interface between propylammonium nitrate and highly ordered pyrolytic graphene. The theoretical study of these 3D effects would add new insights into the behaviour of these systems as it would the analysis of packing fraction effects on the transitions between ordered and disordered phases in ILs in such confined geometries. Work in this direction is now in progress.

5 Acknowledgement

Prof. R. M. Lynden-Bell is greatly acknowledged for her most useful comments, suggestions and careful reading of our manuscript. The authors wish to thank the financial support of the Spanish Ministry of Economy and Competitiveness (Projects MAT2014-57943-C3-1-P, MAT2014-57943-C3-3-P and FIS2012-33126). Moreover, this work was funded by the Xunta de Galicia (AGRUP2015/11 and GRC ED431C 2016/001). All these research projects were partially supported by FEDER. H. M.-C. and J. M. O.-M. thank the Spanish Ministry of Education for their FPU grants. Facilities provided by the Galician Supercomputing Centre (CESGA) are also acknowledged. Funding from the European Union (COST Actions CM1206 and MP1303) is also acknowledged. This project has also received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 734276, and from the National Science Center grant 2015/19/B/ST3/03122.