The standard continuum approach to modelling the interfacial properties in solid electrolytes is commonly to seek solutions to the Poisson-Boltzmann equation. It is possible to calculate the defect concentration in boundary layers due to the equality of electrochemical potentials between the boundary layer and the bulk at spatial equilibrium. In crystals structures, defects occupy defined sites, therefore the chemical potential of a Fermi-Dirac distribution gives a more realistic expression for the defect concentrations.

$$ \mu^o_{i,x} + RT\ln \left( \frac{c_{i,x}}{1-c_{i,x}} \right) + z_i F \Phi_x = \mu^o_{i,\infty} + RT\ln \left(\frac{c_{i,\infty}}{1- c_{i,\infty}}\right) + z_i F \Phi_{\infty}, $$where $\mu^o$ is the standard chemical potential, $c_{\infty}$ is the concentration of defects in the bulk, $c_i,x$ is the concentration of defect $i$ at a given site $x$, $z$ is the valence and $F$ is the Faraday constant, $R$ is the gas constant and $T$ is the temperature in Kelvin. The equality of chemical potentials is rearranged to give an expression for the defect concentration at each site,

$$ c_i = \frac{c_\infty \exp\left(\frac{-z_i\Phi_x + \mu_i}{kT}\right)}{1+ c_{\infty} \left(\exp\left( \frac{-z_i\Phi_x + \mu_i}{kT} \right) -1 \right) }. $$In a system of point charges, the charge density is proportional to the defect concentrations, given by

$$ \rho = \sum_i c_i z_i F, $$where $\rho$ is the charge density, which is linked to the electrostatic potential via Poisson's equation,

$$ \nabla^2\Phi = \frac{-\rho}{\epsilon \epsilon_0}, $$where $\Phi$ is the local electrostatic potential, $\epsilon$ is the relative permittivity of a given material and $\epsilon_0$ is the vacuum permittivity.

For solid electrolytes, the Poisson-Boltzmann equation is commonly solved by setting the boundary condition to the space charge potential. However, the grain boundary can be treated explicitly by solving the Poisson-Boltzmann equation using segregation energies of defects $\Delta E_{i,x}$. The segregation energy is defined as the difference in energy between the energy of a defect at its site and the energy of a defect in the bulk,

$$ \Delta E_{i,x} = E_{i,x} - E_{i,\infty}. $$Including the correct site statistics and the segregation energies, the Poisson-Boltzmann equation becomes

$$ \nabla^2\Phi = - \frac{1}{\epsilon \epsilon_0} \sum_i c_{i,\infty} \frac{ x_{i,\infty} \exp \frac{-z_i F \Phi_x + \Delta E_{i,x} }{k_BT} }{1 + x_{i,\infty} \left( \exp \frac{-z_i F \Phi_x + \Delta E_{i,x} }{k_BT} -1 \right) } z_i F. $$The Poisson-Boltzmann equation can be solved numerically using a finite difference approximation for the second order partial differential equation. To solve the Poisson-Boltzmann equation numerically, the second order partial derivatives need to be replaced with second-order finite difference approximations using Taylor series expansion. The differential equations then become a system of linear algebraic equations, which can be solved using matrix algebra.

In this work, the sites described are the one-dimensional projection of defects in a crystal lattice onto a grid. In real crystals, lattice sites are not necessarily equally spaced and therefore the finite difference equations must be defined to take into account different $\Delta x$ spacings.

For an irregularly space grid, the first derivatives are given as,

$$ f'_{x,a} = \frac{f_{x0} - f_{x-1}}{\Delta x_1}, $$$$ f'_{x,b} = \frac{f_{x1} - f_{x0}}{\Delta x_2}, $$and the second derivative is given as

$$ f''_0 = \frac{ f'_{x,b} - f'_{x,a} }{ \frac{1}{2}(\Delta x_1 + \Delta x_2 )}. $$This can be expanded using the full expression for $f'_{x,a}$, and rearranged,

$$ f''_0 = \frac{2 ( \Delta x_2 ( f_{x+1} - f_{x0} ) - \Delta x_1(f_{x0} - f_{x-1} ) )}{\Delta x_1 \Delta x_2 (\Delta x_1 + \Delta x_2 ) }. $$Taking the second derivative finite difference at every point in $[x_0, x_2]$ then gives the following expression,

$$ f''_0 = \frac{2}{\Delta x_1 \Delta x_2 (\Delta x_1 +\Delta x_2 ) } . \left[\Delta x_2f_{x-1} - ( \Delta x_1 + \Delta x_2 ) f_0 + \Delta x_1 f_{x+1}\right]. $$To solve the Poisson-Boltzmann equation, the function at each position is the electrostatic potential and the second derivative is solved to calculate the charge density at each position. Substituting $\Phi$ and $\rho$ yields the equation

$$-\rho = \frac{2}{\Delta x_1 \Delta x_2 (\Delta x_1 +\Delta x_2 ) } . \left[ \Delta x_2\Phi_{x-1} - ( \Delta x_1 + \Delta x_2 ) \Phi_0 + \Delta x_1 \Phi_{x+1} \right]. $$To find the solution of the discretised Poisson-Boltzmann equation, the above expression must be satisfied at all points on the grid. Therefore, these equations become a system of linear algebraic equations solved using matrix inversion. Matrix inversion is known as the process of finding matrix $b$ which satisfies the equation of an invertible matrix $A$.

Taking prefactor ($p$) as:

$$ p = \frac{2}{\Delta x_1 \Delta x_2 (\Delta x_1 +\Delta x_2 ) } $$$A$ as:

$$ A = \Delta x_2\Phi_{x-1} - ( \Delta x_1 + \Delta x_2 ) \Phi_0 + \Delta x_1 \Phi_{x+1} $$and $b$ as:

$$ b = -\rho $$The Poisson-Boltzmann equation can be solved.

The invertible matrix $A$ becomes:

$A$ = $\begin{pmatrix} (\Delta x_i + \Delta x_{i+1}) & \Delta x_{i+1} & 0 & 0 & 0 \\ \Delta x_i & (\Delta x_i + \Delta x_{i+1}) & \Delta x_{i+1} & 0 & 0 \\ 0 & \Delta x_i & (\Delta x_i + \Delta x_{i+1}) & \Delta x_{i+1} & 0 \\ 0 & 0 & \Delta x_i & (\Delta x_i + \Delta x_{i+1}) & \Delta x_{i+1} \\ 0 & 0 & 0 & \Delta x_i & (\Delta x_i + \Delta x_{i+1}) \end{pmatrix}$

To account for the prefactor, a new matrix $L$ is constructed. $L$ is created by transposing matrix $A$, multiplying it by the prefactor and then transposing the matrix again.

$$ L = ( A^T . p )^T $$Vector $b$ is created from values of $\rho$

$b$ = $\begin{pmatrix} \rho_1 \\ \rho_2 \\ \rho_3 \\ \rho_4 \\ \rho_5 \\ \end{pmatrix}$

$ \Phi = L^{-1} b $ is solved using `numpy.linalg.solve`

which gives $\Phi_{\mathrm{new}}$. This process is iterated until a preset convergence limit between is met.

To control numerical stability in reaching the convergence limit, a damping factor is used to reduce the jump in values of $\Phi$ between iterations.

$$ \Phi = \alpha\Phi_{\mathrm{new}} + ( 1.0 - \alpha)\Phi_{\mathrm{old}} $$Once the Poisson-Boltzmann solver has converged, the potential, charge density and defect distributions have reached their equilibrium values. Using these equilibrium values the resistivity and activation energy can be calculated

Taking the resistivity ($\rho$) as $$ \rho = \frac{1}{\sigma} = \frac{1}{ c \mu z }.$$

Where $\sigma$ is the conductivity, $c$ is the concentration, $\mu$ is the mobility and $z$ is the valence. Treating each site as a resistor in series, the perpendicular grain boundary resistivity can be calculated by taking the ratio of the sum of the resistivity at all sites to the resistivity in the bulk,

$$ r_{\mathrm{GB}} = \frac{\rho_{i,x}}{\rho_{i, \infty }} = \frac{ c_{i, \infty} \mu_i x_i }{\sum_x c_{i,x} \mu_i z_i }. $$Assuming constant mobility for the mobile defect across the space charge regions, the grain boundary resistivity becomes

$$ r_{\mathrm{GB}} = \frac{c_{i,\infty}}{\sum_x c_{i,x} }. $$The resistivity ratio calculation in the code is only suitable when the Poisson-Boltzmann solver is used under the Mott-Schottky approximation as this method assumes that there is only one mobile defect. Once the resisitivity ratio has been calculated, it can be extended into an activation energy from the Arrhenius relationship,

$$ E_a = \frac{\mathrm{d ln}r_{\mathrm{GB}}}{\mathrm{d}\left( \frac{1}{k_BT} \right) }. $$Analysis of grain boundary effects on ionic conductivity typically considers the space charge potential ($\Phi_0$) as a characteristic parameter, however $\Phi_0$ cannot be directly obtained from experimental analyses. $r_{\mathrm{GB}}$ can be obtained directly from experimental impedance spectroscopy and is used to approximate $\Phi_0$ using the Mott-Schottky model.

$$ r_{\mathrm{GB}} = \frac{\rho_{i,x}}{\rho_{i,\infty}} = \frac{\mathrm{exp}(z\Phi_0 / k_BT)}{2z \Phi_0 / k_BT } $$The conventional Mott-Schottky analysis assumes that the grain boundary region is negligibly thin and that oxygen vacancies are fully depleted in the space charge region. This gives an analytical description of the space charge behaviour which can be completely characterised from $\Phi_0$.

By calculating the grain boundary resistivity of the system, the Mott-Schottky model can be used to calculate $\Phi_0$ allowing a comparison between the value calculated using the Poisson-Boltzmann solver and the value calculated using the Mott-Schottky model.

Using `sympy.solvers.solve`

the function required to solve the Mott-Schottky model for $\Phi_0$ was calculated.

```
x = Symbol('x')
y = Symbol('y')
print( solve (exp(x) / ( 2 * x ) - y, x ) )
[-LambertW(-1/(2*y))]
```

A function named `solve_MS_for_phi`

was created in the source code which takes the input `y`

, the grain boundary resistivity, and returns `-mpmath.lambertw(-1/(2*y),k=-1)`

. The output for this is then multiplied by the Boltzmann constant, and the temperature, and divided by the charge to give $\Phi_0$.

The function to calculate the Mott-Schottky approximation for $\Phi_0$ takes one argument, and that is the charge of the mobile defect.

`c_o.solve_MS_approx_for_phi( valence[0] )`

The Poisson-Boltzmann equation can be manipulated to give the Debye length. The Debye length is the "characteristic length scale" from the Debye-Huckel equation, which is the linearised Poisson-Boltzmann equation used for dilute systems.

$$ L_D = \left( \frac{\epsilon \epsilon_0 k_BT}{2 F^2 c_{i, \infty } } \right) ^ \frac{1}{2} $$Where $L_D$ is the Debye length, $\epsilon$ is the dielectric, $\epsilon_0$ is the vacuum permittivity, $k_B$ is the Boltzmann constant, $T$ is the temperature, $F$ is Faraday's constant and $c_{i,\infty}$ is the bulk defect density. $L_D$ is then used to calculate the width of the space charge region.

$$ \delta_{sc} = 2L_D \left(\frac{\Phi_0 z} {k_BT} \right)^ \frac{1}{2}$$The space charge potential can be calculated from the Poisson-Boltzmann equation or from the Mott-Schottky model to find their respective space charge widths.