Migration Barrier Estimation of Carbon in Lead for Lead–Acid Battery Applications: A Density Functional Theory Approach

Recent efforts towards developing novel lead electrodes involving carbon and lead composites have shown potential for increasing the cycle life of lead–acid (LA) batteries used to store energy in various applications. In this study, first-principles calculations are used to examine the structural stability, defect formation energy, and migration barrier of C in Pb for LA batteries. Density functional theory with the GGA-PBE functional performed the best out of various functionals used for structural stability calculations. Furthermore, with the complete incorporation of C in the Pb matrix, the results show that C is energetically preferred to be at the octahedral interstitial (COcta i ) site in the FCC structure of Pb. Additionally, climbing-image nudged elastic band calculations show a minimum energy pathway for C diffusing from a stable octahedral site to the adjacent octahedral site assisted by a tetrahedral intermediate site. Therefore, the minimum energy pathway for C migration is envisioned to be COcta i → C Tetra i → C Octa i , where the total energy barrier is observed to be ~90% and more than 100% lower than the CTetra i → C Tetra i and C Octa i → C Octa i barriers, respectively.


Introduction
A transition of the global energy needs from fossil fuels to renewable resources has increased the need for batteries [1]. The research and development of batteries has increased over the past decade, and the battery market is expected to keep growing at a rate of 14.1% until 2027 [2,3]. The total market shares of different types of batteries are shown in Figure 1. Currently, the dominant battery types are lead-acid (LA), lithium-ion (Li-ion), nickel-metal hydride (Ni-MH), and nickel-cadmium (Ni-Cd) batteries [4,5]. Of all of these types, LA and Li-ion batteries are dominating the market. While Li-ion batteries are getting significant attention in electrical vehicles and electronic devices [6,7], the global LA battery demand is expected to keep growing due to their increased usage in automobiles, uninterrupted power sources (UPSs), and manufacturing methods [3,8,9]. The primary driving force for this sustained growth is attributed to the lower cost and overall recyclability (~99%) of LA batteries. Recent development in LA batteries has focused on reducing their total weight while maintaining/enhancing the overall battery performance [10,11] by adding different forms of carbon (C) to the individual electrodes [12][13][14][15]. The primary location of incorporation of C in the LA battery architecture is at the negative active material (NAM). Figure S1 provides a schematic of LA batteries, highlighting the positive active material (traditionally PbO 2 ) and NAM (traditionally Pb), separated by a separator. The NAM is shown to include carbon, where carbon is completely incorporated into the Pb matrix. In this work we are focusing on enhancing our fundamental understanding of the incorporation of C in a Pb matrix.
Solids 2022, 2, FOR PEER REVIEW 2 fundamental understanding of the incorporation of C in a Pb matrix. . Adapted from [16].
The addition of carbon to NAM mostly improves the battery performance [17][18][19][20], due to (1) increase in electronic conductivity, (2) restriction of lead sulfate (PbSO4) crystal growth during the charging and discharging cycles, and (3) enhancement of capacitive contribution, which improves the discharge performance [21]. Carbon is added to NAM in different forms (e.g., activated carbon [12,18,22,23], carbon black [13,22,24], graphite [23,25,26], graphene [26][27][28], carbon nanotubes [26,[29][30][31], porous carbon [32,33], etc.). While most of these carbons improve the cell performance, microporous carbon with a high surface area has been reported to reduce the capacity of the battery and, hence, is not suitable [20]. Recent studies focusing on the addition of nanostructured carbon to NAM have resulted in higher utilization efficiency and enhancement of charge/discharge performance in LA batteries [29]. There are different approaches to including carbon in NAM, where carbon is either incorporated in the Pb matrix during synthesis or it is adsorbed on the surface of the Pb electrode [22].
Both experimental and theoretical studies are critical to advancing our understanding of different battery technologies. Computational methods at the level of electronic and atomistic methods have been successfully used to enhance our understanding of battery materials [34][35][36][37][38][39][40]. While the impact of carbon on the performance of LA batteries is well established experimentally, it is important to enhance our fundamental understanding of the addition of C to the Pb matrix for further improvement. First-principles calculations are used in this study to address the stability of C in the Pb matrix. Section 2 describes the computational details used in this study. The results on the structural stability of C and diffusion barriers for the movement of C in the Pb matrix are included in Section 3. The conclusions of this work are summarized in Section 4.

Density Functional Theory
Density functional theory (DFT) at the level of plane waves was used for all of the calculations performed on the Pb and Pb-C systems in this work. The DFT calculations were carried out using the Vienna Ab initio Simulation Package (VASP) [41,42] available within the MedeA Platform [43]. The generalized gradient approximation (GGA) functional parametrized by Perdew, Burke, and Ernzerhof (PBE) was used to describe the exchange and correlation energies of the electrons [44]. The projector-augmented wave (PAW) pseudopotentials described the valence electron configurations of 6s 2 5d 10 6p 2 and The addition of carbon to NAM mostly improves the battery performance [17][18][19][20], due to (1) increase in electronic conductivity, (2) restriction of lead sulfate (PbSO 4 ) crystal growth during the charging and discharging cycles, and (3) enhancement of capacitive contribution, which improves the discharge performance [21]. Carbon is added to NAM in different forms (e.g., activated carbon [12,18,22,23], carbon black [13,22,24], graphite [23,25,26], graphene [26][27][28], carbon nanotubes [26,[29][30][31], porous carbon [32,33], etc.). While most of these carbons improve the cell performance, microporous carbon with a high surface area has been reported to reduce the capacity of the battery and, hence, is not suitable [20]. Recent studies focusing on the addition of nanostructured carbon to NAM have resulted in higher utilization efficiency and enhancement of charge/discharge performance in LA batteries [29]. There are different approaches to including carbon in NAM, where carbon is either incorporated in the Pb matrix during synthesis or it is adsorbed on the surface of the Pb electrode [22].
Both experimental and theoretical studies are critical to advancing our understanding of different battery technologies. Computational methods at the level of electronic and atomistic methods have been successfully used to enhance our understanding of battery materials [34][35][36][37][38][39][40]. While the impact of carbon on the performance of LA batteries is well established experimentally, it is important to enhance our fundamental understanding of the addition of C to the Pb matrix for further improvement. First-principles calculations are used in this study to address the stability of C in the Pb matrix. Section 2 describes the computational details used in this study. The results on the structural stability of C and diffusion barriers for the movement of C in the Pb matrix are included in Section 3. The conclusions of this work are summarized in Section 4.

Computational Details Density Functional Theory
Density functional theory (DFT) at the level of plane waves was used for all of the calculations performed on the Pb and Pb-C systems in this work. The DFT calculations were carried out using the Vienna Ab initio Simulation Package (VASP) [41,42] available within the MedeA Platform [43]. The generalized gradient approximation (GGA) functional parametrized by Perdew, Burke, and Ernzerhof (PBE) was used to describe the exchange and correlation energies of the electrons [44]. The projector-augmented wave (PAW) pseudopotentials described the valence electron configurations of 6s 2 5d 10 6p 2 and 2s 2 2p 2 for Pb and C, respectively. All of the calculations performed for this computational examination used a plane-wave cutoff energy of 500 eV. To accommodate varying C concentrations within the Pb matrix, different system sizes along the X-, Y-, and Z-axes were examined.
Due to the variation in supercell size, different k-point Monkhorst-Pack meshes were used for integrations over the Brillouin zone [45]. However, the actual k-spacing was maintained close to 0.25 × 0.25 × 0.25 per Å for all of the structures. For structural optimization, all single-point calculations were set with a convergence value of 1.0 × 10 −6 eV; for ionic relaxation and complete optimization (position, supercell shape, and size), the calculation was performed until the Hellmann-Feynman force on each atom was less than 0.02 eV/Å. Experimentally, the amount of C (C c ) added to NAM of LA batteries was reported to be in the range of 0.2-2.0 wt% [12,13,17,22]. As shown in Supplementary Section S1, assuming a two-element system of Pb and C, the wt%-at% conversion shows that the addition of one carbon to every unit cell of Pb (four Pb atoms for the FCC structure) results in a C c = 1.43 wt% C. By increasing the system size to 2 × 2 × 2 (8 unit cells or 32 Pb atoms), the addition of one C atom results in a C c = 0.18 wt% C, which is closer to the lower range of the experimental limit reported in the literature. Therefore, this work focused on adding one C to Pb supercells of 1 × 1 × 1, 2 × 1 × 1, 2 × 2 × 1, and 2 × 2 × 2 where the numbers refer to the number of unit cells along the X-, Y-, and Z-axes, respectively. The at% of C (C c ) for these cases varied between 20% and~3%, as shown in Table 1, as the supercell size increased from 1 × 1 × 1 to 2 × 2 × 2. There are two different systems simulated in this study: (1) a pure Pb system without any C, and (2) a Pb system with one C. The pure Pb system was simulated to validate the computational approach, while the system with C was simulated to understand the stability of C in the Pb unit cell, the defect formation energy, and the diffusion barrier for movement of C in the Pb matrix.

Validation and Selection of Functional for Pb
The crystal structure of Pb shows an Fm3m symmetry with an experimentally measured unit cell length of 4.9508 Å [46] ( Figure S2). Both LDA and GGA-PBE functionals have been used to examine crystalline Pb in the literature, where LDA calculations underestimated the lattice parameter while GGA calculations overestimated the lattice parameter [47][48][49][50]. Prior to examining the structural prediction by different functionals, the convergence tests on both the plane-wave cutoff energy and k-point mesh were performed with the GGA-PBE (details are included in Supplementary Section S2). After establishing convergence, a full optimization (atom position, cell shape, and cell volume) of the 1 × 1 × 1 supercell with a 9 × 9 × 9 k-point mesh (corresponding to 0.24 × 0.24 × 0.24 per Å actual k-spacing) was performed. Various density functionals (LDA [51,52], GGA-PBE [53], GGA-AM05 [54], GGA-PBESol [55], GGA-rPBE [56], and GGA-BLYP [57]) and hybrid functionals (PBE0 [58,59] and HSE06 [60] with GGA-PBE) were used to compute the ground state lattice parameter of Pb. Table 2 shows the output we obtained from these simulations. The LDA functional underestimated the lattice parameter by 1.4%. This is similar to the previously reported lattice parameters of Pb predicted using LDA and a tight binding [47] approach. Our GGA-PBE functional calculation overestimated the lattice parameter by 1.6%, which is consistent with the previous DFT predictions at the same level of theory [48][49][50]. A much closer prediction of structural details of Pb (lattice parameter within ± 0.5%) was obtained with the GGA-AM05, GGA-PBESol, and both of the hybrid functionals (HSE06 and PBE0), while the GGA-rPBE and GGA-BLYP functionals overestimated the lattice parameter by 3.0% and 4.5%, respectively. Experimental lattice parameters and volumes are taken from [61,62]. Along with the lattice parameter, elastic constants and bulk modulus values were evaluated for each of the potentials using 11 × 11 × 11 k-points with a strain of 0.5%. The mechanical property calculations were performed within the MedeA-MT module [43], where the atomic position in the strained lattice was optimized until the Hellmann-Feynman force on each atom was less than 0.01 eV/Å. Table 2 includes the computed elastic constants and bulk modulus values for each case. The elastic constants from experimental studies were reported to be 55.54 GPa, 45.42 GPa, and 19.42 GPa for C 11 , C 12 , and C 44 , respectively, at 0 K, while they were 47.65 GPa, 40.28 GPa, and 14.41 GPa for C 11 , C 12 , and C 44 , respectively, at room temperature (RT) [62]. The bulk modulus [(C 11 + 2C 12 )/3] for these temperatures was estimated to be 48.79 GPa and 42.73 GPa at 0 K and RT, respectively. From the individual elastic stiffness values evaluated for each functional, GGA-AM05 and the HSE06 and PBE0 hybrid functionals predicted C 44 to be negative; this indicates that the structure of Pb is unstable using these functionals. For the remaining cases, the LDA functional predicted the elastic constants as being very close to the experimental value at 0 K, where the errors in C 11 , C 12 , and C 44 were 4.5%, −1.7%, and −2.6%, respectively. The GGA-PBE functional predictions underestimated these values by 17.1-25.1%. A tighter range of underestimation was observed for GGA-PBESol-between 3.0 and 10.6%. Both GGA-rPBE and GGA-BLYP underestimated the elastic properties by between 26.4 and 47.1% in comparison to experimental values referring to 0 K.
Since most of the experimental measurements were performed at room temperature, we also compared our computational predictions against the experimental elastic constants at RT. The comparison shows that the LDA functional overestimated the elastic constants by 21.8%, 10.9%, and 31.3% for C 11 , C 12 , and C 44 , respectively. The GGA-PBE functional underestimated the C 11 and C 12 stiffness values by 3.4% and 15.5%, respectively, while it overestimated the C 44 constant by 8.2%. The GGA-PBESol functional predicted C 12 to be within 1% of the experimental value; however, it overestimated C 11 and C 44 by 13.1% and 17.4%, respectively. GGA-rPBE and GGA-BLYP underestimated the elastic constants much more compared to the GGA-PBE functional.
Therefore, by comparing the structural and mechanical properties estimated, we neglected GGA-rPBE and GGA-BLYP, as these two overestimated the lattice parameter the most compared to the other functionals. Based on the negative C 44 values, we eliminated GGA-AM05, HSE06, and PBE0 from further consideration. Out of LDA, GGA-PBE, and GGA-PBESol, GGA-PBESol predicted the lattice parameter closest to the experimental value, while the elastic constants were closer to the 0 K values. As our primary focus was to compare our predictions against experimental values at finite temperatures, we considered GGA-PBE for rest of the studies, as the GGA-PBE functional performed the best in predicting the mechanical properties at room temperature compared to the LDA and GGA-PBESol functionals.
After establishing the functional, we confirmed the metallic character of Pb by analyzing the band structure and density of states (DOS) ( Figure S4). A closer look at the partial DOS shows that the metallic characteristics are predominantly dominated by the outermost p-orbitals (valence electrons for Pb are taken as 5d 10 6s 2 6p 2 with the pseudopotential). After examining the pure Pb matrix, the next step was to incorporate C into the Pb matrix.

Carbon Placement in the Pb Matrix
Based on the discussion above, we added one C atom to the Pb matrix by systematically increasing the supercell size from 1 × 1 × 1 to 2 × 1 × 1, 2 × 2 × 1, and 2 × 2 × 2, representing a C concentration of C c = 1.43 wt% − 0.18 wt%, as shown in Table 1. This matches the experimental C range in Pb-C batteries, and represents a complete incorporation of C into the Pb matrix [22]. While the number of C atoms in the Pb matrix is kept in the experimental range, the placement of C in the Pb matrix is not clear from the experimental studies. Therefore, we performed a structural stability test by placing C at different interstitial and substitutional sites. The results on the structural stability of C are presented based on the 2 × 2 × 2 supercell size.
Initially, we examined the total energy of the system for the C atom sitting at an interstitial site, i.e., either at the octahedral site (C Octa i ) or at the tetrahedral site (C Tetra i ) in the Pb matrix ( Figure 2). A self-consistent field (SCF) calculation while keeping the cell shape and volume fixed resulted in C Octa i being more favorable than C Tetra i for the 2 × 2 × 2 supercell with one C atom: where E c Octa i is the total energy of the 2 × 2 × 2 Pb supercell with one C in the octahedral site, and E c Tetra i is the total energy of the 2 × 2 × 2 Pb supercell with one C in the tetrahedral site. The results shown in Equation (1) confirm that the C atoms prefer to sit at the octahedral site, which is similar to the placement of C in the γ-Fe matrix (FCC structure) [63]. After establishing the interstitial site stability for C atoms in the Pb matrix, we also performed a stability check between the interstitial and substitutional sites for C. In order to complete this comparison, we prepared the substitutional site by replacing one Pb atom randomly in the 2 × 2 × 2 supercell with a C atom. Therefore, there is a mismatch between the total number of atoms in the supercells for the interstitial (32 Pb atoms + 1 C atom) and substitutional (31 Pb atom + 1 C atom) cases. The comparison of energy requires the calculation of energy per Pb atom. The energy of each Pb atom ( = −3.533 eV/aom) was obtained from the 2 × 2 × 2 supercell without any C in the structure.
The two energies compared were as follows: After establishing the interstitial site stability for C atoms in the Pb matrix, we also performed a stability check between the interstitial and substitutional sites for C. In order to complete this comparison, we prepared the substitutional site by replacing one Pb atom randomly in the 2 × 2 × 2 supercell with a C atom. Therefore, there is a mismatch between the total number of atoms in the supercells for the interstitial (32 Pb atoms + 1 C atom) and substitutional (31 Pb atom + 1 C atom) cases. The comparison of energy requires the calculation of energy per Pb atom. The energy of each Pb atom (E 1Pb = −3.533 eV/aom) was obtained from the 2 × 2 × 2 supercell without any C in the structure. The two energies compared were as follows: (1) for the interstitial case: (2) and (2) for the substitutional case Total is the total energy of the 2 × 2 × 2 supercell with one C in the octahedral site, and E 31Pb+1C Total is the total energy of the 2 × 2 × 2 supercell with one C in the substitutional site. The difference in energies between Equation (2) and Equation (3) or (E i − E s ) for the 2 × 2 × 2 system size was −2.078 eV (−0.065 eV/atom). This analysis confirms that the interstitial site is preferred over the substitutional site for C in the Pb matrix. It is important to point out that the absolute difference in energy reported here for the two sites depends on the supercell size.

Activation Barrier for C Diffusion in the Pb Matrix
The performance of NAM in LA batteries is improved by the addition of C to the Pb matrix. Since we are considering the incorporation of C in the Pb matrix ( Figure S1), diffusion on C in the Pb electrode will be ideal to tackle sulfation issues at different locations of the NAM. In this section, we attempted to estimate the different migration barriers for C movement in the Pb lattice. All of the calculations in this section were completed using the transition state search (TSS) module within MedeA. The climbing-image nudged elastic band (CI-NEB) [64] approach was used to map the minimum energy pathway between the initial and final configurations. We used a similar computational setting for this part of the examination using the 2 × 2 × 2 supercell, where a 500 eV cutoff energy was used with 3 × 3 × 3 (corresponding to an actual k-spacing of 0.24 × 0.24 × 0.24/Å) k-points. The SCF convergence was set to be 10 −6 eV, and the atomic relaxation was achieved until the force on each atom was 0.02 eV/Å. For all of the combinations simulated, a total of nine images were taken along the unidirectional search path, with a spring constant of 5 eV/Å 2 . The CI-NEB method made sure to estimate the saddle points as the highest points during the migration barrier calculation.
From the structural stability calculations shown in Section 3.2, we performed the CI-NEB transition state search for C moving from one octahedral site to another adjacent octahedral site within the 2 × 2 × 2 supercell (referred to as C Octa i − C Octa i ). Figure 3 shows all of the migration barriers estimated in this study. Figure 3a shows the calculated energies for the initial, final, and other images in between for the C Octa along the transition path. The migration barrier for − was evaluated to be 4.36 eV, which is significantly higher. A similar TSS was completed for the − pathway (Figure 3b). The transition state profile was similar to the − transition with a migration barrier of 3.95 eV. Our calculations suggest that direct movement of C from one interstitial position to the adjacent interstitial position with the same configuration is associated with very large migration barriers. In order to explore other possibilities, we examined the linear transition state search with the initial octahedral to final tetrahedral position, and vice versa. Figure 3c,d provide the CI-NEB energy estimation for − and − diffusion within the Pb matrix. The asymmetry in the initial and final energies was obvious between the octahedral and tetrahedral positions for the two combinations. For the − transition, the energy barrier was estimated to be 1.24 eV, while a relatively smaller barrier of 0.82 eV was predicted for C diffusion from the tetrahedral site to transition, the energy barrier was estimated to be 1.24 eV, while a relatively smaller barrier of 0.82 eV was predicted for C diffusion from the tetrahedral site to the octahedral site, as the absolute energy of the structure with C in the tetrahedral coordinates was higher than for the corresponding octahedral site. Figure S5 in the Supplementary Materials highlights the migration barriers, with data points estimated along the diffusion pathway for each image, and for the climbing image if present.
As established in Section 3.2, C Octa i is the most stable structural configuration of C in the Pb FCC matrix. It is expected that C will diffuse from one octahedral position to other octahedral positions to remain in the minimum energy configuration. Based on our migration barrier calculations, we expect that the C Octa i − C Octa i diffusion pathway should include a tetrahedral site to complete the diffusion in two steps. Figure 4 represents the proposed pathway for C diffusion in the Pb matrix, where the octahedral-to-octahedral site migration is assisted by an intermediate tetrahedral position, i.e., C Octa i → C Tetra i → C Octa i is the preferred migration pathway. Thus, the minimum energy pathway for C migration is proposed to be C Octa i → C Tetra i → C Octa i , where the total energy barrier is 1.24 eV + 0.82 eV; when calculating the percentage difference, it is noted to be~52% of the C Tetra i → C Tetra i migration barrier (3.95 eV) and~47% of the 4.36 eV barrier for the direct C Octa i → C Octa i migration. A similar octahedral-tetrahedral-octahedral diffusion pathway has been reported for Li migration in Li-ion batteries [65,66].

→
→ is the preferred migration pathway. Thus, the minimum energy pathway for C migration is proposed to be → → , where the total energy barrier is 1.24 eV + 0.82 eV; when calculating the percentage difference, it is noted to be ~52% of the → migration barrier (3.95 eV) and ~47% of the 4.36 eV barrier for the direct → migration. A similar octahedral-tetrahedraloctahedral diffusion pathway has been reported for Li migration in Li-ion batteries [65,66].

Conclusions
In this study, we used first-principles calculations to examine the structural stability, defect formation energy, and migration barrier of C in Pb for LA batteries. Our analysis of lattice parameters and elastic constants at room temperature for a pure Pb matrix led to the selection of GGA-PBE functionals in comparison to seven other functionals for examining the Pb-C system within density functional theory. A comparison of system energies of C at different locations within the Pb matrix resulted in the octahedral interstitial site ( ) being more favorable than the remaining interstitial and substitutional sites. In addition to the structural stability, we performed CI-NEB calculations to find the migration barriers between different interstitial sites. Our results predict that C will diffuse from one to another site through an intermediate position

Conclusions
In this study, we used first-principles calculations to examine the structural stability, defect formation energy, and migration barrier of C in Pb for LA batteries. Our analysis of lattice parameters and elastic constants at room temperature for a pure Pb matrix led to the selection of GGA-PBE functionals in comparison to seven other functionals for examining the Pb-C system within density functional theory. A comparison of system energies of C at different locations within the Pb matrix resulted in the octahedral interstitial site (C Octa i ) being more favorable than the remaining interstitial and substitutional sites. In addition to the structural stability, we performed CI-NEB calculations to find the migration barriers between different interstitial sites. Our results predict that C will diffuse from one C Octa The overall activation energy required is less than for a direct C Octa i − C Octa i or similar tetrahedral to tetrahedral migration. With the current understanding of the Pb-C system, future investigation of this system will establish the impact of defects on the structural and electronic properties of Pb-C electrodes. Additionally, the impact of C on the stability of PbSO 4 on Pb surfaces will be part of future work examining sulfation in the presence of C using computational methods.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/solids3020012/s1: Supporting material consisting of additional data from this work has been included in the supplementary document. Reference [67] is cited in the Supplementary Materials.

Data Availability Statement:
The data that support the findings of this study are available within the article and its Supplementary Materials.