Self-interaction and transport of solvated electrons in molten salts

Extremely diluted alkali-metal/alkali-halide solutions feature solvated electrons released by the excess metal atoms that tend to localize in bound states analogous to polarons in dielectric solids.^1–9 Solvated electrons are the simplest anions in nature. They often appear as reaction intermediates in diverse chemical processes such as, e.g., radiolysis, photolysis, and electrolysis of polar materials.¹⁰ Despite having been observed for more than two centuries, since the solutions of potassium in gaseous ammonia were examined by Davy,¹¹ their properties are far from being completely explained, with some important advancements in their full understanding having appeared relatively recently in the literature, aided by the increasing accuracy and affordability of electronic-structure and machine-learning methods.^12–15 

Solvated electrons in molten metal–metal halide solutions have been experimentally investigated, especially since the 1940s, when molten salts were employed in the context of nuclear technologies.^1,16–20 Far from the NonMetal-to-Metal (NM-M) transition, the dynamics of such electronic states is adiabatic;^21,22 therefore, the distribution of the excess electrons at each moment is entirely determined by the instantaneous ionic configuration, and the electronic motion is due to the ionic one.

The adiabatic variation of the potential energy surface determined by the nuclear dynamics is a natural playground for Thouless’ theory of charge pumping;^23,24 in particular, the theorem of charge quantization, together with a recently discovered gauge invariance of transport coefficients,^25–27 provides a theoretical foundation for describing the charge-transport properties of ionic conductors according to the topology of their electronic structure. Notably, topological arguments demonstrate that in non-stoichiometric systems, which feature dissolved electrons, nontrivial charge transport can occur, meaning that adiabatic transport of charge can take place even without a net ionic displacement.^22,28 This happens in non-stoichiometric molten salts such as metal/metal-halide solutions, where the electrical (ionic) conductivity can be recast as the sum of a part due to ions and one due to solvated electrons alone, the two contributions being uncorrelated from one another,^22,28 resulting in a much increased electrical conductivity even before the NM-M transition.² The whole machinery behind these concepts is rooted in the modern theory of polarization;^29,30 the latter provides also a means to rigorously characterize the electronically insulating state by exploiting its defining feature, i.e., the absence of dc conductivity, in terms of the localization of the electronic wavefunction.³¹ 

In practical calculations, electronic localization—and, vice versa, electronic diffusion—is determined by how self-interaction is accounted for in the employed theoretical framework. It is well known that standard Density Functional Theory (DFT) is affected by self-interaction errors due to the interaction of each electron with the total one-body electron density, including its own density.³² The spurious contribution is partially removed by the approximate eXchange and Correlation (XC) functional, but errors are still large, especially for local and semi-local XC functionals.^32,33 The excess electrons in metal–metal halide solutions are effectively few-electron systems and, as such, are particularly affected by spurious self-interactions.³³ Since the fluids are structurally disordered materials, the localization of a small number of electrons is often facilitated with respect to crystalline systems, so even calculations employing semi-local XC functionals result in rather localized electronic states whose dynamics has been studied from Ab Initio Molecular Dynamics (AIMD) simulations based on standard DFT.^3–6,21,22 Nonetheless, the inclusion of a fraction of EXact (Fock) eXchange (EXX) naturally leads to a stronger degree of electronic localization since the EXX partially removes the effect of self-interaction, and helps describing, e.g., polarons in solids,³⁴ and the localization in cavities of excess electrons in fluids.¹⁴ 

It is thus expected that the EXX would quantitatively alter the charge transport properties of an ionic system containing solvated electrons, as the latter tend to be more localized and their interaction with the ionic species becomes more pronounced. Focusing on the paradigmatic case of non-stoichiometric molten salts, we show how the electrical conductivity remains separately determined by a purely ionic contribution and one uniquely due to the excess electrons’ dynamics, with the cross contribution still being vanishingly small. At the same time, we demonstrate that the effect of EXX is observed primarily when examining average structural and electronic properties, rather than instantaneous quantities.

We study the non-stoichiometric molten salt Na_1+xCl_1−x, with x ≈ 0.06, and we compare its properties using two different functionals. The first, chosen as our reference, is the semi-local PBE functional³⁵ that has previously been employed in ab initio simulations of molten salts.²² The second, hybrid functional, aims to enhance the localization of excess electrons. However, determining the appropriate amount of EXX to include in the calculation is often challenging. Advanced techniques have been developed to determine the optimal EXX fraction and provide an accurate description of excess electrons, especially polarons, in materials. These techniques involve approaches such as enforcing piecewise linearity on the DFT energy with respect to electron occupation, addressing in this way the self-interaction errors of DFT.^36,37 Here, we have instead chosen to employ the widely used PBE0 formulation,³⁸ which is parameter-free and features 25% of EXX. We have found this fraction of EXX to be sufficient to induce a suitable level of localization of the excess electrons that allows us to qualitatively compare the differences in properties of non-stoichiometric molten salts under the influence of EXX.

We focus on the simple case of 33 Na atoms and 31 Cl atoms. Simulations employing both functionals are carried out in a cubic cell with a side of 13 Å at a density of 1.40 g/cm³ and a temperature of 1300 K. Further information is provided in  Appendix B. The presence of two extra Na atoms with respect to the stoichiometric formula leads to ionization and the release of one electron each, forming a solvated electronic pair called a bipolaron.^{3–5,21,22,28,39} The Highest Occupied Molecular Orbital (HOMO) corresponds to the bipolaron’s wavefunction.^22,28 To comprehend the bipolaron’s characteristics, we determine the location and spatial extent of the HOMO. This can be achieved by transforming the Kohn–Sham Bloch states into a localized basis, such as the Wannier Functions (WFs). We can thus determine Wannier Centers (WCs) and Spreads (WSs), where the WS quantifies the spatial extent of the orbitals. In particular, the HOMO WS represents the bipolaron’s position, while the HOMO WS provides insights into its degree of localization.^22,31,40 We employ the widely used Maximally Localized Wannier Functions (MLWFs),⁴⁰ where the HOMO WS is defined as the variance of the position operator evaluated over the HOMO WF. In the following, we will loosely speak of spread referring also to its square root, which, having the dimensions of a length, facilitates comparisons with distances.

A typical snapshot taken from a simulation of non-stoichiometric molten NaCl is shown in the left panel of Fig. 1, where several isosurfaces of the HOMO charge density are colored green. In the right panel, we show the distribution of the HOMO WSs computed during the two simulations. Here, it is worth mentioning that the MLWF construction may underestimate the WS when the latter is not small compared to the simulation box size,⁴¹ as it occurs when EXX is not considered. Therefore, the difference in the WS distribution between PBE and PBE0 may even be amplified, if the WS is estimated according to the accurate formulation of Ref. 41. Including a portion of EXX in the XC functional impacts both the structural and dynamical properties of non-stoichiometric NaCl. In the subsequent sections, we will address these aspects.

To understand the structural effects of adding the EXX to the XC functional, we compute the Radial Pair Distribution Function (RPDF), g(r), of the bipolaron (hereon labeled by b) and the atomic species in NaCl that are shown in Fig. 2. We observe several distinct features that reveal the differences between PBE and PBE0. First, the likelihood of locating a Cl ion in close proximity to the bipolaron is significantly lower in the case of PBE0 as compared to PBE. Second, RPDF for the interaction between the bipolaron and any ion exhibits a much more pronounced structure with PBE0. This indicates the existence of a well-defined shell structure surrounding, on average, the bipolaron, which is notably absent when using the PBE functional. In contrast, PBE predicts a comparatively uniform distribution of ions around the bipolaron, lacking any discernible shell-like organization. It must be noted that size effects can affect the results at large distances since the RPDFs have not completely decayed within half the simulation cell’s size. Features at shorter distances are nonetheless significantly different for PBE and PBE0.

This analysis suggests that the effect of EXX can be fully understood only by examining average quantities sampled throughout the entire dynamics, rather than focusing on instantaneous values. This is in accordance with the fact that the distribution of the values of the spread—a statistical quantity—is markedly different in the two simulations, but the allowed values—instantaneous quantities—are partially overlapping. As a result, by evaluating a single snapshot, it is almost impossible to discern whether the simulation was performed using a PBE or PBE0 functional. This is confirmed by a Principal Component Analysis (PCA) performed on Smooth Overlap of Atomic Positions (SOAPs) descriptors^42–46 associated with the bipolaron’s local environments, whose results are presented in  Appendix A.

The expectation values appearing here, denoted by angled brackets, are estimated by block-averaging over trajectory segments and, within each block, via an average over initial times. This is implemented in the software analisi.^55, Figure 4 shows the plot of the Mean Square Displaced Dipole (MSDD), $(6L3kBT)−1⟨Δμ(t)2⟩$⁠, as a function of time for various displaced dipoles related to the ions alone and to the bipolaron. The purely ionic contributions, whose slopes are the respective values of σ_ions, are comparable in both simulations. However, the bipolaron’s contribution, σ_b, in the PBE0 case is significantly smaller than that in the PBE simulation, where it played a leading role in determining the electrical conductivity value. In fact, the bipolaron diffusivity is reduced from 38 × 10⁻⁴ cm² s⁻¹ for PBE to 3 × 10⁻⁴ cm² s⁻¹ for PBE0, while the ionic diffusivity is of the order of 10⁻⁴ cm² s⁻¹ in both cases. Despite the reduced bipolaron contribution due to the inclusion of EXX, the total conductivity is consistent with the sum of the ionic and bipolaronic contributions, as already demonstrated for PBE in non-stoichiometric KCl in Ref. 22 and confirmed here, thus showing the lack of correlation between the two. The total ionic conductivity is significantly reduced after the inclusion of EXX, going from σ ≈ 18 S cm⁻¹ in the PBE case to σ ≈ 5 S cm⁻¹ in the PBE0 case. The latter compares fairly well with experimental conductivity measurements on Na–NaCl melts at similar metal concentration,² while the former appears rather overestimated.

In this work, we have explored the impact of self-interaction on the structural and transport properties of dissolved electrons in non-stoichiometric molten salts through MD simulations of a binary NaCl melt with excess Na. We have found that, when EXX is taken into account, the bipolaron exhibits a tendency to localize within well-defined solvation cells, whereas the RPDF appears structureless when using a semi-local functional. The localization of the bipolaron in the solvation cell is manifested as a statistical property, rather than an instantaneous one, as the RPDFs computed separately on each step of the trajectory are uncorrelated with the bipolaron’s spatial extent, both statically and dynamically. The distribution of the bipolaron’s spread, entailing its spatial extension, testifies the larger average degree of localization in the PBE0 simulation, in accordance with the reduced self-interaction induced by the presence of EXX.

The implications of these observations on the ionic transport properties of the melt are substantial. Notably, the inclusion of EXX significantly decreases the electrical conductivity. On a local level, charge transport due to ions alone correlates with the bipolaron’s motion. This feat notwithstanding, this correlation dissipates on larger scales, resulting in a total electrical conductivity that can be decomposed into a purely ionic contribution, rationalized through integer and constant atomic OSs, and a purely bipolaronic contribution, associated with the motion of the WCs related to the HOMO.

Our study sheds light on the intricate interplay between charge and mass transport in non-stoichiometric molten salts, highlighting the importance of accurately accounting for self-interaction in simulations to capture the underlying mechanisms. By confirming the nontrivial regime where charge and mass transport are effectively uncorrelated, this work represents a first step toward further investigations into the fundamental processes governing the transport properties of complex fluids.

We thank A. Grisafi, K. Rossi, and D. Tisi for fruitful discussions, and G. Fraux for technical help. This work was partially supported by the European Commission through the MaX Center of Excellence for supercomputing applications (Grant No. 101093374) and by the Italian MUR through the PRIN project FERMAT (Grant No. 2017KFY7XF) and the Italian National Center from HPC, Big Data, and Quantum Computing (Grant No. CN00000013). FG acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Action IF-EF-ST, Grant Agreement No. 101018557 (TRANQUIL).

The authors have no conflicts to disclose.

Paolo Pegolo: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Stefano Baroni: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Resources (lead); Supervision (equal); Writing – review & editing (equal). Federico Grasselli: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data, sample input files, and data-analysis scripts that support the plots and relevant results within this paper are available on the Materials Cloud platform.⁵⁸ See DOI: https://doi.org/10.24435/materialscloud:8f-d7.

The distinction between the two functionals is further emphasized by a principal component analysis (PCA) conducted to investigate the presence of clusters within the features characterizing the local atomic environments of the system.⁴⁶ To describe local environments we employed the widely used SOAP descriptors,^42–44 with parameters provided in Table I, as implemented in librascal.⁴⁵ Snapshots were sampled every 100 fs to ensure data uncorrelation. The SOAP descriptors for all the considered local environments and for both functionals were collected into a single feature matrix, [X]_e,f, where rows represent local environments and columns represent SOAP features. The matrix, X, was subsequently centered, and its covariance matrix, C = X^TX, was diagonalized with eigenvalues sorted in descending order. The PCA was then obtained by projecting X onto the first nPC eigenvectors of C, thereby accounting for as much variance in X as possible while simultaneously reducing its dimensionality.

Figure 6 shows the correlation between the first PC, which accounts for 17% of the variance, and the distance of the first peak of the instantaneous Na-b RPDF, r_peak. The Pearson correlation coefficients between these quantities for PBE and PBE0 are −0.43 and −0.69, respectively. Notably, the correlation for the PBE data gets to −0.55 when considering only structures where Na is the closest ion. This suggests that the first PC of the bipolaron’s SOAP descriptors is related to the ionic species and the proximity of its closest neighbor. The correlation is more robust for the PBE0 calculation, further corroborating the fact that the inclusion of EXX influences the first solvation shell of the bipolaron.

MD simulations are performed with the cp2k code,⁵⁹ version 9.1. The PBE and PBE0 functionals are parameterized according to the revised formulation, also known as revPBE.⁶⁰ The electronic density is expanded in the TZV2P-GTH Gaussian basis set. We employed Gödecker–Teter–Hutter (GTH) pseudopotentials⁶¹ encompassing electrons lying in the second shell, allowing for polarization effects that may contribute significantly to the accuracy of the simulations.⁶² Gaussian functions are mapped to a multi-grid with four levels; the plane-wave cutoff for the finest level is set to 400 Ry, with a relative cutoff of 60 Ry. Computations are sped up with the Auxiliary Density Matrix Method (ADMM).⁶³ The hybrid calculations employ a Coulomb operator truncated⁶⁴ at 6 Å to further reduce the computational effort while retaining its accuracy. Long-range Van der Waals interactions are modeled through Grimme’s D3 corrections.⁶⁵ All calculations are spin-polarized in the singlet state, which is energetically favored.^22,39

Both MD simulations sample the canonical ensemble at 1300 K by using a Bussi–Donadio–Parrinello thermostat⁶⁶ with a time constant of 1 ps, i.e., two thousand times the integration time-step of 0.5 fs. The Kohn–Sham wavefunctions are transformed to the basis of MLWFs⁴⁰ through Jacobi rotations every 1 fs to collect WCs and spreads to be used to compute the charge flux of Eq. (5). The PBE simulation has been thermalized for 5 ps. The production run is 50 ps long. The PBE0 simulation has been initialized from a snapshot drawn from the equilibrated PBE simulation and then further thermalized for 2 ps. The production run is 15 ps long. Input files can be found in the Materials Cloud.⁵⁸