MAR_CCT: Marburg program for modelling charge carrier transport

Abstract

A software packet for simulating charge carrier transport in solid electrolytes on the basis of Nernst–Planck Poisson equations is presented. The software is capable of handling a variety of electrochemical ion exchange processes ranging from charge attachment induced transport (CAIT) to alkali proton substitution (APS) and thermal electro-poling with a focus on concentration depth profiles. The software package includes a graphical user interface allowing for simple user input, visualization and output.

---

Introduction

Ion transport in solid matter plays an important role in a wide range of fields in science and technology, including chemistry, physics and biology. From another point of view, this includes aspects of material science and life science. Among the challenges in material science is the advancement of materials, devices and concepts for energy conversion and energy storage on our way towards improved sustainability. Here, the process of transport of protons and oxygen anions is pivotal for fuel cells^1,2 membrane technology³ and sensorics.⁴ Evidently, the transport of Lithium ions is pivotal for the operation of lithium ion batteries (LIBs).^5,6 On the other hand, electrode materials in LIBs also have to be electron conducting, defining the field of charge carrier transport as the umbrella term. Among the challenges in life science is a better understanding of elementary steps in neurotransmission involving aspects of chemistry and biophysics. These can be traced back to the transport of alkali and earth alkali ions through ion channels which is pivotal to sensory perception.^7–9 The field of ion transport is in principle mature.^10–12 However, this does not mean, there is no room for improvement, in particular for knowledge-based advancement.

In order to understand ion transport processes in detail, theoretical and experimental concepts are required that help to understand the full range of processes, from microscopic processes at the atomic level to macroscopically measurable particle flux and concentration profiles. Microscopic transport processes can be described, for example, with the help of programs based on Monte–Carlo simulations.^13–16 In the limit of macroscopic concentration-depth profiles often the coupled set of Nernst–Planck and Poisson (NPP) equations are used. Based on the concepts of Walter Nernst¹⁷ and Max Planck¹⁸ this set of equations provides an appropriate framework for the analysis of transport properties in solids and liquids.^19–22 The coupled equations describe the transport of charge carriers due to a gradient in the electrochemical potential. The electrochemical potential gradient includes possible contributions from concentration gradients and electric potential gradients. Driven by these gradients, the ion translocation takes place which in turn self-consistently affects the electric field in the sample. In this context, NPP calculations in liquids and solids have been performed including calculations for biomolecular applications,^23,24 in neuroscience,^25,26 for diffusion through porous media^27–29 and membranes,^30–33 at electrodes^34–36 and inside battery cells,^37–39 as well as calculations in solid electrolytes.^31,40–42 Aspects like diffusion along grain boundaries,^43–45 across grain boundaries,⁴⁶ finite size effects^47,48 and crowding^49,50 have been included into the NPP transport theory.

If an electrochemical potential gradient is applied across an ion conducting solid electrolyte for an extended period of time, a directed translocation of the ions occurs, such that concentration-depth profiles can arise.^51–56 These profiles then contain information regarding the underlying transport processes and, in addition to simple information concerning the basic DC conductivity and activation energies,^57,58 it is even possible to infer the energy landscape of the conductive ions in the sample.^41,59–62

The energy landscape of native ions in a sample strongly depends on the microscopic structure of the sample.^63–65 While mobile charge carriers in a crystalline structure generally all see the same atomic environment and therefore their properties, i.e. their activation energies, are almost identical, the situation in amorphous materials is significantly different. In these materials, the ions have specific environments, which are also expressed in specific activation energies and site energies. The result is a site energy distribution (SED), the populated part of which (PSED) determines the transport properties. If, under these conditions, a directed potential and concentration gradient is applied for an extended period of time, the ions need to overcome an effective threshold energy for long range transport and their effective activation energy mainly depends on the energy difference between the threshold energy value for long-range transport and the original site energy. If the directed transport continues for an extended time, the PSED can be probed systematically. Ions with lower effective activation energy will move first, before ions with higher effective activation energy contribute to DC transport. The PSED will then be depopulated systematically top-down in the energy domain such that the effective activation energy for the remaining ions becomes larger.^41,62 As a consequence, the effective diffusion coefficient becomes concentration dependent.

To date, several program packets capable of calculating different kind of transport properties have been published. Tu et al. developed for example a parallel finite element simulator for ion transport through ion channels based on NPP simulations.⁶⁶ It is capable to deliver electric potential, ion currents and concentration profiles inside the ion channel. However, the quantities are mostly evaluated in the electrostatic limit. IonMonger2.0 created by Clarke et al. is capable to calculate current, voltage and impedance response of perovskite solar cells.⁶⁷ COMSOL⁶⁸ is an extended but commercial program packet for many physical applications. It also contains an ion exchange module employing NPP calculations to derive concentration profiles and uses Fick's second law to derive ion fluxes.

The ion-exchange processes mentioned above represent, indeed, one important sub-class of charge carrier transport processes. Thermal ion exchange is e.g. the basis for glass strengthening.⁶⁹ Electric field assisted ion exchange can e.g. be considered as a basis for proton exchange membranes.⁷⁰ A complementary approach to electric field assisted ion exchange has been developed in the authors labs over the last years.^40,57 That approach is termed Charge Carrier Attachment Induced Transport (CAIT).⁷¹ A solid sample is brought into contact with a single metal electrode at the backside. A charge carrier beam is shined at the front side giving rise to a well-defined electrochemical surface potential. The gradient of the electrochemical potential towards the backside induces charge carrier transport in the sample, which can either be detected as a current in the outer circuit or as a concentration depth profile inside the sample, provided that the native and the foreign carrier are distinguishable (foreign ion CAIT). Concentration depth profiles are here in general quantified by means of secondary ion mass spectrometry (SIMS). The pivotal advantage of the CAIT approach is that it completely avoids any charge carrier blocking. As a matter of fact, the development and the success of the CAIT/SIMS approach created the necessity to develop a dedicated NPP code, which forms the basis of the program suite reported in this work. Successful applications of the NPP code include quantification of replacement zones in ion conducting solids,^{51–54,58,72} the quantification of site energy distribution,^{41,61,62,73,74} diffusion of ions in polymers,^31,75,76 polyelectrolyte membranes,^33,77 mold compounds,⁴⁵ and disintegration of metal electrodes.⁵⁵

Ultimately, the MAR_CCT program suite is considered to currently be the only published code capable of modelling CAIT/SIMS experiments, but at the same time the code is capable of covering most of the other transport phenomena covered by other codes in the field. More specifically the code also allows to simulate Alkali-Proton Substitution (APS), Plasma-CAIT or electro-thermal poling experiments. Note, that MAR_CCT is made available under public license in contrast to some of the other powerful program suites known.

In the following we present a NPP-based simulator program, which is able to calculate and evaluate time dependent concentration-depth profiles that evolve when concentration or electric potential gradients are applied to a solid electrolyte. Information concerning the space dependent electric potential and effective diffusion coefficients as well as the time dependent ion fluxes are provided. A variety of experimental conditions can be simulated including CAIT experiment, thermal electro-poling and alkali-proton substitution. The diffusion coefficients can either be externally provided or they can be calculated employing a selection of parameterized functions. Alternatively, the diffusion coefficients can be calculated from a user-defined SED. The program suite is available at.⁷⁸

Program structure

Upon opening the program, a graphical user interface (GUI) with the main navigation menu appears. The user interface is written in Python 3.10. It contains three main menu items. The default setting is “Input and Execution”. Here, the users have the option to define experimental conditions and sample parameters with the help of which an input file is created for the actual calculation. The calculation itself employs a Fortran95 code that can then be started to run in the background. The calculation continues to run even after the user interface is closed, as long as the computer is not switched off.

When entering the sheet “Experiment and sample parameters”, default values are already filled into all input fields. In general, these values have to be modified in order to match specific experimental conditions. Using the default values leads to the benchmark calculations presented in this paper.

Calculations already finished can be evaluated under the menu point “Visualization”. The concentration-depth profiles, the space dependent electric potential and the carrier flux through the sample as well as their time development can be displayed. If a SED is provided, the program can calculate the corresponding concentration dependent diffusion coefficients.

The menu item “Evolution Video” makes it possible to create a movie from the snapshots of the concentration-depth profiles, which then shows the temporal evolution of the depth profile as a film in mp4 format. It is also possible to choose an extended movie showing the local occupation of the SED and the associated diffusion coefficient at three different locations within the profile are displayed as a function of time together with the actual concentration profile.

In the following section we will describe the functionality of the Fortran95 code as well as the options of the GUI.

The transport program

Fundamental equations

As a basis for the transport calculation, the program uses the coupled set of one-dimensional Nernst–Planck and Poisson equations. The Nernst–Planck equations (one per chemical species i) describes the ion flux density, J_i, due to a gradient in the electrochemical potential that includes concentration gradients as well as electric potential gradients. The electric potential includes externally applied voltages as well as the electric potential that may occur due to the accumulation of charge inside the sample. The equation is given bywhere D_i(n_i) is the concentration dependent diffusion coefficient, n_i is the respective charge carrier density and ϕ is the electrical potential. The charge carriers have the charge Z_ie and k_BT is Boltzmann's constant times the temperature. The electric potential is determined by the externally applied potential and the charge carrier distribution inside the investigated material. Therefore, the charge carrier translocation inside the sample influences the functional dependence of the electric potential on space such that the electric potential needs to be updated when the charge carrier distribution is modified. Therefore, the Poisson equation is calculated self-consistentlywith the vacuum dielectric permeability ε₀ and the dielectric constant ε_r of the sample. The bulk ion density of the respective species i is given by n_i,0 which defines the state of the electroneutral sample prior of the experiment. The time dependence of the charge carrier motion enters via the continuity equation (Fick's 2nd law).The boundary conditions of the electric potential depend on the experimental conditions. The program code allows to decide between two options: (a) one electrode at the back side of the sample and an open front side where charge carriers may enter via an ion beam (CAIT) or (b) two electrodes, one at the front side and the second on the backside of the sample (thermal electro-poling or APS).

Discretization and boundary conditions

The Fortran code evaluates eqn (1)–(3) on a finite space and time grid with piece wise equi-distant grid elements. The grid elements are usually small in the vicinity of the electrodes or at the open front side and larger deep in the bulk of the sample. An appropriate choice of the space grid is important to guarantee accurate calculations and short computation times. In order to calculate eqn (1)–(3) on the grid, the discretized form of these equations is set up. Therefore, we assign an average charge carrier density n_i(x_a) and an average potential ϕ(x_a) to each of the space grid elements a. We findwhere J^a↔a+1_i is the ion flux between space grid element a and a + 1. The diffusion coefficient at the boundary between these two space grid elements is given by D^a,a+1_i. Additionally, we use the abbreviations and . The time grid element is given by Δt. The Fortran95 code uses a Runge–Kutta routine fourth order to calculate the time evolution of the concentrations.

The electric potential at the front and back side of the sample is determined by the mode of measurement. In the case of thermal electro-poling or APS, both sides of the sample are in contact with an electrode which is set to a fixed electric potential; ϕ_front = constant and ϕ_back = constant. If the CAIT option is applied, the front side of the sample is not in direct contact with a solid metal electrode but an ion beam impinges on the surface. Since the ions are decelerated in front of the sample surface there is always a cloud of charge carriers in front of the sample surface which can quickly adjust to the charging state of the front surface and serves as an infinite ion source. We assume that this cloud of carriers is in equilibrium with the surface population such that effectively the electric potential gradient at the surface becomes zero representing an ideal reversible electrode. Note, that energy barriers which could potentially be operative when charge carriers are injected into a sample are neglected. So far, we do not have evidence of such a barrier. We allow charge carriers to enter the sample as long as the surface potential ϕ_front is lower than the source potential U_R of the charge carriers in the beam. This is consistent with the picture of ions of a given kinetic energy E_kin = eU_R moving toward a charged surface.⁴⁰ In such a situation, the electric potential shows a maximum at or close to the sample front surface such that the gradient of the electric potential vanishes there. The ion flux density that reaches the surface when ϕ_front < U_R is given by with I_i,s the source current of ion species i. If the front potential exceeds the surface potential the ion flux is set to zero. In the case of APS, the source current can directly be given by the (fixed) electrode potential.

Depending on the experimental mode, different approaches for the calculation of the electric potential are beneficial. If two electrodes are in contact with the sample (thermal electro-poling and APS), it is beneficial to accurately calculate the electric potential every time step. Under these conditions, the program uses an analytical formula for the calculation of the potential:

where L is the thickness of the sample.

In the case of a CAIT experiment, the program uses an iterative solution of the eqn (5) to calculate the electric potential since the deposition of charge carriers at the surface combined with the electrically open front surface may ultimately lead to numerical challenges if the analytical formula (7) is used. Within every time step charge carriers are deposited as long as ϕ_front < U_R. The deposition of charge carriers will lead to a finite overshooting of the electric potential beyond U_R at the front side followed by a time period where no charge carriers may reach the surface. Therefore, the potential at the front side will oscillate around the equilibrium value. Since the size of the overshooting is directly related to the size of the time grid element (longer deposition periods mean higher potential changes), the time grid increments must be kept very small which strongly affects the computation time. However, the overshooting is only an artifact of the deposition model. In a real system the cloud of charge carriers in front of the sample surface will ensure that ϕ_front = U_R holds, as long as the source current is sufficiently strong. When the source current is not strong enough a dynamic equilibrium with ϕ_front < U_R will be reached. However, the value of ϕ_front is not known prior to the calculation such that it cannot be set as an external parameter. Additionally, the resistance of the sample can change, when charge carriers are exchanged. This directly leads to a change of ϕ_front such that ϕ_front in general can be time dependent. Since we also want to be able to model experiments where the front potential is not known, we need to use a numeric trick to enhance computation time and at the same time keep the errors small. It shows that an iterative solution of eqn (5) is computationally beneficial.

Eqn (5) shows that the electric potential in grid element a is a function of the electric potential in the two adjacent space grid elements a + 1 and a − 1 as well as the charge concentration in the space grid element a itself. Hence, the modification of the electric potential in a certain space grid element influences the potential in all other space grid elements while a modification of the potential in these space grid elements couples back to the potential in the original space grid element. As a consequence, eqn (5) need to be solved iteratively. In general, it needs millions of iterations to fully converge the potential every time step which again would lead to unacceptably long computation times. However, one may converge the potential during the time evolution of the profile as CAIT. This leads to a situation where the potential quickly approaches a stationary situation and additionally it allows for much larger time steps.

The trick is not to wait for the full convergence of the potential but calculate only few iterations and then continue with the next time step. This way only charge carriers in the direct vicinity of the surface feel the changed potential situation there. Charge carriers deeper in the sample realize the changes only with a little time delay as the iteration of the potential only affects one additional grid element every iteration step. During this time delay the potential in the first increment oscillates around the equilibrium value due to the finite time steps, the charge carriers in the deeper increments then only experience an average potential which is much closer to the equilibrium value. In practice, ten potential iterations per time step lead to a reasonably quick convergence while the computation time remains short and the error to the profile is negligible. If the time grid is nevertheless too coarse, the Runge–Kutta routine will lead to a divergence after some few time steps. In the output files NAN appears instead of a number value.

Please note, whenever potential gradients occur that are exceeding the dielectric breakdown stability of the material, e.g., in thermal electro-poling experiments, this kind of approximation becomes inapplicable. The partially converged potential will lead to artificial dielectric breakdown at positions where there would not be an electric breakdown if the potential is fully converged. The consequence is that the electronic background charge is irreversibly modified and the calculation will yield fundamentally wrong charge carrier distributions. The treatment of dielectric breakdown is the topic of the following section.

Dielectric breakdown

Prior to the start of the calculations, the sample is assumed to be electroneutral. In order to compensate the charge of the mobile ionic charge carriers a corresponding amount of electrons is assigned to each increment.

In the case of thermal electro-poling calculations, the motion of charge carriers leads to the buildup of large electric fields in the vicinity of the electrodes. As no external charge carrier source is present that could compensate the accumulation of charge, the electric field strength will exceed the dielectric beak-down field-strength at a certain point of time.

The program checks every time step whether the electric potential exceeds the dielectric breakdown field strength in any increment. If this is the case, electrons are removed from that increment. The amount of charge carriers removed in one step corresponds to

After removing the electrons, the electric potential is recalculated in the sample. The procedure of removing and recalculating the potential is repeated until the electric field remains smaller than the breakdown field strength everywhere in the sample.

Concentration dependent diffusion coefficient

The diffusion coefficients that enter eqn (1) is in general concentration dependent. There are three options how the diffusion coefficients may be entered into the calculation. They can either be provided in a separate file. The code requires a certain format to handle the data. We come back to that later. Alternatively, the program is capable to calculate the diffusion coefficients from a generic function. One may choose between four functional types of variation of D with carrier density: constant diffusion coefficients, linearly varying diffusion coefficients, exponentially varying diffusion coefficients and sigmoidal dependence.

As a last option, the diffusion coefficient can be calculated from giving a site energy distribution (SED). In this case, one needs to give the density of all sites S₀ which is usually some percent higher than the density of mobile ions and the width of the site energy distribution. The SED is calculated via a sin² function of the form

where Γ is a measure for the energetic width of the distribution. E₀ is calculated such that for bulk ion density the highest occupied site energy is set equal to −E_act, where E_act is the bulk activation energy for the dc transport limit. The parameters S₀, E_act, and Γ need to be provided as an input where the first two are given by experimental observables and the width Γ is a true parameter. We assume that the SED is depopulated from higher energy states toward lower energy states such that the effective activation energy increases. The concentration dependent activation energy is then calculated according toWe assume that the temperature dependence of the diffusion coefficient follows a simple Arrhenius behaviorIn many cases, best agreement between experiment and theory could be found, if the diffusion coefficient of the foreign ion was set to be (almost) concentration independent.^{40,41,53,54,58,61,62} This can be realized in the program by artificially setting the width Γ for the foreign ion to very small values (<μeV). The pre-exponential factor D₀ is a parameter that is determined by matching the calculated conductivity of the sample to the measured conductivity.

The effective diffusion coefficients that enter the Nernst–Planck eqn (1) are directly related to diagonal Onsager transport coefficients. The relation has been elaborated in ref. 62. In a system of two mobile charge carrier species, the effective diffusion coefficient includes the dependence of the chemical potential of one ion species on the concentration of the second ion species. Off-diagonal Onsager coefficients are neglected in the current version of MAR_CCT but may be included in future updates.

It may be conjectured that the concentration dependence of observables in electrochemistry is often modeled through the usage of activities involving concentration dependent thermodynamic factors. In general, diffusion coefficients may either be modelled in the particle density (i.e. concentration) domain or in the activity domain. Here, preference is given to the particle focused picture, motivated by research interest aiming at atomically (i.e. particle resolved) structure and transport properties.⁷⁹

Flow chart of the Fortran code

To illustrate the structure of the Fortran program, the following section contains a flow chart of the code.

When the Fortran code is executed, at first the natural constants and the file input.dat are imported. From the parameters given, a finite grid for a real space coordinate and the time coordinate is created, the concentration dependent diffusion coefficients are calculated and the initial values of charge carrier concentration and electric potential are assigned to the space grid elements. These starting conditions are used to calculate the electric potential inside the sample and the initial ion fluxes. The ion flux enters the continuity equation (Fick's 2nd law) which is evaluated using a 4th order Runge–Kutta routine. The continuity equation yields an update of the ion densities which enter Poisson's equation and yields a new updated electric potential. The potential and the updated ion densities then enter the Nernst–Planck equation which yields the flux for the continuity equation. The cycle continues until the time loop is completed. In between, the space dependent carrier densities and electric potential are stored in snapshots such that the time evolution of both can be recorded.

Input of parameters and program execution

The “Input and Execution” entry in the main menu leads to the part of the program in which the input file for executing the calculation can be created. The input file is created in several steps.

On the start page, the user can choose between experiments with one or two electrodes. The collective term “Experiments with one electrode” also includes CAIT experiments in particular. Thermal electro-poling and field-driven proton exchange experiments fall in the category of experiments using two electrodes.

One electrode. Selecting “Experiments with one electrode”, leads to a page on which the voltage applied to the back of the electrode has to be entered. Further entries are the kinetic energy of the ions in the ion beam and the ion current directed at the sample. The temperature of the sample must also be set here. Pressing “Specify sample” confirms the settings and opens the form in which the sample specifications can be entered.

On the top left the number of the mobile charge carrier species are entered. If the number of charge carrier species is changed, the setting needs to be confirmed by pressing “Set”. Depending on the number of charge carrier species, the table then shows a corresponding number of columns. In the first row of the table, one has the name the charge carrier species, for example with their chemical symbol. The next row is reserved for the bulk carrier densities. If a charge carrier species is not contained in the material prior to the experiment one simply enters “0d0”. The format of the number value corresponds to the format Fortran uses for the number values. The “d” indicates that the number is read in with double precision (16 digits accuracy) and is short for “ten to the power”. Entering 7d27 therefore means 7 × 10²⁷. The next line reads the carrier charge in units of the elementary charge +e. The fourth line defines whether the charge carrier species may be neutralized at the back side electrode. Enter ‘y’ for yes or ‘n’ for no. If the charge carriers neutralize, they are effectively removed from the glass and a corresponding charge enters the calculation of the neutralization current at the back side electrode. Please note, that the neutralization process experimentally leads to the formation of a neutral metal layer between sample and electrode, thus outside the original sample. The neutralization process at the backside electrode is much faster than the flux of ions through any grid element. As a consequence, the transport is solely limited by the flux described by the NPP equations. At this point, an explicit account of the neutralization reaction is not yet included in the program code but is considered to be included in future program versions.

The last line in the table defines whether the ion species enters the sample via the ion beam. Please enter again ‘y’ or ‘n’. Complementing the table four additional parameter have to be given. The first one is the thickness of the sample and the second one the diameter of the by the ion beam irradiated area. Additionally, the dielectric constant of the material and the breakdown voltage of the material need to be provided. If the dielectric breakdown should be artificially excluded, put a very large value here (>10²⁰). Pressing one of the three buttons at the right bottom confirms the settings and opens a window in which the diffusion coefficient can be specified.

The concentration dependent diffusion coefficient can either be determined form a SED that is given in a file or it can be directly read in or in a third case it can be determined employing a generic functional form. The top button leads to the window in which the SED can be read in from a file. This option allows to read complex SEDs – for example bimodal SEDs or even more complicated forms. The file needs to contain two columns. In the first column, the site energies referenced to the threshold for long-range transport are given as negative numbers in eV and the second column contains the site density. The curve should be normalized such that the integral over the curve yields 1. Since Fortran inputs are position sensitive, each column must have the format: 0.12345678 × 10⁺¹². The two columns should be separated by 8 blanks.

The center of the window is given by a table with as many rows as there are ion species. Pressing select file allows to choose a file from the local folders. Additionally, the bulk activation energy needs to be given as well as the pre-exponential factor, D₀, as the code uses eqn (11) to calculate the diffusion coefficient from the concentration dependent activation energies. Pressing “Define Grids” confirms the settings. Reading the diffusion coefficients from a file is analogous to the procedure described above. The sole difference is that the two columns must yield the charge carrier density in 1 m⁻³ in the first column and the diffusion coefficients in m² s⁻¹ in the second.

If the diffusion coefficient shall be computed using a functional form of the diffusion coefficient itself or by a SED that is provided by a sin² function according to eqn (9), press “Diffusion coefficients from parameters”. The next window allows then to decide between five different approximations of the diffusion coefficient. Option 1 is calculating the diffusion coefficient from a SED. Therefore, the bulk activation energy, the Γ Parameter in eqn (9), the overall density of sites and the pre-exponential factor D₀ need to be given for every ion species. The diffusion coefficient can be shown in the preview or in an extra window. Options 2 to 5 are calculating the diffusion coefficient by the formulas given in Table 1. One can change between the options by changing the number in the field “Please choose the approximation”. Confirm the choice by pressing “Done”.

Table 1 Functional forms for defining diffusion coefficients

Constant	Di(ni) = D0,i
Linear dependence
Exponential dependence
Sigmoidal dependence

---

After confirming the input of the diffusion coefficients, one may enter the space and time grid. The grid is piece wise equidistant. It is strongly recommended to use fine space grids wherever concentrations are expected to strongly change (e.g. in regions where a diffusion front is expected). Additionally, one should not vary the size of the grid by much more than a single order of magnitude between two adjacent regions. Large differences can cause numerical errors and artifacts to appear in the calculation. It is possible to set up to 10 space intervals and 5 time intervals. The time grid should be chosen fine when the potential changes rapidly at the front surface. The window provides a thickness check in which the applied grid is compared to the size of the sample. Both numbers: sample thickness and the grid range must agree, otherwise the calculation leads to unphysical results. The number of potential iterations per time interval may be entered here has well. 10 iterations is a value that allows fast convergence at acceptable computation times.

After confirming one reaches the last window, where the output options are defined. If the output for the ion current is activated, one has to provide a number that defines after how many time steps the ion current should be written into a file. The number of time steps is typically on the order of 10⁸. If the number is chosen too small, the output file for the ion current will eventually become very big with potentially millions of entries. Under these conditions, also the code will become slower as the output routine is called very often. Eventually, the number of snapshots can be set. If this option is chosen, the program code will generate a given number of snapshots where quantities as the charge carrier profiles and the electric potential will be written in files. This option is mandatory if the temporal evolution of concentration is of interest.

Eventually, the input file for the Fortran code can be created by pressing “Create input file”. Subsequently, the calculation can be executed by pressing “Start calculation”. The Fortran code is then executed employing the input file created above. The folder in which the results can be found is the one which has been selected as storage path. Please note, that the storage path should be empty to avoid overwriting calculations and ensure that executable, input file and results in the chosen folder are consistent to each other. If the folder is non-empty the program will complain. The code runs in the background such that the GUI can be closed without affecting the calculations. In the end one can decide whether the GUI is closed or a new input file is created.

Two electrodes. If the user decides to calculate an experiment using two metal electrodes in contact with the sample, the procedure is almost identical to the procedure described above. The difference occurs only in the window where the sample setup is defined. For two electrodes, two potentials (one at the front side and one at the back side) as well as the temperature need to be provided. Additionally, one needs to decide whether ions are injected into the sample, for example via an APS process, or whether no ions source is provided, for example when thermal electro-poling should be calculated. Enabling the APS option will lead to the injection of the same amount of charge carriers as is neutralized at the back side of the sample such that the sample remains electro-neutral.

Vizualization

The menu entry “Visualization” allows to analyze simulations that have been run with the Fortran code. At first one needs to specify the folder in which the simulation data are stored. One may also give a run Id which is then added to the title of the created figures. Computed concentration profiles may be compared to experimental profiles. Therefore, the folder in which the experimental data are stored need to be provided. If the path for the experimental data is kept empty, no experimental data will be included in the evaluation.

By pressing “Open data folders”, the simulation is checked and analyzed. The fields “Bulk ion density”, “Sample surface area” and “List of charge carriers will be filled in”. The final concentration profiles appear in the large white field in the center of the window. The experimental profiles should be provided in a comma-separated value file (csv). The names of the experimental data files should be the same as the name of the carrier species and first and second column in these files should be named with “depth” and “density” such that the program can handle the experimental data. There is one data file per mobile species considered.

The toolbar above the graphic allows the user to zoom into certain areas of the graphic or to define the displayed axis range directly. The finished image can then be saved under a selected name using the diskette symbol. The storage format is usually png.

On the right-hand side of the window, one has the option to plot a selection of different quantities. The plots are then shown under the respective menu entry. The following plots can be automatically created:

The “Profile evolution” tab shows how the concentration profiles evolve with ongoing time.

Pressing the “Final potential” tab yields the electric potential across the entire sample at the end of the calculation as well as a zoom into the region of the diffusion profile.

The “Potential evolution” tab shows the time evolution of the electric potential. Note, that the first potential snapshots may not be fully converged when the potential is calculated by iteration. As a consequence, the respective window allows the user to skip the first snapshots.

Pressing the “Diffusion coefficient” tab creates a plot where the diffusion coefficient is shown as function of the fractional abundance of the respective ion species. This plot is usually shown in the format of a semi-logarithmic plot as the diffusion coefficient may vary over many orders of magnitude.

The “Site energy distributions” tab allows checking the potential energy landscape used to calculate the diffusion coefficients (if in use).

The “Current” tab allows checking the neutralization current operative. If the potential is calculated by the iteration method, the neutralization current will first oscillate around the stationary state current and eventually converge to that value. The value of the injected charge is equal to the integral of the neutralization current in very good approximation. The oscillations do not affect this integral as they are sinusoidal such that their integral converges to the value of the stationary state current already after few cycles, even when the oscillation itself still continues.

All plots contain a tool bar where the user can modify the axis and eventually plot the final graphics into a png-file.

Video of time evolution

The program offers the possibility to create a video from the time evolution of the concentration depth profiles. Therefore, the user needs to choose the menu entry “Evolution Video” and select the path where the simulation data are stored. The entries Run Id, bulk ion density, sample surface, list of charge carriers and native carriers are filled in automatically from reading the input file with which the profiles have been created.

Two different kinds of videos are available. On one hand, one can create a video where exclusively the concentration profile is evaluated. The two buttons on the left side belong to this option. Before creating the video by clicking “Make the film” one needs to create png-graphics from the profile snapshots that the program has calculated during the runtime. This can be done by clicking “Create the plots”. If the calculation is performed using a site energy distribution from which the diffusion coefficients are calculated, one can also create a film where the profile is evaluated at three different positions with respect to the charge carrier density of the native ion species. From the charge carrier density, the program calculates the actual population of the SED and the locally active effective diffusion coefficient. With ongoing diffusion, on can follow the depopulation of the SED and the simultaneously decreasing diffusion coefficient.

Benchmark calculations

If the user starts the GUI to create an input file for the calculation, all input fields in the interface are already pre-filled. Using these values, one arrives at one of the three benchmark calculations presented below.

We provide benchmark calculations for CAIT, Poling, and APS experiments. These are generic calculations not based on any actual experiment or sample. However, the parameters are chosen so that they are typical for the respective type of experiments and can be easily adapted to a real experiment. The parameters used are given in Tables 2–4.

Table 2 Parameters of the CAIT benchmark calculation

Cait benchmark
GUI window	Input field	Value
One electrode	Ekin/eV	10
One electrode	Beam current/nA	5
One electrode	Temperature/K	350
One electrode	electrode potential/V	0
Sample	Number of charge carrier species	2
Sample	Name of the carrier species	Na, K
Sample	Bulk carrier density	1d28, 0
Sample	Carrier charge/e	1, 1
Sample	Neutralization at electrodes	y, y
Sample	Carrier is part of the beam	n, y
Sample	Sample thickness/mm	1.5
Sample	Mask diameter/mm	10
Sample	Dielectric constant	10
Sample	Breakdown voltage/(V m^−1)	5d8
Diff-coefficient (from parameters)	Approximation	1
Diff-coefficient (from parameters)	Eact(bulk)/eV	1.00, 1.20
Diff-coefficient (from parameters)	Width parameter (sin^2)/eV	0.18, 0.0001
Diff-coefficient (from parameters)	Density of sites (total)/(1 m^−3)	1.05d28, 1.05d28
Diff-coefficient (from parameters)	D0/(m^2 s^−1)	3d−5, 3d−5
Diff-coefficient (from parameters)	Number of integration steps	5000
Grids	Space grid intervals	5
Grids	Deta x/μm, increments	3d−3, 50 0.135d0, 10 1.35d0, 10 13.5d0, 10 135d0, 10
Grids	Time grid intervals	1
Grids	Delta t/s, increments	3d−2, 70 000 000
Grids	Factor	1
Grids	Potential iterations	10
Output	All outputs enabled
	Current output every × time steps	10 000
	Number of snapshots	100

---

Table 3 Parameters of the poling benchmark calculation

Poling benchmark
GUI Window	Input field	Value
Two electrodes	Electrode potential (left)/V	50
Two electrodes	Temperature/K	350
Two electrodes	Electrode potential(right)/V	0
Sample	Number of charge carrier species	1
Sample	Name of the carrier species	K
Sample	Bulk carrier density	7d27
Sample	Carrier charge/e	1
Sample	Neutralization at electrodes	y
Sample	Sample thickness/mm	0.5
Sample	Mask diameter/mm	8
Sample	Dielectric constant	10
Sample	Breakdown voltage/(V m^−1)	5.5d8
Diff-coefficient (from parameters)	Approximation	1
Diff-coefficient (from parameters)	Eact(bulk)/eV	0.99
Diff-coefficient (from parameters)	Width parameter (sin^2)/eV	0.18
Diff-coefficient (from parameters)	Density of sites (total)/(1 m^−3)	7.35d27
Diff-coefficient (from parameters)	D0/(m^2 s^−1)	1d−4
Diff-coefficient (from parameters)	Number of integration steps	5000
Grids	Space grid intervals	5
Grids	Delta x/μm, increments	2d−3, 200 0.1d0, 30 1.0d0, 10 10.d0, 10 38.66d0, 10
Grids	Time grid intervals	1
Grids	Delta t/s, increments	6d−3, 40 000 000
Grids	Factor	1
Grids	Potential iterations	10
Output	All outputs enabled
	Current output every × time steps	10 000
	Number of snapshots	100

---

Table 4 Parameters of the APS benchmark calculation

APS benchmark
GUI window	Input field	Value
Two electrodes	Electrode potential (left)/V	40
Two electrodes	Temperature/K	380
Two electrodes	Electrode potential(right)/V	0
Sample	Number of charge carrier species	2
Sample	Name of the carrier species	Na, H
Sample	Bulk carrier density	1d28, 0d0
Sample	Carrier charge/e	1, 1
Sample	Neutralization at electrodes	y, y
Sample	Carrier entering the sample	n, y
Sample	Sample thickness/mm	1.0
Sample	Mask diameter/mm	14
Sample	Dielectric constant	10
Diff-coefficient (from parameters)	Approximation	1
Diff-coefficient (from parameters)	Eact(bulk)/eV	0.80, 0
Diff-coefficient (from parameters)	Width parameter (sin^2)/eV	0.066d0, 0.0001d0
Diff-coefficient (from parameters)	Density of sites (total)/(1 m^−3)	1.05d28, 1.05d28
Diff-coefficient (from parameters)	D0/(m^2 s^−1)	3d−8, 1.1d−19
Diff-coefficient (from parameters)	Number of integration steps	1000
Grids	Space grid intervals	3
Grids	Deta x/μm, increments	5d−3, 300 19.925d0, 20 150.0d0, 4
Grids	Time grid intervals	1
Grids	Deta t/s, increments	2.5d−3, 130 000 000
Grids	Factor	1
Grids	Potential iterations	10
Output	All outputs enabled
	Current output every × time steps	1 368 000
	Number of snapshots	100

---

CAIT benchmark

The corresponding calculation yields the results of a K⁺-CAIT on a sample that natively contains mobile Na⁺. The profiles show that K⁺ has replaced the native Na⁺ in the first 70 nm below the sample surface. The diffusion zone however shows that about 20% of the Na⁺ has not yet been replaced, a direct consequence of the concentration dependent Na⁺-Diffusion coefficient. The profile is shown in Fig. 1. The corresponding diffusion coefficients are given in Fig. 2. The diffusion coefficient of the native ion species Na⁺ shows a pronounced concentration dependence while the diffusion coefficient of the foreign ion K⁺ is essentially concentration independent. It has been discussed in literature that a concentration independent foreign ion diffusion coefficient is observed in many CAIT experiments when the foreign ions species has a lower diffusion coefficient than the native ion species.^{40,41,53,54,58,61,62} It has been suggested that the constant foreign diffusion coefficient is a consequence of the replacement process where the site energy distribution of the native ion is depopulated energetically top down while the foreign ion occupies states with low activation energies already at the very beginning of the experiment. With ongoing time in the CAIT experiment only states with high site energies are occupied such that the population of the foreign ion also occurs top down. As a consequence, the effective activation energy for the foreign ion remains constant throughout the experiment.

Fig. 1 Concentration depth profile – CAIT benchmark calculation.

Fig. 2 Diffusion coefficient – CAIT benchmark calculation.

When a constant foreign diffusion coefficient is desired and the SED is used to calculate the concentration dependent diffusion coefficient, the SED width Γ of the foreign ion should be artificially set to very small values. In this situation the code will create an (almost) concentration independent foreign diffusion coefficient. Please note, that choosing zero for the width Γ will lead to infinite number values in eqn (9).

The final electric potential is provided in Fig. 3. The graphic shows that the slope of the electric potential in the diffusion zone is larger than the slope of the electric potential in the unmodified bulk of the sample. The behavior reflects the larger local specific resistance of the glass in the diffusion zone due to the replacement of Na⁺ by the slower K⁺ ions. Since the gradient of the potential is the electric field, there is a larger electric field in the diffusion zone than in the rest of the sample. A direct consequence of this behavior is that there is a small negative excess charge at the diffusion front which stems from unoccupied sites that have already been abandoned from Na⁺ and K⁺ has not yet occupied them.

Fig. 3 Final electric potential – CAIT benchmark calculation.

Poling benchmark

The poling benchmark calculation shows the results of a poling experiment on a generic glass that contains mobile K⁺ ions. The resulting profile shown in Fig. 4: indicates a zone of about 90 nm below the front electrode where 90% of the native K⁺ concentration is depleted. The observation that about 10% of the K⁺ remains in the zone is a result of the strongly concentration dependent K⁺ diffusion coefficient which decreases by up to seven orders (Fig. 5) of magnitude. If the concentration dependence of the diffusion coefficient is less pronounced as it often is the case in realistic glasses, the depletion zone will be completely emptied due to the strongly increased electric field in the depletion zone (Fig. 6). Note, that only native mobile cations have to be entered by the user. The necessary density of electrons is automatically calculated by the program based on initial electro-neutrality.

Fig. 4 Concentration depth profile – poling benchmark calculation.

Fig. 5 Diffusion coefficient – poling benchmark calculation.

Fig. 6 Final electric potential – poling benchmark calculation.

APS benchmark

The APS benchmark calculation simulates a proton substitution in a generic Na⁺ conducting glass. The final profile of the calculations shows that in the first 300 nm below the sample surface a significant amount of Na⁺ has been replaced by H⁺ (Fig. 7). As the proton diffusion coefficient is higher than the typical diffusion coefficient of alkali ions entering the sample like in a CAIT experiment, the H⁺ profiles reach deep into the material (Fig. 8). Profiles as deep as a μm are quite common. The slopes of the two concentration profiles are smaller than in typical CAIT experiments regarding absolute values. Since the diffusion coefficient for low concentrations of Na⁺ drops below the diffusion coefficient of H⁺, the electric field in the diffusion zone is larger than in the bulk (Fig. 9). If the H⁺ diffusion coefficient remains larger than the native diffusion coefficients the electric field in the diffusion zone will be smaller than in the bulk.

Fig. 7 Concentration depth profile – APS benchmark calculation.

Fig. 8 Diffusion coefficient – APS benchmark calculation.

Fig. 9 Final electric potential – APS benchmark calculation.

Conclusions

A new program package, termed MAR_CCT (MARburg program suite for modelling Charge Carrier Transport), capable of describing ion transport in solid-state electrolytes on the basis of the Nernst–Planck–Poisson (NPP) equations is presented. Within this NPP formalism the flux of charge carriers induced by the gradient of the electrochemical potential is considered. The motivation for developing this program suite originated from the charge attachment induces transport (CAIT) experiment developed in the Weitzel group. To date, there is no other program code available for describing CAIT experiments. However, MAR_CCT has already been successfully applied to the description of a broad variety of charge carrier transport experiments, including, electric field assisted ion-exchange, alkali proton substitution (APS), electro-poling but also chemical diffusion experiments.

The MAR_CCT calculations are performed by solving the coupled Nernst–Planck and Poisson equations as a function of time on a given spatial grid. The core of the program is a 4th order Runge–Kutta formalism. The source code is written in Fortran. This Fortran code is embedded in a graphical user interface (GUI) based on a python code that includes the possibility to start the Fortran code right from the GUI. Additional interfaces provide the possibility to visualize, evaluate and export the data. For example, the charge carrier distribution as well as the electric potential within the sample can be visualized as a function of time. The GUI also allows the creation of short movies showing the evolution of the depth profiles.

Within this manuscript, three examples for an application of MAR_CCT are provided, (i) for CAIT, (ii) for electro-poling and (iii) for APS experiments demonstrating the broad applicability of the program suite.

The program has been developed for describing charge carrier transport experiments, where the crucial charge carriers can be mono- and bi-valent ions of either polarity and electrons. Processes driven by forces not dominated by the gradient of the electrochemical potential (e.g., electron–phonon scattering) are currently not covered by the program. Coupling of charge carriers is currently restricted to the Poisson level. Extensions may be required for describing mixed ionic electronic conductors (MIECs) in the context of Lithium-ion batteries in the future. Also, polarization dominated processes, e.g., electrode polarization, are not in the focus of the current version of the program. However, extensions will be made available in the future.

Conflicts of interest

There are no conflicts to declare.

Data availability

The program package described in this work is available at https://www.uni-marburg.de/en/fb15/for5065/MAR_CCT.

Acknowledgements

This research project is funded by the German Science Foundation (DFG) as part of the research unit FOR 5065 (“Energy Landscapes and Structure in Ion Conducting Solids”, ELSICS) project number 428906592 (project P1).

References

1. R. M. Ormerod, Chem. Soc. Rev., 2003, 32(1), 17.
2. A. J. Jacobson, Chem. Mater., 2010, 22(3), 660.
3. Z. Tao, L. Yan, J. Qiao, B. Wang, L. Zhang and J. Zhang, Prog. Mater. Sci., 2015, 74, 1.
4. R. Ramamoorthy, P. K. Dutta and S. A. Akbar, J. Mater. Sci., 2003, 38(21), 4271.
5. M. Li, J. Lu, Z. Chen and K. Amine, Adv. Mater., 2018, e1800561.
6. T. Kim, W. Song, D.-Y. Son, L. K. Ono and Y. Qi, J. Mater. Chem. A, 2019, 7(7), 2942.
7. W. A. Catterall, Neuron, 2010, 67(6), 915.
8. S. Choe, Nat. Rev. Neurosci., 2002, 3(2), 115.
9. E. Bourinet, C. Altier, M. E. Hildebrand, T. Trang, M. W. Salter and G. W. Zamponi, Physiol. Rev., 2014, 94(1), 81.
10. ed. B. Scrosati, Fast ion transport in solids: Proceedings of the NATO Advanced Research Workshop on Fast Ion Transport in Solids, Belgirate, Italy, September 20–26, 1992, Kluwer Acad. Publ, Dordrecht, NATO ASI series Ser. E, Applied sciences, 250, 1993.
11. J. Maier, Physical chemistry of ionic materials: Ions and electrons in solids, Wiley, Hoboken, NJ, 2023.
12. Ion and Molecule Transport in Membrane Systems, ed. V. V. Nikonenko, MDPI – Multidisciplinary Digital Publishing Institute, Basel, Switzerland, 2021.
13. K. Kohary, H. Cordes, S. D. Baranovskii, P. Thomas, S. Yamasaki, F. Hensel and J.-H. Wendorff, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63(9), 094202.
14. K. A. Fichthorn and W. H. Weinberg, J. Chem. Phys., 1991, 95(2), 1090.
15. D. T. Gillespie, J. Comput. Phys., 1978, 28(3), 395.
16. B. L. Hammond, W. A. Lester and P. J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, 1994.
17. W. Nernst, Z. Phys. Chem., 1889, 4U(1), 129.
18. M. Planck, Ann. Phys., 1890, 275(2), 161.
19. G. Gagneux and O. Millet, Appl. Math. Modell., 2016, 40(2), 846.
20. D. Meng, B. Zheng, G. Lin and M. L. Sushko, Commun. Comput. Phys., 2014, 16(5), 1298.
21. M. S. Metti, J. Xu and C. Liu, J. Comput. Phys., 2016, 306, 1.
22. M. Kato, J. Theor. Biol., 1995, 177(3), 299.
23. B. Lu, M. J. Holst, J. A. McCammon and Y. C. Zhou, J. Comput. Phys., 2010, 229(19), 6979.
24. B. Lu and Y. C. Zhou, Biophys. J., 2011, 100(10), 2475.
25. J. J. Jasielec, Electrochem, 2021, 2(2), 197.
26. J. Pods, J. Schönke and P. Bastian, Biophys. J., 2013, 105(1), 242.
27. M. Agnaou, M. A. Sadeghi, T. G. Tranter and J. T. Gostick, Comput. Geosci., 2020, 140, 104505.
28. M. Schmuck and M. Z. Bazant, SIAM J. Appl. Math., 2015, 75(3), 1369.
29. E. Samson, J. Marchand, J.-L. Robert and J.-P. Bournazel, Int. J. Numer. Meth. Eng., 1999, 46(12), 2043.
30. B. Grysakowski, B. Bożek and M. Danielewski, Defect Diffus. Forum, 2008, 273–276, 113.
31. S. Schulze, M. Schäfer, A. Greiner and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2013, 15(5), 1481.
32. T. Sokalski, P. Lingenfelter and A. Lewenstam, J. Phys. Chem. B, 2003, 107(11), 2443.
33. V. Motta, M. Schaefer, J. Huehn, R. Zierold, R. H. Blick, W. J. Parak and K.-M. Weitzel, Adv. Mater. Interfaces, 2020, 7(13), 2000419.
34. T. Sokalski, W. Kucza, M. Danielewski and A. Lewenstam, Anal. Chem., 2009, 81(12), 5016.
35. T. Sokalski and A. Lewenstam, Electrochem. Commun., 2001, 3(3), 107.
36. J. J. Jasielec, G. Lisak, M. Wagner, T. Sokalski and A. Lewenstam, Electroanalysis, 2013, 25(1), 133.
37. A. Subramaniam, J. Chen, T. Jang, N. R. Geise, R. M. Kasse, M. F. Toney and V. R. Subramanian, J. Electrochem. Soc., 2019, 166(15), A3806–A3819.
38. M. Lagnoni, C. Nicolella and A. Bertei, J. Electrochem. Soc., 2022, 169(2), 20570.
39. H. Zhu and R. J. Kee, Electrochim. Acta, 2016, 219, 70.
40. M. Schäfer and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2011, 13(45), 20112.
41. M. Schäfer and K.-M. Weitzel, Mater. Today Phys., 2018, 5(3), 12.
42. X. Jin, R. E. White and K. Huang, J. Electrochem. Soc., 2016, 163(13), A2702–A2719.
43. M. P. Tautschnig, N. M. Harrison and M. W. Finnis, Acta Mater., 2017, 132, 503.
44. Y. Cao, J. Wang and L.-Q. Chen, J. Am. Ceram. Soc., 2019, 102(7), 4136.
45. V. Motta, M. Schaefer, F. Noll, N. Hampp and K.-M. Weitzel, ECS J. Solid State Sci. Technol., 2020, 9(5), 53001.
46. J. Jamnik and J. Maier, Ber. Bunsenges. Phys. Chem., 1997, 101(1), 23.
47. Y. Fu, W. Liu, H. Mofidi and M. Zhang, J. Dyn. Diff. Equat., 2023, 35(2), 1585.
48. M. Ma, Z. Xu and L. Zhang, SIAM J. Appl. Math., 2021, 81(4), 1645.
49. P. W. Bates, J. N. Chen and M. J. Zhang, Math. Biosci. Eng., 2020, 17(4), 3736.
50. M. S. Kilic, M. Z. Bazant and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75(2 Pt 1), 21503.
51. D. Budina, J. Zakel, J. Martin, P. Menezes, M. Schaefer and K.-M. Weitzel, Z. Phys. Chem., 2014, 228(4–5, SI), 609.
52. J. Zakel, P. V. Menezes, M. Schaefer and K.-M. Weitzel, Solid State Ionics, 2013, 242, 20.
53. A. Mager, P. V. Menezes, J. Zakel, M. Schaefer and K.-M. Weitzel, J. Non-Cryst. Solids, 2016, 452, 231.
54. J. Martin, S. Mehrwald, M. Schaefer, T. Kramer, C. Jooss and K.-M. Weitzel, Electrochim. Acta, 2016, 191, 616.
55. Z. Herdegen, M. Schaefer and K.-M. Weitzel, Solid State Ionics, 2021, 359, 115533.
56. J. L. Wiemer, S. Mardeck, C. Zuelch and K.-M. Weitzel, Solid State Ionics, 2020, 357, 115469.
57. P. V. Menezes, J. Martin, M. Schäfer, H. Staesche, B. Roling and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2011, 13, 20123.
58. L. Rossrucker, P. V. Menezes, J. Zakel, M. Schäfer, B. Roling and K.-M. Weitzel, Z. Phys. Chem., 2012, 226, 11083.
59. J. Martin, M. Schäfer and K.-M. Weitzel, J. Non-Cryst. Solids, 2015, 430, 73.
60. J. L. Wiemer, M. Schäfer and K.-M. Weitzel, J. Phys. Chem. C, 2021, 125(9), 4977.
61. M. Schäfer, D. Budina and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2019, 21(47), 26251.
62. V. H. Gunawan, M. Schäfer and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2024, 26(19), 14430.
63. R. Kirchheim and D. Paulmann, J. Non-Cryst. Solids, 2001, 286(3), 210.
64. R. Kirchheim, J. Non-Cryst. Solids, 2000, 272(2), 85.
65. R. Kirchheim, Solid Solutions of Hydrogen in Complex Materials, Solid State Physics: Elsevier, 2004. p. 203–291.
66. B. Tu, M. Chen, Y. Xie, L. Zhang, B. Eisenberg and B. Lu, J. Comput. Chem., 2013, 34(24), 2065.
67. W. Clarke, L. J. Bennett, Y. Grudeva, J. M. Foster, G. Richardson and N. E. Courtier, J. Comput. Electron., 2023, 22, 364–382.
68. COMSOL Multiphysics®, COMSOL AB, Stockholm, Schweden, 2024.
69. A. K. Varshneya, Int. J. Appl. Glass Sci., 2010, 1(2), 131.
70. W. T. Grubb and L. W. Niedrach, J. Electrochem. Soc., 1960, 107(2), 131.
71. J. Martin, M. Gräf, T. Kramer, C. Jooss, M.-J. Choe, K. Thornton and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2017, 19, 9762.
72. A. Hein, J. Martin, M. Schäfer and K.-M. Weitzel, J. Phys. Chem. C, 2017, 121(6), 3203.
73. J. L. Wiemer, M. Schaefer and K.-M. Weitzel, J. Phys. Chem. C, 2021, 125(9, SI), 4977.
74. K. Rein and K.-M. Weitzel, J. Mater. Chem. A, 2024, 12(23), 14117.
75. S. Schulze, J. Zakel, M. Schäfer and K.-M. Weitzel, IEEE Trans. Dielectr. Electr. Insul., 2012, 19(4), 1167.
76. V. Wesp, J. Zakel, M. Schäfer, I. Paulus, A. Greiner and K.-M. Weitzel, Electrochim. Acta, 2015, 170, 122.
77. V. Wesp, M. Hermann, M. Schaefer, J. Huehn, W. J. Parak and K.-M. Weitzel, Phys. Chem. Chem. Phys., 2016, 18(6), 4345.
78. https://www.uni-marburg.de/en/fb15/for5065/MAR_CCT.
79. https://www.uni-marburg.de/en/fb15/for5065 (accessed July 17, 2025).

---