Overscreening and Underscreening in Solid-Electrolyte Grain Boundary Space-Charge Layers

Polycrystalline solids can exhibit material properties that differ significantly from those of equivalent single-crystal samples, in part, because of a spontaneous redistribution of mobile point defects into so-called space-charge regions adjacent to grain boundaries. The general analytical form of these space-charge regions is known only in the dilute limit, where defect-defect correlations can be neglected. Using kinetic Monte Carlo simulations of a three-dimensional Coulomb lattice gas, we show that grain-boundary space-charge regions in non-dilute solid electrolytes exhibit overscreening -- damped oscillatory space-charge profiles -- and underscreening -- decay lengths that are longer than the corresponding Debye length and that increase with increasing defect-defect interaction strength. Overscreening and underscreening are known phenomena in concentrated liquid electrolytes, and the observation of functionally analogous behaviour in solid electrolyte space-charge regions suggests that the same underlying physics drives behaviour in both classes of systems. We therefore expect theoretical approaches developed to study non-dilute liquid electrolytes to be equally applicable to future studies of solid electrolytes.

Predicting the equilibrium distribution of charged ions within crystalline solids at, or near, structural discontinuities such as grain boundaries or heterointerfaces is a long-standing problem in solid-state physics [1][2][3][4][5][6]. Within these near-interface regions, the concentration of individual ionic species can deviate significantly from their formal bulk values, giving rise to so-called "spacecharge" regions. The spatial profiles of space-charge regions is particularly significant in solid electrolytes [7,8], such as those used in solid-oxide fuel cells and all-solidstate lithium ion batteries [5,6,9,10]. In these cases a local decrease in the concentration of the charge-carrying mobile ionic species within a space-charge region is expected to contribute to interfacial resistance and decreased device performance [7,11,12].
The classic treatment of space-charge formation in solid electrolytes considers the distribution of mobile ions in terms of charge-carrying point defects-typically interstitials or vacancies-that behave as an ideal gas interacting only through mean-field electrostatics [7]. For this simple model, the equilibrium defect distribution perpendicular to the grain-boundary interface can be found by solving the one-dimensional Poisson-Boltzmann equation. For low space-charge potentials [13], the approximate linearised Poisson-Boltzmann equation can be used, which has a general analytical solution; the resulting space-charge profiles decay exponentially towards the asymptotic bulk defect concentration with a characteristic decay length equal to the Debye length, λ D ; where ε 0 is the vacuum permittivity, ε r is the relative permittivity of the solid, k B is the Boltzmann constant, T is the temperature, z i are the charged species' valences, e is the elementary charge, and c i are the bulk charged defect concentrations. For high space-charge potentials the full non-linear Poisson-Boltzmann equation must be solved, either numerically or by making further approximations to obtain an approximate analytical result [14]. Space-charge profiles for high space-charge potentials again decay monotonically, exhibiting superexponential decay close to the grain-boundary, converging to exponential decay with decay length λ D at larger distances.
This ideal-gas plus mean-field-electrostatics model of space-charge behaviour is valid only in the dilute limit of low defect density or high relative permittivity, i.e. weak electrostatic coupling. To address this limitation, a number of more complex one-dimensional models have been proposed that aim to describe space-charge behaviour under non-dilute conditions [10,[14][15][16][17][18]. These models extend the simple model described above by including additional non-ideal local or non-local defect chemical potential terms. The equilibrium space-charge profile is then computed by minimising the global free energy as a function of the one-dimensional defect concentration. These extended models can give complex non-exponential or even non-monotonic space-charge profiles, although the exact functional form of the resulting profile depends on the choice of non-ideal chemical potential terms and their parameterisation.
As is the case for the classic space-charge model, these non-dilute space-charge models solve a one-dimensional problem defined in terms of the mobile-defect density perpendicular to the interfacial plane. Each point in this one-dimensional density can be considered a planar average from some corresponding implicit three-dimensional defect distribution; the full three-dimensional distribution of mobile defects within the solid-electrolyte host, however, is never explicitly calculated. As a consequence, defect-defect correlations, which become increasingly important with increasing defect concentration or decreasing relative permittivity, are not explicitly computed. Instead, the effect of these correlations is approximated through the choice of non-ideal defect chemical potential terms used in each model.
Here, we follow a different approach to modelling solidelectrolyte space-charge profiles. Instead of solving an effective one-dimensional problem, we consider an explicit three-dimensional Hamiltonian for point charges in a system with a single grain boundary. We then sample the configuration space of this model, using kinetic Monte Carlo, and construct space-charge profiles as time-averages over specific simulation trajectories. Because this approach treats defect-defect interactions explicitly in three-dimensional space, any defect-defect correlations emerge directly from the simulation trajectories, avoiding the need to describe these correlations through analytical free-energy terms.
We find that in non-dilute solid electrolytes (high defect concentrations) or for strong defect-defect interactions (low relative permittivities) grain-boundaryadjacent space-charge profiles have a qualitatively different functional form to the monotonic decay predicted by classic space-charge theory. We observe damped oscillatory space-charge profiles that can be well described by a simple analytical function. We also observe spacecharge decay lengths that are significantly larger than the corresponding Debye length, and that increase with increasing defect concentration or defect-defect interaction strength (decreasing relative permittivity), giving the opposite behaviour to that predicted by dilute-limit models. The deviation of space-charge decay length from the Debye length is shown to follow a universal scaling law as a function of electrostatic coupling strength. We note that analogous phenomena have been previously reported for concentrated liquid electrolytes, which suggests that the same underlying physics drives emergent behaviour in solid electrolytes as in their liquid counterparts. Finally, we comment on the implications of these results for a possible unified theoretical framework for understanding non-dilute electrolyte behaviour in both liquids and solids.
Methods: Our computational model consists of a 3Dperiodic simple-cubic lattice of m × m × m sites, populated by a fixed number, N , of mobile point defects with +1 charge (Fig. 1). To ensure our model has net zero charge, we include a uniform array of δ − point charges located at the cube-centre-interstitial positions, with δ − = -E gb −N/m 3 . The mobile defects and their counter-charges interact through point-charge electrostatics, which we scale by a relative permittivity ε r . This Coulomb latticegas model has historically been well studied as a bulk system [19], particularly in the context of the classical one-component plasma [20,21]. To model the segregation of mobile defects to a grain boundary, we assign one plane of sites an on-site occupation energy of −E gb ; this term represents the difference in standard chemical potential or "segregation energy" for defects in the grain boundary core relative to the bulk. All sites outside the grain-boundary core have on-site energies of zero.
Our kinetic Monte Carlo simulations were performed using the kmc fmm code, which implements the accurate and optimally scaling electrostatic solver described in Ref. 23 using the ppmd framework [24]. All simulations were performed using a 75 × 75 × 75 simple-cubic lattice of sites, with three-dimensional periodic boundary conditions. The nearest-neighbour mobile-site spacing is 2.5 A, which approximates the typical lattice spacing of a solid electrolyte. For all simulations, we set the onsite occupation energy for the grain-boundary plane as E gb = −0.05 eV. Simulations were run at 300 K. Each individual simulation was initialised with a random distribution of mobile charges and run for a total of 1,500,000 steps. The initial 500,000 steps were discarded to allow for equilibration, and space-charge profiles were computed as a time-average over the subsequent 1,000,000 steps. For each combination of relative permittivity, ε r , and bulk concentration, n ∞ , we performed between 30 and 100 simulations, to obtain statistically converged space-charge profiles. To characterise and quantify the simulated space-charge profiles we performed maximum likelihood sampling, parameter posterior sampling, and Bayesian model selection for competing functional forms using the uravu Python package [25]. The posterior distributions and Bayesian evidence for specific functional models were found by nested sampling [26], implemented in uravu using the dynesty package [27]. Further details are given in the Supporting Information. Results: We first consider simulated space-charge profiles obtained for n ∞ = 0.005, where n ∞ is the bulk per-site number density of mobile defects, and ε r = 100, 10 and 1 (Fig. 2), which spans the range of relative permittivities for typical solid electrolytes [10]. At higher values of ε r the electrostatic interaction between the mobile defects is more effectively screened, and each model will more closely approximate the non-interacting (dilute) ideal lattice-gas limit. Conversely, as ε r is decreased the electrostatic interactions between mobile defects become stronger and non-ideal defect-defect correlations are expected to become more significant.
For the highest permittivity considered here (ε r = 100) we obtain a monotonically decaying space-charge profile ( Fig. 2(a)) that is qualitatively consistent with the dilute-limit behaviour predicted by the classic Poisson-Boltzmann model. More precisely, the simulated spacecharge profile is well described by an exponentially decaying function where n(x) is the planar-average defect number density, n ∞ is the bulk defect number density, and α = 1/λ sys , with λ sys the observed characteristic screening length. At lower permittivities the simulated space-charge profiles are no longer monotonic ( Fig. 2(b)), and are therefore not even qualitatively well-described by Eqn. 2. We instead observe oscillatory behaviour that decays into the bulk. Further decreases in ε r cause both the magnitude of these oscillations and the distance over which they decay to increase (Fig. 2(c)). The appearance of oscillations at higher defect-defect interaction strengths (decreasing ε r ) mirrors the behaviour of non-dilute liquid electrolytes and ionic liquids, which exhibit similar oscillatory charge profiles at electrified interfaces, where this phenomenon is termed "overscreening" [28]. Oscillatory local order has also been observed in previous bulk Coulomb lattice-gas and one-component plasma simulations [20,21,29]. By analogy to the analytically-known behaviour of liquid systems [30,31], we propose a damped oscillatory ansatz for the space-charge profiles in this low-permittivity regime: As shown in Figs. 2(b) and 2(c), this oscillatory functional form gives a significantly better description of these space-charge profiles at low permittivities than the dilute-limit monotonic (purely exponential) model [32]. The simulated space-charge profiles presented in Fig. 2 also show increasing space-charge widths-i.e. increased screening lengths-with decreasing relative permittivities. This is the opposite behaviour to that predicted by the classic space-charge treatment, wherein the Debye length decreases with decreasing ε r (Eqn. 1). The observation of a screening length that is larger than the classical Debye length again mirrors the behaviour of concentrated liquid electrolytes, where this phenomenon is termed "underscreening" [30,31,33].
In ionic liquids and concentrated liquid electrolytes, underscreening has been found to follow scaling laws of the form where λ sys is the observed screening length of the system, d is the diameter of the charged species, and ν is a fixed scaling factor [34][35][36][37]. For our Coulomb lattice-gas model there is no direct analogue for d. We instead analyse our simulation data by plotting ln (λ sys /λ D ) versus ln(Γ), where Γ is a dimensionless parameter that describes the strength of electrostatic coupling in a given system. Γ = λ B /a, where λ B = e 2 /(4πε 0 ε r k B T ) is the Bjerrum length, and a = [3/(4πn ∞ )] 1/3 is a Wigner-Seitz defect-sphere radius [20,21,38]. Further discussion of the choice of an appropriate scaling parameter is provided in the SI. Fig. 3 shows ln (λ sys /λ D ) versus ln(Γ) plotted for simulations at four bulk concentrations: n ∞ = 0.00025, 0.0005, 0.001, and 0.005. For each simulation, we use Bayesian model selection to select which modelpurely exponential or damped oscillatory-is best supported by the simulated space-charge-profile data, and assign λ sys = 1/α. For ln(Γ) > 0.4, we find analogous scaling behaviour to that reported for liquid systems [30,36,37,39], i.e. λ sys /λ D ∝ Γ ν , with ν = 1.087 ± 0.005 (95% CI). For a given bulk concentration of mobile defects, higher values of Γ correspond to stronger defect-defect interactions, i.e. decreased ε r . In this regime, as defect-defect interactions increase in magnitude the observed screening length increasingly positively deviates from the Debye length. For ln(Γ) < 0.4 we observe a region where the observed screening length is smaller than the Debye length. Again, this behaviour has been reported for liquid systems in studies that focus on the application of different closure relationships to the Ornstein-Zernike equation [34,36,[39][40][41][42][43][44][45][46]. At the lowest values of Γ (≈ 0.25, obtained for n ∞ = 0.00025 and ε r = 100) we recover the low-potential dilute-limit result that the simulated screening length is equal to the Debye length; λ sys = λ D . We note that the same scaling behaviour is observed for all four simulated concentrations, suggesting a universal scaling law operates in solid electrolytes that contain a single mobile defect species.
Summary and Discussion: We have performed kinetic Monte Carlo simulations of populations of mobile point defects on a three-dimensional lattice, interacting via point-charge electrostatics, as a simple model for space-charge formation at solid-electrolyte grainboundaries, for a range of bulk defect concentrations and relative permittivities. For dilute, weakly coupled systems (low particle concentration and high ε r ) we recover the behaviour predicted by classic treatments of space-charge behaviour, i.e. the linearised Poisson-Boltzmann model-the space-charge profile decays exponentially to the asymptotic bulk value with a decay length equal to the Debye length λ D . In non-dilute, strongly coupled systems (high particle concentration or low ε r ) we observe damped oscillatory space-charge profiles-overscreening-and space-charge decay lengths that are larger than the corresponding Debye lengthunderscreening.
Overscreening and underscreening are known phenomena for non-dilute liquid systems of charged particles (concentrated liquid electrolytes and ionic liquids) [28,30,31,35,[47][48][49][50][51][52][53][54]. Our simulations predict not only that overscreening and underscreening also occur in nondilute solid electrolytes, but that these phenomena have the same functional behaviour in solids as in analogous liquid systems. Space-charge profiles in systems that exhibit overscreening are well described by a damped oscillatory function and the degree of underscreening, given by the ratio of the observed space-charge decay length to the classic Debye length, follows an apparently universal scaling law above some critical effective defect-defect interaction strength. We also predict an intermediate regime, where the simulated space-charge decay length is smaller than the classic Debye length. Again, this mirrors previous results obtained for the decay of radial distribution functions in non-dilute liquid systems [34,36,39,46]. Our results also suggest an approximate empirical threshold for when the dilute-limit classic space-charge treatment may be reliably used to quantitatively model space-charge behaviour in solid electrolytes of Γ = λ B /a 0.25.
The strong correspondence between our results for space-charge formation at a grain boundary in a model solid electrolyte and those previously reported for liquid electrolytes at an electrified interface suggests that the same underlying physics is responsible for these analogous emergent phenomena in both liquid and solid electrolytes. This suggests that many of the theoretical methods that have been developed previously by researchers seeking to better understand the behaviour of liquid electrolytes, or analogous models such as the one-component plasma [20,21], may be equally applicable to developing an improved understanding of the emergent behaviour of ensembles of mobile charged defects in solid electrolytes. Supporting Information Discussion of Fig. 2 when replotted on a semi-logarithmic scale; discussion of the choice of non-dimensional coupling parameter used in Fig. 3; further details of the Bayesian analysis of the simulated space-charge profiles and subsequent parameter estimation.
Data availability: A dataset containing the full set of time-averaged simulation data that support the findings presented here, and the analysis code used to generate Fig. 2 and Fig. 3 is available as Ref. 22 under the CC BY 4.0 and MIT licences.
In Fig. 2 of the main manuscript we present the onedimensional time-averaged mobile defect distributions for n ∞ = 0.005 and relative permittivities (a) ε r = 100, (b) ε r = 10, and (c) ε r = 1, plotted on a linear scale. For ε r = 10 and ε r = 1 (Figs. 2b and 2c) the damped oscillatory decay model gives qualitatively better agreement with the simulation data than the purely exponential decay model. For ε r = 100 (Fig. 2a) the two models give similar space-charge profiles, which both agree qualitatively with the simulation data. To better show the differences between these two models, we have replotted the data from Fig. 2 on a semi-logarithmic scale at all three relative permittivities (Fig. S1). The oscillatory behaviour of the simulated space-charge profiles is now more easily visible, even at ε r = 100. This preferential agreement with the oscillatory model, even at this high relative permittivity, is also supported by the calculated Bayes factor for this simulation of B = 196 (indicating strong support for the more complex damped oscillatory model). In Figs. S1b and Figs. S1c the replotted data highlights that the long-ranged decay of these oscillatory space-charge profiles is still exponential (following a straight line on this semi-log plot).

DISCUSSION OF THE CHOICE OF NON-DIMENSIONAL COUPLING PARAMETER USED IN FIG. 3 OF THE MAIN MANUSCRIPT.
In Fig. 3 of the main manuscript we present the variation in simulated screening length, normalised with respect to the corresponding Debye length, as a function of the parameter Γ = λ B /a, where λ B = e 2 /(4πε 0 ε r k B T ) is the Bjerrum length, and a = [3/(4πn ∞ )] 1/3 is a Wigner-Seitz defect-sphere radius [2,3]. Γ provides a non-dimensionalised measure of the strength of electrostatic coupling between mobile defects [4]. The choice of Γ for our analysis is inspired by its historical use in studies of one-component plasmas [3,4]; these include a number of studies of one-component Coulomb lattice-gas models, which correspond to the bulk analogue of the grain boundary model we investigate here. Γ also provides an intuitive sense of how changes in model parameters-defect concentration, n ∞ , relative permittivity, ε r , or temperature, T -are expected to change the relative strength of electrostatic coupling. The scaling of Γ with respect to n ∞ , ε r , and T is given by We therefore obtain the intuitive result that increasing defect concentration, or decreasing either the relative permittivity (increasing charge-charge screening) or the temperature will increase the electrostatic coupling strength. The form of Γ also highlights that we can expect to be in the weakly coupled regime if a >> λ B , i.e. the typical defect-defect distance is larger than the distance at which defect-defect interactions are comparable to k B T , and the effect of Coulombic interactions can be treated as a mean-field perturbation [4]. Conversely, we can expect to be in the strongly coupled regime if a << λ B , i.e. the typical defect-defect distance is shorter than the distance at which the electrostatic defect-defect interactions are comparable to k B T . We note that other choices of dimensionless scaling parameter with equivalent scaling behaviour would give the same phenomenological behaviour when plotted on a double-logarithmic scale. For example, one analogous scaling parameter that might be considered to be closer Log-log plot of scaled simulated decay length λsys/λD = 1/(αλD) versus a/λD for simulations performed at four different defect concentrations. The diagonal dashed line shows the maxiumum likelihood estimate for Eqn. 4 in the underscreening regime. Open and closed circles correspond to space-charge profiles that best support the pure exponential decay model or the oscillatory decay model respectively, as determined through Bayesian model selection (see below for details). Error bars show 95% confidence intervals for λsys/λD. Source: Raw simulation data and scripts to generate this figure are available under CC BY 4.0 / MIT licences as Ref. [1].
in spirit to d/λ D , where d is an ion diameter, used to describe underscreening in ionic liquids and liquid electrolytes, is a/λ D . This scaling parameter scales with respect to n ∞ , ε r , and T as a λ D ∝ n 1 6 (ε r T ) 1 2 .

(S.3)
A log-log plot of λ sys /λ D as a function of a/λ D (Fig. S2) therefore shows the same universal scaling behaviour as plotting against Γ. As expected from the scaling analysis above, the scaling parameter ν obtained when using a/λ D is related to the parameter obtained when using Γ by a factor of 2. The gradient of the linear region in Fig. S2, calculated using the same Bayesian analysis procedure used for Fig. 3, is 2.17 ± 0.01 (95% CI), twice that of the gradient calculated in the main text where Γ is used.

FURTHER DETAILS OF THE BAYESIAN ANALYSIS OF THE SIMULATED SPACE-CHARGE PROFILES AND SUBSEQUENT PARAMETER ESTIMATION
The statistical modelling of space charge profiles ( Fig. 2 and Fig. S1) and the scaling law ( Fig. 3 and Fig. S2) in the main paper was performed using Bayesian analysis. For the data presented in Fig. 2 and Fig. S1, we have performed maximum likelihood sampling to obtain posterior distributions for the model parameters (the use of a uniform prior probability means that the maximum likelihood parameters and the maximum a posteriori values are equivalent). We have then performed nested sampling to estimate the Bayesian evidence B for the more complex damped oscillatory model versus the simpler purely exponential model. To estimate the scaling parameters for the data presented in Fig. 3 and Fig. S2 we used maximum likelihood sampling.

Space charge profile modelling
The Bayesian analysis of the exponential and oscillatory models used uniform prior parameter distributions. These bounding values for these distributions are given in Table S1 for each bulk defect concentration. The prior distributions used for both models are identical for α, A, and n ∞ at a given concentration. Additionally, we have calculated Bayes factors for the comparative evidence supporting the two proposed space-charge models (damped oscillatory versus purely exponential), where the prior probability of both models was 0.5.
Because of the size of our simulations is limited by computational cost, the time-average space-charge profiles, n(x) , are not guaranteed to relax to the exact nominal bulk defect concentration, n ∞ . When fitting the exponential and oscillatory model functions to our simulation data, we therefore consider n ∞ as a free fitting parameter. The layers immediately adjacent to the grainboundary xy plane do not follow the same asymptotic behaviour as the space-charge profile. We therefore exclude three layers for the n ∞ = 0.00025, and n ∞ = 0.0005 distributions and two layers for the n ∞ = 0.001, and n ∞ = 0.005 distributions from our fit to get a more accurate description of the long-range behaviour.
The use of nested sampling estimates the Bayesian evidence, Z, for a given model with respect to the data. This evidence can then be used to calculate the Bayes factor, B 01 , between the null hypothesis model 0 and the more complex model 1, 2 ln B 01 = 2(ln Z x − ln Z y ). (S.4) Individual Bayes factors were interpreted using the classification given by Kass and Raftery [5], where 2 ln B 01 > 10 (defined as "very strong" support for the more complex model) was considered sufficient use of the damped oscillatory model for subsquently estimating α. The evidence values and resulting Bayes factors for each of the space charge profiles are provided in the supporting dataset as Ref. 1.

Scaling law modelling
Uniform prior parameter distributions were also used in the maximum likelihood sampling of the scaling law. The constant of proportionality was bounded by (−1, +1) and ν was bounded by (0, +5). S1. Bounds used in Bayesian analysis calculations for each of the four concentrations considered. Curves were fit to distributions of average number of charges in a given plane, not the average fractional site-occupation number n(x) .