Simulation of the liquid–vapor interface of molten LiBeF<sub>3</sub>

Mathieu Salanne; Christian Simon; Pierre Turq; Paul A. Madden

doi:10.1016/j.crci.2007.03.002

Résumés

Anglais
Français

The structural properties of the liquid–vapor interface of an equimolar mixture of LiF and BeF₂ are studied by molecular dynamics. The interaction potential include anion polarization and was constructed from first-principle calculations. The simulation result for the surface tension is in good agreement with the experimental measure. Examination of the electrostatic potential variation across the interface reveals an important charge-separation effect at the interface, which is mostly due to the charge asymmetry between the two cations, Be²⁺ and Li⁺.

Les propriétés structurales de l'interface liquide–vapeur d'un mélange équimolaire de LiF et de BeF₂ sont étudiées par des simulations de dynamique moléculaire. Le potentiel d'interaction inclut la polarisation de l'anion et a été construit à partir de calculs ab initio. Le résultat obtenu pour la tension de surface est en bon accord avec les mesures expérimentales. L'examen de la variation du potentiel électrostatique à travers l'interface révèle d'importants effets de séparation de charge, qui sont notamment dus à l'asymétrie de charge entre les deux cations présents dans le système, Be²⁺ et Li⁺.

1 Introduction

Interfacial properties of ionic liquids have been examined by computer simulation in recent years. These studies concern simple alkali halide electrolytes as well as the more complicated room-temperature ionic liquids. Size asymmetry effects are now well characterized, but no consideration has yet been given to the role of charge asymmetry. Recently, we have developed a first-principle-based simulation model for mixtures of lithium fluoride (LiF) with beryllium fluoride (BeF₂), which accurately reproduced the properties of the real liquid [1]. In particular, the density of the melts was accurately reproduced and the computation of dynamic properties like the conductivity or viscosity gave results in perfect agreement with experimental data. The simulations allow the nature of the various fluoroberyllate ions $({Be}_{n} F_{3 n + 1}^{(n + 1) -})$ and their concentration to be determined for different compositions of the melt. As in this system the beryllium cation carries a +2 formal charge, it presents an opportunity to examine the importance of charge asymmetry on the interfacial properties as well as the way the fluoroberyllate ions are accommodated at the interface.

Furthermore, from a practical point of view, it is crucial to understand interfaces, since their properties can dramatically impact kinetics, or transport in systems involving liquid molten salts. “FLiBe” itself is an important solvent and heat-transfer material in numerous proposed nuclear technologies, including the molten-salt reactor and the fusion-breeder reactor. In the case of Generation IV molten-salt reactors, interfaces will certainly play a role in corrosion, and predicting wettability is also important to assign good boundary conditions in thermo-hydraulic process models. To address feasibility and reliability of such predictions, we simulated a liquid–vapor interface of molten LiBeF₃ (equimolar LiF–BeF₂). The technical details of the model, which accounts for anion polarization, and methods used are first described. The second part of this work is devoted to the study of the structure of the interface, which differs a lot from the charge-symmetric system: in particular an important charge separation phenomenon appears, and the local concentration of beryllium cations close to the surface is enhanced. Finally the profile followed by the electrostatic potential along the liquid is presented, and the surface tension is estimated: the result obtained matches well with the experimental data.

2 Technical details

The interaction potential used in those simulations is a “polarizable-ion model” and consists of a pair potential of Born–Mayer form together with an account of ionic polarization [2]. The pair potential is written as:

V (r_{i j}) = B^{i j} e^{- α^{i j} r_{i j}} - f_{i j}^{6} (r_{i j}) \frac{C_{i j}^{6}}{r_{i j}^{6}} - f_{i j}^{8} (r_{i j}) \frac{C_{i j}^{8}}{r_{i j}^{8}}

(1)

where

C_{i j}^{6}

and

C_{i j}^{8}

are the dispersion coefficients and fⁿ_ij are dispersion damping functions [3] given by:

f_{i j}^{n} (r_{i j}) = 1 - e^{- b_{i j}^{n} r_{i j}} \sum_{k = 0}^{n} \frac{{(b_{i j}^{n} r_{i j})}^{k}}{k!}

(2)

Values for all parameters necessary to simulate mixtures of LiF with BeF₂ are given in Table 1. A full account of the determination of these parameters is given in Ref. [4]. All the parameters are given in atomic units; the atomic unit of length is the Bohr radius, 0.52918 Å, and of energy, the Hartree, 4.3597 × 10^–18 J. The polarization part of the potential includes fluoride ion polarization only. We used a fluoride ion polarizability of 7.09 Bohr³ and applied a damping function:

g_{i j} (r_{i j}) = 1 - c_{i j} e^{- b_{i j} r_{i j}} \sum_{k = 0}^{4} \frac{{(b_{i j} r_{i j})}^{k}}{k!}

(3)

Ion pair	B_ij	α_ij	$C_{i j}^{6}$	$C_{i j}^{8}$	$b_{i j}^{6}$	$b_{i j}^{8}$
F⁻–F⁻	181.84	2.267	15.0	150.0	1.9	1.9
F⁻–Be²⁺	41.72	2.254	0.0	0.0	1.0	1.0
F⁻–Li⁺	20.42	2.052	0.0	0.0	1.0	1.0
Be²⁺–Be²⁺	106.16	3.944	0.0	0.0	1.0	1.0
Li⁺–Li⁺	195.91	4.252	0.0	0.0	1.0	1.0
Be²⁺–Li⁺	0.00	1.000	0.0	0.0	1.0	1.0

to the interaction between the anion dipoles and the cation charges, with b_FBe = 1.78 au⁻¹, b_FLi = 1.81 au⁻¹, c_FBe = 0.99 and c_FLi = 1.40.

The simulation box consists of a slab with square cross-section (of lateral dimension L and perpendicular dimension D), periodically replicated in all spatial directions. A vacuum region is created “above” and “below” this slab. The experimental vapor pressures of these systems are so low that very long simulation times would be needed to observe an evaporation event [5–7]. Thus, in practical terms, our simulated liquid–vapor interface is in reality a liquid–vacuum interface. The simulation supercell therefore has dimensions L × L × L_z, with L_z = D + L_vacuum. The long-range electrostatic interactions are dealt with by employing a three-dimensional (3D) Ewald summation method, with vacuum boundary conditions in the direction perpendicular to the interface (conducting boundary conductions are employed in the two lateral directions). The following terms then have to be added to the energy, forces and pressure tensor resulting from the usual Ewald sum with conducting boundary conditions [8]:

E = \frac{2 π}{V} M_{z}^{2}

(4)

F_{z, i} = - \frac{4 π q_{i}}{V} M_{z}

(5)

Π_{z z} = \frac{2 π}{V^{2}} M_{z} (M_{z} - 2 M_{z, d})

(6)

where V is the volume of the simulation cell, M_z is the total dipole moment in the perpendicular direction and M_z,_d is just its displacement component. We also performed an Ewald summation of dispersion interactions [9], which is necessary to avoid substantial truncation effects [10].

The method employed to enforce canonical (NVT) ensemble sampling is the Nosé–Hoover chain thermostat method [11,12]. The system consisted of 576 F⁻, 192 Be²⁺ and 192 Li⁺ ions, and we studied it at a temperature of 1060 K, which is the temperature for which experimental surface tensions of the mixtures of LiF with BeF₂ were obtained [13]. The lateral length L of the simulation box was obtained by performing a bulk simulation of the system in the NPT ensemble following the method described by Martyna et al. [14], with a pressure fixed at 0 GPa. We obtained a value of 34.56 Bohr for L, and L_z was set to 210 Bohr. We used a time step of 0.5 fs and the system was simulated for 1.35 ns. A snapshot of the simulation box is shown in Fig. 1.

Fig. 1
Example of a snapshot of a simulation box (green: F⁻ ions, white: Be²⁺ ions, and pink: Li⁺ ions).

3 Results

3.1 Structure of the interface

In order to study the structural properties of the interface, we have constructed the density profiles along the direction perpendicular to the interface. The simulation cell is divided into N_bin narrow bins, each of width δz = L_z/N_bin. We have then calculated profiles of the ionic number density:

ρ (z_{n}) = \frac{N_{bin}}{V} 〈 N (z_{n}) 〉

(7)

where 〈N(z_n)〉 is the average number of atoms in the histogram bin z_n. We used a value of N_bin = 200, i.e. a resolution of δz = 0.56 Å. The result for the three different ions is given in Fig. 2. From the inset of this figure, we can deduce that the interface is separated into two zones:

• at the edge of the bulk, where the Li⁺ concentration already begins to decrease, appears a Be-enriched zone in which the Be²⁺ concentration becomes greater than its bulk value;
• then, in the outer region of the interface, the Be²⁺ concentration drops abruptly and reaches values smaller than the Li⁺ concentration. In this zone, the overall density is low compared to the bulk.

Studying the coordination number profiles of the two cations might give an explanation of the particular aspect of the density profiles across the interface. To do so, we considered an anion to be in a cation's coordination shell if the distance between them was smaller than the position of the first minimum of the corresponding radial distribution function (these functions are given for this system in Ref. [1]), namely 2.1 Å for Be–F and 2.7 Å for Li–F. The profiles obtained are given in Fig. 3 (in which the F⁻ ion density profile also appears to help interpretation). This shows that Be²⁺ ions are always tetracoordinated, even through the interface, while the Li⁺ ions coordination number decreases quite continuously. The $Be F_{4}^{2 -}$ units then remain intact at the interface. In our previous study [1] we showed that the system consisted of a corner-linked network of $Be F_{4}^{2 -}$ units, so that two types of fluoride ions can be considered: bridging and terminal fluorides. Of course, at the interface with vacuum all the fluorides are terminal and so there is an excess of negative charges that can only be filled in by lithium ions. This explains the bigger concentration of lithium in the outer part of the interface, which is a consequence of maintaining local charge neutrality.

Fig. 3
Coordination number profiles. The dotted line corresponds to the F⁻ density profile (in arbitrary units).

Still, this is not sufficient to explain the segregation of the Be²⁺ ions in the inner part of the interface. Bresme et al. showed that size asymmetry could cause this type of effects [15]: in their study, the bigger anions have a higher concentration at the interface. Madden and Aguado also observed such effects in equimolar mixtures of KCl with NaCl and LiCl [16], in which the bigger potassium ions were more concentrated at the interface. In LiBeF₃ the segregation is much more important and cannot be explained by size-asymmetric effects only. At this concentration of LiF, an extended network of corner-linked tetrahedral $Be F_{4}^{2 -}$ units has already been formed, and we can only imagine that the higher concentration of the Be²⁺ at the interface can be associated with a stabilization of this network. Indeed, near the vapor phase the density decreases and so the chains present in that region acquire more degrees of freedom for their conformational motion. Such effects are observed in room-temperature ionic liquids: in their study of the surface structure of butylmethylimidazolium ionic liquids, Lynden-Bell and Del Pópolo showed that the surface is structured with a top monolayer containing oriented cations and anions [17], and that the butyl side chains of the cations tend to face the vacuum, while the methyl side chains tend to stay in the bulk of the liquid. An alternative explanation for the segregation is that an excess of polarizable particles at an interface lowers the surface energy [16]. This would imply a tendency of the fluoride ions to occupy the surface and therefore to drag the Be²⁺ along with them.

On passing from the surface into the bulk, one can observe that the Be²⁺ and Li⁺ profiles do not tend to a perfectly constant value and show a number of oscillations, with local regions of alternating Be²⁺ and Li⁺ excess. This type of oscillatory structure has been seen in diffraction studies of the surfaces of liquid metals; recent progress on this topic is reviewed in Ref. [18]. In their study on mixtures of KCl with NaCl, Madden and Aguado observed that potassium segregation towards the surface leads to a local Na⁺ excess in the bulk, and that this excess translates along the direction perpendicular to the interface, bouncing back and forth from the two ends of the slab. They showed that 5 ns of dynamics were necessary to average out these fluctuations, after which no oscillatory structure was observed. In our study, the Be²⁺ segregation at the interface is far more pronounced, so that the possibility of an oscillatory concentration profile seems more likely. To really confirm this possibility would require longer runs than we have so far attempted. We have chosen to study a wide slab but on a relatively short timescale, in order to get completely independent interfaces. The Be²⁺ excess at the interface appears stationary, so the interface itself is well sampled in our simulations, but the oscillatory nature of the Li⁺ profile must be confirmed in longer runs.

3.2 Electrostatic potential

A consequence of the arrangement of the ions close to the interface is the formation of charge ordering. The excess of F⁻ in the outer part leads to an excess of negative charges in that region, which is followed by an excess of positive charges in the Be²⁺ rich region. The charge distribution ρ_q(z) is shown in the upper part of Fig. 4. This leads to a local violation of the electroneutrality condition. The resulting electrostatic field, E, can be computed using Gauss's theorem:

E (z) = \frac{1}{ε_{0}} \int_{- ∞}^{z} ρ_{q} (z^{'}) ⅆ z^{'}

(8)

where ρ_q(z) is the charge density at the position z. Then the electrostatic potential difference across the interface ΔΦ_q(z) is given by:

Δ Φ_{q} (z) = Φ_{q} (z) - Φ_{q} (z_{0}) = - \int_{z_{0}}^{z} E (z^{'}) ⅆ z^{'}

(9)

where z₀ is a point in the vapor region. Bresme et al. already observed such charge density profiles in the case of size-asymmetric ionic liquids [15]. The electrostatic potential provides a quantitative measurement of the magnitude of the charge separation at the interface. These authors obtained electrostatic potentials of the order of 0.1 V. The case of ionic liquids is different: in these systems there is an excess of positive charges [17] due to the migration of the apolar butyl groups towards the vapor. Lynden-Bell and Del Pópolo then obtained negative electrostatic potentials of ≈–0.2 V (when the anion was

P F_{6}^{-}

) to ≈–0.7 V (when the anion was Cl^–). In our case, the charge excess is negative and much more pronounced than in those studies; consequently, the electrostatic potential reaches a very large value of 2.55 V.

Fig. 4
Charge density, induced dipoles' density and electrostatic potential profiles.

However, there is a second contributor to the electrostatic potential due to the induced dipole moment associated with each F⁻ ion. On an average, the induced dipoles of the surface F⁻ ions tend to be oriented so that their negative ends point towards the bulk, as can be seen in the middle part of Fig. 4. The induced dipole contribution to the surface potential can also be computed [19] (ρ_μ(z) is the z-component dipole moment density):

Δ Φ_{μ} (z) = Φ_{μ} (z) - Φ_{μ} (z_{0}) = \frac{1}{ε_{0}} \int_{z_{0}}^{z} ρ_{μ} (z^{'}) ⅆ z^{'}

(10)

The induced dipole contribution to the electrostatic potential is negative and opposes that due to the charge distribution, so that the total electrostatic potential is of ≈0.98 V (the last part of Fig. 4 shows the behaviour of both components of the electrostatic potential). This value is a lot higher than the one obtained in the case of size-asymmetric ionic liquids, which shows that charge asymmetry leads to important charge-separation effects at the interface.

3.3 Surface tension

The surface tension has been obtained through its mechanical definition:

γ = \frac{L_{z}}{2} (〈 Π_{z z} 〉 - \frac{1}{2} 〈 Π_{x x} + Π_{y y} 〉)

(11)

where the factor 1/2 comes from the fact that the total surface area is equal to 2L² in a slab geometry. We obtained a value of 190.7 mN m⁻¹, to be compared with the experimental value [13] of 170.6 mN m⁻¹. The errors for this experimental results are +30% and −10%, so our value seems very reasonable. Furthermore, as shown in Ref. [20] to obtain an estimate of the bulk surface tension we should really perform calculations for several system sizes L and extrapolate to L = ∞. This could easily reduce the estimated bulk surface tension by 10% from the value we have estimated.

An interesting aspect of the experimental results, which were obtained for compositions of the melt ranging from pure LiF to pure BeF₂, is that for a given temperature the surface tension has a minimum at about 40 mol% of BeF₂. The authors attributed it to the formation of complex fluoroberyllate ions. This is in total agreement with the increase in fluoroberyllate concentration at the interface noted in the earlier section. For concentrations between 33.3 and 100 mol%, the variation of the surface tension is small, which might be due to the fact that the fluoroberyllate concentration varies less at the interface than in the bulk. The study of the interface structure for other compositions of the melt would be interesting to address this point.

4 Conclusion

In this paper we have shown that charge asymmetry in LiBeF₃ had important effects on its interfacial properties. Fluoroberyllate species formed in this liquid are stabilized at the interface, leading to a local enhancement of the Be²⁺ concentration. We show that this segregation leads to important charge-separation effects at the interface. An electronic dipole moment is also created at the interface, which softens the electrostatic potential variation across the interface. Still, this electrostatic potential reaches higher values than in the case of size-asymmetric ionic liquids. The surface tension of the liquid has also been computed and the value is in good agreement with experimental measurements. The concentration of fluoroberyllate species at the interface seems to be responsible for the slow variation of the surface tension with composition in BeF₂-rich mixtures of FLiBe. Future works will explore the evolution of all that properties with charge asymmetry by studying systems containing highly charged species (MF₃ and MF₄).

Acknowledgements

The authors would like to acknowledge the support of the European Commission's Research infrastructures activity of the Structuring the European Research Area programme, contract number RII3-CT-2003-506079 (HPC-Europa). Financial support of PCR-RSF (Programme concerté de recherches – Réacteur à sels fondus) is gratefully acknowledged.