Single ion dynamics in molten sodium bromide O. Alcaraz, F. Demmel, and J. Trullas

Citation: The Journal of Chemical Physics 141, 244508 (2014); doi: 10.1063/1.4904821 View online: http://dx.doi.org/10.1063/1.4904821 View Table of Contents: http://aip.scitation.org/toc/jcp/141/24 Published by the American Institute of Physics

THE JOURNAL OF CHEMICAL PHYSICS 141, 244508 (2014)

Single ion dynamics in molten sodium bromide O. Alcaraz,1 F. Demmel,2 and J. Trullas1 1 2

Departament de Física i Enginyeria Nuclear, Universitat Politècnica de Catalunya, 08034 Barcelona, Spain ISIS Facility, Rutherford Appleton Laboratory, Didcot OX11 0QX, United Kingdom

(Received 11 July 2014; accepted 9 December 2014; published online 31 December 2014) We present a study on the single ion dynamics in the molten alkali halide NaBr. Quasielastic neutron scattering was employed to extract the self-diffusion coefficient of the sodium ions at three temperatures. Molecular dynamics simulations using rigid and polarizable ion models have been performed in parallel to extract the sodium and bromide single dynamics and ionic conductivities. Two methods have been employed to derive the ion diffusion, calculating the mean squared displacements and the velocity autocorrelation functions, as well as analysing the increase of the line widths of the self-dynamic structure factors. The sodium diffusion coefficients show a remarkable good agreement between experiment and simulation utilising the polarisable potential. C 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4904821] I. INTRODUCTION

Renewed interest has arisen in the understanding of molten salt properties on the microscopic scale. This research field was initiated in the early 1960s, and much effort has already been invested into these systems. A reason for the regained attractiveness may be the recent activities in the field of room temperature ionic liquids based on organic molecules.1,2 Here, the molten salts serve as the uncomplicated counterparts where the principles of particle motion in interacting Coulomb systems can be studied without the disturbing influence of internal molecular structure. On the other hand, molten salts exhibit highly interesting thermo-physical, chemical, and electrical properties, and a large number of technological applications exist where molten salts are employed as basic working fluids. The large heat capacities of molten salts make them attractive as heat transport media, as an example for heat storage in modern solar power plants or as coolant in next generation nuclear reactors. The most recent and the most popular application of molten salts is perhaps their use as the dielectric medium in high temperature fuel cells.3 The understanding of the structure and dynamics of charged fluids has much wider implications when one considers the apparent importance of Coulomb fluids in macromolecular liquids.4 One of the exciting questions that link application oriented research with fundamental physics is how the bulk properties of these systems are related to the particle behaviour on microscopic time and length scales. For this, however, a detailed knowledge of the microscopic structural and dynamic properties in molten salt systems is indispensable. The study of binary ionic liquids and, as their prototypes, molten alkali halides, has a long tradition in theory, experiment, as well as simulation.5–7 This research was motivated by the obviously simple interatomic interactions in these systems. From the structural point of view, Coulomb liquids exhibit short-range order due to the electrical charges, which are alternating in successive coordination shells around a central ion. Employing isotope substitution techniques in neutron 0021-9606/2014/141(24)/244508/9/$30.00

scattering the partial structure factors were already measured.8 In fact, a comparison of the resulting partial pair distribution functions showed that strong correlations exist up to relatively high distances, i.e., the maxima of the cationic g++(r) and g−−(r) are located at distances where the anionic g+−(r) comprise minima and vice versa. This solid-like structural arrangement is probably responsible for the peculiar dynamic response, which is similar to the solid state. A prominent feature in the collective dynamics of binary ionic liquids is the prediction of optic-type modes as was shown in a pioneering molecular dynamics (MD) simulation of a symmetric molten salt by Hansen and McDonald.9 In the following years, the research was driven by the investigation of collective dynamics in molten salts. The computer-based results motivated many neutron scattering groups to verify the calculated particle dynamics of the molten salts and especially to search for the predicted optic modes and the missing acoustic ones.10–12 No direct indications for collective modes were found in the spectra; only broad peaks could be detected in the current correlation functions. However, the advent of a new technique, inelastic x-ray scattering, paved the way to new high quality experimental results on collective dynamics in molten alkali halides.13 Transport properties, like diffusion, shear viscosity, and ionic conductivity have been studied in parallel through computer simulations.14,15 These quantities have a practical significance in, e.g., electrochemical devices. On a macroscopic scale, diffusion coefficients can be determined through tracer diffusion measurements.16 A restriction for this method is that convection will influence the results. A method that avoids this complication in one-component liquids is quasielastic neutron scattering. It probes the single particle dynamics on a microscopic scale and, from the dependence on momentum transfer, a self-diffusion coefficient can be derived. However, the extraction of useful information in binary systems as molten alkali halides is limited to a small number of ions with appropriate cross sections. In the best case, one ion is scattering incoherently only with a reasonable

141, 244508-1

© 2014 AIP Publishing LLC

244508-2

Alcaraz, Demmel, and Trullas

large cross section. One of the examples where the quasielastic neutron scattering was analysed is molten NaI.12 However, no diffusion was extracted from the data. First simulations of molten alkali halides used classical rigid ion pair-potentials to derive transport coefficients, such as self-diffusion coefficients, ionic conductivity, and shear viscosity,15 which were in quite good agreement with experimental data.17 The Tosi-Fumi pair potentials18 have been well studied for over 40 yr and are still used in contemporary publications.19 Though they neglect polarization effects, deliver a reasonable good description of the structure,20 and thermodynamics,21 of molten alkali halides. However, at the same time inclusion of polarization effects was undertaken for molten potassium iodide (KI) with an increase of the self-diffusion constant for the non-polarizable positive ion.14 Further studies assessed the influence of polarisation on the ion dynamics using shell models,22,23 fluctuating charge models,24 and polarizable point dipole models.25,26 More recently, ab-initio based methods have been applied to study molten alkali halides.27 Hazebroucq et al.28 determined radial distribution functions and selfdiffusion coefficients from a density functional-based MD simulation of molten NaCl and KCl. First principles of Hellmann-Feynman MD (HFMD) of molten NaCl and NaI were reported recently.29 We present here an experimental study on the single ionic dynamics in molten NaBr employing quasielastic neutron scattering (QENS) complemented by MD simulations using rigid and polarizable ion models. The neutron cross section for bromine is mostly coherent and sodium has an equal amount of incoherent and incoherent cross section. Then, at small momentum transfers, the signal will be dominated by the incoherent scattering from sodium, and hence the self-diffusion of the sodium ions only can be measured. The limited insight into the diffusion behaviour of molten NaBr from the experiment is complemented by extensive MD simulations based on rigid ion and polarizable ion models. Since the bromide polarizability is between that of chloride (low) and iodide (high) and polarization cannot be neglected in alkaline iodides,22 it is expected that the anion polarization in NaBr could play a noticeable role. Simulations will deliver diffusivities for both ions, and the influence of ion polarizability can be studied. Furthermore, correlations of ion motions can be elucidated, which may be responsible for deviations from the Nernst-Einstein relation for the ionic conductivity.

II. FORMALISM

The basic dynamical variables to describe the collective dynamic properties of binary charged systems with N+ cations and N− anions are the fluctuating local partial number densities of each species with respect to the averaged number partial density ρa ,δρa (r,t) = {Σiaδ[r − ria(t)]} − ρa , whose components in the reciprocal space (for wave numbers k , 0) are Na Na   ρa (k,t) = ρia(k,t) = exp[−ik · ria(t)], (1) ia=1

ia=1

where a can be either + or −, and ria(t) is the position at time t of the ion i a of species a. Their time correlations are the partial

J. Chem. Phys. 141, 244508 (2014)

intermediate scattering functions Fab(k,t) =

1 ⟨ρa (k,t)ρb (−k,0)⟩ , √ N ca cb

(2)

where N = N+ + N−, ca = Na/N, and the brackets denote the ensemble average over equilibrium configurations. Because of the liquids isotropy, the k vector dependence disappears in Fab(k,t). We define Fab(k,t) as Price and Copley,10 and Rovere and Tosi,6 whose initial time values are the AshcroftLangreth partial static structure factors Sab(k). Other authors, e.g., Hansen and McDonald,7 and McGreevy,11 do not include the factor (ca cb )−1/2. By using this factor, Sab(k) approaches unity as k → ∞ for a = b, while it approaches zero for a , b. The self-contribution to the autocorrelation Faa (k,t) is the selfintermediate scattering function, Fas (k,t) = ⟨ρia(k,t)ρia(k,0)⟩  Na  1  exp{−ik · [ria(t) − ria(0)]} , = Na ia

(3)

that is the autocorrelation of the reciprocal space components of the single density δ [r − ria(t)]. The partial density fluctuations spectra, namely the time Fourier transform ofFab(k,t), are the partial dynamic structure ∞ factors Sab(k,ω) = π1 0 Fab(k,t)cos(ωt)dt, and those of Fas (k,t) are the self-dynamic structure factors Sas (k,ω). The initial values of the intermediate scattering functions, which are the areas under the corresponding dynamic structure factors, are Fab(k,0) = Sab(k) and Fas (k,t) = 1. Concerning the behaviour of Sas (k,ω) in liquids, it decays monotonically and its halfwidth at half-maximum Γa (k), increases with k. This increase indicates that shorter wavelength fluctuations in the single particle densities die out more rapidly.30,31 Moreover, its height Sas (k,ω = 0) decreases as k increases because the area under Sas (k,ω) is the same at all ks. In the hydrodynamic limit (k → 0),Sas (k,ω) is a single Lorentzian centred at ω = 0,30 1 Γa (k) , π ω2 + Γa (k)2 with a half width at half maximum (HWHM), Sas (k,ω) =

Γa (k) = Da k 2

(for k → 0),

(4)

(5)

where Da is self-diffusion coefficient for translational diffusion on long distances of single a-type particles. The diffusion coefficient Da is related to the msd of single cations and anions, 

2  1

∆r ia(t) = |ria(t) − ria(0)|2 3 Na  1 1 

= |ria(t) − ria(0)|2 , (6) 3 Na ia=1 and their self-velocity autocorrelation functions (vcf), N

a 1 1 1  Λa (t) = ⟨via(t)·via(0)⟩ = ⟨via(t)·via(0)⟩, 3 3 Na ia

(7)

through the Einstein relation and the Kubo integral formula,

2   ∞ r (t) Da = lim ia = Λa (t)dt. (8) t→ ∞ 2t 0

244508-3

Alcaraz, Demmel, and Trullas

J. Chem. Phys. 141, 244508 (2014)

The second equality is a particular case of the general relation between any variable A(t) and its time 

dynamical  7,30 ˙ lim |∆A(t)|2 /2t = ∞⟨ A(t) ˙ A(0)⟩dt. ˙ derivative A(t), Any 0 t→ ∞ empirical transport constant is related to a general mean square ˙ deviation of A(t) and the autocorrelation function of A(t) through Einstein- and Kubo-like formulae. For the ionic conductivity κ, they are the mean square displacement (MSD)  of the charge-centre position rz = i z i ri /N, where z i is the ionic charges in units of the fundamental charge e,

  N  rz (t) − rz (0)|2 , ∆r z2 (t) = 3

(9)

and the autocorrelation function of the charge-centre velocity  vz = i z i vi /N, Λz (t) =

N ⟨vz (t) · vz (0)⟩, 3

(10)

and κ=

⟨r 2 (t)⟩ ρ N e2 ρ N e2 lim Z = k BT t→ ∞ 2t k BT



∞

Λ Z (t)dt,

(11)

0

where ρ N = N/V is the ionic density and k B the Boltzmann’s constant. For molten alkali halides, it can be written as κ=

ρ N e2 1 (D+ + D−)(1 − ∆) = κ NE(1 − ∆), kBT 2

(12)

where κ NE is the Nernst-Einstein approximation, and the deviation factor ∆ is usually a small positive quantity that accounts for correlations between distinct ions.7 The origin of this deviation can be understood, in an averaged sense, as a motion of both ions in the same direction, which contributes to the diffusion coefficient but not to the electrical conductivity. III. EXPERIMENTAL DETAILS

Inelastic neutron scattering measures the total dynamic structure factor S(k,ω) of the density fluctuations. Nonmagnetic neutron scattering interacts with the nuclei through two cross sections, the incoherent and coherent one. The incoherent cross section connects with the self-correlations of a single particle, and the coherent cross section interferes with correlations between the particles. In general, the total double differential cross section per particle is a cross section weighted sum of these two contributions, Scoh(k,ω) and Sinc(k,ω),10 d 2σ d 2σ =~ dΩdω dΩdE f  k σ k f  σcoh σinc s f T Scoh (k,ω) + S (k,ω) = S(k,ω), = k0 4π 4π k 0 4π (13) where k f and k0 are the modulus of the scattered and incident neutron wave vectors k f and  k0, the momentum and energy
transfers are ~k = ~ k f − k0 and ~ω = E f − E0 = ~ k12 − k 02 / 2m n with m n the neutron mass. In a multi-component system, the cross sections are averaged values, which for alkali halides are σcoh = (σcoh+ + σcoh−)/2 and σinc = (σinc+ + σinc−)/2, with

σT = σcoh + σinc, while 4π 1  2 b S++ (k,ω) + 2b+b− S+− (k,ω) σcoh 2 +  + b2− S++ (k,ω) , (14)  1 1 σinc+ S+s (k,ω) + σinc− S−s (k,ω) , (15) S s (k,ω) = σinc 2

Scoh (k,ω) =

where b2a = σcoh,a/(4π) is the squared coherent scattering length of the a-type ions. Replacing Sab(k,ω) and Sas (k,ω) by Fab(k,t) and Fas (k,t), we get the coherent and selfintermediate scattering functions, whose initial values are S s (k) = 1, and Scoh(k) = 21 [b2+ S++(k)+2b+b− S+−(k)+ b2− S−−(k)] is the coherent static structure factor normalized in such a way it approaches to unity as k→ ∞. Under conditions of quasielastic neutron scattering k f ≈ k0 and, thus, we are exploring long wavelengths near the hydrodynamic regime (k → 0 and ω → 0). Furthermore, the factor k f /k0 ≈ 1 can be neglected in the above equations. Table I provides the neutron cross sections and the scattering lengths of both ions.32 Since σinc,Na is quite large in contrast to σinc,Br, in molten NaBr S s (k,ω) is practically determined s by SNa (k,ω). Furthermore, since Scoh(k) at small wave vectors is much lower than S s (k) = 1, the measured intensity at small momentum transfers is dominated by the incoherent scattering from sodium, and hence the self-diffusion of the sodium ions can be extracted from the relation Γ ≈ DNa k 2. NaBr powder was filled into a flat niobium can and then electron beam welded. The niobium wall thickness was 0.5 mm and the sample thickness 3.8 mm, which provided a scattering power of about 12%. Niobium is a nearly perfect coherent scatterer, and hence will not contribute to the elastic line except where Debye-Scherrer lines appear. The first reflection of niobium is at k = 2.7 Å−1, which is beyond the first structure factor maximum of molten NaBr (k ≈ 1.7 Å−1). The cell was installed under a 45◦ orientation in a transmission geometry into a standard furnace with niobium shields. The measured temperatures were 1043 K, 1123 K, and 1223 K for molten NaBr (Tmelt = 1000 K). The temperature uncertainty was smaller than ±1.5 K during all measurements. An identical cell was used for empty cell runs. QENS measurements were performed at the OSIRIS spectrometer at the ISIS Facility, UK. This instrument is an indirect time of flight backscattering spectrometer that was operated with a final energy E f = 1.845 meV. The energy resolution deduced from a vanadium measurement was 0.025 meV. The data analysis included monitor normalisation, empty cell subtraction, and conversion into constant k-spectra. Absorption coefficients for the empty cell subtraction have been calculated according to a method of Paalman and Pings.33 TABLE I. Neutron scattering lengths, cross sections, and absorption cross sections.32

Na Br

ba (fm)

σ coh, a (barn)

σ inc.a (barn)

σ abs, a (barn)

3.63 6.795

1.66 5.80

1.62 0.1

0.53 0.9

244508-4

Alcaraz, Demmel, and Trullas

J. Chem. Phys. 141, 244508 (2014)

Fig. 1 shows an intensity map of molten NaBr at 1043 K. Towards small k-vectors, a strong increase in the amplitude can be observed, which decreases with increasing k. This shape demonstrates that the measured signal stems from sodium diffusion within this range of k-vectors. At energy transfers away from the quasielastic signal, a small increase in intensity towards larger momentum transfers can be noticed, which is the result of the smoothly increasing coherent structure factor. Measured spectra are shown in Fig. 2 for a few selected momentum transfers and the temperature dependence. The lines are fits with a single Lorentzian convoluted with the measured resolution function, which describes quite well the line shape in the considered k-range. Included into the fit was a linear sloping background and, if necessary, a delta function for not completely subtracted elastic contributions.

IV. MOLECULAR DYNAMICS SIMULATIONS

Using our own software, we carried out MD simulations of a rigid and a polarizable ion model for molten NaBr. The rigid ion model (RIM) is that given by the Fumi-Tosi effective pair potential.18 It is in the Born-Mayer-Huggins form,34 Cab Dab z a z b e2 (16) + βabb exp[−γ(r − Rab)] − 6 − 8 , r r r where Cab and Dab are the van der Waals coefficients with the Mayer values35 given in Table II. The Pauling factors in the overlap repulsive term are β++ = 1.25, β+− = 1, and β−− = 0.75, and the strength parameter b = 0.221 eV is the same for all alkali halides. The hardness parameter is γ = 3.155 for NaBr, and Rab = Ra + Rb with RNa = 1.17 Å and RBr = 1.716 Å closely related to the ionic radii. The polarizable ion model (PIM) is constructed by adding the many-body induced polarization interactions to the FumiTosi pair potential φab(r). We assume that a point dipole is induced on an ion placed at position ri as result of two contributions. The first one is the dipole induced by the local electric field Ei due to all other ions, whose moment is given by the electronic polarizability α i in the linear approximation φab(r) =

FIG. 2. Spectra S(k, ω) of molten NaBr for (a) different small momentum transfers at 1043 K and (b) different temperatures for k = 0.35 Å−1. Included in plot (a) is the energy resolution from a vanadium scan and the empty cell contribution. The lines in plot (b) depict the fits with Lorentzian functions.

α i Ei . The second is the short-range damping contribution due to the neighbouring ions. According to Wilson and Madden,24 the constitutive relation for the resulting dipole moment can be written as pi = α i Ei − α i

N 

z je

f (r ij)

r ij3

j,i

rij,

(17)

where rij = ri − r j and r ij is its modulus, Ei = Eqi + Eip =

N  z je j,i

r ij3

rij +

N 

*3 (p j · rij) rij − 1 p j + , (18) r ij5 r ij3 j,i ,

is the local field at ri with Eqi , the field due to all point charges except qi = z i e and Eip that due to all dipole moments except pi . A convenient form for f (r) is the Tang and Toennies dispersion

TABLE II. Mayer values of van der Waals coefficients.35 ++ 6

FIG. 1. Intensity map against momentum and energy transfer for molten NaBr at 1043 K. Towards small k-vectors, a strong increase in amplitude occurs from the incoherent scattering of the sodium ions.

C ab (eVÅ ) 8

Dab (eVÅ )

+−

−−

1.05

8.74

122.33

0.50

11.86

280.87

244508-5

Alcaraz, Demmel, and Trullas

J. Chem. Phys. 141, 244508 (2014)

damping function,36 f (r) = exp(−r/ρ)

4  (r/ρ)k k=0

k!

.

(19)

By choosing this form, it is ensured that the short-range terms cancel out te z j e/(r ij)3 contributions of Ei when two ions approach at too short distances. The potential energy of this polarizable ion model may be conveniently written as U=

N N N N  1 1  φij(r ij) − pi · Eqi − pi · Eip 2 i=1 j,i 2 i=1 i=1

+

N N  N  pi2  z je f (r ij) 3 rij · pi . + 2α i i=1 j,i r ij i=1

(20)

In carrying out the MD simulations, we actually have to calculate the forces. The force acting in an ion is that given in Ref. 37 for simple PIM with an extra term due to the shortrange damping interactions as given in Ref. 38. For NaBr, we assume that only the anions are polarizable with αNa = 0 and αBr = 4.16 Å3 because the former is much lower than the latter.39 The value of the damping parameter ρ has been estimated according to the relationship ρ = (σ+ + σ−)/c proposed by Wilson and Madden, with c = 4.65 and 7.42 for divalent and trivalent metal halides, respectively.40 Taking into account the ρ value used in their simulations of molten LiF,24 we get c = 4.1 for alkali halides and ρ = 0.71 Å for NaBr. In MD simulations of copper and silver halides assuming α+ = 0 (and effective charges |z a | < 1), the short-range damping polarization it is not needed.37,38,41 However, in molten alkali halides (with |z a | = 1), it cannot be neglected because two unlike ions can approach to unphysical separations and become over-polarized in such a way polarization catastrophe occurs.25 Furthermore, if α+ is not neglected in copper and silver halides, a more sophisticated short-range damping polarization must be introduced to avoid unphysical over-polarizations.42 By using the potential models described above, we simulated molten NaBr at the measured at temperatures 1043 K, 1123 K, and 1223 K, with densities 0.02716, 0.02645, and 3 0.02551 ions/Å , respectively.43 MD simulations were carried out at constant energy, following the simple procedure described in Ref. 25 to reach the desired temperatures, considering 1000 ions placed in a cubic box of side L with periodic boundary conditions, and Beeman’s algorithm34 was used with a time step of 1 fs = 10−15 s. Once equilibrium was achieved, the properties were averaged over 450 000 time steps in three runs of 150 000 (around 20 h each in a single processor) with energy keeping practically constant (it changes less than 0.1%). This large number of time steps is necessary to reduce the statistical fluctuations in the partial intermediate scattering functions, which decay very slowly at low wavenumbers and because the ionic conductivity is a collective transport property that needs to be averaged over long simulations. The electric fields and the corresponding long-range interactions were evaluated by the Ewald method34 with the Ewald parameter equal to 6.5/L. The real space terms were truncated at distances longer than L/2, and the reciprocal space contributions of wave vectors k = n(2π/L), where n being a vector of integer

components, were truncated for n2 > 12. The dipole moments have been evaluated by using the prediction-correction iterative method proposed by Vesely,44 with the same convergence limit used in Ref. 25. Other computational details are those described in Refs. 25 and 37. In Fig. 3 a comparison is made between the PIM and RIM structural properties at 1043 K. The partial structure factors Sab(k) were computed by a hybrid method described in Ref. 25. The differences between PIM and RIM are very small, as higher the temperature, smaller the differences, and they can only be distinguish in the Na-Na correlations, as it was also observed in MD simulations of molten NaI,25 but less marked because the bromide polarizability is lower than the iodide one. Due to the induced dipoles on bromides, their negative ends attract the sodium neighbours in such a way the repulsion between them is screened. Then, their separation can be smaller than it would be if the bromides were not polarized, and the first peak of the PIM gNaNa(r) is slightly shifted to lower distances than that of gBrBr(r). Moreover, when polarization effects are added, the gNaNa(r) exhibits a more pronounced penetration into the first coordination (the distances d Na at which it becomes significantly different from zero is slightly lower), and its maximum becomes slightly lower. Since, to a certain extent, d Na is a measure of the effective ionic diameter,

FIG. 3. (a) Pair correlation functions g ab(r ) at 1043 K for both PIM and RIM. (b) Partial static structure factors Sab(k) for PIM at 1043 K. The inset shows a comparison of SNaNa(k) between PIM and RIM.

244508-6

Alcaraz, Demmel, and Trullas

J. Chem. Phys. 141, 244508 (2014)

the sodium ions can behave as they were slightly smaller. This phenomena were already observed in MD simulations of PIM for molten silver and copper halides,37,38,41,42 but much more intense because the short-range damping polarization were not needed. In reciprocal space, the main difference is a small reduction in the peak height of the first sharp diffraction peak of SNaNa(k) with inclusion of polarization. It appears that the sodium ion correlations are slightly decreased. This increased disorder in the cation arrangement might be due to the increased accessible space from the polarisable bromide. 2 To study the ion dynamics, we computed the MSDs ⟨r ia (t)⟩ 2 and ⟨r Z(t)⟩, as well as their velocity counterparts Λa (t) and Λ Z (t), and by using the Kubo and Einstein relations, we calculated Da , κ, and ∆, as described in Ref. 25. The estimated uncertainties for Da are lower than 5%, whereas those for κ are about 10% despite the long simulations. The differences between the values calculated using the Kubo and Einstein relations are within the error intervals. Furthermore, we also computed the self-intermediate scattering functions Fas (k,t), and their spectra Sas (k,ω), as described in Ref. 19, for the 13 lowest values of k allowed by the periodic boundary conditions −1 up to 0.71 Å to cover the same range as QENS measurements. Despite the inclusion of the anion polarization leads to very small changes in the structure, the effects on the transport properties are noticeable, mainly in the sodium diffusivity. In Fig. 4 we present the normalised vcf Ca (t) = Λa (t)/ Λa (0) for sodium and bromide. There is a clear difference between both ions due to the mass differences. The sodium CNa(t) shows a quite well expressed negative part and oscillate beyond the minimum, whereas the bromide CBr(t) decays slower and only shows a small negative part, as it is usual for the light and heavy ions of molten alkali halides.45 The nega tive part of CNa(t) can be understood as a backscattering effect for the specific particle. The sodium ion cannot immediately escape from the cage of the surrounding bromine ions and is scattered back until it finally escapes. The polarization of the surrounding bromine ions facilitate the escape chances of the sodium ions, and hence the negative part of CNa(t) is reduced and sodium diffusivity is increased. However, for the bromide, there is only a small difference for the two potentials. Similar

At first, we will discuss the results for the sodium diffusion. Fig. 5 shows the extracted HWHM against k 2 from the measured S(k,ω) at the three temperatures. The linear dependence demonstrates that the observed dynamics stems from translational diffusion of the sodium ions. Included is a linear fit to extract the diffusion coefficient according to Γ = Dk2. In Fig. 6 we show the HWHM extracted from the sodium s spectra SNa (k,ω) of both RIM and PIM simulations at 1043 K. Since the introduction of polarization facilitates the diffusivity of the positive ions, as we see above, the slope for PIM results is higher than that for RIM. Included are the values from the neutron experiment. There is a remarkable good agreement between the experiment and the results from the PIM potential.

FIG. 4. Normalised vcf C a (t) = Λa (t)/Λa (0) for the sodium and bromine ions for both PIM and RIM. The main change occurs for the sodium ions.

s (k, ω) for both RIM and PIM FIG. 6. The HWHM against k 2 from the SNa simulations at 1043 K. Included are the QENS results.

FIG. 5. The HWHM against k 2 extracted from the measured S(k, ω) for all three different temperatures. The lines are linear fits to extract the sodium diffusion constant.

trends were observed for molten KI14 and NaI25 utilising a polarizable potential.

V. RESULTS AND DISCUSSION

244508-7

Alcaraz, Demmel, and Trullas

J. Chem. Phys. 141, 244508 (2014)

TABLE III. Sodium self-diffusion coefficient values, D Na in units of 10−5 cm2/s, from QENS data and PIM (and RIM) results.

TABLE IV. Bromide self-diffusion coefficient values, DBr in units of 10−5 cm2/s, from PIM (and RIM) results.

T (K)

QENS

MD HWHM

MD msd-vcf

T (K)

MD HWHM

MD msd-vcf

1043 1123 1223

8.41 ± 0.02 10.70 ± 0.02 13.41 ± 0.02

8.5 (6.3) 10.2 (7.7) 12.7 (9.5)

9.3 (6.7) 11.1 (8.2) 13.6 (10.2)

1043 1123 1223

5.8 (5.1) 7.2 (6.3) 8.8 (7.8)

6.2 (5.5) 7.5 (6.7) 9.3 (8.3)

In Table III, we present the extracted diffusion coefficients DNa for all three temperatures, and in Fig. 7 the results are plotted against the inverse temperature. The errors in the DNa values from the QENS data are stemming from the fit errors. Note the good agreement between the DNa from QENS and DNa from HWHM of the PIM simulated spectra. The diffusivities s obtained through the HWHM of SNa (k,ω) are systematically smaller than the ones calculated from the msd or integrating the vcf. This difference is due to the approximat nature of a single Lorentzian to describe the line shape of a diffusive process in the non-hydrodynamic regime. Deviations from a simple Lorentzian in space-dependent diffusion are becoming increasingly important with increasing k-vectors.31 Thus, from msd-vcf results, we suspect that the actual sodium diffusion may be slightly higher than the values extracted from QENS. The fits in Fig. 2 demonstrate that a Lorentzian is a very good description for the line shape but also show some deviations, in particular, at the peak. The simulated DNa from RIM are about 30% smaller than the ones from PIM. A similar difference for the sodium movements in a polarisable anion environment has been observed in a molten NaI simulation,22,25 even though DNa from the polarisable model still underestimated the experimental value at that time. Included in Fig. 7 are measured sodium diffusion coefficient of molten NaCl and NaI obtained through tracer diffusion.16 The sodium ions show a very similar diffusivity in all alkali halides. Since the bromide self-diffusion coefficient, DBr, cannot be obtained from QENS, we rely on the simulation results provided in Table IV. The changes between the RIM and PIM are smaller compared to the ones for the sodium ion. There is an increase of DBr with inclusion of polarization, however, only

in the range of 10%. From tracer measurements, no diffusion value for the bromide is available to our knowledge. In Fig. 8 the simulated and measured anion diffusion coefficient of molten NaCl and NaI obtained through tracer diffusion14 are plotted against the inverse temperature. The microscopic steps leading to the diffusion process can be regarded as thermally activated jump processes. These processes result in an Arrhenius-type temperature dependence for the diffusion constant: D = D0 exp(−E/k BT), where E denotes the activation energy. From the temperature dependence of DNa, the activation energy can be derived. Analysing the temperature dependence accordingly, we get from the QENS data an activation energy of E = 285 meV. From the MD simulation results of the HWHM (and msd-vcf), we get E = 257 meV (252 meV) for RIM and E = 235 meV (243 meV) for PIM. For bromine ions, we obtain E = 244 meV (252 meV) for RIM and E = 252 meV (263 meV) for the PIM. Young and O’Connell17 fitted the experimental diffusion values assuming the same activation energy for both cation and anions. They get 289 meV for NaCl and 282 meV for NaI, very close to the value obtained for sodium ions from QENS data of molten NaBr. In Table V, we present the ionic conductivity (κ) and deviation factor (∆) values calculated for PIM (and RIM). The former are compared with the available experimental data.16,17 Although the DNa values for PIM are in better agreement with QENS measurements, the RIM conductivity results are closer to the experimental data. Since the PIM conductivity values are higher than the experimental ones, we conclude that polarizable effects, as well as details of the Fumi-Tosi effective pair potential used in PIM, lead to so high values of

FIG. 7. Sodium self-diffusion coefficients derived from QENS and MD results for RIM and PIM against the inverse temperature. Included are values for molten NaCl and NaI obtained from tracer diffusion experiments.16,17

FIG. 8. Bromide self-diffusion coefficients from MD results for RIM and PIM against the inverse temperature. Experimental values from tracer diffusion experiments for anion diffusion in molten NaCl and NaI are included.16,17

244508-8

Alcaraz, Demmel, and Trullas

J. Chem. Phys. 141, 244508 (2014)

TABLE V. Experimental data of ionic conductivity, κ in units of (Ω cm)−1, from Bockris and Hooper,16 also given by Young and O’Connel,17 as well as the PIM (and RIM) results for κ and the deviation factor ∆. T (K)

κ exp

κ PIM (RIM)

∆ PIM (RIM)

1043 1123 1223

2.96 (2.82) 3.2 (3.06) 3.4 (3.39)

3.44 (3.00) 3.57 (3.09) 3.61 (3.54)

0.09 (0.00) 0.14 (0.06) 0.20 (0.03)

DBr. The ∆ values for RIM are close to zero and do not change with temperature, while the PIM results show an increase of ∆ with temperature. From the experimental data collected by Janz et al.,46 it is found that the ∆ value for NaCl increases with temperature from 0.08 at 1123 K to 0.22 at 1323 K. One might argue that with rising temperature, a more stochastic motion of the particles, and hence a loss of correlation is expected. However, these Coulomb liquids provide here the opposite picture. This trend point to a behaviour that with reduced density, both ions are more correlated. From the Young and O’Connel data,17 ∆ is practically constant around 0.2. An experimental approach to extract the bromine diffusion coefficient can be made from the measured κ and DNa values, and using Eq. (12) where the deviation factor ∆ is a priori unknown value that can only be derived from simulation or when diffusion measurements of both ions are available. Experimentally, for molten alkali halides, ∆ use to be in the range 0.1–0.2.14 Using these values in Eq. (12) with the measured κ in Table V by Bockris and Hooper16 and the measured DNa of this paper, the range values of DBr at 1043 K is 5.0–5.7·10−5 cm2/s, at 1123 K is 6.8–9.0·10−5 cm2/s and at 1223 K is 5.8–8.2·10−5 cm2/s. The DBr values from PIM simulations exceed these ranges as it has been predicted above from the so height values for the PIM conductivities.

VI. CONCLUSIONS

In this work, we extracted a single ion diffusion constant in a binary Coulomb liquid through QENS. Even though the total dynamic structure factor consists of five terms, with a sensible choice of the cross sections and wave vector region, the scattered intensity is dominated by a single s one SNa (k,ω). In this vein, analysing the increase of the line widths of the total spectra at low wavenumbers, we extracted the sodium diffusion coefficient in molten NaBr at three temperatures (1043, 1123 and 1223 K). In parallel, we performed a comprehensive MD study on molten NaBr at the same temperatures. Two methods were employed to derive the ion diffusion, calculating the msd and the vcf, as well as analysing the line widths of the self-dynamic structure factors. A rigid ion model potential (RIM) and a polarizable ion model potential (PIM) were simulated. RIM underestimated the sodium diffusion constants, which are up to 30% smaller than the experimental ones. Introducing anion polarization, although structural changes are not significant, ionic diffusion increases and a remarkable good agreement with the experimental data was achieved. Since the selfdiffusion values obtained from the msd-vcf are slightly higher

than from the line widths, we believe that the actual sodium self-diffusion may be slightly higher than the values extracted from QENS. From MD simulations, the bromine diffusion coefficient, for which no experimental value is available, was obtained, as well as the ionic conductivity and the Nernst-Einstein deviation factor that accounts for correlations between distinct ions. The PIM conductivities are slightly higher than the experimental ones, and the value of the deviation factor increases with temperature from 0.09 at 1043 K to 0.20 at 1223 K, as it was found from experimental data of molten NaCl. From the experimental ionic conductivity and sodium diffusion constant and assuming that the deviation factor lies between 0.1 and 0.2, we estimated the range values for the actual bromine diffusion constant. Since PIM values exceed these ranges and PIM conductivities are higher than the experimental ones, we conclude that the actual bromine diffusion coefficients may be slightly lower than the PIM ones. In spite of these small inaccuracies, the good agreement between experiment and simulation for the sodium movements makes us confident to obtain further quantities from simulations with this potential, in particular, properties not easily accessible from experiment, e.g., shear and bulk viscosities or thermal conductivities, which might be relevant for practical applications.

ACKNOWLEDGMENTS

O. Alcaraz and J. Trullàs thank Generalitat de Catalunya for Grant No. 2009SGR-1003 and MINECO of Spain for Grant No. FIS20012-39443-C02-01. We are grateful to the ISIS furnace section for the excellent support. 1Ionic Liquids in Synthesis, edited by P. Wasserscheid and T. Welton

(WileyVCH, Weinheim, 2007). 2N. V. Plechkova and K. R. Seddon, Chem. Soc. Rev. 37, 123 (2008). 3J. R. Selman, J. Power Sources 160, 852 (2006). 4A. Naji, M. Kanduc, J. Forsman, and R. Podgornik, J. Chem. Phys. 139, 150901 (2013). 5N. H. March and M. P. Tosi, Coulomb Liquids (Academic Press, San Diego, 1984). 6M. Rovere and M. P. Tosi, Rep. Prog. Phys. 49, 1001 (1986). 7J. P. Hansen and I. McDonald, Theory of Simple Liquids (Academic Press, London, 2006). 8F. G. Edwards, J. E. Enderby, R. A. Howe, and D. I. Page, J. Phys. C 8, 3483 (1975); S. Biggin and J. E. Enderby, ibid. 15, L305 (1982); J. Y. Derrien and J. Dupuy, J. Phys. Arch. 36, 191 (1975); J. Locke, S. Messoloras, R. J. Stewart, R. L. McGreevy, and E. W. J. Mitchell, Philos. Mag. Part B 51, 301 (1985); E. W. J. Mitchell, P. F. J. Poncet, and R. J. Stewart, Philos. Mag. 34, 721 (1976); M. A. Howe, R. L. McGreevy, and W. S. Howells, J. Phys.: Condens. Matter 1, 3433 (1989). 9J. P. Hansen and I. R. McDonald, Phys. Rev. A 11, 2111 (1975). 10D. L. Price and J. R. D. Copley, Phys. Rev. A 11, 2124 (1975). 11R. L. McGreevy, Solid State Phys. 40, 247 (1987). 12R. L. McGreevy, E. W. J. Michell, and F. M. A. Margaca, J. Phys. C 17, 775 (1984). 13F. Demmel, S. Hosokawa, M. Lorenzen, and W. C. Pilgrim, Phys. Rev. B 69, 012203 (2004); F. Demmel, S. Hosokawa, and W. C. Pilgrim, J. Alloys Compd. 452, 143 (2008); S. Hosokawa, F. Demmel, W. C. Pilgrim, M. Inui, S. Tsutsui, and A. Baron, Electrochemistry 77, 608 (2009). 14G. Jacucci, I. R. McDonald, and A. Rahman, Phys. Rev. A 13, 1581 (1976). 15G. Ciccotti, G. Jacucci, and I. R. McDonald, Phys. Rev. A 13, 426 (1976). 16J. O. Bockris and G. W. Hooper, Discuss. Faraday Soc. 32, 218 (1961). 17R. E. Young and J. P. O’connell, Ind. Eng. Chem. Fundam. 10, 418 (1971). 18M. P. Tosi and F. G. Fumi, J. Phys. Chem. Solids 25, 45 (1964). 19O. Alcaraz and J. Trullas, J. Chem. Phys. 132, 054503 (2010).

244508-9 20M.

Alcaraz, Demmel, and Trullas

Dixon and M. J. Gillan, Philos. Mag. Part B 43, 1099 (1981). Guissani and B. Guillot, J. Chem. Phys. 101, 490 (1984). 22M. Dixon, Philos. Mag. Part B 47, 509 (1983); 47, 531 (1983); 48, 13 (1983). 23P. J. Mitchell and D. Finchman, J. Phys.: Condens. Matter 5, 1031 (1993). 24M. Wilson and P. A. Madden, J. Phys.: Condens. Matter 5, 2687 (1993). 25O. Alcaraz, V. Bitrian, and J. Trullas, J. Chem. Phys. 127, 154508 (2007). 26V. Bitrian, O. Alcaraz, and J. Trullas, J. Chem. Phys. 130, 234504 (2009). 27G. Tabacchi, C. J. Mundy, J. Hutter, and M. Parrinello, J. Chem. Phys. 117, 1433 (2002). 28S. Hazebroucq, G. S. Picard, C. Adamo, T. Heine, S. Gemming, and G. Seifert, J. Chem. Phys. 123, 134510 (2005). 29N. Galamba and B. J. Costa Cabral, J. Chem. Phys. 126, 124502 (2007); 127, 094506 (2007). 30U. Balucani and M. Zoppi, Dynamics of the Liquid State (Clarendon Press, Oxford, 1994). 31J. P. Boon and S. Yip, Molecular Hydrodynamics (McGraw Hill, 1980). 32V. F. Sears, Neutron News 3(3), 26–37 (1992). 33H. H. Paalman and C. J. Pings, J. Appl. Phys. 33, 2635 (1962). 34M. J. L. Sangster and M. Dixon, Adv. Phys. 25, 247 (1976). 21Y.

J. Chem. Phys. 141, 244508 (2014) 35J.

E. Mayer, J. Chem. Phys 1, 270 (1933). T. Tang and J. P. Toennies, J. Chem. Phys. 80, 3726 (1984). 37J. Trullas, O. Alcaraz, L. E. Gonzalez, and M. Silbert, J. Phys. Chem. B 107, 282 (2003). 38V. Bitrian and J. Trullas, J. Phys. Chem. B 110, 7490 (2006). 39J. R. Tessman, A. H. Kahn, and W. Shockley, Phys. Rev. 92, 890 (1953). 40M. Wilson and P. A. Madden, Mol. Phys. 92, 197 (1997); P. A. Madden and M. Wilson, J. Phys.: Condens. Matter 12, A95 (2000); F. Hutchinson, M. Wilson, and P. A. Madden, Mol. Phys. 99, 811 (2001). 41V. Bitrian, J. Trullas, and M. Silbert, J. Chem. Phys. 126, 021105 (2007); O. Alcaraz, V. Bitrian, and J. Trullas, ibid. 134, 014505 (2011); V. Bitrian, O. Alcaraz, and J. Trullas, ibid. 134, 044501 (2011). 42V. Bitrian and J. Trullas, J. Phys. Chem. 112, 1718 (2008). 43G. J. Janz, Molten Salts Handbook (Academic Press, New York, 1967). 44F. J. Vesely, J. Comput. Phys. 24, 361 (1977). 45O. Alcaraz and J. Trullas, J. Chem. Phys. 113, 10635 (2000). 46G. J. Janz, C. B. Allen, N. P. Bansal, R. M. Murphy, and R. P. T. Tomkins, Physical Properties Data Compilations Relevant to Energy Storage II. Molten Salts: Data on Single and Multi-Component Salt Systems, NSRDSNBS 61, Part II (U.S. GPO, Washington, 1979). 36K.