We report the results of constant-potential molecular dynamics simulations of the double layer interface between molten 2LiF–BeF2 (FLiBe) and 23LiF–6NaF–21KF (FLiNaK) fluoride mixtures and idealized solid electrodes. Employing methods similar to those used in studies of chloride double layers, we compute the structure and differential capacitance of molten fluoride electric double layers as a function of applied voltage. The role of molten salt structure is probed through comparisons between FLiBe and FLiNaK, which serve as models for strong and weak associate-forming salts, respectively. In FLiBe, screening involves changes in Be–F–Be angles and alignment of the oligomers parallel to the electrode, while in FLiNaK, the electric field is screened mainly by rearrangement of individual ions, predominantly the polarizable potassium cation.

Molten salts are effective high-temperature solvents and thermal transport media that find use in a variety of applications. Nuclear energy is one application of interest that utilizes high-temperature molten halide salts. The fluoride salt-cooled high-temperature reactor (FHR) is a nuclear reactor under development that uses 2LiF–BeF2 (FLiBe) as a primary coolant,1,2 the molten salt reactor (MSR) is a nuclear reactor that uses a fluoride (such as FLiBe + actinide fluorides) or chloride molten salt mixture as liquid fuel,3 and the affordable, robust, compact (ARC) fusion reactor uses FLiBe as a tritium breeding blanket.4 The thermophysical and thermochemical properties of molten salts enable reactor designers to produce heat at high temperature and near-atmospheric pressure. These characteristics help boost power conversion efficiency for electricity production, diversify the energy products produced by nuclear reactors, simplify implementation of high standards for reactor safety, ease the demonstration of safety for reactors in development, and improve economic competitiveness.5–8 

During the Molten Salt Reactor Experiment (MSRE), which ended in the 1970s, corrosion of metal alloys by molten salt was extensively studied and corrosion control strategies were proposed.9 As nuclear reactors today are being developed for commercialization, corrosion control continues to be a topic of relevance to the economic performance and design of the reactors. In order to understand and, hence, predict and control corrosion processes in molten salt,10–13 a detailed understanding of the electrochemical double layer that forms at the metal/molten salt interface is essential. The structure of the metal/salt interface informs the charge transfer model for corrosion, the reaction mechanisms that occur at the metal/salt interface, and the rates and mechanisms of transport of oxidants and corrosion products. By understanding the potential dependence of the double layer capacitance, we can begin to develop models for the molecular structure of the double layer and the relationship between melt chemistry, fluoroacidity, and the role of the double layer in corrosion mechanisms. This fundamental knowledge then enables the generation of corrosion models that can predict the sensitivity of corrosion performance to changes in salt chemistry, to the presence of ionizing irradiation, and to the occurrence of nuclear transmutation reactions and nuclear decay reactions at the metal/salt interface.

Experimental investigations of the double layer in molten salts are complicated by high temperatures, material compatibility, air sensitivity, toxicity, and limited availability of reference electrodes.14–17 Mean-field treatments of electric double layers common to other solvents typically assume dilute concentrations of electrolytes, a framework strictly inapplicable to molten salts.18 As a result, molecular simulation has emerged as a powerful tool for understanding molten salt double layers, beginning with the pioneering work of Heyes and Clarke19 on molten alkali chlorides. Subsequent studies have extended this approach to account for constant-potential boundary conditions20 and ionic polarizability,21 which allow for more faithful comparison with experimental results. For example, experimental differential capacitance measurements of the molten chloride interfaces with Mg have shown good agreement with simulation results of molten chloride interfaces with idealized electrodes.20,22

Recently, experimental studies of molten fluoride salt/electrode interfaces have reported capacitance values that are orders of magnitude larger than those reported in molten chloride interfaces.23–25 Molecular simulation of molten fluoride interfaces may help the understanding of these results. We are unaware of previously published computational studies of fluoride double layers similar to those undertaken for molten chlorides.21,26 Furthermore, we note that the existing studies of molten salt double layers tend to focus on weakly associated melts. Weakly associated ionic melts such as eutectic LiF–NaF–KF (FLiNaK) are not found to form long-lived associates between the cations and fluoride anions, with their structure being best described by weakly coordinated complexes.27 Strongly associated salts possess associates stable on picosecond timescales, such as FLiBe in which BeF42 has a cage correlation time near 20 ps.28,29 Furthermore, these complex ions can form extended structures of interconnected fluoride coordination environments such as BexFy dimers and trimers and longer polymers.28,30–33 It is unclear to what extent the strength of the associates and degree of polymerization will impact the structure of the double layer in molten fluorides.

In this work, we use the recently developed SEM-Drude polarizable ion model34 to conduct constant-potential molecular dynamics simulations of FLiNaK and FLiBe double layers. Sampling these trajectories allows us to directly compare the FLiNaK and FLiBe double layers in terms of ionic density profiles, differential capacitance, and in-plane structure as a function of distance normal to the interface. Importantly, the current results, employing model electrodes, yield differential capacitance values similar to those reported for other halide melts.35–37 We discuss possible experimental effects that could lead to the much larger values previously reported for fluoride salts.

Molecular dynamics simulations of molten 2LiF–BeF2 (FLiBe) and 23LiF–6NaF–21KF (FLiNaK)38 confined between immobile model electrodes were executed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package.39 The electrodes were constrained to a constant-potential difference via the method of Siepmann and Sprik,40 initially implemented into LAMMPS by Wang et al.41 and later extended by Tee and Searles.42 Electrode charges were treated as Gaussian distributions with width 1.979 Å−1 and were updated every ten time steps. In order to model ion–ion interactions within the salt, we employed the SEM-Drude model developed by Sharma et al.,34 which adapts the Drude oscillator model proposed by Lamoureux43 to include damping terms necessary for accurate representation of molten salt interactions. Details of the salt interaction potential, fitting procedures, and used parameters are given in the supplementary material. Long-range electrostatic interactions were evaluated using a standard Ewald summation with an accuracy of 10−6. Short-range interactions were evaluated with a cutoff radius of 15 Å. The ion-dipole damping potential was evaluated via tabulation. Short-range interactions between the salt ions and the electrode atoms were chosen such that the electrode interacts with the salt ions as if it were fluorine. This is analogous to the convention used for simulations of molten chloride interfaces,20,44 allowing for more direct comparison between the presently studied molten fluoride systems and previous published results for molten chlorides.

Bulk simulations (periodic in all directions) of FLiBe and FLiNaK in the NPT ensemble at 973 K and zero pressure were equilibrated using a three-chain Nose–Hoover barostat and thermostat until a steady volume was reached. Following bulk equilibration, three FCC (100) planes of model electrode atoms were placed 2 Å in the z direction away from each end of the simulation cell. The electrode has associated lattice constant of 4.093Å and 4.067 Å for the FLiNaK and FLiBe simulations, respectively, such that eight unit cells matched the length of the x and y dimensions. Periodic boundary conditions in the x and y directions were maintained but turned off in the z direction. The final dimensions of the FLiBe cell was 32.746 Å in the x and y dimensions with 64.45 Å separating the planes of electrode atoms closest to the melt. The final dimensions of the FLiNaK cell was 32.50 Å in the x and y dimensions with 63.965 Å separating the planes of electrode atoms closest to the melt. The FLiBe and FLiNaK simulation cells contained a total of 5166 and 3600 salt ions, respectively. Thus, the cell volumes are similar but the ion count is smaller for FLiNaK due to its lower ionic density. A schematic of the final setup for FLiBe is illustrated in Fig. 1.

FIG. 1.

(a) FLiBe simulation cell setup. Visualization generated by OVITO.45 (b) Schematic of electrochemical notation used throughout the text. Dashed vertical lines correspond to the position of electrode layers nearest the melt.

FIG. 1.

(a) FLiBe simulation cell setup. Visualization generated by OVITO.45 (b) Schematic of electrochemical notation used throughout the text. Dashed vertical lines correspond to the position of electrode layers nearest the melt.

Close modal

Production runs were executed in the NVT ensemble at 973 K using a three-chain Nose–Hoover thermostat with a time step of 0.5 fs. The interfacial simulations were allowed to equilibrate for 250 ps at an applied potential difference of 0.0 V. The potential was then increased incrementally, with 50 ps of equilibration time for each subsequent incremental value of the potential, except for FLiNaK at low target potentials, where 100 ps intervals were used in the equilibration. In the equilibration, the magnitude of the increments in voltage ranged between 0.2 and 0.7 V in magnitude, with smaller voltage increases at lower applied potentials. The snapshot after each 100 ps equilibration was used to seed production runs at each voltage, during which data were collected for 1.5 ns. Production runs were executed for electrode potential differences of 0.0, 0.2, 0.5, 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, and 3.0 V. For further clarification, plots of the electrode charge density as a function of time during the equilibration process for the 3.0 V production runs have been included in Sec. IV of the supplementary material.

Statistically independent runs for FLiBe and FLiNaK at 0.0 and 3.0 V were executed with the same procedure as the production runs, but with initial coordinates reflected about all three axes. No statistically significant differences between the results reported below for the initial production runs and those derived from these independent runs are observed, as demonstrated for the ion density profiles and surface charge density data included in Sec. III of the supplementary material. In order to determine whether statistically sufficient sampling was performed, production runs were divided into five sequential time bins and surface charge density as a function of interfacial potential difference was calculated for each bin. The spread in the data between each bin is low, indicating that statistically sufficient sampling was measured. Details on this analysis are included in Sec. II of the supplementary material. In addition, the equilibrium electrode atom charge time autocorrelation functions for an applied potential of 0.0 V were calculated and are included in Sec. IV of the supplementary material. These results demonstrate that although the timescale for decay of charge fluctuations was slightly longer for FLiBe relative to FLiNaK, for both systems, sampling was averaged over timescales multiple orders of magnitude larger than the timescale for electrode charges to decorrelate.

Ionic density profiles were constructed by averaging a histogram of ionic positions over 1.5 ns of dynamics, sampled every 100 fs. Figure 2 shows the density profile for each ionic species as well as the net charge density profile of FLiBe and FLiNaK at electrode potential differences of ΔV = 0.0, 1.0, and 3.0 V. Positive (negative) potentials were applied to the electrode with a positive (negative) z coordinate. The net charge density profiles were obtained via a charge weighted sum of all ions and dipoles. All net charge density profiles display the oscillating behavior reported by Heyes and Clarke.19 In FLiBe, these charge density oscillations exponentially decay with a screening length of ∼3.4 Å, this length exhibits little dependence on applied potential at both the negative and the positive electrode. In FLiNaK, the charge density oscillations decay with a screening length of ∼5.4 Å at the positive (and neutral) electrode and 2.0 Å at the negative electrode. While FLiNaK’s screening length varied significantly depending on electrode polarity, it displayed little dependence on the applied potential. This insensitivity of the screening length to applied potential is in agreement with studies on binary chloride interfaces.26,44

FIG. 2.

Ionic density and net charge profiles of FLiBe and FLiNaK at electrode potential differences of ΔV = 0.0, 1.0, and 3.0 V. Net charge obtained from a charge weighted sum of all ions and dipoles. Ionic density in units of number of ions per Å3; net charge in units of e3. Vertical dashed lines correspond to the plane of electrode atoms closest to the melt.

FIG. 2.

Ionic density and net charge profiles of FLiBe and FLiNaK at electrode potential differences of ΔV = 0.0, 1.0, and 3.0 V. Net charge obtained from a charge weighted sum of all ions and dipoles. Ionic density in units of number of ions per Å3; net charge in units of e3. Vertical dashed lines correspond to the plane of electrode atoms closest to the melt.

Close modal

While the screening length varied little with applied potential in both FLiBe and FLiNaK, the magnitude of the charge density oscillations increased with applied potential. The behavior of the ion density profiles can be attributed to properties of the individual ions and the salt systems as a whole. For instance, lithium (nonpolarizable) displays below-stoichiometric density at distances 3–7 Å from electrodes of either polarity in both salt systems. In the same region of the FLiNaK double layer, excess potassium (polarizable) is found. Nonpolarizable species are ill-suited to screen electric field in comparison to their polarizable counterparts; therefore, we find that the double layer is depleted in the less polarizable cation. It may then come as a surprise to find beryllium (nonpolarizable) enriched at 3–7 Å from electrodes of both polarizations in FLiBe. However, beryllium ions are tightly bound inside of tetrahedral BeF4 coordination complexes,28,30,46 and these complexes are polarizable by tetrahedral orientation, bond-angle changes, and polarization of individual fluoride ions in the potential model.

The net charge density (ρq) can be related to the electrostatic potential profile using Poisson’s equation,

(1)

The net charge density profiles in Fig. 2 were averaged over the x and y dimensions, so Poisson’s equation can be collapsed into one dimension,

(2)

The net charge density profile was numerically integrated using the one-sided Green’s function method presented by Wang et al.,47 

(3)

where L is the simulation cell box length in the z dimension and C1, C2 are constants solved for by applying the constant electrode potential boundary conditions

(4)

with ΔV corresponding to the potential difference between the electrodes.

We note that while the constant-potential method enforces the applied potential only between atoms of the respective electrodes, the potential becomes invariant after three layers of electrode atoms.42 Therefore, applying constant-electrode potential boundary conditions to the simulation cell edges and the electrode atoms should be equivalent. This procedure was applied to the net charge density profiles of FLiBe and FLiNaK at all simulated electrode potential differences in order to generate Fig. 3. All potential profiles reach a constant value far (>10Å) from the interface, indicating that the electric double layer screens electric fields at each simulated potential difference. However, it appears that the inner layer of electrode atoms in Fig. 3 are held at a potential difference slightly less than the applied value. This is an artifact of the finite-size bins used for averaging the charge density data to calculate the Poisson potential. The plotted value corresponds to the average potential within a bin of width 0.1 Å. The empty space within that bin that does not contain electrode atoms is not constrained to the applied potential difference. This effect is not typically noticeable in constant-potential simulations of ionic liquids,42,47 but it has been observed in constant-potential simulations of molten salts.48 A possible reason for this discrepancy is high temperatures allowing charged species to approach the electrode more closely, causing potential variations within the bin containing the inner electrode plane to be more significant.

FIG. 3.

Poisson potential profiles for FLiBe and FLiNaK at all simulated electrode potential differences.

FIG. 3.

Poisson potential profiles for FLiBe and FLiNaK at all simulated electrode potential differences.

Close modal

A finite potential difference between the electrode and the bulk salt in the middle of the simulation cell exists even at zero applied potential for FLiBe and (to a lesser extent) FLiNaK. Based on previously published simulations of double layers for chloride salts,21,26 this is caused by differences among the salt constituents in the ion–electrode short-range interactions, allowing for some charged species to approach the electrode more closely than others even in the absence of applied potential. On average, this results in a buildup of charge near a neutrally polarized electrode and, therefore, creates finite potential difference. Vatamanu et al.26 tested this explanation by simulating the double layer of a hypothetical Cl+–Cl system where cations and anions had identical size and identical short-range interactions with the electrode. While simulating this system at zero applied potential, they found no potential difference between the electrode and bulk. Modeling results by Painter et. al.37 also illustrate an effect of ion–cation differences on the potential at zero charge (PZC). In our simulated FLiBe and FLiNaK, the absolute value of the bulk potentials is determined by the chosen convention for treating electrode–ion short-range repulsions, and the larger variability among the ionic sizes among the ions in FLiBe than among the ions in FLiNaK results in a larger magnitude, nonzero bulk potential at zero applied voltage across the electrodes: −0.7 V for FLiBe and close to zero for FLiNaK (Figs. 4 and 3).

FIG. 4.

Surface charge density (top) and differential capacitance (bottom) of FLiBe and FLiNaK as a function of ΔΨ.

FIG. 4.

Surface charge density (top) and differential capacitance (bottom) of FLiBe and FLiNaK as a function of ΔΨ.

Close modal

Each interface between the electrode and the molten fluoride can be treated as a parallel capacitor with potential drop

(5)

where ΨE is the electrode potential (±ΔV/2) and ΨB is the average potential of the bulk molten fluoride (the center of the profiles given by Fig. 3). ΔΨ is shown schematically in Fig. 1. The recorded values of the electrode atom charges can be summed and averaged over time in order to find the mean surface charge density (σ) associated with the capacitor. The differential capacitance is given by

(6)

We solved for CD as a function of ΔΨ using a procedure similar to Vatamanu et al.26 involving the fitting of σ(ΔΨ) to a third-order polynomial. CD(ΔΨ) was found by taking the first derivative of this polynomial. The result of this procedure is shown in Fig. 4.

To determine the optimal order of a polynomial to minimize the residuals errors without overfitting, a k-fold cross-validation test was used, as discussed in Sec. II of the supplementary material. The differential capacitance that employs a third-order fit closely agrees with the results from employing a fourth-order fit (as in Vatamanu et al.) near the PZC where data are densely populated, but it disagrees by up to 15% further from the PZC where data are sparsely populated. Based on the cross-validation analysis, this difference is interpreted to be a result of the fourth-order model overfitting the data in the range of potentials with a denser sampling, and it is for this reason that we based the results here on a third-order polynomial fit. Although the approach was not pursued in the present work, we note that an alternative method for computing capacitance is the linear-response framework presented by Scalfi et al.,49 which would avoid the ambiguities of the polynomial fitting.

The calculated values of differential capacitance range from 5.8 to 7.8 µF/cm2 for FLiNaK and from 3.8 to 6.0 µF/cm2 for FLiBe, and they have a dependence on ΔΨ. The differential capacitance values are on the same order of magnitude as reported in previous simulation studies of chloride molten salts20,26,44 and also similar to though somewhat smaller than values experimentally measured for single-component molten halide melts.35–37 

FLiNaK displays a higher differential capacitance than FLiBe at all simulated potentials. This can be attributed to a greater ability to screen electric fields due to the higher polarizability of individual cations as well as weaker tendency to form strong associates in the internal structure, implying greater freedom to rearrange. It is unclear how much of the differential capacitance is governed by the structure of the double layer as opposed to the chemistry between the salts and the electrode. A study by Pounds et al.21 used density functional theory to fit the salt–electrode interactions of a molten LiCl interface with Al. Doing so resulted in a differential capacitance of 10.2–10.4 μF cm−2 at 1200 K. Vatamanu et al. modeled the same LiCl interface while treating salt–electrode interactions as if the electrode interacted as the system anion (the same convention used in this paper), and reported a capacitance of 4.1 μF cm−2 at 900 K.26 Although the disparate temperatures complicate directly comparing the values, these results indicate the need for future work predicting the differential capacitance of molten fluorides with a chemically accurate description of salt–electrode interaction.

The modeled differential capacitance minimum falls around ΔΨ = −0.6 V for FLiNAK and 0 V for FLiBe. In FLiNaK, the polarizable cation, K, contributes predominantly to shielding, and above ±0.6 V polarization, Li (nonpolarizable) also begins to contribute. In FLiBe, the double layer consists of a similarly sized population of oligomers as in the bulk, but the Be–F–Be angles are larger (i.e., the oligomers are more stretched out), and the oligomers preferentially align parallel to the electrode. The structure in the double layer in FLiBe is further discussed in Sec. III E with regard to changes to the in-plane radial distribution of ions and the changes that occur in the oligomer geometry in the FLiBe double layer.

The calculated capacitance can be compared with differential capacitances measured through electrochemical experiments. Differential double layer capacitance can be measured50 through Electrochemical Impedance Spectroscopy (EIS),51 Cyclic Voltammetry (CV),52,53 and Galvanostatic Charge/Discharge Methods (GCDMs).52 Experimental work aimed at estimating double layer capacitance needs to (i) use inert electrodes to prevent reactions at the electrode–electrolyte interface; (ii) verify by current–potential scans that no reactions are occurring in the window of voltage used for the measurement; (iii) use a thermodynamic reference electrode to allow for the quantification of the electrode potential; and (iv) measure the dependence of the capacitance on the voltage applied between the reference electrode and the working electrode. Examples of studies following these criteria are seen for Room-Temperature Ionic Liquids (RTILs)54 but to the knowledge of the authors are absent for molten fluorides.

While for molten fluoride salts, there are no CV nor GCDM studies attempting to measure differential capacitance, there are three corrosion studies23–25 that report EIS data that can be used to estimate a magnitude for differential capacitance. The differential capacitance inferred from these EIS studies on FLiNaK is summarized in Table I. Unfortunately, none of the studies report the potential dependence or the potential. The main purpose of these studies was to explore corrosion behavior in fluoride melts and not to measure the capacitance of the double layer. The experimental design of these studies is, therefore, not optimized to calculate the double layer capacitance, and it does not have the experimental features highlighted above. Most notably, all three studies use non-inert working electrodes (Cr25 and stainless steel16,24) that undergo corrosion reactions at the interface with the salt. As a result, the capacitance measured in the studies is affected by the reactions at the electrode surface and is not representative of the double layer capacitance considered in the current simulations; thus, a comparison to the modeling results presented here is not of relevance.

TABLE I.

Differential capacitance of molten fluoride salts from experimental studies. In all studies, the capacitance is extracted through fitting of EIS with an R-(CPE-R) fitting circuit.55,56 In all studies, data are collected at the open circuit voltage, which likely varies and is not quantified relative to a thermodynamic reference potential. None of the studies performed cyclic voltammetry to verify that reactions are not contributing to the measured capacitance.

Molten saltWorking electrodeTemperature (K)CD (μF/cm2)
FLiNaK23  Stainless Steel 316 973 0.2 × 106 
FLiNaK24  Stainless Steel 316 923 16.4a 
FLiNaK + 0.2 wt. % CrF325  Cr 873 205 × 106b 
FLiBe + 0.2 wt. % CrF325  Cr 873 0.8 × 106b 
Molten saltWorking electrodeTemperature (K)CD (μF/cm2)
FLiNaK23  Stainless Steel 316 973 0.2 × 106 
FLiNaK24  Stainless Steel 316 923 16.4a 
FLiNaK + 0.2 wt. % CrF325  Cr 873 205 × 106b 
FLiBe + 0.2 wt. % CrF325  Cr 873 0.8 × 106b 
a

Inconsistencies in the reported electrolyte resistance and electrode surface area preclude recalculation of the result.

b

CD not reported in the study25 and calculated here based on the reported fitting circuit parameters.

The in-plane ordering of ions was quantified using a radial distribution function (RDF) modified to measure density correlations in the directions parallel to the electrode, as proposed by Heyes and Clarke.19 While a standard RDF corresponds to relative density within a spherical shell in the three-dimensional interval r + dr, the tangential RDF corresponds to the relative density of particles within a cylindrical shell in the two-dimensional interval r + dr. The height of this cylindrical shell is given by the parameter Δz, which is the width of the interval of z coordinates under consideration,

(7)

Here, ρ corresponds to the local density of the radial ion and dnr is a function that counts the number of ions within a cylindrical shell of thickness dr and within the considered z coordinates. Equation (7) was calculated near and far from the interface every 100 fs, averaged over 1.5 ns, in order to generate Fig. 5. Our choice of Δz balanced two competing factors. A value of Δz comparable to the distance between peaks in the charge density profile (roughly 2.2 Å for FLiBe and 2.9 Å for FLiNaK) allows for sampling of an entire ionic layer. However, r in the tangential RDF measures distance in the x, y plane. Therefore, a Δz greater than the nearest neighbor distance will result in a nonzero tangential RDF at low r because it would count more than one ionic layer in z, making the tangential RDF difficult to interpret. We, therefore, choose a Δz = 1.5 Å in order to sample the majority of an ionic layer in both FLiBe and FLiNaK without having the value of Δz exceed the nearest-neighbor separation of Be or Li with F.

FIG. 5.

Tangential radial distribution functions of FLiBe and FLiNaK near and far from the interface at zero applied potential.

FIG. 5.

Tangential radial distribution functions of FLiBe and FLiNaK near and far from the interface at zero applied potential.

Close modal

The tangential RDF shape and peak heights display little dependence on applied potential or electrode polarization. The shape of the FLiBe tangential RDF far from the interface shows qualitative agreement with the bulk radial distribution function.28 Nearest neighbor peaks corresponding to weakly associated cation–ion pairs in both FLiBe and FLiNaK (i.e., Li–F, Na–F and K–F) increase near the interface, indicating tighter coordination within the double layer of the typically weak associates. In FLiNaK, second-nearest neighbor pairings do not display changes in peak positions.

In contrast, peak heights corresponding to second-nearest neighbors in FLiBe (such as Be–Li and Be–Be) increase in height, and the Be–Be tangential RDF develops a second peak at roughly 2.3 Å separation (Fig. 6). In the bulk, there is no Be–Be peak until 2.8 Å. These modifications to the second-nearest neighbor peaks demonstrate FLiBe rearranging its medium-range inter-oligomer order in the double layer (second-nearest neighbors and further). However, the intra-oligomer order does not appear to change, characterized by the Be–Be peak centered at roughly 5 Å. Thus, the screening electric field is mediated by rearrangement of the individual oligomers and modification of oligomer shape, which is further discussed in Sec. III E.

FIG. 6.

Beryllium–beryllium tangential radial distribution functions near and far from interface at zero applied potential.

FIG. 6.

Beryllium–beryllium tangential radial distribution functions near and far from interface at zero applied potential.

Close modal

In order to further assess the changes in FLiBe medium-range order within the double layer, we analyze the length and geometry of chains of BeF42 associates (i.e., “oligomers”). Oligomer geometry statistics were collected every 100 fs by constructing graphs of all beryllium and fluorine ions. Connections were added between each fluorine–beryllium pair within 2.0 Å of each other. No changes in the distribution of oligomer chain length near or far from the interface were observed; however, we report changes in the geometry of molecules larger than monomers, as shown in Fig. 7. The relative frequencies of Be–F–Be connection angles were measured in order to quantify the tendency of BeF4 chains to occupy “jagged” or “stretched” states. In both the bulk and double layer regions, Be–F–Be connections with angles below ∼25° were highly unfavorable, corresponding to the volume exclusion of ions within the BeF4 molecules. In the bulk region, the most likely connection angle was ∼125°, agreeing with ab initio calculations, but angles lower than 125° until ∼75° remained relatively favorable in disagreement with ab initio calculations.28 In the double layer Be–F–Be, connections become biased toward wider angles, indicating a tendency for oligomers to stretch out in comparison to their bulk counterparts.

FIG. 7.

Probability distributions of Be–F–Be angle connecting dimers and orientation of dimers with respect to electrode. All graphs sampled beryllium ions from a range of z coordinates 4 Å wide. Distributions labeled “Bulk” correspond to graphs centered 30 Å from the electrode and distributions labeled “Double Layer” correspond to graphs centered 4 Å from the electrode.

FIG. 7.

Probability distributions of Be–F–Be angle connecting dimers and orientation of dimers with respect to electrode. All graphs sampled beryllium ions from a range of z coordinates 4 Å wide. Distributions labeled “Bulk” correspond to graphs centered 30 Å from the electrode and distributions labeled “Double Layer” correspond to graphs centered 4 Å from the electrode.

Close modal

In addition, we measure the orientation of dimers relative to the electrode by measuring the angle between a vector drawn between two connected beryllium ions relative to the (x, y) plane. A dimer with ideal “bulk” behavior should be isotropic, with no bias toward any particular orientation between parallel and perpendicular. We find that this is true for the bulk region of our system up until ∼45°, where the distribution begins to decay, indicating a lasting effect of the interface on oligomer orientation beyond the lengthscale of our simulation size. In the double layer, parallel orientations become heavily biased, in agreement with the intuition of similarly charged ions layering parallel to the interface in the double layer. The changes to the Be–Be tangential RDF in the double layer can be qualitatively explained by these angular distributions. The increase in peak height can be attributed to the bias in parallel orientations of oligomers. The tangential RDFs are calculated only considering ions separated within a cutoff z separation, as defined by Eq. (7). Biasing the orientations of oligomers to be parallel to the electrode increases the probability that two beryllium ions will have a difference in z coordinate that does not exceed this cutoff. Because the Be–F first nearest neighbor peak is so tight, the absolute value of the separation between beryllium ions within an oligomer is controlled by the Be–F–Be connection angle. The segregation of the bulk Be–Be peak into two separate peaks can be attributed to the biasing of the Be–F–Be toward higher angles and resultant bimodal distribution shown in Fig. 7.

Employing constant-potential molecular dynamics approaches demonstrated previously for chloride salts, we present simulation results for the structure and capacitance of the double layers between fluoride molten salts and model electrodes. Specifically, we compare results for FLiBe and FLiNaK fluoride mixtures, chosen as model systems that are known to form strong and relatively weak associates in the molten state, respectively.28,57 The results show that the polarizability of the cations and the formation of associates may dictate the mechanism by which charge is screened in the electric double layer. While weakly associating melts are free to screen electric field by rearranging ionic positions, strongly associating melts that form oligomers may screen field by modification of the oligomer ordering and oligomer structure.

We calculated the double layer differential capacitance as a function of potential for FLiBe and FLiNaK. It was observed that FLiNAK (a weakly associated melt) has a higher differential double layer capacitance than FLiBe (a former of strong associates and oligomers): 5.8–7.8 µF/cm2 for FLiNaK and 3.8–6.0 µF/cm2 for FLiBe. The screening length in FLiNaK is on the order of 2.0 Å at negative polarity and 5.4 Å at positive polarity. The screening length in FLiBe is on the order of 3.4 Å regardless of electrode polarity.

In FLiNaK, the polarizable cation, K, contributes predominantly to shielding, and above ΔΨ = ±0.6, Li also begins to contribute. In FLiBe, the double layer consists of a similarly sized population of oligomers as in the bulk, but the Be–F–Be angles are larger (i.e., the oligomers are more stretched out), and the oligomers preferentially align parallel to the electrode.

This study provides initial results for fluoride molten salt double layers using an idealized electrode interaction. Comparison to experimental results shows slightly lower differential capacitance than 20 to 50 µF/cm2 that has been previously measured for chloride, bromide, and iodides melts.35 This is expected, given that double layer differential capacitance is dependent on the salt–electrode interaction, and an idealized interaction is treated here. Future modeling efforts should include the development of more chemically realistic models of the fluoride/electrode interactions and the use of different metallicity models extending the constant-potential framework such as a Thomas–Fermi semiclassical model58 or chemical-potential-equalization approach.59 

It would be highly valuable to perform future experimental measurements,50 by GCDM, CV, and EIS, of double layer differential and integral capacitance in molten fluoride salts that generate potential-dependent capacitance and that carefully verify the absence of chemical reactions at electrode surfaces that otherwise can mask the double layer capacitance. Unfortunately, three EIS studies on FLiNaK and FLiBe are the only measurements to date that we could find for fluoride salts, and because they do not meet these requirements, they cannot be used for the purpose of extracting the double layer capacitance. Finally, experimental measurements of double layer capacitance for different electrodes that vary in the nature of their reactivity with the salt would also be useful in informing (or validating) the development of the salt–electrode interaction models.

Further descriptions of the potential model, simulation details, and results and analysis can be found in the supplementary material.

This research was supported as part of FUTURE (Fundamental Understanding of Transport Under Reactor Extremes), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES). L.L. also acknowledges the support provided by the Rose Hills Foundation during the early parts of this research through the UC Berkeley Summer Undergraduate Research Fellowship. The authors also thank Shern Tee for helpful discussions in interpreting the behavior of the electrode.

The authors have no conflicts to disclose.

Luke Langford: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Nicholas Winner: Conceptualization (equal); Data curation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal). Andrea Hwang: Data curation (supporting); Methodology (supporting); Software (supporting); Writing – original draft (equal). Haley Williams: Writing – original draft (equal); Writing – review & editing (equal). Lorenzo Vergari: Conceptualization (equal); Data curation (equal); Writing – original draft (equal); Writing – review & editing (equal). Raluca Olga Scarlat: Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (equal). Mark Asta: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
C.
Andreades
,
A. T.
Cisneros
,
J. K.
Choi
,
A. Y. K.
Chong
,
M.
Fratoni
,
S.
Hong
,
L. R.
Huddar
,
K. D.
Huff
,
J.
Kendrick
,
D. L.
Krumwiede
 et al.,
Nucl. Technol.
195
,
223
(
2016
).
2.
R. O.
Scarlat
,
M. R.
Laufer
,
E. D.
Blandford
,
N.
Zweibaum
,
D. L.
Krumwiede
,
A. T.
Cisneros
,
C.
Andreades
,
C. W.
Forsberg
,
E.
Greenspan
,
L.-W.
Hu
, and
P. F.
Peterson
,
Prog. Nucl. Energy
77
,
406
(
2014
).
3.
J.
Serp
,
M.
Allibert
,
O.
Beneš
,
S.
Delpech
,
O.
Feynberg
,
V.
Ghetta
,
D.
Heuer
,
D.
Holcomb
,
V.
Ignatiev
,
J. L.
Kloosterman
,
L.
Luzzi
,
E.
Merle-Lucotte
,
J.
Uhlíř
,
R.
Yoshioka
, and
D.
Zhimin
,
Prog. Nucl. Energy
77
,
308
(
2014
).
4.
B. N.
Sorbom
,
J.
Ball
,
T. R.
Palmer
,
F. J.
Mangiarotti
,
J. M.
Sierchio
,
P.
Bonoli
,
C.
Kasten
,
D. A.
Sutherland
,
H. S.
Barnard
,
C. B.
Haakonsen
 et al.,
Fusion Eng. Des.
100
,
378
(
2015
).
5.
B.
Mignacca
and
G.
Locatelli
,
Prog. Nucl. Energy
129
,
103503
(
2020
).
6.
R.
Moir
, in
ICENES’2007, 13th International Conference on Emerging Nuclear Energy Systems
,
İstanbul, Turkiye
,
3–8 June 2007
[
Energy Convers. Manage.
49
,
1849
(
2008
)].
7.
C.
Andreades
,
R. O.
Scarlat
,
L.
Dempsey
, and
P.
Peterson
,
J. Eng. Gas Turbines Power
136
,
062001
(
2014
).
8.
C. W.
Forsberg
,
Nucl. Technol.
206
,
1659
(
2020
).
9.
P. N.
Haubenreich
and
J. R.
Engel
,
Nucl. Appl. Technol.
8
,
118
(
1970
).
10.
K.
Sridharan
and
T. R.
Allen
,
Molten Salts Chemistry
(
Elsevier
,
Oxford
,
2013
), Chap. 12 p.
241
.
11.
V.
Ignatiev
and
A.
Surenkov
,
Structural Materials for Generation IV Nuclear Reactors
(
Elsevier
,
2017
), Chap. 5 p.
153
.
12.
F.
Schmidt
,
P.
Hosemann
,
R. O.
Scarlat
,
D. K.
Schreiber
,
J. R.
Scully
, and
B. P.
Uberuaga
,
Annu. Rev. Mater. Res.
51
,
293
(
2021
).
13.
H.
Baes
, Jr.
and
W.
Tanghe
,
J. Nucl. Mater.
149
,
149
(
1974
).
14.
F.
Carotti
,
H.
Wu
, and
R. O.
Scarlat
,
J. Electrochem. Soc.
164
,
H854
(
2017
).
15.
J.
Zhang
,
C. W.
Forsberg
,
M. F.
Simpson
,
S.
Guo
,
S. T.
Lam
,
R. O.
Scarlat
,
F.
Carotti
,
K. J.
Chan
,
P. M.
Singh
,
W.
Doniger
,
K.
Sridharan
, and
J. R.
Keiser
,
Corros. Sci.
144
,
44
(
2018
).
16.
J.
Qiu
,
A.
Wu
,
Y.
Li
,
Y.
Xu
,
R.
Scarlat
, and
D. D.
Macdonald
,
Corros. Sci.
170
,
108677
(
2020
).
17.
F.
Carotti
,
E.
Liu
,
D. D.
Macdonald
, and
R. O.
Scarlat
,
Electrochim. Acta
367
,
137114
(
2021
).
18.
C.
Merlet
,
D. T.
Limmer
,
M.
Salanne
,
R.
van Roij
,
P. A.
Madden
,
D.
Chandler
, and
B.
Rotenberg
,
J. Phys. Chem. C
118
,
18291
(
2014
).
19.
D. M.
Heyes
and
J. H. R.
Clarke
,
J. Chem. Soc., Faraday Trans. 2
77
,
1089
(
1981
).
20.
S. K.
Reed
,
O. J.
Lanning
, and
P. A.
Madden
,
J. Phys. Chem.
126
,
084704
(
2007
).
21.
M.
Pounds
,
S.
Tazi
,
M.
Salanne
, and
P. A.
Madden
,
J. Phys.: Condens. Matter
21
,
424109
(
2009
).
22.
23.
J.
Qiu
,
D. D.
Macdonald
,
R.
Schoell
,
J.
Han
,
S.
Mastromarino
,
J. R.
Scully
,
D.
Kaoumi
, and
P.
Hosemann
,
Corros. Sci.
186
,
109457
(
2021
).
24.
H.
Ai
,
M.
Shen
,
H.
Sun
,
K. P.
Dolan
,
C.
Wang
,
M.
Ge
,
H.
Yin
,
X.
Li
,
H.
Peng
,
N.
Li
, and
L.
Xie
,
Corros. Sci.
150
,
175
(
2019
).
25.
Y.
Liu
,
Y.
Song
,
H.
Ai
,
M.
Shen
,
H.
Liu
,
S.
Zhao
,
Y.
Liu
,
Z.
Fei
,
X.
Fu
, and
J.
Cheng
,
Corros. Sci.
169
,
108636
(
2020
).
26.
J.
Vatamanu
,
O.
Borodin
, and
G. D.
Smith
,
Phys. Chem. Chem. Phys.
12
,
170
(
2010
).
27.
C.
Bessada
,
A.-L.
Rollet
,
A.
Rakhmatullin
,
I.
Nuta
,
P.
Florian
, and
D.
Massiot
,
C. R. Chim.
9
,
374
(
2006
), part of Special Issue: GERM 2005.
28.
N.
Winner
,
H.
Williams
,
R. O.
Scarlat
, and
M.
Asta
,
J. Mol. Liq.
335
,
116351
(
2021
).
29.
M.
Salanne
,
C.
Simon
,
P.
Turq
,
R. J.
Heaton
, and
P. A.
Madden
,
J. Phys. Chem. B
110
,
11461
11467
(
2006
).
30.
C. F.
Baes
,
J. Solid State Chem.
1
,
159
(
1970
).
31.
L. M.
Toth
,
J. B.
Bates
, and
G. E.
Boyd
,
J. Phys. Chem.
77
,
216
(
1973
).
32.
A. S.
Quist
,
J. B.
Bates
, and
G. E.
Boyd
,
J. Phys. Chem.
76
,
78
(
1972
).
33.
S.
Delpech
,
S.
Jaskierowicz
, and
D.
Rodrigues
,
Electrochim. Acta
144
,
383
(
2014
).
34.
S.
Sharma
,
M. S.
Emerson
,
F.
Wu
,
H.
Wang
,
E. J.
Maginn
, and
C. J.
Margulis
,
J. Phys. Chem. A
124
,
7832
(
2020
).
35.
A. D.
Graves
and
D.
Inman
,
J. Electroanal. Chem. Interfacial Electrochem.
25
,
357
372
(
1970
).
36.
E. V.
Kirillova
and
V. P.
Stepanov
,
Russ. Metall.
2016
,
691
697
.
37.
K. R.
Painter
,
P.
Ballone
,
M. P.
Tosi
,
P. J.
Grout
, and
N. H.
March
,
Surf. Sci.
133
,
89
100
(
1983
).
38.

Note, FLiNaK’s reported eutectic ratio is 23.25-5.75-21; so, we have used a composition slightly off-stoichiometry. We do not expect this small difference to meaningfully affect our results.

39.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
40.
J. I.
Siepmann
and
M.
Sprik
,
J. Chem. Phys.
102
,
511
(
1995
).
41.
Z.
Wang
,
Y.
Yang
,
D. L.
Olmsted
,
M.
Asta
, and
B. B.
Laird
,
J. Chem. Phys.
141
,
184102
(
2014
).
42.
S. R.
Tee
and
D. J.
Searles
,
J. Chem. Phys.
156
,
184101
(
2022
).
43.
G.
Lamoureux
and
B.
Roux
,
J. Chem. Phys.
119
,
3025
(
2003
).
44.
O. J.
Lanning
and
P. A.
Madden
,
J. Phys. Chem. B
108
,
11069
(
2004
).
45.
A.
Stukowski
,
Modell. Simul. Mater. Sci. Eng.
18
,
015012
(
2009
).
46.
R. J.
Heaton
,
R.
Brooks
,
P. A.
Madden
,
M.
Salanne
,
C.
Simon
, and
P.
Turq
,
J. Phys. Chem. B
110
,
11454
(
2006
).
47.
Z.
Wang
,
D. L.
Olmsted
,
M.
Asta
, and
B. B.
Laird
,
J. Phys.: Condens. Matter
28
,
464006
(
2016
).
48.
R.
Sundararaman
,
D.
Vigil-Fowler
, and
K.
Schwarz
,
Chem. Rev.
122
,
10651
10674
(
2022
).
49.
L.
Scalfi
,
D. T.
Limmer
,
A.
Coretti
,
S.
Bonella
,
P. A.
Madden
,
M.
Salanne
, and
B.
Rotenberg
,
Phys. Chem. Chem. Phys.
22
,
10480
10489
(
2020
).
50.
H.
Wang
and
L.
Pilon
,
Electrochim. Acta
76
,
529
531
(
2012
).
51.
D.
Macdonald
,
Transient Techniques in Electrochemistry
(
Springer Science & Business Media
,
2012
).
52.
M.
Arulepp
,
L.
Permann
,
J.
Leis
,
A.
Perkson
,
K.
Rumma
,
A.
Jänes
, and
E.
Lust
,
J. Power Sources
133
,
320
328
(
2004
).
53.
R. B.
Rakhi
,
D.
Cha
,
W.
Chen
, and
H. N.
Alshareef
,
J. Phys. Chem. C
115
,
14392
14399
(
2011
).
54.
M.
Figueiredo
,
C.
Gomes
,
R.
Costa
,
A.
Martins
,
C. M.
Pereira
, and
F.
Silva
,
Electrochim. Acta
54
,
2630
(
2009
).
55.
G. J.
Brug
,
A. L. G.
van den Eeden
,
M.
Sluyters-Rehbach
, and
J. H.
Sluyters
,
J. Electroanal. Chem. Interfacial Electrochem.
176
,
275
(
1984
).
56.
C.
Hsu
and
F.
Mansfeld
,
Corrosion
57
,
747
(
2001
).
57.
E. W.
Dewing
,
Metall. Trans. B
20
,
675
(
1989
).
58.
L.
Scalfi
,
T.
Dufils
,
K. G.
Reeves
,
B.
Rotenberg
, and
M.
Salanne
,
J. Chem. Phys.
153
,
174704
(
2020
).
59.
H.
Nakano
and
H.
Sato
,
J. Chem. Phys.
151
,
164123
(
2019
).

Supplementary Material