Voltage-Relaxation GITT and Reverse Monte Carlo to Determine Lithium Diffusion and Distribution in TiO₂ and Highly-Ordered Nanoporous Hard Carbons

In this study two-electrode coin cells were used for the electrochemical analyses. All carbon materials and polyvinylidene fluoride (PVDF) (301F) were baked using custom built Sievert type apparatus at 300°C and at 110°C respectively overnight. The purpose was to remove moisture from the electrode materials. The whole assembly process was done inside an argon filled glove box with oxygen content below 1 ppm. Four different slurries were made using synthesized carbon samples as active material, PVDF as binder and N-methyl-2-pyrrolidone (NMP) as solvent. The ratio between active material and binder was 9:1. No carbon black was added in the process. Electrodes were prepared by coating thin layer of the slurry on a battery grade Cu foil. These electrodes were baked at 70°C for 10 minutes and then at 120°C for overnight inside the glove box. Coin cells were assembled with metallic lithium chip as reference/counter electrode and a carbon electrode as the working electrode. 1 M LiPF₆ in a mixture of 1:1 ethylene carbonate (EC) and diethylene carbonate (DEC) was used as the electrolyte. Two layers of Celgard brand separators were sandwiched between Li metallic chip and working carbon electrode which was eventually sealed under pressure of 50 kg/cm² using hydraulic crimping machine (MTI; MSK 110). All cells were rested for two hours to saturate before taking any measurements.

First, cells were cycled for 7 cycles between 0.00 V and 3.00 V vs Li/Li⁺ with constant current density of 30 mA/g using Maccor 4300 model to avoid effects of solid electrolyte interphase formation (SEI) for the diffusion measurements. Then all the cells were rested for one day and subsequently subjected to galvanostatic intermittent titration technique (GITT). This method is composed of constant current pluses followed by a relaxation period. In the experiment, 15 minutes current pulses with an approximately C/2 current and 4 hours of relaxation were used. The C value was calculated based on the cycling measurement results for each battery. More details can be found in the supporting information. The GITT process was continued until the carbon electrode reached a voltage of zero with respect to Li/Li⁺. Before the de-insertion current pulses began, the sample was allowed to soak at zero volts for 10 hours to ensure full insertion. De-insertion continued with current pulses of C/2 for 15 minutes and 4 hour relaxation period. All the voltage vs time measurements were recorded at room temperature using a Maccor model 4300 instrument.

Powder X-ray diffraction data were fitted using fityk⁴⁹ software in order to obtain Bragg peak widths from the experimental results. All the recovery voltage curves were fitted using xmgrace software⁵⁰ to extract the VR-GITT parameters, α and τ, as discussed below.

According to the BET results in Figure 1a, all of the carbon types except GNPC exhibit a type IV isotherm based on the IUPAC classification.⁵¹ This indicates the presence of mesoporosity and microporosity in the samples resulting in large pore volume and surface area. GNPC has a type I isotherm, exhibiting very low surface area and a micro-porous nature in contrast to the other samples. GNPC was made with phenolic resins without polymer (F127), it therefore has no ordered morphology at the nanoscale, yet retains a disordered carbon structure.⁵² Figure 1b shows the pore size distribution of the samples from NLDFT models. NLDFT method provides a better understanding about the lower pore size regions (<2 nm) compared to the BJH method.⁵³ Nonetheless, poor fitting in the micropore region caused higher standard deviations (Table S1) for CMK3 and NPC-1 compared to NPC-2. We use 2D Cylindrical model for CMK3 even though it has more slit like pores as evident from the reported H3 type hysteresis loop.³⁶ This is probably the reason we obtain higher standard deviations for CMK3 carbon. Pore size distributions from the BJH method are also included in supplementary documents for comparison (Figure S2). All of the carbon types except GNPC contain micropores and mesopores according to the pore size distributions. Summary of N₂ Physisorption results are given in Table I. The SEM images (Figure S3) show the particle size variation between different carbon samples. The particle size of the GNPC was assumed to be similar to NPC-1 since these samples were sieved through the same 74 μm size sieve. The largest particle size for the carbons CMK3, NPC-1, NPC-2, and GNPC are shown in Table I.

Zoom In Zoom Out Reset image size

Figure 2 shows wide-angle XRD patterns for all the carbon samples. These XRD patterns have been used in Scherrer equations to calculate the (002) interlayer spacing and coherence length which are shown in Table I. The broad peaks corresponding to graphitic (002) and (100) at around 26° and 43° two theta, respectively, indicates no significant long-range order; the carbon network itself is glassy, with very small graphitic nanodomains. However, the GNPC structure shows higher intensity and somewhat narrow peaks with slightly lower d spacing and higher coherence lengths compared to the other carbon samples. This suggests that the GNPC structure may be slightly more graphitic on the nanoscale, compared to the other ordered carbons.

Zoom In Zoom Out Reset image size

In addition to analyzing the voltage response during the current pulse, it is useful to examine the relaxation of the open cell voltage (OCV) following current pulses. The depolarization technique, developed for planar interfaces,^54,55 can be used to analyze the OCV relaxation curve in order to get diffusion coefficients in thick electrode pellets. In the depolarization method, the diffusion coefficient is calculated from the gradient obtained from a plot of the logarithm of [OCV(t) - OCV(t = ∞)] vs time. In contrast we employ a stretched exponential to analyze the relaxation curves and refer to our approach as voltage relaxation (VR-GITT), and we calculate diffusion constants assuming a spherical particle geometry. It is important to note that during the diffusive relaxation process the dominant contribution comes from the largest particles in the material. SEM indicates these sizes are in the micron range.

It is known that the stretched exponential is more appropriate to model relaxation and decay processes than the standard exponential function.^56–58 Prominent examples include polymer dynamics, and diffusive processes such as rehydration. The above mentioned different systems can be explained with stretched exponential function, exp ⌊ − (t − t₀/τ)^α⌋, having either fixed or variable alpha (α) between zero and one and tau (τ) larger than zero depending on the system.⁵⁶ In diffusive processes the α parameter is shape-related and is important behavior to understand in our material where the pores and channels are ordered. Marabi et al. has shown that the stretched exponential shape factor α can have distinct values depending on the geometry and the mechanism of water uptake such as 0.67, 0.72 0.81 for spheres, cylinders and slabs respectively when considering rehydration.⁵⁹

A simple stretched exponential for a typical voltage relaxation in our electrochemical data is captured by the expression in Eq. 1

Where, V_max is the relaxed open cell voltage obtained from the fitting, V_min is the minimum voltage of the recovery curve, (t−t₀) is the relaxation time (difference between initial and final relaxation time), τ is the time constant, and α is the shape-related factor. Fitting of the stretched exponential to a typical VR curve is shown in Figure S4. The diffusion constants following every current pulse were calculated using

where L is the radius of the largest particle. To validate the diffusion constants calculated from our extracted time constants, diffusion coefficients were also calculated based on the standard GITT equation which is based on a one-dimensional diffusion process:⁶⁰

Here r is the radius of the particle, t is the duration of the current pulse, ΔE_s is the change in equilibrium voltage during the current pulse, ΔE_t is the total transient voltage change after the IR drop. This equation holds when t ≪ (r²/D). The voltages of the experimental GITT curves are shown in Figure S5. The VR-GITT data show that NPC-1 has the lowest recovery voltage and fastest relaxation process. This suggests that NPC-1 may not have as deep trap sites discussed in Mochida et al. as compared to CMK3 and NPC-2.¹² A possible explanation may be the following. When some portion of Li ions intercalate into deep trap sites on the first insertion they will have diffusion barrier activation energies larger than the other intercalated Li ions during the voltage recovery. Therefore, the recovery voltages may provide some information about the trap sites. Li resident in a deep trap site will increase the diffusion rate for diffusing Li that cannot enter the site. A competing mechanism, however involves the blocking of a diffusive path by the resident Li. These competing effects may be difficult to determine from the VR-GITT data.

We note one drawback of the standard GITT technique is the ΔE_s term in the numerator of Equation 3 which represents the change in equilibrium voltage of the material upon, e.g. insertion of Li. In a two-phase region where the equilibrium voltage does not change ΔE_s would be zero, the diffusion constant may not be calculated in this manner. See for example Zhu and Wang.⁶¹ However, if one uses the voltage relaxation curve and extracts the time constant, a calculation of the diffusion constant is possible with knowledge of the relevant length scale for the diffusion. As it is sometimes difficult to measure diffusion constants, especially over a large concentration region, the VR-GITT technique provides one more characterization tool. For these purposes, and also to provide a more well-studied standard for our analysis that follows, we show the VR-GITT technique applied to nanoparticle anatase TiO₂. This material has been investigated for lithium-ion anodes and its nanoparticle phase diagram is understood as a function of Li concentration.^40–46 We present diffusion measurement results across the entire two-phase region of anatase TiO₂ + beta Li titanate during Li insertion and de-insertion. The VR-GITT results for anatase TiO₂ demonstrate the ease of use and power of the method.

In the Li-Ti-O system there are three phases of relevance for lithium-ion battery applications: (1) alpha-phase anatase with stoichiometry Li_δTiO₂, (2) beta lithium titanate, Li_0.5Ti₁O₂, and (3) gamma phase Li₁Ti₁O₂. These phases and the nanoparticle phase diagram have been investigated previously.^40–46 (The preparation of our nanoparticle anatase TiO₂ is described in the supplementary.)

The VR-GITT results for 35 nm anatase TiO₂ is shown in Figure 3. Figure 3a shows the fitted final voltages from the stretched exponential yielding the plateau voltage. First, it is clear that there is a plateau region for the 35 nm particles. Insertion of lithium into anatase TiO₂ proceeds via a two-phase region corresponding to the transformation of alpha-TiO₂ to the beta lithium titanate phase. The length of the plateau becomes shorter with smaller particle size (See the supplementary Figure S6 for 11 nm particle size), in agreement with the reported nanoparticle phase diagram in the literature.⁴⁰ The drop in potential near full insertion corresponds to the formation of the gamma phase with 1:1 ratio of Li:Ti. The ionic conductivity in this phase is lower, and there is a corresponding change in the extracted VR-GITT parameters.

Zoom In Zoom Out Reset image size

The extracted tau parameters and corresponding diffusion coefficients are shown in Figures 3b, and 3d. Most interestingly, small changes in the extracted tau parameter indicate a diffusion coefficient that increases as the concentration of lithium increases during the phase transition from alpha to beta. The reason that the diffusion coefficient of Li⁺ increases is the following.⁴² Density functional theory calculations show that at low lithium concentration in anatase TiO₂ there is strong localization of the electron e⁻ from the intercalated Li to a 3d state of the closest Ti center, forming a bound polaronic state of the Li⁺ to the Ti center. This binding is roughly 56 meV. As the lithium concentration increases in the anatase TiO₂ and these gap states begin to fill and form a delocalized band, the binding of Li⁺ to these states falls, manifesting as a higher diffusion rate of Li⁺. VR-GITT measurements of lithium intercalation into anatase TiO₂ illustrate the above explanation in detail. Further, Figure 3c shows the extracted shape parameters from the VR-GITT fits and indicate that the diffusion proceeds largely with the same character (spherical) except at the point where the formation of the gamma Li₁Ti₁ phase appears. There the alpha parameter increases noticeably.

Figure 4 contains the extracted time constants and α parameters vs Li concentration using Equation 1 for different carbon samples. Voltages that relax sharply at the beginning of the relaxation are given by smaller values of α. The extracted α values generally lie between 0.2–0.5. As evident in Figure 4a the α shape parameter clearly distinguishes carbons synthesized with and without silicate in the preparation process. The shape parameter data fall into two distinct groups. A strong variation of α for carbons NPC-1 and GNPC without silicate and a non-varying α for carbons NPC-2 and CMK3 prepared with silicate. The time constants are generally increasing and shown in Figure 4b, with the exception of the one non-porous carbon GNPC, which shows a rapid drop in τ. The low capacity in these measurements is due to large current density used during the VR-GITT pulse.

Zoom In Zoom Out Reset image size

Figure 5 shows the diffusion coefficients for carbon samples from both the VR-GITT time constant (Equation 2) and the standard GITT (Equation 3). These results show very good agreement with each other. All the obtained diffusion coefficients from VR-GITT and GITT are in the range of 10⁻⁷ and 10⁻⁹ cm²s. This shows an agreement with the range of reported diffusion coefficients for carbonaceous materials.^39,62–64 However, due to the lack of the Li diffusion measurements for the exact materials discussed in this study, we cannot compare the values directly.

Zoom In Zoom Out Reset image size

GNPC shows the highest diffusion coefficients. This may be due to more graphitic nature of the GNPC sample. Persson et al. have shown that Li diffusion values in highly oriented pyrolytic graphite (HOPG) in the direction parallel to graphene layers is about ∼10⁻⁷–10⁻⁶ cm²s⁻¹ whereas the diffusion along the grain boundaries is about ∼10⁻¹¹ cm²s⁻¹ through computational calculations.⁶⁵ This suggests if our GNPC has more graphitic nano-domains it should have a higher diffusion coefficient.

Li plating is a well-known phenomenon which occurs specially under low temperatures, high current pulse and high State of Charge (SOC). There have been several studies reported in literature about the investigation of Li deposition inside electrode material under different conditions. One common finding is that uneven Li plating leads to form dendrites and may even pierces the separator causing internal short circuit of the battery.^66–72 Li plating inside the electrode material without leading to dendritic growth depending on the morphology of the electrode material is also found to be possible. This phenomenon can lead to very high reversible capacities for the materials. In these cases, formation of metallic Li occurs inside the interior of the pores or in the cavities/defect sites of the material and it will not lead to form destructive dendrites on the electrode surface.^73–76

The lithium capacities in excess of LiC₆ in the cycling data strongly suggest that lithium must be plating inside the pores of the carbon. Interestingly, NPC-2 type carbon with very high specific area and pore volume can give reversible capacities above 1000 mAh/g which is more than two times of the standard graphite capacity (Figure S7). One hypothesis is that the lithium plates during the relatively high current pulse on insertion in the GITT measurement, and then diffuses through the carbon slowly. If this is the case, we expect a non-uniform distribution of lithium in the carbon particles during insertion. After full insertion and during the de-insertion VR-GITT measurements, we expect that this non-uniformity will be smaller.

In order to investigate the distribution profile of lithium in the carbon particles we employed computational modeling of the diffusion process during the voltage relaxation on each insertion and de-insertion. We extracted density profiles numerically by integrating the diffusion equation and using Monte Carlo methods to determine the lithium density profile following every VR-GITT current pulse. The details are as follows.

The outer edge of the particle that is in contact with the electrolyte is the surface that provides the measured voltage response. Our objective was to use this voltage response to calculate the concentration profile of lithium in the particles immediately following the current pulse in the VR-GITT experiments. We deduced the initial lithium profile by assuming lithium was incorporated into the (mostly surface region) of the carbon particles during the (short) current pulse, and that this high-concentration diffused into the carbon particles over the relaxation period of the measurement (Figure S8). The resulting voltage response is given by the concentration of the lithium on the outer edge of the particle (Figure S9). In order to make the computations tractable we assumed spherical particles. To obtain the time evolution of the lithium concentration we used a standard Crank-Nicolson integration routine with Dufort-Frankel near r = 0. Our computations allowed for a diffusion constant that varied as a function of distance from the center of the particle, D(r). This assumption allows for the difference in diffusion rates along the pore channels and/or through the carbon matrix itself. The simplification that D is only a function of r and not a complicated function of the carbon pore morphology is a reasonable one on based on the length scale ratio of the pore size (10⁻⁹ m range) to the particle size (10⁻⁶ m range). The carbon particles are not monolithic ordered structures, but are comprised of small domains of ordered pores. The size of the domains is much smaller than the micrometer size of the largest carbon particles.

For each voltage response curve, collected during the VR-GITT measurement, we determined the initial lithium density profile in the particle using standard Monte Carlo simulated annealing (MCSA) to fit the experimental voltage response. The cost function was therefore the difference between the experimental voltage relaxation and the computed voltage response, given an initial density profile. An example MCSA computation is shown in Figure S10. For all computations the integrated density was normalized to unity and the resulting density profiles f(r), and diffusion constant profiles D(r) are dimensionless.

We used N = 20 spherically concentric shells or bins in a particle for the simulations. In each shell the density and diffusion constant profiles were allowed to vary according to three different configuration changes, (i) scaling of the overall profile, f→ (1 +/− δ)^*f (ii) increase or decrease of profile, f→ (f +/− δ), (iii) increase or decrease profile of single bin, f_i → (f_i +/− δ). An exponential annealing schedule was used with T(n) = exp (−n/n₀), with n₀ = 2, and a total of 10 temperature steps. The number of configuration changes for each temperature was N_c = 10³. The value of δ for each of the configuration changes was chosen so that rejection percentages were in the range of 65% to 75%.

The relaxation curve at each VR-GITT data point represents the diffusion of lithium in the particles at either an increasing or decreasing total Li concentration if the measurement was taken as an insertion or de-insertion. The initial lithium density profiles f(r), are plotted as a function of the cycle number and the particle radius in Figure 6. Peaks appearing along the density profiles f(r) vs radius graph denote high Li concentration positions in the particle. In our model we consider outer surface as the surface where electrolyte makes contact with the carbon particle. This includes both the outer edge of the particle and the inner surface of the nanopore channel. High Li concentration points appearing along the radius of the particle should be due to agglomerations/plating of Li inside the nanopore channels. In previous studies Li plating in defects and cavities have been observed for graphene and mesopore carbons.^74,75 Moreover, intercalation of plated Li into graphite during relaxation period have been reported.⁷⁷ Uhlmann et al. have reported that the shape of the voltage curve during charging and during relaxation for graphite in the presence of plating. However, we have no two-phase region and we did not observe any voltage plateau in voltage relaxation curves for our amorphous hard carbon materials.⁷⁷

Zoom In Zoom Out Reset image size

There are no reports in the literature of direct observation of Li plating or intercalation of plated Li for our carbons. However, plating and intercalation occur in graphite and graphene support the possibility of a similar phenomenon happening inside the amorphous carbons used in this paper. According to BET analysis our carbon samples have hexagonal nanopore structures. These pore channels are not perfectly smooth and have cavities and sharp protrusions along the pore channels. So, during high current pulse, it is likely that Li may plate in these areas and diffuse into the particle during the relaxation.

The results of the MCSA computations are revealing. First, we note that the variation of the inferred diffusion coefficients was small (within a factor of 5) across the VR-GITT measurement of a sample and a varying diffusion constant was not necessary to obtain the density profiles. Second, the density profile becomes flat as the concentration of lithium reaches its maximum value and returns to a more featured landscape upon de-insertion (Figure 6). A flat voltage profile indicates evenly distributed Li inside the pores. The zero voltage soaking step following the VR-GITT insertion pulses evidently increases the Li concentration and smooths its density profile. After the voltage soaking period, the density profiles stay flat until few de-insertion current pulses for GNPC and NPC-1. This may suggest that a higher amount of Li plating occurs in these samples compared to NPC-2 and CMK3. This is in contrast to a well-studied material such as anatase TiO₂ nanoparticles. An example of the VR-GITT applied to this material shows a very flat landscape of lithium concentration profile across insertion and de-insertion and only an increase in the density right at the edge of the particle upon full insertion (Figure S11).

To support our density profile analysis, we quantified the degree of fluctuation in the density plots. The results indicate that the computed lithium profiles of the carbon samples are significantly more non-uniform than the standard provided by the TiO₂ anatase nanoparticles. Figure 7 shows the standard deviation for the calculated density profiles at each step of the VR-GITT process, i.e., the fluctuations in the profile from r = 0 to r = 1, for each relaxation curve that was fitted using MCSA. The similarity between the samples prepared with silicate (CMK3, NPC-2) or without (NPC-1, GNPC) is again very apparent in the de-insertion standard deviations. A direct comparison of the carbon data with the standard TiO₂ is not possible because of the difference in particle size and due to the differences in the experimental procedure. In our TiO₂ experimental procedure, no soaking step was involved between insertion and de-insertion. However, the time scale for relaxation of all of the samples is of the order required by the approximate known diffusion constants for carbon and TiO₂, and so the standard deviation profiles do provide some indication of the initial lithium density profiles inside the particles.

Zoom In Zoom Out Reset image size

Our results show that the density profiles in the carbons become more uniform (smaller standard deviation) as the Li concentration increases and the most uniform following the zero voltage soaking step. This indicates that the Li is most evenly distributed at full insertion. A calculation of the pore volume of the materials and the Li concentration from the VR-GITT measurements shown in the SI indicate that if plating occurs in the pores, it fills no more than about 10–20% of the pore volume. Interestingly for the GNPC-type carbon, which has the lowest pore volume, this number is calculated to be around 100%. This suggests that if plating initiates on the 'sharp' edges of the pore walls during a high current pulse, it may evenly distribute itself over the surface of the pore wall and/or intercalate over time.