Mean field approach to first and second order phase transitions in ionic liquids

TL;DR

This paper uses a mean-field model to analyze phase transitions in ionic liquids, revealing how nano-structures form and influence the electrical double layer, with implications for their self-assembly and electrochemical properties.

Contribution

It introduces a mean-field framework based on Onsager's relations to describe phase transitions and nano-structure formation in ionic liquids, linking them to experimental observations.

Findings

Bulk nano-structures arise via type-I and II phase transitions.

Type-I phase transition is induced by thermal fluctuations and depends on ion ratio.

Insights relate to decaying charge layers observed in experiments.

Abstract

Ionic liquids are solvent-free electrolytes, some of which possess an intriguing self-assembly property. Using a mean-field framework (based on Onsager's relations) we show that bulk nano-structures arise via type-I and II phase transitions (PT), which directly affect the electrical double layer (EDL) structure. Ginzburg-Landau equation is derived and PT are related to temperature, potential and interactions. The type-I PT occurs for geometrically dissimilar anion/cation ratio and, surprisingly, is induced by perturbations on order of thermal fluctuations. Finally, we compare the insights with the decaying charge layers within the EDL, as widely observed in experiments.

Equations

$$\mathcal{E}=\int \text{d}\mathbf{x}\{k_{B}T\left[p\ln \frac{p}{\bar{c}}+n\ln \frac{n}{\bar{c}}\right]+\frac{\bar{c}\beta }{p_{\max }{n}_{\max }}np+$$

$$\dfrac{\partial}{\partial t}\left(\begin{array}[]{c}p\\ n\end{array}\right)=\frac{M}{\bar{c}}\nabla\cdot\left[pn\left(\begin{array}[]{rr}\gamma&-1\\ -1&1/\gamma\end{array}\right)\nabla\left(\begin{array}[]{c}\delta{\mathcal{E}}/\delta p\\ \delta{\mathcal{E}}/\delta n\end{array}\right)\right],$$

$${\partial }_{t}p$$

$${\nabla }^2\varphi$$

$$\left(\begin{array}{c}p\\ \varphi \end{array}\right)=\left(\begin{array}{c}{p}_0+{ϵ}^{1/2}{p}_1+ϵp_2+{ϵ}^{3/2}{p}_3+...\\ {ϵ}^{1/2}{\varphi }_1+ϵ{\varphi }_2+{ϵ}^{3/2}{\varphi }_3+...\end{array}\right),$$

$$p_1=A(ϵt,\sqrt{ϵ}x)\exp (i{k}_{c}x)+\text{complex conjugate},$$

$${\partial }_{t}q=$$

$$+\frac{1}{1+\gamma }{\partial }_x^2\varphi ,$$

$${\partial }_x^2\varphi =$$

$$s=-1+\left[\chi \frac{2\gamma }{{\left(1+\gamma \right)}^3}-\frac{1+{\gamma }^2}{\gamma \left(1+\gamma \right)}\right]k^2-\sigma \frac{2\gamma }{{\left(1+\gamma \right)}^3}{k}^4.$$

$${\sigma }_c=\frac{2{\gamma }^2\chi -(1+{\gamma }^2)(1+\gamma {)}^2{)}^2}{8{\gamma }^3(1+\gamma {)}^3},k_c^2=\sqrt{\frac{(1+\gamma {)}^3}{2{\sigma }_c\gamma }}.$$

$$\chi >{\chi }_k=\frac{\left(1+{\gamma }^2\right){\left(1+\gamma \right)}^2}{2{\gamma }^2},$$

$$\frac{\partial A}{\partial t}=\frac{{\sigma }_c-\sigma }{{\sigma }_c}A+\alpha |A{|}^{2}A+\frac{4}{k_c^2}\frac{{\partial }^{2}A}{\partial x^2}$$

$$p\simeq \frac{1}{1+\gamma }+2\sqrt{\frac{\sigma -{\sigma }_c}{{\sigma }_c\alpha }}\cos k_{c}x,$$

$$\varphi \simeq 2\frac{1}{p_0{k}_c^2}\sqrt{\frac{\sigma -{\sigma }_c}{{\sigma }_c\alpha }}\cos k_{c}x.$$

$${\chi }_s$$

$$\sigma -{\sigma }_c>\left|\frac{\alpha {\sigma }_c{p}_0^2{k}_c^4}{4}\right|\sim \frac{32}{\chi -4}+\mathcal{O}\left(\gamma -1\right).$$

Full text

Molten salts are comprised of large and asymmetric anions and cations, with their molecular structure consisting of a charged ion attached to a hydrophilic or hydrophobic functional group. With significant charge delocalization and irregular geometries, the ions do not readily form a tightly-bound lattice and remain liquid even at room temperatures and in the absence of any solvent Dupont (2011), hence the name “ionic liquids” (ILs). Their tunable molecular structure enables the tailoring of ILs to a large number of applications Plechkova and Seddon (2008); Silvester and Compton (2006); Wishart (2009); Zhou et al. (2009); Brandt et al. (2013); Li et al. (2014); MacFarlane et al. (2014); Armand et al. (2009), e.g., batteries, supercapacitors, dye-sensitized solar cells, lubricants and nanoparticle syntheses, where they are advantageous due to their high charge density, low-volatility, and high chemical, thermal and electrochemical stability.

It is the amphiphilic-type structure, however, that gives ILs another intriguing property- the ability to self-assemble, see Yochelis et al. (2015) and references therein. IL molecules spontaneously form bicontinuous, hexagonal or lamellar phases, see Hayes et al. (2015) and references therein, in a fashion similar to the morphologies of block copolymers and liquid crystals. The bulk nano-structure effects not only the mechanical and transport properties Erbas and de la Cruz (2016), but also the electrical double layer (EDL) structure and thus charge transfer properties Endres et al. (2010); Uysal et al. (2013); Rotenberg and Salanne (2015); Fedorov and Kornyshev (2014). Obtaining insight into the emergence of nanostructure is therefore essential for the integration and control of ILs in scientific and industrial applications Plechkova and Seddon (2008); Silvester and Compton (2006); Wishart (2009); Zhou et al. (2009); MacFarlane et al. (2014); Armand et al. (2009); Fedorov and Kornyshev (2014).

It was recently proposed that the coupling of bulk and EDL morphologies may be portrayed as the result of competition between short-range intermolecular interactions (e.g., hydrogen bonds, solvation interactions, steric effects) and the long-range Coulombic interactions Gavish and Yochelis (2016), where the contribution of the former may be incorporated via the Cahn–Hilliard theory of phase separation. This nonlocal Cahn–Hilliard framework has successfully reproduced the emergence of bulk nano-morphologies that have been observed experimentally Hayes et al. (2015); Yochelis et al. (2015), while providing further insights into electrokinetic phenomena in ILs, namely in the form of transient currents. However, the methodology was based on geometrically identical ions, which is clearly an abstraction of any real-world system Hayes et al. (2015). Furthermore, recent experiments have implied that phase transitions at the EDL should be attributed to electrode polarization Rotenberg and Salanne (2015). However, it is difficult to determine solely on the basis of experiments and atomistic simulations whether phase transitions in EDL structure depend only on the specific electrode polarization or whether there is an additional competing mechanism originating in the bulk Yochelis (2014). Isolating the separate influences of the bulk and interface properties is challenging as the coupling is inherent in the system Gavish and Yochelis (2016), and measurements of the EDL alone may lead to erroneous conclusions about the bulk, namely, that decaying charge layering in the EDL immediately implies an unstructured bulk Bazant et al. (2011).

In this Letter, we show that nano–patterning in ILs is a result of first and second order phase transitions. To this end, we extend the symmetric nonlocal Cahn-Hilliard mean field theory Gavish and Yochelis (2016) to account for ion size asymmetry; the new formulation also relies on Onsager theory. Then, using weakly nonlinear analysis, we derive a Ginzburg–Landau equation and perform numerical continuation to distinguish between the types of phase transitions that account for the emergence of nano-morphology both in the bulk and in the EDL. Particularly, we show that the formation of ordered layers near the solid surfaces is induced by naturally occurring energy fluctuations, while the observed decaying spatial charge oscillations may be attributed to an isotropic bulk morphology (independent of the phase transition type), rather than an unstructured bulk.

We start by considering a fully dissociated IL of monovalent anions ($n$) and cations ($p$), where their molar concentrations range as $0\leq n\leq n_{\max}$ and $0\leq p\leq p_{\max}$, respectively. The free energy of such an asymmetric system in the mean-field framework is generally comprised of the phase separation and Coulombic components Gavish and Yochelis (2016):

[TABLE]

where, $\bar{c}$ is the harmonic average density of ions, $k_{B}$ is the Boltzmann constant, $T$ is the temperature, $\beta$ is the interaction parameter for the anion/cation mixture and has units of energy, $E_{0}\kappa^{2}/4$ is the interfacial energy coefficient where $E_{0}$ has units of energy and $\kappa$ has units of length, $\phi$ is the electric potential, $q$ is the elementary charge, and $\epsilon$ is the permittivity. The first three terms in (Mean field approach to first and second order phase transitions in ionic liquids) correspond to the Flory–Huggins functional, which for $\beta>\beta_{c}=4k_{B}T$ takes the form of a double well potential, thus driving phase separation Huggins (1942); Flory (1953); Cahn and Hilliard (1958), while the last two terms of (Mean field approach to first and second order phase transitions in ionic liquids) are related to electrostatics,i.e., Poissons’ equation.

The phase transitions to bulk morphology are related to spatiotemporal symmetry–breaking of a uniform bulk concentration. To derive the dynamic equations for the charge densities, we exploit the Onsager framework Onsager and Fuoss (1932); Hilliard (1970)

[TABLE]

where $M$ is the diffusion coefficient De Gennes (1980); Brochard et al. (1983); Psaltis and Farrell (2011) and $\gamma=p_{\max}/n_{\max}>0$ is the packing ratio of the ions affected by their geometry; $\gamma=1$ corresponds to the symmetric (1:1) case where $p_{\max}=n_{\max}=1/2$ Gavish and Yochelis (2016). For the sake of analysis, we introduce non-dimensional variables and parameters $\tilde{p}=p/p_{\max}$, $\tilde{n}=n/n_{\max}$, $\bf{\tilde{x}}=\bf{x}/\lambda$, $\tilde{t}=t/\tau$, $\tilde{\phi}=q/\left(k_{B}T\right)\phi$, $\lambda=\sqrt{\gamma\epsilon k_{B}T/\left[\left(1+\gamma\right)q^{2}{c}_{\max}\right]}$, $\tau=\lambda^{2}/\left[(1+\gamma)K_{B}T\right]$, and obtain (after omitting the tildes) the final dimensionless equations for the asymmetric ILs:

[TABLE]

with $\partial_{t}n=-\gamma^{-1}\partial_{t}p$. The parameter $\sigma=\frac{E_{0}}{2k_{B}T}\frac{\kappa^{2}}{\lambda^{2}}$ depends on the ratio of the strengths of short and long-range interactions and reflects the competition between them, and $\chi=\beta/(k_{B}T)$ is the Flory parameter. For consistency with traditional dimensionless analysis, we have chosen to scale space by the Debye-like scale $\lambda$ while noting that this choice does not reflect a typical electric screening length, as for dilute electrolytes.

The emergence of spatially inhomogeneous solutions for (3), see, e.g., Figure 1(b), is associated with the finite wavenumber instability about a uniform solution as in the symmetric ($\gamma=1$) case Gavish and Yochelis (2016). For the sake of analysis, we first consider an infinite domain in one space dimension, using $\sigma$ as a control parameter. We then proceed to demonstrate the validity of the results on a finite domain $x\in(-L,L)$ as well, using no-flux (i.e., Neumann) boundary conditions for $p$, and fixed potential (i.e., Dirichlet) boundary conditions $\phi(x=\pm L)=\pm V/2$, where $V$ is the applied potential difference.

Employing the weakly nonlinear theory (multiple time scale analysis), we expand

[TABLE]

where $p_{0}=1/(1+\gamma)$ is the uniform cation concentration, $\epsilon\ll 1$ measures the distance from the instability onset, and following Shiwa (1997); Gavish et al. (2016)

[TABLE]

and $\phi_{1}=p_{1}/(p_{0}k_{c}^{2})$. The complex amplitude $A$ varies slowly in time, and $k_{c}$ is the critical wavenumber at the instability onset $\epsilon:=\left(\sigma-\sigma_{c}\right)/\sigma_{c}=0$. Both quantities are obtained via the linear stability analysis about the uniform state $(p,\phi)=(p_{0},0)$ to periodic perturbations Cross and Hohenberg (1993). Linearizing (3) around $(p_{0},0)$, i.e., setting $p=p_{0}+q$ and retaining linear terms in $(q,\phi)$ yields

[TABLE]

and substituting the second into the first equation yields the dispersion relation of finite wavenumber instability type Gavish and Yochelis (2016)

[TABLE]

$s$ is the temporal growth rate of periodic perturbations associated with non vanishing wavenumbers $k$, i.e., $s(k)>0$ implies linear instability. Varying $\sigma$, the instability onset is obtained by seeking $\sigma=\sigma_{c}$ (the phase transition point) and $k_{c}$ such that $s(k=k_{c},\sigma=\sigma_{c})=0$, $s(k\neq k_{c})<0$ and $\frac{\mathrm{d}s}{\mathrm{d}k}(k_{c})=0$, yielding

[TABLE]

Notably, although the model equations are gradient the coupling to Poisson equation damps the large scale mode Shiwa (1997); Gavish et al. (2016); Gavish and Yochelis (2016), i.e., $s(0)<0$. Since $k^{2}_{c},\gamma>0$, the instability may exist only for

[TABLE]

as identified by the shaded region in Fig. 1. For the symmetric case ($\gamma=1$) we retain the result $\chi_{k}=\chi_{c}\equiv\beta_{c}/(k_{B}T)=4$ from Gavish and Yochelis (2016).

Slightly above the instability onset, $\epsilon\ll 1$, the dispersion relation admits a finite band of wavenumbers around $k_{c}$ for which $s>0$. However, direct numerical integration shows that at $\epsilon=-0.001$ and $\gamma=2$, we find two distinct cases: for $\chi=8.3$ the solutions correspond to small amplitude periodic solutions (see dark line in Fig. 1(b)), while for $\chi=7$ the solutions are highly nonlinear (as evidenced by the spike-like shape of the light line in 1(b)) and of large amplitude (relative to $\epsilon$). This behavior indicates two distinct distinct phase transitions mechanisms, which can be analyzed using the amplitude equation for $A$. Thus we substitute (4) into (3), and after imposing the solvability condition at order $\epsilon^{3/2}$ obtain

[TABLE]

where $\alpha=(1+\gamma)^{2}/\gamma+k_{c}^{2}(2\sigma_{c}k_{c}^{2}-2\chi)/(1+\gamma)+\frac{2(\gamma^{3}-1)k_{c}^{2}}{9\gamma}\left[\frac{\gamma^{2}-1}{2\gamma}-\left(\frac{1-\gamma}{\gamma}+\frac{2\chi(1-\gamma)}{(1+\gamma)^{2}}\right)k_{c}^{2}+7\frac{2\sigma_{c}(1-\gamma)}{(1+\gamma)^{2}}k_{c}^{4}\right].$ The Ginzburg–Landau partial differential amplitude equation (7) can be used to approximate large scale modulations of the basic pattern $\cos(k_{c}x)$ near the (phase transition) onset Shiwa (1997); Gavish et al. (2016) but here for simplicity we assume $A$ to be independent of space. Then

[TABLE]

From the right hand sides of (8) it is evident that the signs of $\sigma-\sigma_{c}(\gamma,\chi)$ and of $\alpha(\gamma,\chi)$ differentiate between two types of periodic solutions: (i) The super–critical (phase transition of type-II) bifurcation when $\sigma<\sigma_{c}$ and $\alpha<0$, with the amplitude scaled as $\sqrt{|\sigma-\sigma_{c}|}$ (see Fig. 2(a)), and (ii) The sub–critical (phase transition of type-I) bifurcation when $\sigma>\sigma_{c}$ and $\alpha>0$, with the amplitude exhibiting a discontinuous jump to amplitudes $\mathcal{O}(1)$ for $\sigma<\sigma_{c}$. Typical solutions of both cases for $\sigma\lesssim\sigma_{c}$ obtained via direct numerical integration are depicted in Fig. 1(b), while the transition between the two types is obtained by solving for $\alpha(\chi,\gamma)=0$, which yields

[TABLE]

as depicted by the dashed line in Fig. 1(a).

The branches of bifurcating solutions beyond the weakly nonlinear limit are computed in both cases with numerical continuation method Uecker et al. (2014), as depicted in Fig. 2. While in the super-critical case (type II) the bifurcating solutions are stable (Fig. 2(a)), in the sub-critical case (type I) these solutions are linearly unstable until reaching a fold, $\sigma_{c}^{SN}$, whereon they change direction and grow toward the instability onset of the homogeneous state. However, as in the Ohta-Kawasaki case Gavish et al. (2016), the large amplitude periodic solutions belonging to a branch portion after the fold ($\sigma<\sigma_{SN}$, see figure) are stable on $2L=\lambda_{c}=2\pi/k_{c}$ (exhibiting hysteresis) but become Eckhaus unstable on large domains; due to applications our interest here is in domains that are larger than the wavelength of the bifurcating states, i.e., $L\gg\pi/k_{c}$. Consequently, the hysteresis region $\sigma_{c}<\sigma<\sigma^{SN}_{c}$, between the $k_{c}$ solutions and the uniform state is destroyed on large domains. Nevertheless, direct numerical integration shows that there are additional periodic solutions that emerge and form hysteresis for $\sigma_{c}<\sigma<\sigma_{L}^{SN}$ (Fig. 1(b)). These solutions belong to one of the secondary branches of periodic solutions (denoted by $k_{L}$) that emerge below the onset, i.e., at $\sigma<\sigma_{c}$ as shown by thick line in Fig. 2(b).

From a mathematical point of view, bistability of the uniform and periodic state persists on finite domains. However, in the physical context, perturbations are related to energy fluctuations which operate at the scale of $k_{B}T$ or $\phi\sim 1$ in dimensionless units. The size of fluctuations should be compared to the amplitude of the unstable parts of the subcritical branches. Setting $\phi\approx 1$ in (8b) yields

[TABLE]

The numerical continuation in Fig. 2(b) shows that in the parameter regimes considered here the folds are too close to $\sigma_{c}$ for (10) to be fulfilled. Indeed, even a weak perturbation $\phi=V$ with $|V|\ll 1$ at a boundary will result in a large–amplitude pattern after a temporal transient, as shown in Fig. 3. In other words, in a physicochemical system with the required temperature and anion/cation asymmetry for first order phase transitions, energy fluctuations will drive the formation of nano-structures for $\sigma_{c}<\sigma<\sigma_{L}^{SN}$, even though the homogeneous state is mathematically stable in this domain.

Finally, the large amplitude periodic patterns in 1D take on a much richer form in 3D, such as laminar, bicontinuous, and sponge-like morphologies Cross and Hohenberg (1993); Choksi and Ren (2005); Nyrkova et al. (1994). Thus, the isotropic neutral bulk region is averaged upon projection on the parallel axis in between the electrodes to $(p,\phi)=(p_{0},0)$, which is also consistent with the electroneutrality. Near the electrodes the isotropy is broken, and the periodic states align with the surface Gavish and Yochelis (2016), with the alignment more pronounced closer to the electrode. Consequently, experimental measurements such as atomic force microscopy (AFM) Atkin et al. (2014), atomic force apparatus Perkin (2012) and X-ray reflectivity Uysal et al. (2013), which are essentially 1D methods, observe only charge layering. We emphasize that while charge layering is an aspect phase transition of type-I in the bulk and loss of isotropy near the electrodes, it does not automatically imply bulk nano-structure. The distinction between structured and unstructured bulk may be further advanced, for example, using time dependent measurement which shows different behavior in each case Gavish and Yochelis (2016). In higher spatial dimensions the 1D periodic bulk structure in takes the universal forms of labyrinthine and hexagonal patterns. The former may be expected in the region of type-II phase transitions for nearly–symmetric cases ($\gamma\approx 1$) and emerges via a zigzag instability Gavish and Yochelis (2016), while the latter emerges via a sub-critical bifurcation Golovin and Nepomnyashchy (2006). In both cases, similarly to the first order phase transition region, perturbations on the order of $\phi=1$ will cause the emergence of periodic morphology that will average to electroneutrality in the bulk.

To conclude, we have demonstrated the impact of cation/anion size asymmetry on ionic liquid self-assembly in both the bulk and the electrical double layer regions. Importantly, the framework is consistent with the second law of thermodynamics through the nonlocal Cahn-Hilliard approach that is based on Onsager’s relations Gavish and Yochelis (2016). The nano-morphology is shown to arise first in the bulk via phase transitions of types I and II under certain conditions, following the derived universal Ginzburg-Landau equation (7). In both types, near the transition threshold, $\sigma_{c}$, the bulk is sensitive to energy fluctuations of about $k_{B}T$ or respectively $25mV$, with the resulting formation of large amplitude nano-morphologies even in the type I case. Since these phase transitions are of universal nature Cross and Hohenberg (1993); Golovin and Nepomnyashchy (2006), in higher space dimensions the morphologies that form are of isotropic nature (e.g., labyrinths and hexagons) and thus, upon averaging about the parallel plane to the electrodes, correspond to a uniform density with a vanishing electrical field, i.e., electroneutrality in 1D context. The second effect of an isotropic structured bulk is its subtle impact on the charge layering within the EDL region Yochelis et al. (2015); Hayes et al. (2015). Since the bulk isotropy breaks down near the solid surface, the patterns tend to align with electrode orientation. Upon averaging about the plane normal to the electrode one naturally observes charge layers near the electrode that gradually lose their orientation toward the bulk. Consequently, spatially decaying oscillations are observed in the direction normal to the electrode, yet they should not be misinterpreted as evidence for absent bulk structure. As many applications Hayes et al. (2015); Fedorov and Kornyshev (2014), e.g., energy- and lubrication-related, exploit and depend on mass transport and charge transfer properties, the framework developed here offers a new perspective regarding the interpretation of empirical observations, as well as opportunities for the enhancement of device efficiency Plechkova and Seddon (2008); Silvester and Compton (2006); Wishart (2009); Zhou et al. (2009); Brandt et al. (2013); Li et al. (2014); MacFarlane et al. (2014); Armand et al. (2009) in terms of conductivity, structural integrity (rheology), and electrochemical reactions.