Electrical Double Layer Properties in Diameter Asymmetric Molten Salt

Results of the Monte Carlo simulation of the electrode | molten salt interface are reported. The system investigated is represented by a primitive mod...
0 downloads 0 Views 789KB Size
Download PDF
J. Phys. Chem. C 2010, 114, 13329–13333

13329

Electrical Double Layer Properties in Diameter Asymmetric Molten Salt Investigated by Grand Canonical Monte Carlo Method Jacek Kłos and Stanisław Lamperski* Department of Physical Chemistry, Faculty of Chemistry, A. Mickiewicz UniVersity, Grunwaldzka 6, 60-780 Poznan´, Poland ReceiVed: May 14, 2010; ReVised Manuscript ReceiVed: June 29, 2010

Results of the Monte Carlo simulation of the electrode | molten salt interface are reported. The system investigated is represented by a primitive model of electrolyte in contact with a planar charged hard wall. The density of an electrolyte is comparable with the density of real molten salts. The calculations are performed for T ) 1500 K, which is above the phase transition temperature determined from Monte Carlo analysis of the heat capacity CV and the radial distribution function. The density and charge profiles reveal a multilayer structure of the electrolyte in the vicinity of the electrode surface. The differential capacitance has a distorted camel-shape similar to that predicted by theory and observed for some ionic liquids. Introduction Molten salts are dense systems that exist at high temperatures. They can replace traditional aqueous electrolytes in some applications, for example, batteries and fuel cells. In molten form, the range of temperatures in which they can exist is substantial, while they also have good electrochemical and thermal stability. The electrochemical properties of molten salts have long been of keen research interest both experimentally and through numerical simulations. Recently, Lamperski and Kłos1 have presented grand-canonical Monte Carlo (GCMC) simulation results for a molten salt, treated as a restricted primitive model electrolyte near a planar electrode surface. The authors found that the structure of an electrolyte in the vicinity of the electrode surface depended on the electrode charge but not so much on temperature. The capacitance curves had a parabola-like shape with a maximum at the electrode charge σ ) 0. This result is not consistent with the Gouy-Chapman theory,2,3 but it is in agreement with the predictions of the modified Poisson-Boltzmann theory (MPB),1 which includes corrections for the correlation and exclusion volume effects.4 Lamperski et al.5 have determined the values of the reduced temperature and density at which the differential capacitance, Cdiff, curves change their shape from having a minimum at zero electrode charge to having a maximum there. The results obtained from the MPB theory and the Poisson-Boltzmann equation with an exclusion volume term (PB + EVT) were compared with those obtained by Monte Carlo (MC) simulation. The exclusion volume and interionic correlations were found responsible for the resultant shape of the differential capacitance curve. Stafiej et al.6 were perhaps the first who showed that, for a dense and asymmetrical system, the Cdiff curve has a distorted bell-like shape. Fawcett et al.7 have also observed the change from a minimum to a maximum in the Cdiff curve with increasing concentration of the solution. Our previous work1 concerned symmetrical electrolytes. It is of interest to know how the asymmetry of ion diameters influences the electric double layer properties. Investigation of this problem is the subject of the present paper. * To whom correspondence should be addressed. Phone: +48 61 8291454. Fax: +48 61 829 1505. E-mail: [email protected].

Some insight into the properties of molten salt double layer can be drawn from a comparison with the theoretical and experimental results obtained for ionic liquids. Besides the main similarities, which are the high density and solventless liquid, there are differences. Ionic liquids are usually composed of organic ions and much smaller inorganic counterions. They occur as fluids at evidently lower temperatures than molten salts and have limited thermal stability of the organic ion. Some theoretical explanation of the properties of electrical double layer has been provided by Kornyshev.8 On the basis of the theory of the Poisson-Boltzmann (PB) lattice-gas analysis, he found that the Cdiff as a function of electrode potential was a bell shaped curve with a maximum at the zero potential when all ions were of the same size, while for ions of different size, the Cdiff curve had a distorted camel-shape. Next, Fedorov and Kornyshev9,10 have carried out molecular dynamics simulations of an ionic system with the nonelectrostatic interactions described by the Lennard-Jones potential. They found that, for ions of the same diameters, the Cdiff curves have a bell-like shape, which is in agreement with the mean-field theory.9,10 For ions of different diameters, the shape changes, and the Cdiff curve is also shifted away from zero charge. Recently, Fedorov et al.11 proposed a more realistic model of ionic liquids by adding neutral parts to cations (linear two- or three-bead rigid conformation). Monte Carlo simulations for this nonsymmetrical system gave highly distorted Cdiff curves that resemble some experimental results. Lockett et al.12 have determined the differential capacitance for several ionic liquids with a chlorine anion on the glassy carbon electrode at different temperatures. Their results show the Cdiff versus electrode potential curves to be generally distorted camel-shape with two maxima, higher on the side of positive potentials. The minimum between the maxima is in the vicinity of zero potential. The same authors reported increase in the Cdiff with increasing temperature. Islam et al.13 carried out the capacitance measurements for several ionic liquids on the gold, platinum, glassy carbon, and mercury electrodes. They observed different shapes of Cdiff curves: belllike shape, camel-shape, and U-shape. Generally, the capacitance curves measured on Pt and Au have a convex-shape while those measured on Hg and glassy carbon have a U-like shape. These

10.1021/jp104402u  2010 American Chemical Society Published on Web 07/16/2010

13330

J. Phys. Chem. C, Vol. 114, No. 31, 2010

Kłos and Lamperski

results show that the capacitance depends not only on the physical interactions but also on the chemical ones. Model and Simulation Details Monte Carlo calculations were carried out in the grand canonical ensemble. The method involves three moves in the simulation box: insertion and removal of an electrically neutral group of ions and displacement of an ion. The expressions for the probability of acceptance of these moves are given elsewhere14 and in our previous papers.1,15 The simulation box was a rectangular prism of the dimensions W × W × L. Two planar hard walls, impenetrable to ions, represented the electrode surfaces. Both of them were charged with uniform electric surface charge density, σ. The interactions of ions with the images in the electrodes were not considered. The simulation box was replicated in the directions parallel to the electrode surface. The ions were treated as hard spheres of diameter d+ (cation) and d- (anion) with the point electric charge z+e or z-e at the center (e is the elementary charge, z is the charge number). They were immersed in the continuous medium characterized by relative permittivity εr. This is the primitive model (PM) of an electrolyte. When the density is high and εr describes the electrical polarizability of ions, this model can be applied to molten salts or ionic liquids. The hard-sphere and electrostatic interactions between the ions present in the simulation box were calculated explicitly with the use of the minimum image convention. Similarly, the interactions of these ions with the electrode surface were calculated exactly. The electrostatic interactions of an ion with the ions from outside the box were estimated by the method introduced by Torrie and Valleau,16 in which infinitely large charged planes with a square hole for the simulation box were applied. The distribution of ion concentration against the electrode surface was described by the singlet distribution function, g. The volume charge density Fv and the mean electrostatic potential Ψ were calculated from the singlet distribution function

FV(x) ) eNA

∑

ziniogi(x)

(1)

i)-,+

ψ(x) ) -

∑

e z no εrε0 i)-,+ i i

L/2 g (x )(x1 - x)dx1 ∫max(x,d /2) i 1 i

(2) where noi is the bulk number density of ions of kind i, NA is the Avogadro constant, εr is the relative permittivity, and ε0 is the vacuum permittivity. Equation 2 for Ψ can be found in ref 17. The interpolation polynomials method5,18 was applied to calculate the differential capacitance. The second-order polynomials were fitted to five, six, or seven successive (Ψ(0), σ) points, and thus, we obtained three sets of Cdiff versus σ dependences. All of them appeared similar in shape (the MC data proved to be numerically stable) and were smoothed with the LOWESS transformation procedure.19 Numerical Details of the Simulation The simulation box parameters were as follows: W ) 5063.085 pm and L ) 5800.000 pm. The difference in the number of ions N+ and N- present in the box had to be constant during the simulation to give the desired electrode surface charge σ:

σ ) -e(z-N- + z+N+)/2W2

(3)

A 1:1 molten salt with the ion diameter d- ) 400 pm (average of the crystallographic radii of Cl-, Br- and I- anions) and d+ ) 320 pm (average of the crystallographic radii of Rb+ and Cs+)20 was investigated. The density of a system is usually given by the packing fraction η: o 3 o 3 η ) π(n+ d+ + nd-)/6

(4)

When using the value of the packing fraction, η ) 0.35, suggested for molten salts by Larsen,21 one obtains the bulk density of ions no- ) no+ ) 6.91 × 1027 m-3. The simulations were carried out for the following electrode charges (in σ/Cm-2) -0.5, -0.4, -0.3, -0.25, -0.20, -0.15, -0.1, -0.05, -0.025, 0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5. The individual activity coefficients γ- and γ+, needed for the GCMC simulations, were obtained by the recently developed inverse GCMC method (IGCMC).15 This method allows evaluation of individual activity coefficients for a specified electrolyte concentration. The number of configurations used to equilibrate the system was 40 × 106. The averages were calculated from the next 100 × 106 configurations. Simulation for each σ was repeated 10 times to improve the statistics of the results. Temperature is one of the physical factors affecting the degree of molecular order. Hence, before starting the simulations of the electrical double layer properties, we wanted to determine the melting temperature of our salt to be able to carry out simulations above this temperature. The melting point can be estimated from the temperature dependence of heat capacity. The heat capacity at constant volume, CV, was calculated from energy fluctuations during the constant volume MC simulation for the bulk electrolyte. The dependence of the reduced heat capacity C*V (C*V ) CV/NkB, where N is the number of ions and kB is the Boltzmann constant) on the reduced temperature T* [T*) 4πε0εrdkT/e2 with d ) (d+ + d-)/2 ] is shown in Figure 1. There are two high peaks at T* ) 0.21 and 0.22, respectively, and two low peaks, decaying up to T* ) 0.25. Their presence suggests structural transformations in the temperature range 0.21-0.22. Earlier, we have determined that, for the equal ion size system, the heat capacity maximum was at about T* ) 0.251 while Boda et al.22 reported that, below the critical temperature T*c ) 0.282, the capacitance increased with increasing temperature, which was an unusual behavior, while above T*, c the capacitance decreased. To confirm our results, we also analyzed the shape of the cation-cation radial distribution at different temperatures. As shown in Figure 2, the curves at T* ) 0.17, 0.19 are evidently different from those at T* ) 0.26, 0.28. This suggests again the structural transformations from a solidlike structure to a liquidlike one. The two curves at T* ) 0.17 and 0.19 with well formed maxima indicate a solidlike phase. Projections of ions present in the simulation box on xy, yz, or zx planes at the end of simulations show the bcc crystallographic structure. Two curves at T* ) 0.26 and 0.28, above the temperature of transformation, are characteristics of liquidlike system. Because the investigation described was carried out under constant volume, not constant pressure, it was hard to say if we observed a typical phase transition. So, all we can say is that the crystal-like structure occurs in temperatures below T* ) 0.21, while above T* ) 0.25, we have a liquid phase converting probably to a supercritical fluid with further increase in temperature. For all these reasons, we decided to carry out the double layer investigation at a temperature higher

Electrical Properties in Molten Salt

Figure 1. Temperature dependence of reduced heat capacity, C*V, for a 1:1 electrolyte with ion diameters d- ) 400 pm and d+ ) 320 pm.

Figure 2. Radial distribution functions g(x) of cations for a 1:1 electrolyte with ions of diameters d- ) 400 pm and d+ ) 320 pm at different temperatures T* ) 0.15 (squares), 0.19 (triangles), 0.26 (circles), and 0.28 (diamonds).

J. Phys. Chem. C, Vol. 114, No. 31, 2010 13331

Figure 3. Singlet distribution functions of anions, g-(x), for a 1:1 (d) 400 pm and d+ ) 320 pm) electrolyte at σ/Cm-2 ) -0.5 (circles), -0.2 (squares), -0.05 (triangles up), 0.05 (hexagons), 0.2 (diamonds), and 0.5 (triangles down); the other parameters are as in the text.

Figure 4. Singlet distribution functions of cations, g+(x), for a 1:1 (d- ) 400 pm and d+ ) 320 pm) electrolyte. The symbols have the same meaning as in Figure 3.

than that of the phase transition, namely, at T*) 0.32. This temperature can be obtained for different values of T and εr (the ion diameters were given earlier). These values should have physical meaning. When assuming that εr ) 29-11 and d ) 380 pm, one obtains T ≈ 7400 K, so we have chosen εr ) 10 as we did earlier.1 This value is the average permittivity of common salt crystals.23 Using it, we approximately obtain T ) 1500 K, and at this temperature, we decided to carry out our investigations. The individual activity coefficients calculated using the IGCMC method for 1500 K were γ- ) 277.051 and γ+ ) 34.050. Results and Discussion The singlet distribution functions of anions and cations are shown in Figures 3 and 4, respectively. They reveal damped oscillatory behavior with increasing distance from the electrode. At the contact distance, the value of the anion distribution increases from 0.24 at σ ) -0.5 Cm-2 up to 17.3 at σ ) 0.5 Cm-2. At the same time, the second maximum diminishes in intensity and moves slightly away from the electrode. The cation distribution functions show similar behavior, but the maximum at the contact distance decreases from 17.5 at σ ) -0.5 Cm-2 to 0.02 at σ ) 0.5 Cm-2. Figure 5 presents the charge density profile at positive electrode charges. Because of the smaller diameter of cations, in the region x/d- ) 0.4 to x/d- ) 0.5, there appears a thin

Figure 5. Volume charge density, Fv(x)/Cm-3 of the electrolyte at σ/Cm-2 ) 0.025 (circles), 0.1 (squares), 0.25 (triangles), and 0.5 (diamonds). The other parameters are as in the text.

positively charged layer (TPCL), despite the fact that the electrode is positively charged. The existence of TPLC at positive electrode charges is due to a geometric packing effect (the drift of smaller size cations into region 0 < x < d-/2, where anions cannot go). The layer of the cations closest to the electrode causes repulsion of other cations located farther from the electrode surface. The TPLC has a small negative slope, which diminishes with increasing electrode charge. The total number of cations in the TPLC diminishes with increasing

13332

J. Phys. Chem. C, Vol. 114, No. 31, 2010

Figure 6. Volume charge density, Fv(x)/Cm-3 of the electrolyte at σ/Cm-2 ) -0.025 (circles), -0.1 (squares), -0.25 (triangles), and -0.5 (diamonds). The other parameters are as in the text.

Figure 7. Profile of mean electrostatic potential for a 1:1 (d- ) 400 pm and d+ ) 320 pm) electrolyte at σ/Cm-2 ) -0.5 (circles), -0.2 (squares), -0.05 (triangles up), 0.05 (hexagons), 0.2 (diamonds), and 0.5 (triangles down). The other parameters are as in the text.

electrode charge due to the electrostatic repulsion. With increasing electrode charge, the TPLC disappears, and its influence on the ion adsorption vanishes. In the range x/d- ) 0.5 up to around 1.0, there is a layer with the predominance of the negative charges while the next layer is again charged positively. Figure 6 shows a reverse situation to that in Figure 5 at the negative electrode charges. At small electrode charges (from σ ) -0.025 Cm-2 to σ ) -0.1 Cm-2), the TPCL evidently screens the electrode charge, so at the layer adjacent to TPLC, the unexpected anion adsorption is observed. With increasing negative electrode charges, the positive charge of TPCL increases, which intensifies the screening effect. Figure 7 shows the dependence of the mean electrostatic potential, Ψ, on the distance from the electrode for selected positive and negative electrode charges. The curves again exhibit damped oscillations with increasing distance from the electrode. In the vicinity of the electrode, the potential increases or decreases rapidly, but with increasing distance, such changes are more moderate. In general, the potential curves have three distinct extremes at x/d- ≈ 0.5, 1.4, and 2.2 (their magnitudes depend on the electrode charge). It is worth noting here that Valisko´ et al.24 carried out a series of MC and density functional theory calculations for a PM electrolyte at a charged interface. They presented potential profiles, which in contrast to our results, do not have the oscillations perhaps due to the lower density of electrolyte and weaker electrostatic interactions. Our potential curves are not symmetrical for the same absolute values

Kłos and Lamperski

Figure 8. Differential capacitance of the electrical double layer, Cdiff, as a function of the surface charge density σ for a 1:1 (d- ) 400 pm and d+ ) 320 pm) electrolyte. Data obtained from the interpolation polynomials method for five points (circles), six points (triangles), and seven points (squares). A solid line and diamonds are the results smoothed by LOWESS transformation procedure. The other parameters are as in the text.

of electrode charges, which is clearly visible in the range of distances x/d- ) 0 to x/d- ) 1.2. In the case of negative electrode charges, the extremum is shifted toward the electrode as compared to the case of positive electrode charges that is due to the smaller size of cations than anions. Also, due to the ion size asymmetry, the smaller cations can come closer to the electrode surface, and this is responsible for greater absolute values of the electrode potential for positive electrode charges as compared to the values for the negative ones. Thus, the potential of zero charge (pzc) obtained from the linear interpolation of two smallest electrode charge potentials amounts to Ψ ) 0.0332 V. This result clearly shows that the nonzero potential at σ ) 0 can also occur due to asymmetry without recourse to specific adsorption. Figure 8 shows a plot of Cdiff versus σ. Generally, it is a nonsymmetrical curve with two maxima of different heights (distorted camel-shape). The higher maximum is located at σ ≈ -0.15 Cm-2 while the lower is at σ ≈ 0.15 Cm-2. Also, the capacitance wing at negative electrode charges has higher values than that at the positive charges. The greater values of Cdiff on the side of negative electrode charges are due to closer contact distance of cations with electrode than with anions. Similar distorted bell-like Cdiff curve for an asymmetrical, high density system was obtained by Stafiej et al.6 Here, it is insightful to compare our results with those obtained for some ionic liquids. The differential capacitance curves of molten salts are expected to show a qualitatively similar shape to that of ionic liquids. Our curve has a shape similar to that obtained by Kornyshev8 for the lattice saturation factors γ+ ) 0.2 and γ- ) 1.0 (do not confuse the Kornyshev’s γ with the activity coefficient); however, the sequence of heights of the maxima is reversed. Also, the curve obtained by Kornyshev is slightly shifted toward the negative potentials as compared to our simulation results. It is interesting to compare our capacitance results with those obtained by Fedorov and Kornyshev10 who considered also different ion diameters. To make the comparison possible, we introduce the reduced capacitance C* (C* ) Cd/4πε0εr) and the reduced potential Ψ* (Ψ* ) Ψ4πε0εrd/e). To be consistent with the Fedorov and Kornyshev model, we swapped over our diameters of anions and cations, so d- ) 320 pm and d+ ) 400 pm. Thus, we obtained a reduced capacitance curve that is the analogue of the reflection of the capacitance in Figure 8

Electrical Properties in Molten Salt

J. Phys. Chem. C, Vol. 114, No. 31, 2010 13333 (primitive model of molten salt). The main achievement of this work is the determination of the Cdiff versus σ dependence, which has a distorted camel-shape. The capacitance curve was compared with the available theoretical and experimental data. The theoretical results show qualitative similarity with those presented in this paper. Also, some similarities were found with the experimental data. However, a full agreement cannot be expected because in the real systems the material of the electrode plays an important role with the chemical bonding (whose nature can only be described adequately by quantum theory) of the surface with the ions. The full description of the properties of the electrical double layer is a complex problem on both physical and chemical levels. Thus, further investigation is needed to recognize precisely the influence of all interactions on its properties.

Figure 9. Comparison of curves of reduced differential capacitance C*diff as a function of reduced potential Ψ* (details in text): circles represent data presented in Figure 8, and squares represent data presented in the Figure 1 in ref 10 (due to the kindness and permission of the authors).

about σ ) 0. A comparison of the differential capacitance results is shown in Figure 9. As seen, both curves are qualitatively similar. A nearly quantitative similarity is observed at the positive electrode potentials at which the adsorption of small ions (anions) takes place. The difference in C*diff at the negative potentials is greater, but it seems that both curves tend to the same value at large negative potentials. In the vicinity of zero potential, the curves have a distorted camel-shape (two nonsymmetrical maxima). The curve obtained by Fedorov and Kornyshev exhibits a stronger distortion perhaps due to the ratio of the cation/anion diameters. This ratio for the Fedorov and Kornyshev diameters is 2 while for our diameters the ratio is 1.25. Also, different shapes of the capacitance maxima can result from the difference in the reduced temperatures and in the packing fractions. Fedorov and Kornyshev used T* ) 0.040 and η ≈ 0.13 while our values were evidently greater and amounted to 0.32 and 0.35, respectively. The experimentally determined differential capacitance curves of ionic liquids are still very rare. Moreover, they strongly depend on the electrode material applied. We have found that a distorted camel-shape, similar to our MC predictions, was observed for 1-methyl-3-ethyimidazolium chloride, 1-methyl3-hexyimidazolium chloride on glassy carbon12 and for N,N,Ntrimethyl-N-propylammonium bis(trifluoromethanesulfonyl)imide on platinum.13 The corresponding ion diameter ratio of these salts varied from for 1.1 to about 1.8, while in our investigation this ratio was 1.25. The distorted camel-like shape of the Cdiff curve can probably be changed when including the specific (nonelectrostatic) interactions with the electrode material. Summary In this paper, we have presented MC simulation results of high density ionic system with ions of different diameters

Acknowledgment. The authors are very grateful to Professors M.V. Fedorov and A.A. Kornyshev for providing the MD capacitance data used in Figure 9. Financial support from Adam Mickiewicz University, Faculty of Chemistry, is appreciated. References and Notes (1) Lamperski, S.; Kłos, J. J. Chem. Phys. 2008, 129, 164503. (2) Gouy, G. J. Phys. (Paris) 1910, 9, 457. (3) Chapman, D. L. Philos. Mag. 1913, 25, 475. (4) Outhwaite, C. W.; Lamperski, S. Condens. Matter Phys. 2001, 4, 739. (5) Lamperski, S.; Outhwaite, C. W.; Bhuiyan, L. B. J. Phys. Chem. B 2009, 113, 8925. (6) Stafiej, J.; Ekoka, A.; Borkowska, Z.; Badiali, J. P. J. Chem. Soc. Faraday Trans. 1996, 92, 3677. (7) Fawcett, W. R.; Ryan, P. J.; Smagala, T. G. J. Phys. Chem. B 2009, 113, 14310. (8) Kornyshev, A. A. J. Phys. Chem. B 2007, 111, 5545. (9) Fedorov, M. V.; Kornyshev, A. A. Electrochim. Acta 2008, 53, 6835. (10) Fedorov, M. V.; Kornyshev, A. A. J. Phys. Chem. B 2008, 112, 11868. (11) Fedorov, M. V.; Georgi, N.; Kornyshev, A. A. Electrochem. Commun. 2010, 12, 296. (12) Lockett, V.; Sedev, R.; Ralston, J.; Horne, M.; Rodopoulos, T. J. Phys. Chem. B 2008, 112, 7486. (13) Islam, Md. M.; Alam, M. T.; Ohsaka, T. J. Phys. Chem. C 2008, 112, 16568. (14) Allen, M. P.; Tildsley, D. J. Computer Simulation of Liquids; Calderon Press: Oxford, U.K., 1987. (15) Lamperski, S. Mol. Simul. 2007, 33, 1193. (16) Torrie, G. M.; Valleau, J. P. J. Chem. Phys. 1980, 73, 5807. (17) Outhwaite, C. W. Can. J. Chem. 1981, 59, 1854. (18) Lamperski, S.; Zydor, A. Electrochim. Acta 2007, 52, 2429. (19) Cleveland, W. S. J. Am. Stat. Assoc. 1979, 74, 829. (20) Shannon, R. D.; Prewitt, C. T. Acta Crystallogr. 1969, B25, 925. (21) Larsen, B. J. Chem. Phys. 1976, 65, 3431. (22) Boda, D.; Henderson, D.; Chan, K. -Y. J. Chem. Phys. 1999, 110, 5346. (23) CRC Handbook of Chemistry and Physics, 86th ed.; Lide, D. R., Ed.; CRC: Boca Raton, FL, 2005; pp12-44. (24) Valisko´, M.; Boda, D.; Gillespie, D. J. Phys. Chem. C 2007, 111, 15575.

JP104402U