Back-exchange: A novel approach to quantifying oxygen diffusion and surface exchange in ambient atmo-spheres

A novel two-step Isotopic Exchange (IE) technique has been developed to investigate the inﬂuence of oxygen containing components of ambient air (such as H 2 O and CO 2 ) on the effective surface exchange coefﬁcient ( k ∗ ) of a common mixed ionic electronic conductor material. The two step ’back-exchange’ technique was used to introduce a tracer diffusion proﬁle, which was sub-sequently measured using Time-of-Flight Secondary Ion Mass Spectrometry (ToF-SIMS). The isotopic fraction of oxygen in a dense sample as a function of distance from the surface, before and after the second exchange step, could then be used to determine the surface exchange coefﬁcient in each atmosphere. A new analytical solution was found to the diffusion equation in a semi-inﬁnite domain with a variable surface exchange boundary, for the special case where D ∗ and k ∗ are constant for all exchange steps. This solution validated the results of a numerical, Crank-Nicolson type ﬁnite-difference simulation, which was used to extract the parameters from the experimental data. When modelling electrodes, D ∗ and k ∗ are important input parameters, which signiﬁcantly impact performance. In this study La 0 . 6 Sr 0 . 4 Co 0 . 2 Fe 0 . 8 O 3 − δ (LSCF6428) was investigated and it was found that the rate of exchange was increased by around 250% in ambient air compared to high purity oxygen at the same p O 2 . The three experiments performed in this study were used to validate the back-exchange approach and show its utility.


Introduction
The self-diffusion coefficient, D * , and the effective surface exchange coefficient, k * , of oxygen conducting materials are important parameters when selecting materials for applications such as Solid Oxide Fuel Cells (SOFC) electrodes, oxygen separation membranes and sensors. The measurement of D * and k * by Isotopic Exchange Depth Profiling (IEDP) was first described by Kilner et al. in 1984 1 . The technique has since been greatly refined to allow for the investigation of higher diffusivity materials using linescans 2 and parallel detection of secondary ions using ToF-SIMS 3 . The IEDP method described in 1 has a key limitation, which the approach presented in this work aims to resolve. The use of high purity isotopically enriched oxygen gas is significantly restricted by its cost and as such it is common to reclaim the gas used after each experiment. This precludes the use of enriched multicomponent gas mixtures, which are needed to model typical operating environments for many common applications. Furthermore, there is some discrepancy in the literature over values of D * and k * for widely studied materials, such as the mixed conductor La 0.6 Sr 0.4 Co 0.2 Fe 0.8 O 3−δ (LSCF6428). This may be the result of trace gasses present during the exchange, which are not always reported.
To extract D * and k * from experimental isotopic fraction data, the appropriate solution to the diffusion equation, eq. 1, is fitted by a non-linear least squares method. This expression was originally derived by Carslaw and Jaeger 4 as the solution to sys. 3, but was more famously communicated by Crank in his book "The Mathematics of Diffusion" (hence "Crank's solution") 5 .
This function is the solution to the 1-dimensional (1D) diffusion equation in a semi-infinite domain, with a surface exchange boundary and uniform initial conditions, where C is the nor-malised isotopic fraction found using eq. 2 (N.B. this normalisation notation is consistent throughout this paper and the dash should not be confused with a derivative); t > 0 is the duration of the exchange; and x ≥ 0 is the depth of the profile into the sample from the exchange surface.
The background isotopic fraction and that of the enriched exchange gas, C bg and C gas respectively, are required for the normalisation in eq. 2. System 3 uses the notation ∂ t to represent a partial derivative with respect to t.

Back-exchange
The back-exchange technique presented in this study starts with a conventional exchange in an 18 O enriched gas, following the protocol described in 3 , but then includes a second step where the sample is exposed to any unenriched gas environment. In order to draw a meaningful comparison, the temperature and partial pressure of oxygen should be the same for all exchange steps. To the authors' knowledge, there is no analytical solution available in the literature for this kind of two-step exchange process; however, a solution has been derived here for the special case where the values of D * and k * are the same for both exchanges.

Analytical solution
Many diffusion problems are amenable to solution by superposition, but due to the exchange type (Robin) boundary condition, this approach was non-trivial. The system describing the second exchange step is contained in sys. 4, where C Crank (φ 1 , x) is the result of eq. 1, keeping the values of D * and k * the same for the first and second exchange steps, after φ 1 seconds. Although variable t is used for exchange time in eq. 1, the non-ideal nature of real exchange experiments (i.e. finite heating and cooling rates) requires an "effective exchange time" to be calculated from temperature time series data. The variable φ i is used to describe this adjusted time, in exchange step i, and is based on the linearisation described by Killoran 6 . In the following discussion, time t is assumed to start at the beginning of the final exchange step, with the results of previous exchanges effecting only the initial conditions.
System 4, has two key differences from the single step system. Firstly, by considering the system time to start at the beginning of the second exchange, the initial condition is the Crank solution resulting from the first exchange, rather than a uniform distribution at background level. Secondly, the boundary condition no longer has the C gas stimulation and now simply sets the surface gradient proportional to the current surface concentration (Neumann boundary condition).
The existence and uniqueness of a weak solution, in the sense of distributions, to this system follows from the Lax-Milgram Theorem once it has been converted into Laplace space 7 . The solution was found by considering the following two physical constraints. Firstly, that the diffusion length, which is the convenient measure of the degree of penetration of a profile, must be approximately 2 D * (t + φ 1 ); secondly, that C must be positive for all values of t and φ 1 . The solution found is given by, which can be expressed more succinctly as the difference between two Crank solutions: By simple differentiation and substitution, it can be shown that eq. 5 is the unique solution to sys. 4. Furthermore, this approach can be extended to n exchanges, eq. 7, where D * and k * are constant; φ j is the duration of exchange j; and the normalised gas concentration at the boundary alternates between 1 and 0 at each step.
Although useful for validation purposes, these analytical solutions are only suitable for fitting profiles where D * and k * are the same in each exchange step. The following section describes a  Fig. 1 The space of possible back-exchange profiles under the constraint that D * was constant between the two exchanges. At each exchange time ratio, θ , a family of profiles was generated by varying the value of λ , which is the ratio of k * 2 to k * 1 . This graph also shows the sigmoid nature of the maximum surface value of the isotopic fraction as well as the maximum depth of the profiles' peak at each value of θ . numerical scheme constructed for fitting all other cases.

Numerical Methods
In order to generate profiles to use for fitting the experimental back-exchange data with multiple D * and k * values, a 1D Crank-Nicolson (CN) 8 type finite-difference simulation was used to numerically solve the system. The CN approach was chosen for its numerical stability and O(∆x 2 , ∆t 2 ) accuracy. For these line-scan type datasets, the spatial step, ∆x, was fixed to correspond to the experimental data raster spacing, while the time step, ∆t, was adjusted to preserve accuracy. The ghost node concept, described in more detail in 9 , was employed at both ends of the simulation domain. A ghost node is a fictitious extra node located outside the boundary of the system, which is not evaluated directly or stored in memory, but instead used to preserve the boundary condition by extrapolation, which maintains the O(∆x 2 , ∆t 2 ) accuracy throughout the domain.
To ensure that the semi-infinite boundary was modelled accurately, a hybridised boundary condition was utilised, which involved running the simulation twice (once with a mirror boundary (Neumann, ∂ x C = 0) and once with a fixed boundary (Dirichlet, C = 0)) and taking the mean (validated using the analytical form described above). However, as all the profiles measured in this study are very close to C = 0 at the edge of the analysis area, the choice of this boundary condition was of minimal significance.
The simulation was used to explore the space of possible profiles, fig. 1, under the constraint that D * was constant between the two exchanges, which is a property we expect from the systems investigated here, i.e. constant temperature and pO 2 . Seven families of curves were generated at equally spaced values of the exchange time ratio, θ , defined as, where φ 2 is the effective exchange time for exchange step 2.
Each family was produced by varying λ , defined as, where k * 1 and k * 2 are the values of k * for first and second exchange steps respectively. The highest profile in each family occurs when λ = 0, which is the case where there is no exchange between the solid and the gas phase during the second exchange step. The lowest profile occurs as λ → ∞, where the surface exchange is relatively very fast in the second exchange step and the process becomes diffusion limited. At both extremes, the families of curves exhibit asymptotic behaviour, which is discussed further in the context of uncertainty quantification later in this work.
In order to allow fig. 1 to generalise to all possible backexchange scenarios (with D * 2 = D * 1 ), the isotopic fraction required an additional normalisation step defined by, and the depth, x was also normalised with respect the to the standard diffusion length.
All the required normalisation equations, described above, have been inset into fig. 1 for convenience.

Fitting
A Nelder-Mead simplex algorithm 10 , which is a non-linear, least squares approach, implemented in the MatLab fminsearch function 11 , was used to fit the CN generated profiles to the experimental isotopic fraction data. The curves were fitted within a region of interest (ROI) between the surface and x = 2 for all profiles, to ensure that the reported values of the coefficient of determination, R 2 , were comparable (N.B. using data beyond x = 2 for fitting will force R 2 to be misleadingly high for this type of problem). This required a certain degree of iteration of the fitting process as the new values of D * altered the derived location of x = 2.
In order to efficiently fit the back-exchange data with the CN derived profiles, it is beneficial to select initial values of D * and k * with the correct order of magnitude. This is achieved by first fitting the data with the relevant analytical solution described earlier in this section and then using the results of this much faster method to initialise a second fitting using the CN profiles.

Experimental
The exchange experiments were conducted using the protocol described by DeSouza et al. for oxygen isotopic exchange in dense ceramic pellets 3 . LCSF powder was pressed and sintered at 1250 • C for 8 h (using a heating rate of 300 • C h −1 ) into cuboid pellets with dimensions 2 mm × 2 mm × 5 mm (H×W×D). This geometry allowed us to assume a 1D system when modelling, as the expected depth of x = 2 (where C < 0.005 for all k * ) was <500 µm. The pellets were then ground and polished to a mirror finish on one face using sequentially finer grades of polishing media. The grinding process was carried out using P800 grade (c. 22 µm particle size), P1200 (c. 15 µm particle size) and P2400 (c. 10 µm particle size) silicon carbide paper with water as a lubricant. After the grinding process, the samples were then polished sequentially using water based diamond suspensions with particle sizes of 6 µm, 3 µm, 1 µm and 0.25 µm. Care was taken not to skip grades to ensure not only that the surface was smooth, but also that the mechanical damage induced by each step was removed by the following step.
Before the first exchange, each pellet was annealed in research grade 99.999% pure dry oxygen ( 16 O) at the desired experimental temperature and oxygen partial pressure. By doing so, the oxygen vacancy concentration of the sample is fixed and ensures chemical equilibrium with the following exchanges. The pre-anneal time is calculated to be ten times that of the first exchange to provide a sufficiently large region with a constant oxygen nonstoichiometry. In order to determine accurately the effect of the atmosphere in the second exchange step, strict precautions were taken to ensure that the pre-anneal step and first exchange were undertaken in a dry, very high purity oxygen atmospheres. All experiments were performed in a bakeable ultra high vacuum exchange apparatus fabricated from stainless components and sealed with "conflat" flanges. Before each experimental step, the exchange volume was evacuated to a pressure of < 4 × 10 −7 mbar using a vacuum system employing a turbomolecular pump backed by a dry membrane pump to reduce the chance of back-streamed oil vapour. The 18 O enriched oxygen for the exchanges was derived from a 5 Å molecular sieve reservoir, which ensured that the gas was dry (c. 1 vppm water). The samples were not re-polished between exchanges.
Three experiments, illustrated in fig. 2, were performed using an initial static volume partial pressure of 200 mbar and a nominal annealing temperature of 785 • C. Experiments A and B were used to validate the back-exchange approach; experiment C was used to demonstrate its utility by investigating the effect of ambient air  Due to the high oxygen diffusivity of LSCF6428 it was necessary to prepare the sample for line-scan SIMS profiling 2 after the final exchange. To do so, the same grinding and polishing regime was applied to a face perpendicular to the exchange surface, prior to imaging. The sample was polished to a depth greater than a predicted value of x = 4 to ensure analysis occurs beyond the region where the diffusion front from this face will affect the isotopic ratio.
The exposed cross-section was then imaged using a ToF-SIMS IV machine (ION-TOF GmbH, Münster, Germany), as illustrated in fig. 3. The technique collects a full mass spectrum at each pixel and rasters over the surface to generate a 2D image. A burst analysis mode was utilised to measure the 16 O/ 18 O ratio due to high oxygen sputter yields 3 . This mode prevents detector saturation by using a series of short primary ion pulses rather than a larger, single ion pulse.
The value of C bg was measured from an unexchanged, but similarly prepared sample using the same imaging method and C gas was measured during the exchange using the silicon wafer technique described in 12 .

Data processing
Before the 2D images were collapsed into 1D line-scans, some data processing was required to correct for several established sources of error. A Poisson correction is applied to each mass peak to correct for counts missed due to the "dead time" of the detector, as discussed in 13 . This correction has a very significant effect in the detection of low isotopic fractions, which is relevant to the accurate measurement of the background concentration of 18 O (C bg ≈ 0.2% = 1/500). Bilinear interpolation was used to produce an image rotation series at 0.1 • increments. After generating a profile at each step, it was possible to align the exchange surface by choosing the angle with the maximum gradient of total counts ( 16 O+ 18 O), this same process was also used to find the location of the "surface". Typically, however, the degree of misalignment was small (< ±2 • ), which has only a modest effect on the fitted values of D * and k * as long as the surface location is specified as the mean surface, rather than the leading point. Figure 4 shows the mean oxygen counts ( 16 O+ 18 O) against distance from the surface, superimposed onto the 2D ToF-SIM image for sample C, which required 1 • of clockwise rotational alignment.
Following alignment, the isotopic fraction was calculated in each pixel and then the mean of each column was used to construct a 1D depth profile. Figure 5 shows the isotopic fraction ( 18 O/( 16 O+ 18 O)) against distance from the surface, superimposed onto the normalised 2D ToF-SIM data.
For each exchange, the time-temperature profile inside the furnace was recorded using a shielded thermocouple. The nominal duration of the first and second exchange steps were 1.5 h and 0.75 h respectively. By analysing the thermocouple data, the effective temperature was set to 785 • C and, using the method described by Killoran 6 , the effective exchange times were calculated to be, φ 1 = 1.4 h and φ 2 = 0.66 h.
In order to make this correction, the Killoran method requires an activation energy of diffusion for the LSCF6428. The value used here was E a =1.917 eV and although this value was itself derived from a similar experiment 14 , a sensitivity study suggests the resulting effective anneal times were highly insensitive to its value (<1 min eV −1 ). A separate activation energy is associated with the surface exchange coefficient, which is more difficult to account for as the surface exchange boundary condition is a function of both D * and k * . However, the rapid heating/cooling rates (c. 170 • C min −1 ) combined with suitable anneal times (>40 min in all cases) should minimise any resulting discrepancies. Figure 6 shows the profiles generated from the line-scans for experiments A-C. In each case the vertical axis was normalised using eq. 2.

Results
The values of k * 1 extracted from the fits was essentially identical to within reasonable experimental accuracy for all three experiments, as can be seen in Table 1. All three fits had values of R 2 above 0.999 measured for the region between the surface and x = 2. The value of D * was assumed to be constant for the two exchange steps.

B .
A . Fig. 6 Graph showing the measured isotopic profiles for the three samples (A, B and C), as a function of depth from the sample surface. Table 1 also shows that k * 2 of sample B, was identical to within experimental accuracy to k * 1 of sample B; whereas k * 2 of sample C was c. 250% larger than k * 1 of sample C. Although it is possible to approximate the uncertainty associated with the fitting process, it is expected that this will be small in comparison to other sources, such as sample to sample variation. As such, no uncertainty values are reported. The surface position and the size of the raster step are determined through calibration and contribute both a systematic and a statistical error.

Discussion
The results of the three experiments appear to show the validity and the utility of the back-exchange approach both through the consistency of k * 1 across all samples, as well as the significant variation between the two values of k * 2 . In all cases, the surface exchange coefficient in pure, dry, isotopically enriched O 2 was in excellent agreement both between the experiments and with the literature 15 . The surface exchange coefficient for the second exchange in experiment B was also in excellent agreement (to within experimental accuracy) with all values of k * 1 , which validates the back-exchange procedure. In all of the dry oxygen exchanges, the oxygen used was very high purity (99.999% pure) and measured to be very dry (< 1vpm H 2 O), which we believe to be crucial in the consistency of the result.
There are several possible mechanism for the augmentation of k * in the presence of air compared to dry oxygen. Ambient air contains a wide range of oxygen containing components and it is known that H 2 O and CO 2 16-18 can both significantly influence the overall oxygen exchange rate. Literature data suggests that both the temperature and overall composition of the anneal environment effects the measured oxygen exchange rate. It was beyond the scope of this paper to quantify the direct effect of each gaseous component; however, it should be clear that the backexchange technique is capable of making these measurements in a systematic study.

Experimental Error
As can be seen in Table 1 the k * 1 values have a small degree of variation. The source of this error may potentially be due to microstructural differences between samples. For this experiment, each sample was cut from the same sintered pellet and so assumed to be identical; however, subsurface defects (such as pores) could not be measured. Closed porosity will affect the local transport properties and thus adjust the measured diffusion profile.
For a number of common perovskite-based MIEC materials, exsolution of A-site cations to form secondary phase surface particles is a common occurrence [19][20][21][22] . These particles are typically non-conductive and thus reduce the catalytically active area for the reduction and incorporation of oxygen. It has been shown that significant 'strontium segregation' can occur in LSCF6428 at relatively low temperatures and short anneal times 20,22 . However, high temperature environmental scanning electron microscopy was used to inspect the surface of the LSCF6428 before and after annealing and no surface particles were observed.
The k * 1 value of sample B and C show excellent accordance, however, sample A appears to be lower. As mentioned previously, the statistical dead-time corrections are an additional source of error. Using the line-scan method and ToF-SIMS analysis, a constant total ion count is necessary to ensure accurate isotopic fraction measurement; however, the geometry and surface roughness of the sample can alter secondary ion extraction rates. Polishing the surface to a 0.25 µm grade is sufficient to negate these errors across the bulk surface, but at the edge of the sample it is very difficult to prevent curvature during to the polishing process. The radius of curvature (and thus number of affected pixels in the 2D ion image) determines surface uncertainty of the profile. In this work samples were adhered together prior to the second polishing step using a commercially available epoxy resin. This minimised the curvature to approximately 1 µm to 2 µm, greatly reducing the uncertainty of the line-scan profile near the surface. Figure 7 is essentially a section through fig. 1 at the time ratio measured in experiment C (θ = 0.47), although the vertical axis was left as C rather than C * . A family of profiles were generate by varying λ by a factor of 2 between each.

Profile uncertainty
The red curve highlights the case where λ = 1 (i.e. k * 2 = k * 1 ), which was close to the measured result of experiment B. The green curve is where k * 2 = 2.5 × k * 1 , similar to experiment C. For very small or very large value of λ , the profiles approach an asymptote, which would undermine the techniques ability to accurately measure the corresponding values of k * 2 . When D * 1 and D * 2 were allowed to float independently during the fitting, a family of similar curves could be generated, which all closely fit the profile. However, all samples tested were sections cut from the same pellet and also exchanged for identical times and temperatures, thus making it reasonable to assume that the profiles generated in the first exchange was the same for experiments A-C. This suggests that, for the back-exchange fitting, profiles where D * 1 and k * 1 are close to that resulting from experi- ment A should be preferred.
All of the data analysis, curve fitting and graphs were generated using a custom software package, TraceX 23 , which is a GUI developed using the MatLab guide platform. The latest version of TraceX is freely available from the supplementary materials section of this publication and those interested are encouraged to contact the author.

Conclusions
A novel back-exchange technique has been developed to quantify the effect of contaminants on the effective surface exchange coefficient. This involves a two-step exchange process, followed by line-scan SIMS imaging of the isotopic fraction. A new analytical solution was derived for a special case of this system and a Crank-Nicholson type finite difference simulation was used where it could not be applied (i.e. where k * 2 = k * 1 ). The technique was validated by comparison with values derived from a conventional exchange profile as well as values from the literature. This work found that ambient air had caused a c. 250% increase in the value of k * compared to pure, dry O 2 , which demonstrates the significance of the technique. The back-exchange approach can also be extended to any atmosphere with the same pO 2 as that of the first exchange step; however, care must be taken when approaching the asymptotic regions of the profiles as uncertainty can become large and asymmetric. Further work must be done to isolate the effect of each gas component in air.