Classical Density Functional Theory of Ionic Liquids

The current technological revolution based on ionic liquids (ILs) is driven by their unique properties. Though ILs have been known to scientists for close to a century (Walden , 1914), they are now the focus of intense activity mainly because of their promise for “environmentally friendly” applications. Some actual and potential uses of ILs include: specific solvents for heterogeneous and homogeneous catalysis; selective solvents for removal of heavy metal contaminants; electrolytes in various electrochemical processes and devices; and as dispersive agents for stabilization of nanoparticles. In all of these applications, the structural properties of ILs and their mixtures in the bulk phase and at interfaces are crucial to their performance (Maier et al. , 2010). In electrowinning processes, such as aluminium refining, metal ion speciation is a crucial consideration in determining the population of electroactive components (Rocher et al. , 2009) and consequently influences the pertinent dynamical processes that lead to oxidation and reduction. These populations are ultimately governed by the free energy of the IL solution, which in turn depends upon the mutual interactions between the IL molecules and the added components. It is now clear that fundamental theoretical studies of structure-property relationships in ILs will provide important new insights, which will assist this kind of research. Not only in the choice and design of ILs for a specific application, but also for scenario based modelling of IL processes. For example, theoretical studies of the structure of ILs at charged surfaces are important for understanding the properties that they impart to electrodeposition processes, including factors such as ionic speciation and mass transport of electroactive species (MacFarlane et al. , 2010). As emphasized by Kornyshev in a recent review article (Kornyshev , 2007), applications of ILs at electrified interfaces to energy-storage systems, electrowetting devices or nanojunction gating media will be strongly promoted by a deeper understanding of the structure and properties of the interfacial double layer. Classical Density Functional Theory of Ionic Liquids


Introduction
The current technological revolution based on ionic liquids (ILs) is driven by their unique properties.Though ILs have been known to scientists for close to a century (Walden , 1914), they are now the focus of intense activity mainly because of their promise for "environmentally friendly" applications.Some actual and potential uses of ILs include: specific solvents for heterogeneous and homogeneous catalysis; selective solvents for removal of heavy metal contaminants; electrolytes in various electrochemical processes and devices; and as dispersive agents for stabilization of nanoparticles.In all of these applications, the structural properties of ILs and their mixtures in the bulk phase and at interfaces are crucial to their performance (Maier et al. , 2010).In electrowinning processes, such as aluminium refining, metal ion speciation is a crucial consideration in determining the population of electroactive components (Rocher et al. , 2009) and consequently influences the pertinent dynamical processes that lead to oxidation and reduction.These populations are ultimately governed by the free energy of the IL solution, which in turn depends upon the mutual interactions between the IL molecules and the added components.It is now clear that fundamental theoretical studies of structure-property relationships in ILs will provide important new insights, which will assist this kind of research.Not only in the choice and design of ILs for a specific application, but also for scenario based modelling of IL processes.For example, theoretical studies of the structure of ILs at charged surfaces are important for understanding the properties that they impart to electrodeposition processes, including factors such as ionic speciation and mass transport of electroactive species (MacFarlane et al. , 2010).As emphasized by Kornyshev in a recent review article (Kornyshev , 2007), applications of ILs at electrified interfaces to energy-storage systems, electrowetting devices or nanojunction gating media will be strongly promoted by a deeper understanding of the structure and properties of the interfacial double layer.

Other approaches
ILs are fascinating from a theoretical standpoint.In order to treat ILs, a theory must account for both electrostatic and steric correlations.Attempts to develop such a theory have appeared in the literature.These include the lattice based model by Kornyshev and co-workers (Fedorov et al. , 2008;Kornyshev , 2007).In another lattice approach by Lauw and coworkers (Lauw et al. , 2009), both cations and anions were modelled by a star-like architecture, with each monomer subject to an independently specified dielectric response.Under certain circumstances, these models were able to qualitatively reproduce the so-called "camel-shaped" differential capacitance, C D , profile observed in imidazolium ILs at carbon electrodes (Islam et al. , 2009;Lockett et al. , 2008a).This behaviour is an indication of subtle structuring effects at charged surfaces and is a stern test for theory.Although camel-shaped C D profiles are predicted, qualitative discrepancy against experimental results remains.This is not surprising since these crude theories lack a proper description of electrostatic correlations despite accounting for steric contributions.Charge-charge correlations have been shown to be important in interfacial phenomena of highly coupled electrostatic fluids.As with dispersion forces, ionic correlations are always attractive and lower the total energy.Thus, we expect that they will also play a role in establishing the camel shape of the C D .That this is the case, is seen from both simulations (Lamperski & Zydor , 2007;Vatamanu et al. , 2010) and modified Poisson-Boltzmann treatments (Lamperski et al. , 2009) of a restricted primitive model (RPM) electrolyte.These investigations show that the RPM displays a camel-shaped C D at low temperatures and densities.However, electrostatic correlations do not tell the whole story for ILs.That is, for the RPM the C D becomes bell-shaped and exhibits a single maximum at densities and temperatures typical for ILs.This suggests that the simple RPM is an inadequate model for ILs and other factors, such as dispersion forces and molecular structure (Fedorov et al. , 2010;Trulsson et al. , 2010), play an important role in determining their properties.In this work, we will use a model for ILs wherein the molecular ions interact via dispersion, as well as electrostatic, forces and the cation has an oligomeric structure consistent with imidazolium ions.This model has also been adopted by us in previous recent publications (Forsman et al. , 2011;Trulsson et al. , 2010), and employs the methods of polymer free energy density functional theory (DFT) (Woodward , 1991).The result will be a simple, albeit approximate DFT for ILs.We have already introduced the salient features of the theory in a previous publication (Forsman et al. , 2011).Here, we will present a much more thorough treatment, emphasizing the relation to theories of simpler fluids.We will also extend the scope of the study with some emphasis on temperature dependence and the effects of ion-specific adsorption.The IL-DFT that we will discuss here is an off-lattice theory, whose accuracy can be straightforwardly tested against computer simulations of the same molecular model.We will begin with a discussion of relevant earlier work, in order to introduce the essential mechanistic components that will be combined in our density functional theory for ionic liquids.

Classical density functional theory (DFT) of simple fluids
The ideal gas law is known and used by most physicists and chemists but a more complete description of a fluid including its ability to condense into a liquid requires a treatment of the interactions between its particles.The simplest realistic model of a molecular fluid assumes pairwise spherically symmetric interactions by a pair-potential φ(r) which has a repulsive (positive) inner region characterised by a length scale σ and an attractive tail that defines an (1) We shall let this familiar pair-potential illustrate our simple fluid in the discussion below.
Even for a simple fluid, the the analysis quickly becomes quite complicated if we attempt to describe configurational states by explicitly locating all the particles as in simulation methods, or if we try to resolve the coupled hierarchy of correlation functions as in rigorous liquid state theory.Here we shall instead consider an approach based on deliberate and efficient approximations seeking to represent the essential features while hiding the less important details of the fluid properties.The most noteworthy early attempt in this direction was taken by J. van der Waals (van der Waals , 1873) who 140 years ago realized that the interactions in molecular fluids introduced two mechanisms: i) an excluded volume effect due to the size of the molecules and ii) a binding energy due to the attraction between molecules.His equation of state relating pressure P to volume per particle v and absolute temperature T reads It incorporates an excluded volume v 0 per particle and a binding energy per particle e b = a/v due to the attractive interactions.Following his reasoning in the simplest way for the Lennard-Jones potential one would assume that particles do not approach closer than r = σ and therefore exclude on the average the volume v 0 = 1 2 × 4πσ 3 /3 = 2πσ 3 /3 per particle as evaluated in the low density limit.The attractive part of the potential is in van der Waals (vdW) theory assumed to be sampled without inducing any correlations, i.e. in the mean-field approximation so that the binding energy parameter becomes: for the Lennard-Jones interaction.The corresponding configurational free energy per particle can be taken to be , where v free is the free volume available to each particle after the excluded volume effect has been accounted for, i.e.
The origin of density functional theory can be traced back to the famous failure of the vdW equation of state in its prediction of unphysical "wiggles" in P(v, T) where the pressure could increase with increasing volume and even become negative.The Maxwell reconstruction (Maxwell , 1875) eliminates these wiggles and is based on interpreting the vdW equation of state as the description of a uniform fluid with density n(r)=n b = 1/v everywhere in the available volume.The total configurational free energy can then be obtained as Now the Maxwell reconstruction can be understood by recognising that the uniform fluid is not the stable state.Instead the fluid separates into a rarefied gas phase (of density n volume V (g) ) coexisting with a condensed liquid phase of density n b and volume V (l) < V (g) .Since the total volume and particle number are conserved we have the two relations: (l) = N, which constrain the search for the most stable form of the fluid.At fixed T, minimizing the configurational free energy with respect to the choice of n

The generalized van der Waals (GvdW) theory
The analysis of the traditional van der Waals theory above shows that: i) equations of state can be understood as approximate local descriptions of energy and entropy in a fluid of given density and temperature; ii) spatially variable densities allow the state of a fluid to be conveniently described.
In the treatment of bulk multiphase fluids as above one can assume that the variation of the density is slow on the length scale of either the repulsive or attractive range of interaction.This is certainly true well inside each phase but not neccessarily true in the interface between phases.Such interfaces are generally of minor extent (N 2/3 ) compared to the bulk phases (N) so that their effect can be neglected in the thermodynamic limit.We may also have slow density variations due to external fields.An example is our atmosphere which is held in place by the gravitational field.We could describe the atmosphere very well with a vdW functional as above parametrized by excluded volume and binding energy.We now consider two improvements of the vdW functional.First we would like to improve the accuracy of the vdW equation of state for simple fluids over their full density range.In liquids, it turns out that the excluded volume per particle v 0 is severely overestimated by its low density limit value 2πσ 3 /3 as used above.Consideration of v 0 as an estimate of the volume per particle in the close packing limit shows that this volume would be just σ 3 in the simple cubic packing arrangement where each hard sphere of diameter σ is assigned a cubic volume in which it fits touching all sides.Thus, setting v 0 = σ 3 will improve the ability of the vdW equation of state to describe the full range of fluid densities.Second, we shall restrict the assumption of local interactions to the short range repulsions (i.e., the excluded volume effect) but give the attractive interactions a nonlocal representation.The GvdW functional for our Lennard-Jones fluid becomes (Nordholm & Haymet , 1980): (5) Here the attractive part of the pair-potential is defined as and we have added an energy due to an external field V ext (r).This functional can treat bulk phase equilibria as well as the main mechanisms of adsorption at a wall and gas-liquid surface tension in a simple fluid of this type (Johnson & Nordholm , 1981;Nordholm & Gibson , 1981). 131 Classical Density Functional Theory of Ionic Liquids

The fine-grained GvdW theory
The key to many phenomena in nonuniform fluids is the nonlocality of interactions, that is the nonvanishing range of interaction on the length scale of density variations.From complete neglect of this nonlocality in the vdW analysis of phase separation we have proceeded above to account for the nonlocality of the pairwise attractive interaction leaving only the shorter range repulsive interactions to local approximation.There are, however, many applications where the nonlocality of the repulsive interactions is also important.We know, for example, that at higher densities hard spheres pack into distinct layers at hard walls.Such packing structures are not predicted by theories based on local representation of excluded volume effects.Moreover, it is clear that repulsive nonlocality will enter also into an accurate prediction of gas-liquid surface tension.For this reason fine-grained GvdW theory has been developed to allow both attractive and repulsive interactions to be treated in a nonlocal representation.The key step is to realize that the excluded volume effect is related to an average particle density over a spatial domain determined by the range of the repulsive interactions.If we, as above, treat the particles as interacting hard spheres of diameter σ then the simplest coarse-graining average reflecting the domain of the repulsive interactions is Thus we assume that the excluded volume effect is related to n(r) while the density itself is n(r) which may be considerably more finely structured.The free energy density functional is then (8) This functional can now describe hard sphere structuring at hard walls and related phenomena (Freasier et al. , 1989;Nordholm et al. , 1980).Other equations of state, and coarse-graining approaches, have been suggested.For instance, Tarazona and Evans (Tarazona & Evans , 1984) based their free energy functional on the Carnahan-Starling (Carnahan & Starling , 1969) equation of state for hard spheres which yields a configurational free energy per particle in terms of the volume fraction η = πσ 3 n b /6.The "weighting" of the particle density was done by Tarazona (Tarazona , 1984) according to where the weighting functions w i are determined to fit the direct correlation function and the radial distribution function for a uniform hard sphere fluid as described by the Percus-Yevick integral equation theory.

Density functional theory (DFT) of simple ionic systems
The simple fluids considered above are all characterized by short range interactions where the second virial coefficient B 2 (T) is finite.In the case of ionic systems the Coulomb energy of two point charges q i and q j separated by a distance r reads where ε 0 is the vacuum permeability and ε is the relative dielectric constant.Hence, the range of interaction is infinite and B 2 (T) diverges.This means that interactions have a more profound effect on these systems and the correlations between the ions play a more fundamental role than in the simple fluids discussed above.

The one-component plasma (OCP)
The simplest IL is the one-component plasma (OCP) where one of two components is treated as an inactive neutralizing background charge while the other is an active ionic component.An example could be a molten metal where the metal ions are taken to move in a background electron gas which is in its ground state.For active ions of magnitude q and density n b the bulk charge density is qn b while the neutralizing background charge density is −qn b .I n the absence of correlations, it follows that the system would behave as an ideal gas with zero potential energy on average.The active ions do repel each other, however, so that compact configurations are suppressed, the average potential energy is negative and the OCP is stable.The long range nature of the interaction means each ion is surrounded by a like-charge-depletion zone (a hole) with a shape that exactly neutralizes its charge q.Thus correlations play a fundamental role in the OCP.The GvdW free energy density functional for the ion density around a central ion can be written as Here ∆n(r)=n(r) − n b is the deviation of the ion density from its bulk value.On introducing the Lagrange multiplier α to accommodate the canonical constraint of fixed particle number, the variational derivative of this functional yields 14) where we recognise the mean-field electrostatic potential ψ(r).Setting this derivative to zero and writing β = 1/(k B T) yields the equation This is the well known Poisson-Boltzmann equation: a mean-field approximation which can be solved by iteration.Linearizing this equation we obtain the Debye-Hückel equation 133 Classical Density Functional Theory of Ionic Liquids www.intechopen.comNow, after defining the dimensionless coupling strength Γ by where a = 4πn b /3 −1/3 is the ion-sphere radius, we get the internal potential energy per particle The electrostatic contribution to the free energy is then The Debye-Hückel theory is both a mean-field and a linear response theory.The radial distribution function is Since physically g(r) is greater or equal to zero the linear response cannot be applicable for small r where g(r) becomes negative.The Debye-Hückel hole (DHH) theory (Nordholm , 1984a) rephrases the Debye-Hückel theory to assign g(r) the value zero at these small r and then adds the exponential part continuously to give Thus the κ-value is retained but the hole is reshaped to be physically realizable and consistent with the asymptotic DH screening of the central charge.At the density n b the "hole" should contain one ion.Thus we have so that In the DHH theory the charge screening distribution can be integrated over the Coulomb energy (12) to yield As before, the electrostatic free energy can be obtained by a charging integration where The DHH theory (Nordholm , 1984a;Penfold et al. , 1991) produces thermodynamic properties in good agreement with simulation results for the full fluid range of coupling constants 1 < Γ < 160.We note that when the ions are imbedded in hard spheres, as in the RPM of electrolytes, the DHH theory above should be applicable when the hard sphere diameter σ is smaller than the electrostatic hole size h above.The theory has been applied in the study of polyelectrolyte screening (Penfold et al. , 1990) and we shall find use for the DHH model of coion correlations in our treatment of ILs below.While the DHH provides an analytic solution to the a charged fluid problem, the more general DFT's discussed later will require numerical iterative solutions.With this in mind, it is useful to consider the convergence properties of the DHH theory, where an analytical solution is available, in order to elucidate some of the issues that may confront the numerical approaches.

An iterative scheme
By considering a simple statistical mechanical theory of bulk point plasmas, we aim to elucidate some general convergence characteristics of successive approximation schemes that are often deployed in free energy density functional calculations for Coulomb fluids.
Although only a very crude level of theory is discussed here, we suggest that the conclusions are relevant for more sophisticated analyses of ILs.
Although the explicit DHH solution for the OCP is very valuable, the corresponding generalisation for arbitrary plasmas presents an ill-posed inverse problem (Penfold et al. , 2005) since the number of parameters grows faster with mixture components than does the number of constraint equations.Nevertheless, a simple iterative scheme can be formulated that resembles the successive approximation algorithms often used for more sophisticated free energy density functional calculations of Coulomb fluids.To study the convergence properties of this scheme in detail, it is useful to re-cast the procedure in terms of the DHH theory for the OCP.Recall the essential idea here is to "renormalise" the Debye-Hückel radial distribution function where g DH is given by ( 21) and ξ = ξ(h) is an asymptotic charge scaling factor fixed to ensure that g remains continuous and non-negative everywhere.By virtue of the previous analysis the appropriate DHH value is ξ (∞) = h exp(κh)/(aΓ), but the following successive approximation scheme immediately suggests itself.Set the iteration counter p = 0 and specify an initial guess ξ (0) .Now, invoke the continuity constraint to fix the current value for electrostatic exclusion hole size as the positive root of where l = κh.Next, appeal to the local electroneutrality condition ( 23) and obtain a new estimate for the renormalisation factor Increment the counter p → p + 1 and repeat the nested cycle of equations ( 29) and (30) until convergence.This iterative procedure is related to the normal method of solving the nonlinear PB equation for screening charge densities.Indeed, the DHH analysis can be regarded as 135 Classical Density Functional Theory of Ionic Liquids www.intechopen.com where l = l ξ (p) is given implicitly by ( 29) and the result ( 24) is clearly manifest as fixed points of the iteration function G.
For the OCP, charge reversal phenomena can be ignored so that ξ must remain positive.Consequently, by virtue of ( 29) and (31), the renormalisation factor is bounded From Figure 1, it is evident that ξ ∈ 0, ξ max implies G(ξ) ∈ 0, ξ max and we are encouraged to expect that G will be a contraction map on a suitable sub-interval that also satisfies the gradient criterion (ii).Differentiation yields where ξ > 0 demands (3Γ) 3/2 > l 3 and the inequality is confirmed by Figure 1.The condition G ′ (ξ) < 1 is satisfied if and only if  0) .The corresponding solution l * given by ( 24) is also indicated (solid line) as well as an upper bound (dotted line) presented by the non-negativity condition on the renormalisation factor ξ.
Figure 2 illustrates the behaviour of L(l).Straightforward analysis establishes a single minimum value of Hence, for 0 < Γ < Γ cut ≈ 2.64 only the non-negatitivity constraint on ξ is active.It follows from the Banach Fixed Point Theorem (Sutherland , 2006) that the iteration sequence generated by G(ξ) will converge to a unique solution from all initial values in the range For Γ ≥ Γ cut , a positive lower bound on ξ (0) arises from the contraction map requirement (ii) Thus, the requirement (i) will also affect the upper bound on ξ (0) (see Figure 1) so that, The results (36), ( 37) and ( 38) are plotted in Figure 3 to visualise the "basin" of attraction as a function of coupling stength.In accord with our intuition, it is evident that the success of the iteration scheme generated by G becomes increasingly more sensitive to the initial guess ξ (0) with increasing Γ.Since G ′ ξ (∞) = 0 (equation ( 33) and Figure 31), we also observe that the fixed point ξ (∞) is super-attracting so that the rate of convergence is at least second order in the neighbourhood of the solution. 137 Classical Density Functional Theory of Ionic Liquids www.intechopen.com

Polymer density functional theory
The polymer density functional theory (DFT) is designed to treat flexible (Woodward , 1991) as well as semi-flexible (Forsman & Woodward , 2003) chains, and has the added advantage that solvent particles can be included explicitly.For this reason it is ideal for describing the underlying skeletal structure of the IL.In some cases, notably when solvent particles and monomers differ considerably in size (Forsman et al. , 2002;Forsman & Woodward , 2005) interesting phenomena arise.The polymer DFT has also been generalized to infinite, as well as polydisperse chains (Woodward & Forsman , 2008;2009).The earliest expression of modern polymer density functional theory was described by Woodward (Woodward , 1991).That work built on the earlier formalism of Chandler and coworkers (Chandler et al. , 1986).Since then, most presentations of polymer density functional theory have used the formalism of Woodward, which separates the ideal and excess components of the free energy.Different versions of the polymer DFT generally amount to different approximations to the excess component.For a reasonably comprehensive review of poymer density functional methods we refer to the review article by Wu and Zhidong (Wu & Zhidong , 2007).
Here we will describe the polymer DFT for finite chains.We denote the position of the i-th monomer in a chain by r i .The complete polymer configuration of an r-mer can then be written as R =(r 1 , ..., r r ) with the corresponding volume element dR = dr 1 ... dr r .Neighbouring monomers along a chain are connected by a bond of a fixed length σ.The Boltzmann factor containing the bonding potential is thus, e −βV B (R) ∝ ∏ r−1 i=1 δ(|r i+1 − r i |−σ), where V B (R) is the bond potential.The r-point density distribution, N(R) is defined such that N(R)dR is the number of polymer molecules having configurations between R and R + dR.We can then extract the monomer density as n(r)= ∑ r i=1 δ(r − r i )N(R)dR.The free energy functional is decomposed into ideal and excess parts so that The ideal chain term is given by where V ext (R) is the externally applied potential.A similar expression is used for the solvent molecules.Note that for ideal chains in an ideal solvent, the theory is exact.The excess part has to be approximated.Following Woodward (Woodward , 1991), we formulate it in terms of the monomer and solvent densities, so that F ex [n m (r); n s (r); nm (r),; ns (r)], where n(r) is a weighted density (Nordholm et al. , 1980) as defined by ( 7).An explicit expression can be obtained by integrating an appropriate equation of state for polymer chains (Forsman et al. , 2002;Wichert et al. , 1996;Woodward & Yethiraj , 1994).
The grand potential is obtained as where μ p and μ s are the chemical potentials of the polymer and solvent fluids in the bulk.By minimising Ω we obtain the corresponding equilibrium density distributions for the polymer and solvent.

Coarse-grained models and dielectric response
In the last few years there has been a move away from full atomistic modelling, toward more coarse-grained descriptions (Fedorov et al. , 2008;Kornyshev , 2007;Lauw et al. , 2009).Coarse-grained models allow us to make some progress in our theoretical analysis.
The models we will use must be simple enough to allow this analytical treatment, but sufficiently detailed to at least discriminate between broad classes of ILs and reproduce observed interfacial structures and behaviours.Steric interactions alone will usually give rise to interfacial structures on the length scale of the molecular size in a dense liquid.However, electrostatic and dispersion forces are also important.Indeed, monolayer enrichment of alkyl chains in imidazolium based ILs has been observed in angle-resolved x-ray photoelectron studies of the vapour/liquid interface (Lockett et al. , 2008b;Maier et al. , 2010), indicating the important role of attractive forces at the interface.Related to this are the observations of small-angle X-ray scattering experiments that appear to indicate structure on the mesoscale due to the incompatibility of the ionic and non-polar parts of cationic molecules, with sufficiently long carbon chains (Qiu & Texter , 2008).It has been conjectured that these structures are responsible for the unique solvent properties of some ILs, and to their ability to impart stability to dispersions of metallic nanoparticles.When mesoscopic structuring is present, added ions will most likely concentrate in the polar regions of the IL and affect their structure and stability.Other types of molecules added to ILs will also affect structure.Polymers and nanoparticles may stabilize existing mesostructures, or else may act as templates for others.More specific interactions can also occur in ILs, especially those based on dialkyl-imidazolium cations, which are known to form hydrogen-bonded and/or π − π stacked supramolecular structures (Mele et al. , 2003).These structures are possibly responsible for the anomalous temperature dependence in DC measurements of imidazolium ILs (Lockett et al. , 2008a  be obtained by using a partial charge distribution, rather than a single point charge.Partial charge models are widely available in the current literature.Further variability can be introduced via dispersion and steric interactions, which will be modelled by Lennard-Jones potentials centred on each sphere.
• the anion in many ILs is usually simpler than the cation and can often be represented as a negatively charged sphere.More complex anionic species can be modelled using similar approaches to that used for the cations.
• electrostatic interactions will be scaled by a relative dielectric constant, which mimics the low frequency polarizability of the molecular constituents.
Within the framework of this general model, other refinements are possible.For example, bond stiffness in the cation can be introduced via next nearest-neighbour interactions.More specific interactions can also be modelled.For example, some imidazoliums are believed to undergo π − π stacking interactions.One area of concern in modelling ILs, especially using approximate theory, is the appropriate choice of dielectric constant.Typical experimental values for the low frequency dielectric constant of ILs range between 10 and 20.As our model is coarse-grained, dielectric effects due to electronic, orientational and distortional polarization are neglected.The electronic contribution to ǫ r is around 2, while the latter two contributions are due to charge distribution and molecular flexibility, respectively.In our model the molecules possess no permanent dipoles, as the unit charges are placed at a single center.Therefore, compared with "real" ILs, the electrostatic coupling may be overestimated.Recent work by Izgorodina and coworkers (Izgorodina et al. , 2009) shows that the orientational contribution to ǫ r , can range between approximately 1 and 20, depending on the size of the molecular dipole moments, while the distortion contribution remains small (approximately unity).
The appropriate value for ǫ r to be used in computer modeling is the subject of considerable debate in the literature (Izgorodina et al. , 2009;Kobrak & Li , 2010;Krossing et al. , 2006).This is because the simulations in principle capture a molecular low-frequency response, but not the electronic response.This leads to some uncertainty regarding the appropriate value with which to scale down the Coulomb interactions in a simulation.Paradoxically, the more approximate DFT approach is less ambiguous in this respect.The electrostatic correlation functional that we will use in this work is based on a one-component plasma (OCP) model.Long-ranged molecular responses of the kind just discussed, are not captured by the theory.This warrants the use of typical experimental values for the DFT calculations, so the value of ǫ r can be chosen in a relatively unambiguous manner.
correction used here accounts only for short-ranged effects, expected to dominate in ILs, and do not account for the ionic polarization contribution to the dielectric screening of electrostatic interactions.Thus the mean-field terms, that is the first two terms in (44), remain unaffected by the correlational term and the ionic polarization contribution must be included in the value of ǫ r .The total grand potential, Ω DFT tot is obtained as The functionals can be simplified by integrating out the (x, y) dimensions, parallel with the "electrode walls".Minimization of the chosen functional, with respect to the fluid particle densities then follows in the usual manner, wherein electroneutrality is ensured via a Donnan potential.

Performance of the DFT for the full IL model
In our previous publication (Forsman et al. , 2011), we demonstrated that our DFT is able to reproduce simulated differential capacitance curves qualitatively, and with a semi-quantitative accuracy.Furthermore, the level of the curves, and their overall shape, agrees well with experimental results (Islam et al. , 2009;Lockett et al. , 2008a), which tends to support our simple model.Still, since we are armed with a theory and a microscopic model, we can go further and make scrutinizing analyses of the microscopic structure at the surfaces.These are not easily accessible by experimental methods, but are nevertheless of fundamental interest.One of the major strengths of a theoretical approach is the ability to correlate microscopic detail to macroscopic observations.We start by analyzing the ability of our DFT to predict density distributions at the electrodes.Comparisons with simulations are given in Figure 5. Details of the simulation parameters and methods are provided elsewhere (Forsman et al. , 2011).We note that the DFT tends to underestimate the height of the primary peak; presumably reflecting a slight underestimation of ion correlation effects.At a negative surface, the DFT predicts a too long-ranged oscillatory structure, and with a period that also seems a bit overestimated.Still, overall accuracy is quite reasonable.At a positive surface, we do observe some larger discrepancies, as seen in Figure 5 (b).In particular, the simulated data display a distinct primary peak of the charged monomers, outside the first anion "layer".A corresponding peak is detectable also with DFT (see inset), but it is much less pronounced.We believe this primarily results from our neglect of Coulomb truncation at short range, due to the hard cores of the particles.We are working on an improved DFT to remedy this problem.The total monomer density (of charged as well as neutral monomers) is given in Figure 5 (c).
We believe the agreement is satisfactory, although the DFT in this case slightly overestimates the primary adsorption peak.
When making these comparisons, we should keep in mind the problem of ambiguity regarding the choice of ǫ r in the simulations.We aim to describe an IL with ǫ r = 12 (used in DFT), but since the simulations capture a part of this (as discussed earlier), we have set ǫ r = 7.6 for the MC runs.This is probably a reasonable value, but there is still some degree of ambiguity.

A model system with only dispersion
Having demonstrated that our DFT is reasonably accurate, albeit not perfect, it might be a good idea to investigate the performance of the DFT that would result if we remove all .Density profiles at "electrodes", as predicted by DFT, and simulated ("MC"), respectively.The "underlying" soft electrode-particle interaction is non-adsorbing, i.e. a α w = 0 (with α = m,c) (a) Simple anions ("s") and (positively) charged monomers ("m") at a strongly negative surface, σ −1 el = −50 Å 2 /e.The densities are normalized by their value at the mid plane.The inset is a zoom-in.(b) Simple anions and charged monomers, at a strongly positive surface, σ −1 el = 50 Å 2 /e.The inset is a zoom-in.(c) Total monomer density, outside at a strongly negative surface, σ −1 el = −50 Å 2 /e.charges, on ions as well as on the surfaces.We then end up with a binary Lennard-Jones mixture, where one component has a uniform oligomeric structure (a "pearl-necklace" of L-J monomers).We use the same densities as used for the fully charged system (in simulations and DFT, respectively).One difference from an ordinary binary L-J mixture is that we retain a "Donnan potential", ensuring that the number of simple particles (anions in the IL) in the slit region equals the number of oligomers.This approach, we believe, ensures that we make the system as similar as possible to the original IL, but without the charges.This should allow us to establish which of our approximations in the full DFT that are likely to be most problematic, i.e., where we should focus our efforts to improve the DFT.Structural comparisons are given in Figure 6, for adsorbing (a α w = 1) as well as non-adsorbing (a α w = 0) surfaces.All qualitative features (humps and oscillations) are well reproduced by the DFT, although it does tend to overestimate the primary solvent peak somewhat.Given this satisfactory (a) Non-adsorbing surfaces, a α w = 0 (b) Adsorbing surfaces, a α w = 1 performance for neutral species, we can almost with certainty state that the discrepancies found for ILs are primarily due to approximations of the Coulomb interactions; specifically the way in which these are altered by correlations in the system.As already mentioned, these problems are currently being addressed.After having evaluated the performance of the DFT (see also (Forsman et al. , 2011)), with the conclusion that it faithfully reproduces the true behaviour of the model, we will henceforth display DFT predictions without the corresponding time-consuming simulations.

Effects of the range of the soft surface repulsion
Let us now consider the role of the width of the depletion region in the case of non-adsorbing surfaces.In Figure 7, we have compared differential capacitance curves with σ = 4Å and a α w = 0, but with two different choices of σ w .In our previous simulations/calculations, we have set σ w = σ, but here we also consider the case σ w = 1.5σ.The differential capacitance curve responds rather dramatically to a change in the distance of closest approach of the fluid!Within the exclusion region the potential varies strongly and linearly, giving rise to the so called "compact layer" contribution to the capacitance (Fedorov et al. , 2008).The "linear potential effect" is illustrated in Figure 8, where the corresponding potential profiles are shown.We see how these profiles are very similar, but "z-shifted" relative to the plane of charge.For the surface with a more long-ranged soft repulsion, σ w = 1.5σ, the linear regime (where there are no responding ions) extends further out, which leads to a larger potential difference relative to the bulk value.While the "exclusion effect" in principle is trivial, we believe its relevance has been underestimated in the literature.

Effects of temperature and dielectric constant
In a recent experimental work, Lockett and coworkers (Lockett et al. , 2008b) claimed to establish the temperature dependence of the differential capacitance for imidazolium chloride ILs at around 393 K.We are not experimentalists, but have learned that it is quite difficult to ascertain the amplitude of C D .Still, we will henceforth proceed under the assumption that the reported findings are correct.Specifically, they found that the C D curve increases with temperature.Can we capture such a behaviour with our simple model, and if so, under what circumstances?Let us first consider the process of merely increasing the temperature, at a constant pressure.The corresponding DFT predictions, for non-adsorbing electrodes, are collected in Figure 9.We see that the differential capacitance in this case drops with temperature, except at low surface potentials.This is in contrast with the experimental findings (Lockett et al. , 2008b).But what if the dielectric constant of the IL also changes with temperature?Specifically, let us assume that the low-frequency dielectric constant increases with temperature.In Figure 10, we highlight the consequences to the response of C D to an increased temperature.Under these circumstances, the differential capacitance is indeed predicted to increase with temperature, in accordance with the observations (Lockett et al. , 2008b)!Establishing the temperature dependence of ǫ r for ILs is not an easy experimental task, but Hunger and coworkers (Hunger et al. , 2009) have recently obtained some estimates.It turns out that the static dielectric constant in most cases seems to descrease slightly with temperature.Still, there are exceptions, and at least one case with a positive value of ∂ǫ r /∂T was detected.We refrain from a discussion about the reliability of the various experimental measurements.Instead, we simply state that our simple theoretical model will predict an overall increase of C D with temperature, provided that ǫ r also increases.

Electrodes with adsorption or specific adsorption
Let us furthermore consider the case of adsorbing, or partially adsorbing electrodes.We recall the wall potential defined by ( 42) that allows us to treat the case of specific adsorption of the charged species, which presumably are more polarizable than CH 2 groups.Furthermore, there may be important contributions from image charge interactions, as discussed below.As we have reported in previous publications (Forsman et al. , 2011;Trulsson et al. , 2010), adsorption has a dramatic influence on the differential capacitance curve.Here we also consider the case of specific adsorption, with results presented in Figure 11.As might be anticipated, ion specific adsorption is quite efficient when it comes to "reducing the camel hump".Recall that the "hump" is the result of ion depletion at low surface potentials.This depletion is of course less pronounced when there is a non-electrostatic attraction to the electrodes.Note that these theoretical findings are commensurate with experimental In this case, neutral monomers adsorb weakly (a m w = 0.5), while anions (simple) and cations (charged monomers) display a strong adsorption (a c w = 1).differential capacitance data on metallic and non-metallic electrodes (Islam et al. , 2009).Specifically, one will in general find a bell-shaped C D curve with metal electrodes (no camel-hump), and U-shaped C D with glassy carbon electrodes.The latter scenario implies a camel-hump shaped complete C D curve, although experimental limitations (redox-reactions etc.) may not reach the saturation limit.These findings can be rationalized if we consider that metallic electrodes are likely to be more strongly adsorbing at low surface potentials, on account of the expected attractive image charge forces, and strong dispersion interactions.As we have seen here, adsorption will tend to diminish "camel-hump" tendencies, ultimately resulting in a purely bell-shaped differential capacitance curve.Such adsorbing interactions are likely to be weaker with glassy carbon electrodes, which naturally will promote the occurrence of a camel-shaped C D curve.Finally, we investigate the temperature and dielectric response dependence for an "ion-specific adsorption" electrode, with a m w = 0.5; a c w = 1.The DFT predictions are provided in Figure 12.Here we quite clearly see how the differential capacitance will drop upon a

Conclusions
We have tried to illustrate how classical DFT applied to simple coarse-grained models, can be used to gain an improved understanding of fundamental physical mechanisms underlying macroscopic behaviours of ILs.Furthermore, the DFT calculations are noise-free and exceptionally rapid, permitting more elaborate models to be studied.Examples of simple extensions/improvements are an intrinsic stiffness of the oligomer chain (Forsman & Woodward , 2003), distributed charges and/or a structurally more detailed model.Such extensions can easily be implemented.

129
Classical Density Functional Theory of Ionic Liquids www.intechopen.comenergy scale ǫ.Such models are referred to as simple fluids.The most popular form of simple fluid pair-potential is the Lennard-Jones potential φ

b
will then produce the familiar Maxwell reconstructed isotherms with flat segments (tielines) replacing the wiggles.

Fig. 1 .
Fig. 1.A representative family of iteration function G for generating the sequence of asymptotic charge renormalisation factors.Coupling constant Γ values corresponding to each curve are indicated and span the range consistent with an OCP in the fluid state (0 < Γ < 140).Lines of slope ±1 are indicated by the symbols.a very simplified PB theory.It will be convenient to express the iteration sequence by the recurrence relation

Fig. 2 .
Fig. 2. Bounds on the Coulomb hole size l required for the iteration function G to be a contraction map are shown.The gradient condition (ii) leads to the result (34) where L(l) is plotted here (broken line).For sufficiently large coupling parameters Γ ≥ Γ cut , the largest root l min of L = Γ sets a non-trival lower bound ξ lo on the initial sequence term ξ(0) .The corresponding solution l * given by (24) is also indicated (solid line) as well as an upper bound (dotted line) presented by the non-negativity condition on the renormalisation factor ξ.

Fig. 3 .
Fig. 3.The range of initial term values ξ (0) is shown as a function of couplping strength Γ that will ensure convergence of the iteration scheme defined by G to a unique fixed point.

Fig. 6 .
Fig.6.Calculated and simulated density profiles, for our model fluids, but with all charges being removed.(a) Non-adsorbing surfaces, a α w = 0 (b) Adsorbing surfaces, a α w = 1 performance for neutral species, we can almost with certainty state that the discrepancies found for ILs are primarily due to approximations of the Coulomb interactions; specifically the way in which these are altered by correlations in the system.As already mentioned, these problems are currently being addressed.After having evaluated the performance of the DFT (see also(Forsman et al. , 2011)), with the conclusion that it faithfully reproduces the true behaviour of the model, we will henceforth display DFT predictions without the corresponding time-consuming simulations.

Fig. 8 .
Fig. 8.The variation of the mean electrostatic potential, Ψ, with distance, z, from a negatively charged surface.The surfaces are non-adsorbing, with two different choices of σ w ,asin Figure 7.

Fig. 9 .Fig. 10 .
Fig.9.Differential capacitance curves, with non-adsorbing electrodes, at two different temperatures.The dielectric constant is set to be temperature independent.

Fig. 11 .
Fig. 11.Effects of non-electrostatic electrode adsorption.The top differential capacitance curve illustrates effects from ion specific adsorption.In this case, neutral monomers adsorb weakly (a m w = 0.5), while anions (simple) and cations (charged monomers) display a strong adsorption (a c w = 1).differential capacitance data on metallic and non-metallic electrodes(Islam et al. , 2009).Specifically, one will in general find a bell-shaped C D curve with metal electrodes (no camel-hump), and U-shaped C D with glassy carbon electrodes.The latter scenario implies a camel-hump shaped complete C D curve, although experimental limitations (redox-reactions etc.) may not reach the saturation limit.These findings can be rationalized if we consider that metallic electrodes are likely to be more strongly adsorbing at low surface potentials, on account of the expected attractive image charge forces, and strong dispersion interactions.As we have seen here, adsorption will tend to diminish "camel-hump" tendencies, ultimately resulting in a purely bell-shaped differential capacitance curve.Such adsorbing interactions are likely to be weaker with glassy carbon electrodes, which naturally will promote the occurrence of a camel-shaped C D curve.Finally, we investigate the temperature and dielectric response dependence for an "ion-specific adsorption" electrode, with a m w = 0.5; a c w = 1.The DFT predictions are provided in Figure12.Here we quite clearly see how the differential capacitance will drop upon a

Fig. 12 .
Fig. 12. Temperature response, for the "ion-specific adsorption" system in Figure11.Cases of a constant as well as temperature-increasing low frequency dielectric response are given.

147
Classical Density Functional Theory of Ionic Liquids www.intechopen.comtemperature increase, if ǫ r remains unchanged.An increasing dielectric constant may on the other hand facilitate a C D curve that increases with temperature.
).The model template for ILs that we will use in our work embodies the essential features of ILs and possesses a number of adjustable parameters to allow us the flexibility of representing different classes of ILs.It is schematically represented in Fig 4 for the case of imidazolium ILs and has the following general features: • the cation is represented as a charged sphere with attached non-polar oligomeric chains (neutral spheres).This is adequate to model a wide variety of ILs including imidazoliums, pyridiums and ammoniums.Further discrimination between these cationic types can 139Classical Density Functional Theory of Ionic Liquids www.intechopen.com