Structure of molten NaCl and the decay of the pair-correlations

Molten salts have enjoyed a resurgence of interest because of their potential applications as the thermal energy storage material in concentrated solar power plants,^1–4 as the electrolyte in the so-called liquid-metal batteries for grid-energy storage, where the electrodes are made from liquid metals,^5,6 and as the fuel and coolant system for next-generation nuclear-fission molten salt reactors.⁷ They also find applications as the flux in the refinement of molten aluminum⁸ and as the electrolyte in the processing of nuclear fuel^9–11 or in the recovery of metals such as the rare-earth elements from scrap waste.¹¹ 

In order to predict the structure-related properties of a material, it is necessary to have an accurate model for describing the interactions between the atoms. Here, molten salts have long been used as a test bed for the development of these models because they are relatively simple ionic systems.^12,13 Often, it is possible to measure the full set of partial structure factors via the method of neutron diffraction with isotope substitution,^14–17 which delivers maximal information on the pair-correlation functions for testing the theoretical predictions. Such work has led to the incorporation of the energy terms associated with the ion polarization, leading to the so-called polarizable ion model (PIM).^18–20 These polarization terms play an important role in accurately predicting the structure and properties of the different phases.^21,22 More recently, molten NaCl has been used as a test-bed system for the development of interatomic potentials via machine learning techniques.²³ 

The nature of the interatomic interactions also makes molten salts an attractive proposition for exploring the asymptotic behavior of the pair-correlation functions, and thereby testing the predictions of simple theories that should also apply to inhomogeneous fluids at large distances from an interface.^24,25 In the case of concentrated electrolytes and ionic liquids, such theories have been called into question: The decay of the pair-correlation functions found from experiment^26–28 is much slower than that expected from simple theories based, e.g., on the restricted primitive model (RPM) in which equally sized hard-sphere ions are embedded in a uniform dielectric medium. In concentrated electrolytes or ionic liquids, the discrepancy between the measured and predicted decay lengths is typically an order of magnitude or more.²⁹ 

We have therefore been motivated to re-examine the structure of the archetypal molten salt NaCl by combining state-of-the-art neutron and x-ray diffraction experiments with molecular dynamics simulations. The objective is to obtain accurate information at the partial pair-correlation function level to provide a comparison with previous neutron¹⁶ and x-ray^30,31 diffraction results and to test, inter alia, the long-ranged decay of the pair-correlation functions and how it relates to simple theory.^24,25

This paper is organized as follows: The diffraction theory is outlined in Sec. II for a molten salt of arbitrary composition in order to retain a general notation for the corrections that are required to calculate the asymptotic behavior of the pair-correlation functions for simulations on a system of finite size. The diffraction experiments are described in Sec. III. The model potentials are outlined in Sec. IV, where a PIM and two rigid-ion models (RIMs) are considered in order to probe the effect of different energy terms on the melt structure. The molecular dynamics simulations are described in Sec. V. The results are presented in Sec. VI using both the Faber–Ziman (FZ)³² and Bhatia–Thornton (BT)³³ representations of the pair-correlation functions. The FZ representation emphasizes the contributions to the structure from the distribution of ion–ion pairs. The BT representation emphasizes the contributions to the structure from the number density and concentration fluctuations, which relate more readily to the thermodynamic properties of the melt. The results are discussed in Sec. VII by reference to the thermodynamic properties of the melt, the dependence of the structure on the interatomic potentials, and the decay of the pair-correlation functions at large-r values, where r is a distance in real-space. Conclusions are drawn in Sec. IX.

The total structure factor F(k) measured in a neutron or x-ray diffraction experiment on a molten salt, such as MX or MX₂, where M and X denote a cation and anion, respectively, can be decomposed into its contributions from the pair-correlation functions describing the atomic species according to³⁴ 

where c_α denotes the atomic fraction of chemical species α (M or X), w_α(k) represents either the k-dependent x-ray form factor f_α(k) or the k-independent coherent neutron scattering length b_α of chemical species α, and $SαβFZ(k)$ is the FZ partial structure factor for chemical species α and β. In an x-ray diffraction experiment, Eq. (1) is usually re-written as

where $w(k)=cMwM(k)+cXwX(k)$⁠.

In experiment, the partial pair-distribution functions g_αβ(r) are obtained from the Fourier transform,

where n₀ = N/V is the ion number density, N is the total number of ions, and V is the volume. In simulation, the g_αβ(r) functions are calculated from direct knowledge of the atomic positions. The mean coordination number of atoms of type β, contained within a volume defined by two concentric spheres of radii r_min and r_max centered on an atom of type α, is given by

The total structure factor can also be decomposed into its contributions from the number density and concentration fluctuations³³ such that Eq. (1) is re-written as

where $SNNBT(k)$⁠, $SCCBT(k)$⁠, and $SNCBT(k)$ represent the BT number–number, concentration–concentration, and number–concentration partial structure factors, respectively. The FZ and BT partial structure factors are related by the expressions

and

The low-k limits of the BT partial structure factors are related to key thermodynamic properties according to the equations

where T is the absolute temperature, p is the pressure, k_B is the Boltzmann constant, χ_T is the isothermal compressibility, G is the Gibbs free energy, and the dilation factor $δ=n0vM−vX$⁠, where v_α is the partial molar volume per particle of chemical species α.

For an ionic system, the charge–charge partial structure factor $SZZ(k)=SCCBT(k)/cMcX$ satisfies the Stillinger–Lovett³⁵ conditions (i) S_ZZ(0) = 0, which expresses charge neutrality, and (ii) $SZZ(k)=ΛD2k2$ in the small-k region, where Λ_D is the Debye screening length given by

where ɛ ≡ 4πɛ_rɛ₀, ɛ_r is the dimensionless relative dielectric constant, ɛ₀ is the vacuum permittivity, and Z_αe is the charge on ion α.³⁶ The Debye screening length is a consequence of requiring perfect screening in the long wavelength limit as k → 0.^37,38 Hence, $SCCBT(0)=SNCBT(0)=0$ and $SNNBT(0)=n0kBTχT$⁠.

The number–number partial pair-distribution function is given by the Fourier transform

Equivalent expressions relate the concentration–concentration partial pair-distribution function $gCCBT(r)$ to $SCCBT(k)/cMcX−1$ and the number–concentration partial pair-distribution function $gNCBT(r)$ to $SNCBT(k)/cMcX$⁠.

The samples investigated by neutron diffraction were Na³⁵Cl (99% enrichment, Sigma-Aldrich), Na^natCl (99.999%, Alfa Aesar), Na^mixCl, and Na³⁷Cl (95% enrichment, Sigma-Aldrich), where nat denotes the natural isotopic abundance of chlorine and mix denotes a 50:50 mixture of ³⁵Cl and ³⁷Cl. The samples were loaded into cylindrical silica ampoules of inner diameter 5.0(2) mm and wall thickness 1.00(5) mm that had been cleaned using a 48 wt % solution of hydrofluoric acid, rinsed using water and then acetone, and baked dry under vacuum. The ampoules were sealed under a vacuum of 0.8 × 10⁻⁵ Torr. The coherent neutron scattering lengths, taking into account the isotopic enrichment, are b_Na = 3.63(2) fm, $b35Cl$ = 11.56(2) fm, $bnatCl$ = 9.5770(8) fm, $bmixCl$ = 7.36(2) fm, and b_37Cl = 3.51(6) fm.³⁹ 

The neutron diffraction experiment was performed using the instrument D4c at the Institut Laue-Langevin (ILL) with an incident wavelength of 0.6956(1) Å.⁴⁰ A vanadium furnace was used with a heating element of diameter 17 mm and wall thickness 0.1 mm, together with a heat shield of diameter 25 mm and wall thickness 0.04 mm. Diffraction patterns were measured at 1093(5) K for each sample in its container in the furnace, an empty container in the furnace, and the empty furnace. Diffraction patterns were also measured at room temperature for the empty furnace, a cylindrical vanadium rod of diameter 6.37(1) mm in the furnace for normalization purposes, and a rod of neutron absorbing cadmium (of dimensions comparable to the sample) in the furnace to account for the effect of the sample self-shielding on the background count-rate at small scattering angles.⁴¹ 

Each diffraction pattern was built up from the intensities measured for the different detector groups. These intensities were saved at regular intervals to check for the stability of the experimental setup.⁴² For the Na^mixCl and Na³⁷Cl samples, the diffraction patterns showed the emergence of Bragg peaks, the most prominent at k ≃ 1.5 Å⁻¹, which grew slowly with time as the sample began to react with its silica container.⁴³ This reaction resulted from contamination of the Na³⁷Cl sample and led to cracked containers when the melt was cooled to room temperature.⁴⁴ Water is a possible impurity,⁴⁵ but the characteristic shape associated with inelastic scattering from light hydrogen, to which neutron diffraction is particularly sensitive because of the large incoherent scattering cross section of light hydrogen,³⁴ was not observed in the diffraction patterns.

The data corrections followed a standard procedure.^46,47 The number density recommended by Janz,⁴⁸ i.e., n₀ = 0.031 86(1) Å⁻³ at 1093 K,⁴⁹ was used in the data analysis. The Bragg peaks in the structure factors measured for the Na^mixCl and Na³⁷Cl samples were eliminated by taking the difference function

with x = 0.9.⁴³ 

The x-ray diffraction experiment was performed using the beamline BL04B2 at SPring-8 with an incident energy of 61.65 keV corresponding to the (220) planes of a silicon monochromator.⁵⁰ The sample of Na^natCl (99.999%, Alfa Aesar) was loaded into a cylindrical silica capillary tube of inner diameter 1.00(1) mm and wall thickness 1.00(5) mm, which had been cleaned using a 48 wt % solution of hydrofluoric acid, rinsed using water and then acetone, and baked dry under vacuum. The tube was sealed under a vacuum of 1.5 × 10⁻⁵ Torr. A vacuum bell jar was used to avoid air scattering, and the diffracted x rays were collected using a Ge solid-state detector. Diffraction patterns were measured at 1091(1) K for the sample in its container and an empty container. The x-ray data analysis procedure followed the method described elsewhere⁵⁰ using ion form factors⁵¹ and the Compton scattering correction taken from Refs. 52 and 53.

The measured functions are related to the FZ partial structure factors according to the equation

where the matrix elements are given by $a11=(1−x)cNa2bNa2$⁠, $a12=2cNacClbNa(bmixCl−xb37Cl)$⁠, $a13=cCl2(bmixCl2−xb37Cl2)$⁠, $a21=a31=cNa2bNa2$⁠, $a22=2cNacClbNabnatCl$⁠, $a23=cCl2bnatCl2$⁠, a₃₂ = 2c_Nac_Clb_Nab_35Cl, $a33=cCl2b35Cl2$⁠, $a41=cNa2fNa(k)2/f(k)2$⁠, $a42=2cNacClfNa(k)fCl(k)/f(k)2$⁠, and $a43=cCl2fCl(k)2/f(k)2$⁠. Equation (13) can be re-written as

where F is the column vector for the total structure factors, S is the column vector for the partial structure factors, and W is the m × n weighting factor matrix (m rows and n columns) with m > n and rank equal to n. W contains k-dependent matrix elements that originate from the x-ray form factors. The notation of Eq. (14) can also be used to relate F to a column vector S representing the BT partial structure factors $SIJBT(k)$ (I, J = N, C).

A unique S can be found by using singular value decomposition (SVD) to solve the least squares problem in which the two-norm ‖r‖₂ of the residual vector r = F − WS is minimized.^54,55 The weighting factor matrix can be written as W = UΣV^T, where U is an m × n orthogonal matrix, Σ is an n × n diagonal matrix with main diagonal elements σ₁ ≥ σ₂ ≥⋯≥ σ_n > 0 known as the singular values, and V^T is the transpose of an n × n orthogonal matrix V. The least squares solution to Eq. (14) is then given by

where the vectors U_l and V_l represent the lth columns of the matrices U and V, respectively. If F is subjected to an uncertainty δF, then, the corresponding uncertainty δS on S is given by

Figure 1 shows the k-dependence of the singular values σ_l (l = 1, 2 or 3) for the FZ and BT representations of the partial structure factors. In each case, σ₃ is much smaller than the other terms, which indicates that S will be most sensitive to the uncertainty associated with the l = 3 term. In the FZ representation, successive components in the column vector $V3=0.877−0.4370.199T$ correspond to the NaNa, NaCl, and ClCl partial structure factors, respectively, where the values are typical and correspond to k = 0. Thus, an uncertainty δF will lead to errors in the NaNa and ClCl partial structure factors that have an opposite sign to the error in the NaCl partial structure factor. In addition, the NaNa and ClCl partial structure factors will be the most and least sensitive to δF, respectively. Similarly, in the Bhatia–Thornton representation, σ₃ is again much smaller than the other singular values. Successive components in the column vector $V3=−0.100−0.980−0.169T$ correspond to the NN, CC, and NC partial structure factors, respectively, where the values are typical and correspond to k = 0. Thus, an uncertainty δF will lead to an error that has the same sign for all the partial structure factors. The CC and NN partial structures factor will be the most and least sensitive to δF, respectively.

The two-norm condition (or Turing) number of the weighting factor matrix is given by κ₂ = σ₁/σ_n, where σ₁ and σ_n are the largest and smallest singular values of W, respectively. It follows that κ₂ ≥ 1, where a value near to unity indicates a well conditioned matrix. Sometimes, the two-norm condition number $κ2′$ for the normalized matrix W′ is also quoted.^54,55 Figure 1 shows the k-dependence of these condition numbers for the FZ and BT partial structure factors. The measured datasets are more favorable for extracting the $SαβFZ(k)$ than the $SIJBT(k)$ functions.

We note that the partial structure factors of this work were obtained from direct inversion of the measured k-space functions, without the application of any constraints, and therefore reflect the statistical error on those functions. It was checked that the extracted BT partial structure factors satisfy the inequality relations $SNNBT(k)≥0$⁠, $SCCBT(k)≥0$⁠, and $SNNBT(k)SCCBT(k)≥SNCBT(k)2$ within the experimental error, which originate from the need for the intensity measured in a diffraction experiment to be either zero or positive. Unlike previous work on molten NaCl,^15,16 the partial structure factors were not obtained from an iterative procedure in which smoothed functions were extracted from the measured datasets via the application of physical constraints. These smoothed functions do not reflect a substantial statistical error on the measured diffraction patterns.¹⁵ 

The potential energy between two ions labeled i and j with charges q_i and q_j, separated by the distance r_ij, can be expressed as

The short-ranged part was represented by a modified Born–Mayer⁵⁶ potential (see Ref. 12 and references therein) where, in atomic units,

B_ij and a_ij are parameters controlling the ionic radius and the rate of decay of the repulsive wall, respectively, $Cijn$ are dispersion parameters, and $f̃n(rij)$ are damping functions. The latter account for the loss of asymptotic behavior in the dispersion interactions at short range,⁵⁷ are of the form suggested by Tang and Toennies,⁵⁸ and are controlled by the parameters $b̃6ij$ and $b̃8ij$⁠. Traditionally, many-body polarization effects are incorporated using a shell model,⁵⁹ which is a mechanical representation of the dipolar charge displacement, composed of a charged shell connected by an harmonic spring to a charged core.^12,60 In the PIM, the charge–charge Coulomb potential is also augmented with a description of ion polarization that takes a many-body character, but at the simplest level of a dipolar model, the induced dipoles μ_i are incorporated as point moments.¹⁹ In atomic units, the Coulomb (electrostatic) interaction is then given by

where the charge q_i = +1 for Na⁺ or q_i = −1 for Cl⁻, T = 1/r_ij, and the interaction tensor $Tαβ⋯=∇α∇β⋯(rij)−1$⁠, where {α, β} = {x, y, z}.⁵⁷ The first three terms represent the charge–charge, charge–dipole and dipole–dipole interactions. The fourth term represents the polarization energy, where k_i = 1/2α_i and α_i is the dipole polarizability of species i. The functions $g̃ij(rij)$ control the short-ranged contributions to the induced dipoles (see, for example, Refs. 19 and 61) and take the form

similar to the form used to model the short-ranged damping of the dispersion interactions. The charge-ordered nature of an ionic liquid means that the short-range terms are controlled entirely by the cation–anion interactions. Hence, $c̃NaNa=c̃ClCl=0$⁠.

A dipolar PIM requires two sets of parameters: the dipole polarizabilities, α_i, and the short-range damping parameters, $b̃ij$ and $c̃ij$⁠. The former characterize the response of an ion’s electron density to the internal electric field that arises from the presence of both the formal ion charges and the induced dipole moments on the other ions. The latter control the effect of overlap of the nearest-neighbor electron density on the induced dipole moments.^62–64 The dipole polarizabilities are obtained from ab initio electronic structure calculations in which ions are studied in the condensed phase.^61,65 The short-range damping parameters are obtained from ab initio calculations by applying specific distortions to the nearest-neighbor shell of counter-ions.^62–64 The remaining parameters controlling the short-range interactions are obtained by “force-fitting” methods (see Ref. 66 for full details). System properties, here the ion forces and cell stresses, are fitted to those extracted from the density-functional-based electronic structure calculations. Dispersion terms are parameterized from the ion polarizability.

The parameters describing U_sr(r_ij) are listed in Table I. In addition, the damping parameters for the dispersion interactions were fixed at $b̃6ij=b̃8ij=1.7$ au; the polarizability of the Na⁺ and Cl⁻ ions was 0.89 and 20.0 au, respectively; and the short-range damping parameters [Eq. (20)] were $b̃NaCl=b̃ClNa=1.76$ au, $c̃ClNa=3.0$ au, and $c̃NaCl=0.697$ au.

Simulations were also performed using two RIMs. In a RIM, there are no induction effects, so the Coulombic interaction energy reduces to U_Coul(r_ij) = q_iTq_j. In the first model, termed RIM1, the ion polarization terms were simply omitted, so the model retains the asymmetry of the short-range interactions between the Na–Na and Cl–Cl ion pairs resulting from differences in both their respective sizes and dispersion interactions. In the second model, termed RIM2, the like-ion short-ranged interactions were also set to zero and all dispersion interactions were removed, so the model maps onto an effective RPM, which shows symmetry in the interactions between the like ions that is broken in the more complete model. In the RPM, the ions are treated as hard-spheres, whereas in the RIM2, the repulsive wall is steep but finite-ranged. This difference of detail should not affect the cation–cation or anion–anion interactions because Coulomb ordering ensures the absence of like-ion nearest-neighbor contacts. There will be an effect, however, on the short-ranged cation–anion interactions, but this is not expected to alter substantially the long-ranged asymptotic behavior of the associated pair-correlation function.

The PIM molecular dynamics simulation was performed in the NpT ensemble in order to assess the ability of the potential model to reproduce the equation of state for the liquid. The fluctuations in the pressure were around the atmospheric value. RIM simulations were performed in the NVT ensemble using the density obtained at each temperature from the NpT PIM simulations. The system contained 1728 NaCl molecules in a tetragonal cell with side lengths L_x = 2L_y = 2L_z. This cell geometry was chosen in order to maximize the accessible r-space distance range at an affordable computational cost. A constant temperature was maintained using a Nosé–Hoover thermostat,^67,68 and a constant pressure was maintained using a barostat.⁶⁹ Crystalline configurations were melted at a temperature of T = 5000 K, well above the estimated melting point, and the temperature was then reduced to the required value and equilibrated for 50 000 time steps (≈30 ps). Production runs were of 1 700 000 time steps (≈1 ns) for the RIM and 600 000 time steps (≈360 ps) for the PIM. The PIM simulations are approximately one order of magnitude slower than their RIM counterparts because of the need to evaluate many-body interactions.

In the calculation of the large-r behavior of the pair-correlation functions, it is important to consider the finite size of the simulation cell. In the large-r limit,^70,71

where N_α and N_β are the number of particles of chemical species α and β, respectively. Equation (21) provides the correction that must be applied to the calculated g_αβ(r) function, where $SαβFZ(0)$ is found from simulation by extrapolation at low-k.

The number–number partial pair-distribution function is given by

In the case of an MX system for which N_M = N_X and c_M = c_X = 1/2, the large-r limit is given by

Similarly, the large-r limit of the concentration–concentration partial pair-distribution function

is given by

and the large-r limit of the number–concentration partial pair-distribution function

is given by

Figure 2 shows the temperature dependence of the number density obtained from experiment^49,72–79 and the PIM simulations. The latter are in the range found by experiment and are larger than the values obtained from previous simulation work,⁸⁰ reflecting the effect of including the ion polarization.

Figure 3 shows the total structure factors ³⁵F(k) and $Fnat(k)$ and the difference function ΔF(k) measured by neutron diffraction. The x-ray total structure factor S_X(k) is also shown and is compared to the function measured in the work of Ohno and Furukawa.^30,81 Discrepancies can be attributed to the use, in the present work, of high-energy x rays to reduce the attenuation corrections and a diffractometer with a vacuum chamber to minimize the background scattering.

Figure 4 compares the measured FZ partial structure factors with those obtained from the neutron diffraction work of Biggin and Enderby¹⁶ and the PIM simulation. The PIM results are in good agreement with the functions measured in the present work, reproducing the principal peak position at k ≃ 1.74 Å⁻¹ and the form of the high-k oscillations.

Figure 5 compares the measured BT partial structure factors with those obtained from the PIM simulation. Also shown are the contributions to the simulated $SIJBT(k)$ functions from the weighted FZ partial structure factors. The insets highlight the low-k regions, plotted as a function of k² in order to fit a straight line and obtain both the gradient and k → 0 limit.

For $SNNBT(k)$⁠, the cancellation at k ∼ 1.74 Å⁻¹ of the principal peaks in the NaNa and ClCl partial structure factors with the principal trough in the NaCl partial structure factor leads to a broad first peak at k ∼ 2.5 Å⁻¹. The high-k oscillations are dominated by the fluctuations in $SNaClFZ(k)$ [Fig. 5(a)]. In contrast, there is a sharp principal peak in $SCCBT(k)$ at k ∼1.74 Å⁻¹ that has a significant contribution from all three of the FZ partial structure factors [Fig. 5(b)]. The high-k oscillations are again dominated by the contribution from $SNaClFZ(k)$⁠. For $SNCBT(k)$⁠, there is no contribution from the NaCl partial structure factor and the similarity between the NaNa and ClCl partial structure factors leads to a weak principal peak at k ∼1.74 Å⁻¹ and to a near-cancellation of the high-k oscillations [Fig. 5(c)].

The prominent principal peak in $SCCBT(k)=cMcXSZZ(k)$ corresponds to the Coulomb ordering of the ions.¹³ The lack of significant structure in $SNNBT(k)$ indicates only weak topological ordering, in stark contrast to molten glass-forming MX₂ systems, such as ZnCl₂.^55,82,83

Figure 6 compares the measured partial pair-distribution functions g_αβ(r) with those obtained from the neutron diffraction work of Biggin and Enderby¹⁶ and the PIM simulation. The measured functions were obtained by Fourier transforming the partial structure factors, which accounts for the low-r oscillations about zero at distances smaller than the distance of closest approach between two ions. In the present diffraction work, these features do not have a significant contribution from the truncation effects associated with terminating $SαβFZ(k)$ at the maximum accessible value k_max ≃ 16 Å⁻¹ because the measured functions do not show oscillations that extend beyond k_max. Low-r oscillations are absent from the simulated g_αβ(r) functions because they are generated from direct knowledge of the ion positions. The nearest-neighbor distances and coordination numbers are summarized in Table II. The simulated functions are in excellent agreement with those obtained from the present diffraction work in terms of the first peak position, coordination number, and longer-range oscillations although the first peak in g_NaNa(r) is more intense.

Figure 7 compares the measured BT partial pair-distribution functions with those obtained from the PIM simulation. Also shown are the contributions to $gIJBT(r)$ from the weighted g_αβ(r) functions obtained from the simulation. For $gNNBT(r)$⁠, the first peak at r ≃ 2.68 Å has a strong contribution from g_NaCl(r) and therefore originates from the packing of ions with unlike charges. The oscillations that extend to larger distances have a periodicity of the order half that of the underlying g_αβ(r) functions. For $gCCBT(r)$⁠, there is a sharp principal trough at r ≃2.68 Å followed by oscillations that extend to large r-values, which reflects the Coulomb ordering. In contrast, $gNCBT(r)$ is almost featureless, reflecting the similarity between g_NaNa(r) and g_ClCl(r) and the weak coupling between the charge and number density fluctuations.¹³ 

The Fourier transform relations between the BT partial structure factors and corresponding partial pair-distribution functions can be written as

The low-k behavior of the BT partial structure factors can be found by making a Taylor expansion of sin(kr) such that⁸⁵ 

The moments $MIJ(2m)$⁠, where m = 0, 1, 2, …, are obtained from the r_max → ∞ limit of the running moments defined by

In the case of the number density fluctuations,

Equivalent expressions define the $ρIJ(2m)(r)$ functions corresponding to $gCCBT(r)$ and $gNCBT(r)$⁠. Figure 8 shows the running moments for m = 0 obtained from the PIM simulations at T = 1100 K.

As discussed in Sec. II, $SCCBT(0)=SNCBT(0)=0$ for an ionic system. It follows that in the limit as r_max → ∞, the m = 0 moments are given by

These moments deliver therefore the k → 0 limiting behavior of the BT partial structure factors. The other moments deliver the low-k behavior. For example, in the case of the concentration fluctuations,

such that in the low-k region,

Figure 9 compares the χ_T value found from the diffraction results by fitting a straight line to $SNNBT(k)$ vs k² in the low-k region and extrapolating to the k → 0 limit [Eq. (31)] with the values found from sound velocity experiments.⁸⁴ Figure 9 also shows the results obtained from the PIM simulations, where χ_T was found from either (i) the low-k behavior of $SNNBT(k)$ in the same way as for the diffraction work or (ii) by calculating $MNN(0)$ from the data shown in Fig. 8 [see Eq. (36)]. The difference in the values calculated from $SNNBT(k)$ and $MNN(0)$ reflects the errors associated with obtaining both the low-k and high-r limits, respectively. The simulated results show good agreement with experiment over the measured range of temperatures up to ∼1273 K. The deviation of the PIM compressibility from the linearly extrapolated experimental values at high temperatures may therefore reflect genuine behavior.

Table III lists the values of $ΛD2$ extracted from the molecular dynamics results by fitting a straight line to $SCCBT(k)/cMcX$ vs k² in the low-k region [Eq. (40)]. The fitted functions are shown in the insets to Figs. 5(b) and 11(b). In comparison, it is often assumed that ɛ_r can be set to the high-frequency dielectric constant ɛ_∞ obtained from the refractive index n measured at optical wavelengths via the relationship ɛ_∞ = n². For molten NaCl at 1093 K, the measured value is between n = 1.413 (Na lamp, wavelength 5893 Å) and n = 1.424 (Hg lamp, wavelength 5461 Å).⁴⁵ It follows from Eq. (10) that $ΛD2$ = 0.0326–0.0331 Ų, close to the value obtained from the PIM simulation. In the case of the diffraction work, a straight line fit to $SCCBT(k)/cMcX$ vs k² in the low-k region [see the inset to Fig. 5(b)] gives $SCCBT(0)/cMcX$ = 0.016(16) or $SCCBT(0)$ = 0.004(4) with $ΛD2$ = 0.17(2) Ų. The latter is therefore larger than expected, which reflects the experimental error.

In the case of $SNCBT(k)$⁠, the k → 0 limit found from the diffraction work by fitting a straight line to $SNCBT(k)/cMcX$ vs k² in the low-k region was 0.144(12) [see the inset to Fig. 5(c)]. Hence, $SNCBT(0)$ = 0.036(3) instead of the theoretically expected value of zero. This discrepancy reflects the small magnitude of the NC partial structure factor and its sensitivity to the experimental error (Sec. III C).

Figure 10 compares the FZ partial structure factors from the PIM simulations with those obtained from the rigid-ion models described in Sec. IV. As expected, the $SNaNaFZ(k)$ and $SClClFZ(k)$ functions are identical for the RIM2, which is equivalent to an effective RPM. The effect of including ion polarization (RIM1 → PIM) is to reduce the height of the principal peak or trough at k ∼ 1.74 Å⁻¹, corresponding to a reduction in the magnitude of the oscillations in the partial pair-distribution functions. This effect is most pronounced for $SNaNaFZ(k)$ and corresponds to an increase in the mobility of the Na⁺ cations. At T = 1200 K, for example, the self-diffusion coefficients for the Na⁺ and Cl⁻ ions are D_Na = 8.8 × 10⁻⁵ cm² s⁻¹ and D_Cl = 7.2 × 10⁻⁵ cm² s⁻¹, respectively, for the PIM as compared to D_Na = 6.2 × 10⁻⁵ cm² s⁻¹ and D_Cl = 5.6 × 10⁻⁵ cm² s⁻¹, respectively, for the RIM1. The inclusion of polarization terms tends to lower the activation barrier associated with the motion of an ion by facilitating charge screening during each step in a diffusion pathway, e.g., as an ion “squeezes” through the nearest-neighbor coordination shell of its counter-ions.¹² 

Figure 11 shows the BT partial structure factors obtained from the different models at 1100 K and compares them to the functions obtained from the diffraction work. For $SNNBT(k)$⁠, the PIM results show the best agreement with experiment in terms of the main peak structure in the k-range of ∼2–3.5 Å⁻¹, the oscillations at high-k, and the low-k behavior as highlighted in the inset to Fig. 11(a). Both RIMs give smaller low-k limits than either the PIM or experiment, corresponding to systems that are less compressible. The isothermal compressibility values predicted by the PIM, RIM1, and RIM2 are shown in Fig. 9. For $SCCBT(k)$⁠, the principal peak is more intense for both of the RIMs, which accompanies a more intense principal peak or trough in the associated $SαβFZ(k)$ functions. For $SNCBT(k)$⁠, the near-cancellation of the contributions from the underlying FZ functions leads to more dramatic differences between the different models. For the RIM2, the symmetric nature of the like-ion interactions leads to near perfect cancellation, so $SNCBT(k)≃$ 0. For the RIM1, the balance of the principal peak intensities in the underlying FZ functions results in a principal peak that is inverted with respect to that extracted from experiment or the PIM.

According to simple theory for the RPM in which the ion interactions are described by pairwise additive potentials having short-ranged repulsive and Coulomb terms with no induction of dispersion effects, the asymptotic behavior of the total pair-correlation functions rh_αβ(r) is given by

In Eq. (41), h_αβ(r) = g_αβ(r) − 1, all three rh_αβ(r) functions share a common decay length $[a0αβ]−1$ and a common wavelength of oscillation $2π/a1αβ$⁠, the amplitudes are related by |A_MM‖A_XX| = |A_MX|², and the phases are related by θ_MM + θ_XX = 2θ_MX.²⁴ Similarly, the asymptotic behavior of the BT total pair-correlation functions $rhNNBT(r)=r[gNNBT(r)−1]$⁠, $rhCCBT(r)=rgCCBT(r)$⁠, and $rhNCBT(r)=rgNCBT(r)$ can be written as³⁶ 

In the RPM, the M⁺ and X⁻ ions have the same size, and g_MM(r) = g_XX(r) such that $gNCBT(r)=0$⁠, i.e., the concentration (or charge) and number density fluctuations are decoupled. In this case, there is no requirement for $rhNNBT(r)$ and $rhCCBT(r)$ to share the same asymptotic behavior.³⁶ If $rhCCBT(r)$ decays more slowly than $rhNNBT(r)$⁠, it is expected that the amplitudes of rh_MM(r) = rh_XX(r) and rh_MX will be equal but opposite in sign and that θ_MM = θ_XX = θ_MX, i.e., rh_MM(r) = rh_XX(r) = −rh_MX(r) at the longest range.²⁴ For molten NaCl, $gNCBT(r)≠0$⁠, but it has, nevertheless, a small magnitude. The large-r behavior of the total pair-correlation functions should therefore resemble that of the RPM.