Lithium and Sodium Ion Binding Mechanisms and Diffusion Rates in Lignin-Based Hard Carbon Models.

Hard carbons are the primary candidate for the anode of next-generation sodium-ion batteries for large-scale energy storage, as they are sustainable and can possess high charge capacity and long cycle life. These properties along with diffusion rates and ion storage mechanisms are highly dependent on nanostructures. This work uses reactive molecular dynamics simulations to examine lithium and sodium ion storage mechanisms and diffusion in lignin-based hard carbon model systems with varying nanostructures. It was found that sodium will preferentially localize on the surface of curved graphene fragments, while lithium will preferentially bind to the hydrogen dense interfaces of crystalline and amorphous carbon domains. The ion storage mechanisms are explained through ion charge and energy distributions in coordination with snapshots of the simulated systems. It was also revealed that hard carbons with small crystalline volume fractions and moderately sized sheets of curved graphene will yield the highest sodium-ion diffusion rates at ∼10-7 cm2/s. Self-diffusion coefficients were determined by mean square displacement of ions in the models with extension through a confined random walk theory.

Introduction

Efficient, sustainable, and low-cost energy storage is a global necessity. For the past 30 years, Li-ion batteries have been the gold standard and workhorse of energy storage needs for mobile electronics, electric vehicles, medical devices, and so forth; however, lithium is not an infinite resource and its storage in earth’s crust is localized to a few countries. Since this is the case, researchers have been exploring options for the replacement of lithium as the charge carrying ion in energy storage devices. Sodium has been identified as one of the most promising options as it is inexpensive, widely globally available, and can be used in cost- and weight-prohibitive situations such as large-scale grid support and stationary energy storage for renewable energy sources.[1−3]

One of the primary challenges of replacing lithium with sodium in current energy storage devices is the inability for sodium to intercalate within graphite and forms binary graphite intercalation compounds with any reasonable charge density.[4,5] It has been shown previously that sodium will only form NaC64 when inserted into graphite.[6] This has led researchers to explore hard carbons as anode materials. Depending on their nanostructure, hard carbons have the potential to possess a greater charge density, higher resistance to degrade from electrolyte interactions, low working voltage, longer cycle life, and a higher degree of sustainability when compared to the current commercial flake-graphite and spherical graphite anodes.[4,7,8]

Recent research has suggested lignin as a sustainable and domestic source for nanostructured hard carbons with far reaching applications in energy storage.[9−12] Lignin is a highly abundant and renewable resource that possesses high carbon content and an amorphous, cross-linked three-dimensional structure of aromatic polymers.[13] Defining a complete processing–structure–property–performance (PSPP) relationship between lignin and carbonaceous products is difficult because lignin is derived from woody plants and grasses, and the relative fractions of the constituent organic compounds are highly variable by feedstock, which in turn influences the nanostructures and properties of the final carbon composites.[14] Research on the PSPP relationships of lignin reveals that pyrolyzing and reducing lignin produces carbon–carbon composites composed of crystalline (graphitic) and amorphous (disordered graphene sheets) domains. This work refers to lignin-based carbon composites (LBCC) as LBCCs and lignin-based hard carbons to distinguish them from other hard carbon materials. The crystalline volume fraction (CVF), crystallite size, and crystallite form (spheres, fullerenes, onion-fullerenes, nanotubes, multiwalled nanotubes, graphite, and so forth) of lignin-based hard carbons can be tuned via the choice of lignin feedstock, processing, and carbonization temperature.[14−16] Hard carbons can be generated from a variety of natural sources, including, for example, coconut husks, walnut shells, mangosteen shells, coffee grounds, and so forth.[8,17,18] These natural materials differ from the lignin used in this work because they contain multiple components, including cellulose and hemicellulose, as well as other impurities, inorganics, proteins, tannins, and possibly other extractive compounds. A recent study on the structure of hard carbons synthesized from coconut husks, walnut shells, and corn silk showed a range of structures from randomly oriented graphene sheets to stacked sheets depending on feed stock. In general, these materials demonstrate larger graphitic domains, at comparable carbonization temperatures, relative to the LBCCs. The authors attribute this to the catalysis caused by the presence of impurities.[17]

The work of García-Negrón et al. demonstrates that pyrolyzing, reducing at 1050 °C, and ball milling of kraft softwood lignin produce a carbon composite material composed of spherical nanocrystallites embedded in an amorphous graphene matrix which, when processed into an anode and tested in a Li-ion coin cell battery, possesses a specific capacity of 444 mAh/g with 98% Coulombic efficiency over extended galvanostatic cycles.[11] This shows that LBCCs can achieve at least a 20% increase in specific capacity over traditional graphitic anodes (372 mAh/g) and can be considered as a high-efficiency, sustainable, and low-cost option for battery electrodes.

Present challenges faced by researchers with hard carbon electrodes lie in understanding the ion storage mechanisms, preferential ion localization, volume change (swelling) during (de)sodiation and (de)lithiation, as well as the optimal nanostructure–porosity–CVF combination to achieve the highest performance.[4,19] To investigate solutions to some of these challenges for LBCCs, McNutt et al. created large-scale models of the LBCCs with varying crystallite sizes, CVFs, and densities to emulate the LBCCs synthesized at different reduction temperatures from hardwood lignin.[20] Molecular dynamics simulations of the LBCC models charged with lithium revealed that the carbon-edge-terminating hydrogen plays a critical role in the ion storage mechanism for LBCCs as Li ions preferentially localize in the hydrogen dense interfacial region between crystallites and amorphous graphene fragments and allows for Li ions to be stored at a greater density than when intercalated between planes of graphite as LiC6.[21,22] McNutt et al. also explain that as the crystallite size decreases, interfacial volume and hydrogen content increase, leading to larger Li-ion storage capacity.[22] To further explain the ion storage mechanism in LBCCs, Kizzire et al. used a small subsystem of the composites reported by McNutt et al. that consisted of a single nanocrystallite embedded in a matrix of amorphous graphene fragments and simulated with lithium and sodium loading configurations using ReaxFF potentials.[23] Reactive potentials consume more computational resources than nonreactive potentials; however, they allow for modeling of the formation and dissociation of chemical bonds and include both the Coulombic interactions and van der Waals forces necessary for accurate modeling of charged graphitic anodes.[23−26] The ReaxFF potentials were deemed necessary, as accurately capturing the charge transfer between ions and host structure is critical to understanding ion migration and preferential ion localization.[23] Kizzire et al. revealed that sodium, if not initially placed in an intercalated site, will preferentially localize in the amorphous graphene region, whereas lithium will migrate from both intercalated and amorphous graphene initial positions to the hydrogen dense interfacial regions and attempt to form a lithium hydride-like structure but are incapable, as the hydrogen is tethered to the relatively immobile carbon matrix.[23] Results from this previous study prompted interest into investigating lithium and sodium in large-scale LBCC models with ReaxFF potentials.

This work builds upon the previous work of McNutt et al. and Kizzire et al. and investigates lithium and sodium in large-scale LBCC models with reactive potentials to determine preferential localization, composite swelling, mesoscale interactions, and lithium/sodium diffusion rates. We accomplish this by analyzing the resulting radial distribution functions (RDFs), charge and energy distributions, mean square displacement (MSD) of lithium and sodium ions extended by a confined random walk (CRW) theory, and snapshots of charged composites. This work is propelled by interest in using LBCCs as sustainable, domestic, and low-cost electrodes for sodium and lithium-ion batteries.

Results and Discussion

Ion Charge and Binding Energy Analysis

In the following section, we compare the energy and charge distributions for the LBCC models with lithium and sodium loading configurations. Figure shows the binding energy and charge distributions for lithium and sodium ions in the intercalated and amorphous initial loading configurations for the 50% CVF system. Examining the Li-ion binding energy and charge distributions in Figure , we can see that after simulating for 1 ns, the respective distributions are nearly identical for both the amorphous and crystalline intercalated initial loading configurations. This result informs us that the Li ions will migrate to the same regions irrelevant of the initial position, which is in good agreement with previous works.[22,23]

Figure 1

Binding energy and charge distributions for lithium (a,b) and sodium (c,d) in the 50% CVF system for ions initialized in the amorphous and crystalline domains.

Examining the Na-ion binding energy and charge distributions in Figure for the 50% CVF system simulated for 1 ns, we can see a single mode distribution for Na ions intercalated in the crystallites and a distinct bimodal distribution for Na ions initialized in the amorphous domain. Through examination of individual ions in snapshots of the simulation frames and identifying their charges and binding energies, we found that Na ions sandwiched between neighboring planes of amorphous graphene fragments had similar binding energies and charges to those Na ions that were intercalated within the crystalline domain. These “doubly bound” Na ions had deeper binding energies and higher charges compared to the Na ions that adsorbed onto the planar surfaces of amorphous graphene fragments and crystallites.

Figure a,b shows the binding energy and charge distributions after 1 ns of simulation for Na ions initialized in the amorphous graphene domain for the 10, 50, and 90% CVF systems. Inspection of Figure a,b shows a large percentage of Na ions having deeper binding energy and greater charge in the 90% CVF system compared to the 10 and 50% CVF systems. Na ions with binding energies that average −37 kcal/mol in the 90% CVF system correlates to Na ions that are sandwiched between adjacent graphene planes or Na ions at intercalation positions at the edge of nanocrystallites with high amounts of disorder in interplanar spacing and angles. Na ions with binding energies near −20 kcal/mol are found adsorbed onto a graphene surface or a basal plane of a nanocrystallite. In Figure a,d, the energy of a sodium ion interacting with one or two carbon planes is evident. The binding energy is essentially doubled when the sodium interacts with two planes. The distribution between these two limits reflects the disordered nature of the composite. The greater percentage of Na ions with deeper binding energy in the 90% CVF system results from the high fraction of graphene planes directly adjacent to crystallites or each other which decreases the amount of adsorption sites. The lower CVF systems allow for a more even distribution between these two Na-ion localizations.

Figure 2

(a,b) Binding energy and charge distribution for sodium initialized in the amorphous domain for the 10, 50, and 90% CVF systems. Colored dots correspond to binding energy and charges for colored sodium ions in (d). (c) Front facing view of the sodiated 10% CVF system with crystalline carbon (red), amorphous graphene fragments (blue), sodium (white), and hydrogen (removed for clarity). (d) Enlarged section of the 10% CVF system with sodium color coded to represent charge and binding location. Na ions bound to the surface of graphene and crystallites (green), Na ions intercalated between neighboring sheets of graphene (light blue), Na ions intercalated within edges of nanocrystallites (purple), and Na ions bound to other Na-ions in a semi-metallic–like state (orange). [Some surface-adsorbed Na ions and edge-terminating hydrogen have been removed for figure clarity in (d)].

Figures a,d and 3a show that most Na ions in the 10% CVF system are adsorbed onto the face of a graphene fragment. Further, though the sodium ions were initialized randomly throughout the composite, there are obvious regions in the amorphous graphene domain with higher and lower concentrations of sodium, as seen in Figure a, suggesting that in these low CVF composite systems, sodium will aggregate.

Figure 3

Snapshot slices of LBCC systems with carbon nanocrystallites (gray) and sodium (red) initialized in the amorphous graphene (blue) and porous domains for (a) 10% CVF, (b) 50% CVF, and (c) 90% CVF.

Interestingly, the charge distribution for Na ions in the 10% CVF system shows a third state of Na-ion charge, centered at 0.06 e, not present in other systems. To identify the source of this third state of Na-ion charge, see Figure d, which presents a zoomed section of Figure c with Na ions color coded to correspond to a charge value. Light blue and purple represent doubly bound Na ions in the amorphous (blue) and crystalline (red) domains, respectively, with an average charge value of 0.36 e. Light green represents the Na ions adsorbed (or singly bound) to the surface of an amorphous or crystalline carbon plane with an average charge value of 0.225 e, while orange represents the third localization only found in the 10% CVF system with an average charge value of 0.06 e and low average binding energy of −14 kcal/mol. These orange Na ions are bound to each other, and the low charge represents a quasi-metallic like state. Higher loadings of Na ions in these moderately porous composites would create more Na-ion clustering within the pores, similar to the orange-colored ions in Figure d. Na-ion clustering inside pores has been reported by others in the literature as stable configurations that have been shown to be highly reversible and enable charge densities near 300 mAh/g in hard carbon anodes.[19,27]

While the binding energy distributions in Figure show that intercalation positions are more energetically favorable for sodium, the barrier for Na-ion intercalation is very high, as reported in the literature.[4,5,28] This is true except for the case where nanocrystallite planes have shifted, and local interplanar distance is larger than 3.6 Å. The shifting of planes in graphitic nanocrystallites has been predicted and analyzed in the previous work.[29]

From these results, we find that lithium and sodium storage mechanisms are fundamentally different in lignin-based hard carbons. In these simulations, we observe that lithium migrates out of intercalation/adsorption sites and into the hydrogen dense interfacial region, as shown in Figure a,b, whereas a portion of the sodium initialized in the amorphous domain are found to migrate past the interfacial region to find intercalation positions in the nanocrystallites with expanded interplanar spacing, as seen in Figure d. This shows that sodium will preferentially bind to the carbon matrix, while lithium will preferentially bind to the hydrogen.

Figure 5

Snapshot slices of the 50% CVF systems after simulation for 1 ns with lithium (yellow), sodium (red), crystalline carbon (gray), amorphous carbon (blue), and hydrogen (removed for clarity); (a) lithium initialized within the amorphous domain, (b) lithium initialized as intercalated within the crystalline domains, (c) sodium initialized within the amorphous domain, and (d) sodium initialized as intercalated within the crystalline domains.

The energetic origin of the different binding sites for lithium and sodium can be explained in part by the energies of the respective hydrides. For this potential, lithium hydride has an energy of −38.65 kcal/mol,[23] which is much more favorable than the per ion energy for lithium in the carbon composite, which drives lithium ions to the hydrogen terminated interfaces. In the case of sodium, the hydride has an energy of −33.57 kcal/mol, which is less favorable than the per ion energy for sodium in the carbon composite. Therefore, the low energy state for sodium ions in the composite interacts with two carbon planes. The lithium-ion binding mechanism is limited to hard carbons with a similar nanostructure and high hydrogen content in the interfacial region between crystalline graphite particles and graphene fragments, while the sodium-ion binding mechanism is applicable to hard carbon anodes that possess similar nanostructures, namely, a porous nanostructure with small nanocrystalline carbon domains embedded in a graphene fragment matrix.

Analysis of the energy and charge distributions in conjunction with the snapshots shown in Figures d, 3a–c, and 5a–c suggests that in application, sodium insertion into LBCC anodes would result in Na ions preferentially adsorbing to the surface of amorphous graphene fragments and the surface planes of nanocrystallites with a smaller fraction intercalating along the edges of nanocrystallites where local interplanar spacing is above 3.6 Å due to an inherent disorder in the system. Figure b,d implies that after the preferential filling of adsorption and intercalation storage sites, sodium will fill porous regions in the composite giving rise to an adsorption–intercalation-pore filling sodiation scheme for LBCC anodes. For this work, pores should be defined as an open space between graphene planes or nanocrystallites with spacing greater than 6.5 Å, in agreement with the density functional theory (DFT) work on hard carbon by Olsson et al.[6]

Previous work using lignin carbon fiber mats (LCFs) has given us guiding information on the nature of solid-electrolyte interphase (SEI) formation in lignin-based carbons. The work conducted by Tenhaeff et al. has stated that for LCFs in electrodes using electrolytes with propylene carbonate (PC), little SEI formation was observed for LCFs carbonized at 1000 °C, while an appreciable SEI was formed upon lithiation of the LCFs carbonized at 2000 °C.[9] The authors state that the difference in SEI formation is due to the amount of disorder present in the 1000 °C LCFs, which limits SEI formation and that PC electrolytes exfoliate the graphitic domains in the 2000 °C LCFs, forming an appreciable SEI.[9] We posit that LBCCs would behave in a similar manner, with lower reduction temperatures and small CVFs delivering the best resistance to large SEI formation.

For glucose-based hard carbons, Au et al. found that pores were highly interconnected at carbonization temperatures of 1000 °C, and while pores were larger for carbonization at 2000 °C, the increasing size of the graphitic regions closed off parts the interconnected pore structures, limiting ion diffusion.[30] It is reasonable that porosity in LBCCs would progress in a similar manner, suggesting that there exists an optimal lignin reduction (carbonization) temperature above 1000 °C which would maximize the sodium storage capacity in the adsorption–intercalation-pore filling scheme, while retaining interconnected pores for fast ion diffusion.

Qualitatively speaking, from these results, it is reasonable that lignin-based hard carbons with lower CVF, smaller nanocrystallites, and moderate porosity would allow for the highest energy density for sodiated LBCC anodes. The same can be said for lithium, as smaller crystallites increase the surface area to volume ratio and thus increase the amount of hydrogen present in the interfacial region allowing a higher density of lithium storage than graphitic intercalation.[11,21−23] In general, the anode should possess a sufficiently strong binding to achieve high capacity but avoid binding that is too strong that it leads to irreversible adsorption. The distributions in sodium binding energy shown in Figure a demonstrate that the binding energy can be controlled through the composite nanostructure generated through controlled processing conditions. Because carbon edges in lignin-based hard carbons are terminated with hydrogen, side reactions between carbon and lithium or sodium should be minimal and not decrease the Coulombic efficiency. García-Negrón et al. has shown that Coulombic efficiencies of over 98% can be achieved after extended galvanostatic cycling of Li-ion coin cells with lignin-based hard carbons as the active anode material.[11]

Anode Swelling

In applications, knowledge of the volume change that occurs in an anode during ion (de)loading is vitally important as excessive volume change can damage battery structure leading to failure with safety concerns. In general, the volume change between empty and fully intercalated graphitic anodes in commercial Li-ion batteries is ≤10–14%.[31,32] The swelling for each of the LBCC simulated systems can be found in Table . We can see that lithium initialized in the amorphous domain produces the least amount of swelling, which is to be expected because lithium preferentially localizes in the interfacial regions, bound to hydrogen at a greater density than when intercalated in graphite.[21−23] LBCCs loaded with sodium exhibit roughly 50% greater swelling than composites loaded with lithium. This is also expected as sodium has a greater ionic radius and does not exhibit the same high-density binding with hydrogen as lithium. We note that these swelling values were obtained from simulating at atmospheric pressure and anode structure could isotropically expand, whereas in application, the anode structure is constrained within the battery housing. Additionally, the Li-ion charge density in these simulated systems is approximately one-third that of fully Li-intercalated graphite because the charge density was chosen to correspond to charge density in previous works as stated in the methods section. Reporting of these swelling values are meant to provide reference for future experimental endeavors in the creation and characterization of Li and Na-ion batteries with LBCC anodes.

Table 2

Collection of Simulated Systems with Relevant Parameters

system	crystalline volume percentage	number of atoms	initial ion placement	charge density (mA h/g)	composite density (g/cm3)	crystallite radius (Å)	swelling percentage
1	10	658,510	none	0	1.68	17
2	10	689,788	Na-amorphous	111.82	1.60	17	15.6
3	50	108,148	none	0	1.62	7
4	50	113,160	Li-intercalated	137.44	1.53	7	10.0
5	50	113,160	Li-amorphous	137.44	1.54	7	8.7
6	50	113,160	Na-intercalated	126.99	1.60	7	13.3
7	50	113,160	Na-amorphous	126.99	1.57	7	15.8
8	90	150,950	none	0	1.54	5
9	90	155,964	Na-amorphous	100.45	1.49	5	12.9

Local Structure Analysis

In Figure a–d, the ion–ion and ion–hydrogen RDFs are shown for the 50% CVF system with amorphous and crystalline initial loading states. The Li–Li and Li–H RDFs found in Figure a,b are highly similar as both initial loading conditions result in Li-ions migrating to the hydrogen dense interfacial region as can be seen in the simulation cell slices in Figure a,b. One would expect there to be more long-range structures in the Li–H RDF due to the Li-ions affinity for bonding to the hydrogen; however, because the hydrogen is essentially tethered to the relatively immobile carbon, no long-range Li–H structures can exist. The Li–Li RDFs found in Figure a,b are similar for both initial loading conditions (either in the crystalline or amorphous domains) because the simulation in which the Li-ions began in the crystalline domain resulted in Li-ions migrating to the hydrogen dense interfacial region as can be seen in the simulation cell slices in Figure a,b. Thus, the two different initial conditions led to the same equilibrium distribution. One would expect there to be more long-range structure in the Li–H RDF due to the Li-ions affinity for bonding to the hydrogen; however, because the hydrogen is essentially tethered to the relatively immobile carbon, no long-range Li–H structure can exist. In the case of sodium, when the sodium ions are initially intercalated in the crystalline domain, they do not migrate out of the crystallites. This energetically favorable initial condition does not lead to the same state as when the Na ions are initially placed in the amorphous domain, leading to different RDFs. The dip occurring in the Na–Na RDF for Na ions intercalated within the crystalline domain shown in Figure c near 9 Å denotes the average distance of a Na ion to the interfacial region where no ions are present, and the subsequent rise near 11 Å is the average distance between Na ions found between separate nanocrystallites as observable in Figure d.

Figure 4

Component RDFs for ions initialized in the amorphous graphene and crystalline intercalation domains for the 50% CVF system. (a) Li–Li RDF, (b) Li–H RDF, (c) Na–Na RDF, and (d) Na–H RDF.

The Na ion component RDFs for the various composites can be seen in Figure a–d along with visual representations of the ion–atom pairs that constitute each peak. The most notable among these RDFs is shown in Figure a, where the increased intensity in Na–Na pairs for the 10% CVF system denotes a greater local density of Na ions suggesting an increased amount of agglomeration, as can be seen in Figure a. Figures a and 3a–c reveal an inverse relationship between CVF and local Na-ion density, with low CVF and moderate porosity displaying the highest degree of Na-ion agglomeration.

Figure 6

Na-atom component RDFs for each of the amorphous sodiated LDCC systems with corresponding snapshots of the general Na-atom pairs representing each peak in the RDFs. (a) Na–Na RDFs, (b) Na–H RDFs, (c) Na-amorphous graphene RDFs, and (d) Na-crystalline carbon RDFs.

Ion Diffusion

To calculate the self-diffusion coefficients for lithium and sodium in the LBCC anodes, we recorded the unwrapped coordinates of sodium and lithium ions during simulations and calculated the MSD of ions through the composites. The MD-generated MSDs were then fit with the CRW simulation at room temperature and extended to 100 ns. The cage radius and cage-to-cage hopping probability reported in Table represent a characteristic length scale of confinement and a probability proportional to the activation barrier to ion diffusion respectively.[33] Where the cage radius is less than the diameter of an atom, this describes the relative volume explored by the point at the center of the ion. The exponent value details the linear proportionality of MSD to observation time, which is required by the Einstein relation. Table reports the MSD values of MD simulation alone and with extension to the long-time limit (represented with an exponent value near 1.0) with CRW theory. The MSD from MD simulation is plotted with their corresponding CRW extensions up to 1 ns in Figure . We note that the CRW was simulated out to 100 ns but plotted to 1 ns for clarity in comparing with the MD simulations. The MSD data from MD simulation are plotted to 0.5 ns because auto correlation functions become noisy near the end because there is a decreasing amount of data in each subsequent point. Likewise, the calculations of diffusion coefficients from MD simulation only used data up to 0.5 ns. The self-diffusion coefficients were calculated using MSD with extensions through CRW theory to reach the long-time limit required by the Einstein relation. This diffusion model is universal for confined diffusion systems (aside from bulk metallic glasses)[33] and applicable to other hard carbon materials.

Table 1

Self-Diffusion Coefficient Values for Sodium and Lithium in Lignin-Based Hard Carbons from MD Simulation and CRW Extension

system	ion-initial-domain-CVF system	MD exponent	CRW exponent	MD diffusion coefficient (cm2/s)	CRW diffusion coefficient (cm2/s)	cage radius (Å)	cage hopping probability
2	Na-amorphous-10	0.77	0.97	2.70 × 10–7	2.76 × 10–7	0.90	0.00141
4	Li-intercalated-50	0.66	0.97	1.53 × 10–7	1.63 × 10–7	1.05	0.00067
5	Li-amorphous-50	0.64	0.97	1.24 × 10–7	1.27 × 10–7	0.92	0.000617
6	Na-intercalated-50	0.45	1.10	1.04 × 10–8	1.43 × 10–8	0.44	0.000134
7	Na-amorphous-50	0.67	1.04	2.00 × 10–8	2.29 × 10–8	0.398	0.000275
9	Na-amorphous-90	0.65	0.95	2.51 × 10–8	2.62 × 10–8	0.51	0.000252

Figure 7

MSD generated from MD simulations (color) with their corresponding CRW extensions plotted to 1 ns. Legend is read as a migrating ion—initial domain (A, amorphous; C, intercalated within nanocrystallite)—LBCC model.

We find that the CRW values for the self-diffusion coefficients for lithium in the 50% CVF system and sodium in the 10% CVF system are on par with the experimentally found and ab initio-calculated diffusion rate of lithium in pristine graphite in the planar direction, 4.4 × 10–7 cm2/s.[34] The CRW values of the diffusion rate for sodium in the 50 and 90% CVF systems are slightly smaller with values ∼10–8 cm2/s. Sodium in the 10% CVF system was found to have the highest diffusion rate of all simulated systems with a value of 2.8 × 10–7 cm2/s, while sodium in the 90% CVF system was found to have the lowest diffusion rate among the systems studied.

Recent DFT studies of alkali metals in hard carbon anodes state that the curved graphene morphology decreases the alkali metal ion migration barrier and that, combined with porosity, majorly contributes to high self-diffusion coefficients and higher cycling performance in hard carbon anodes.[6,35] Because the 10% CVF system possesses a higher degree of porosity and graphene curvature compared to the 50 and 90% CVF systems, the high diffusion rate of sodium in the low CVF system is substantiated.

Tian et al. state that the factors that increase rate performance in anodes are well known and include, among other things, reducing active particle size and increasing solid-state diffusivity and electrode porosity.[36] For specific application where power density or fast charging is paramount, lignin-based hard carbon anodes should possess smaller crystallites (active particles), a lower crystallite volume fraction, and an interconnected pore structure, all of which would increase the diffusion rate of lithium and sodium through LBCC anodes and are achievable through lower reduction temperature of lignin. Additionally, low reduction temperatures (in the neighborhood of 1000 °C) provide the highest surface area to volume ratio and thus the highest hydrogen content, maximizing the lithium storage capacity in lignin-based hard carbon anodes.

Conclusions

Reactive molecular dynamics simulations were carried out for lithium and sodium loaded in three large lignin-based-carbon-composite systems with 10, 50, and 90% CVFs. The reactive potentials used for this work were deemed necessary to accurately capture the ion binding mechanisms, diffusion properties, and the complex mesoscale structure intrinsic to lignin-based hard carbons. Analysis of energy and charge distributions in conjunction with snapshots of the lithiated systems shows that lithium will preferentially localize in the hydrogen dense interfacial region between crystallites and amorphous graphene fragments regardless of initial localization.

Snapshots of the sodiated systems in conjunction with charge and energy distributions reveal that sodium will preferentially bind to the surface of graphene and basal surfaces of nanocrystallites with a small fraction intercalating at the edges of nanocrystallites that have local d-spacing above 3.6 Å due to the inherent disorder in the nanocrystallites. Once the adsorption and intercalation positions have been filled, sodium will agglomerate in pores. This adsorption–intercalation-pore filling sodiation scheme leads to high charge capacity in hard carbon anodes. The lower binding energies found for the adsorption and pore filling sodium ions also suggest these storage mechanisms to be largely reversible.

It was found that the LBCC system with the lowest CVF and curved graphene fragments along pores produces the largest sodium ion diffusion rate among the composites studied in this work. The results of this study indicate that a porous LBCC with low CVF and sheets of curved graphene will produce an anode with high diffusion rate and large charge capacity for both lithium and sodium-ion batteries. We also posit that a lignin-based anode with these structural characteristics would lead to faster charging and more power delivery in ion batteries. The correlation of features in an experimental sodiation process and the atomic level binding mechanisms in this paper are of keen interest. At this time, we show that there is a range of binding energies associated with different types of binding sites. Sodiation of coin cell batteries with LBCC anodes will allow us to explore this correlation of features in future work.

Theoretical Methods

The hard carbon models in this work were designed by McNutt et al. to emulate the nanostructure of hardwood lignin pyrolyzed and reduced at 1000, 1500, and 2000 °C as synthesized and characterized by Tenhaeff et al.[9,20] The hard carbon models possess spherical AB stacked graphite crystallites with radii of 5, 7, and 17 Å embedded in an amorphous graphene fragment matrix at 90, 50, and 10% CVFs, respectively. All crystalline and amorphous edge carbons were terminated with hydrogen. Relaxation of the model resulted in slight bending of the graphene fragments in the amorphous domain and shifts in crystalline planes such that the equilibrium interplanar spacing became 3.4 Å, representative of the disorder in the real LBCC system and verified as accurate by comparison of the simulated and experimental RDFs.[20]

A total of nine reactive simulations (three without ion loading, six with ion loading) were performed using LAMMPS and with ReaxFF potentials developed by Hjertenaes et al. and Raju et al. for the sodiated and lithiated systems, respectively.[24,37,38] Previous works have verified that the two reactive potentials are the same in their handling of carbon–carbon and carbon–hydrogen interactions, and thus, the potential reported by Raju et al. was used for the systems without ions.[23] The nine systems were relaxed at 1 atm in the NPT ensemble at 298 K with 0.25 fs timestep until potential energy was equilibrated. The six systems with lithium/sodium loading were then simulated for 1 ns in the NVT ensemble at 298 K with 0.25 fs timestep. The trajectory files were saved in both wrapped and unwrapped configurations for the RDF and MSD analysis, respectively, and the volume of each system was recorded for swelling calculations. The charge densities for Na-ion systems were set between 100 and 125 mA h/g, consistent with values used in previous work for these composite systems.[21] The differing charge density between sodium and lithium systems is due to the difference in ion mass, as all 50% CVF systems have the same number of ions.

Ideally, the results of a simulation are independent of initial ion placement when the simulation is run a sufficiently long time to drive the system to thermodynamic equilibrium. However, the finite simulation time and kinetic barriers result in systems with distinct initial conditions, such as ions initially placed in the graphitic versus amorphous domains, not arriving at the same state. This was investigated by McNutt for lithium.[22] Because the energy was lower for the amorphous system, they judged that it was the more energetically probable state. Based on this result, in the simulation matrix implemented in this present work, some of the composites are investigated with initial placement of ions in both the crystalline and amorphous domains, while others are investigated exclusively with ions initially placed in the amorphous domain.

The 90 and 10% CVF systems were simulated uncharged and with sodium initialized in the amorphous carbon domain. The 50% CVF system was simulated uncharged, with sodium and lithium initialized in intercalated positions within the crystalline carbon domain, and with sodium and lithium initialized in the amorphous carbon domain. The 90% CVF system contained 155,964 atoms (88,447 crystalline carbon, 8,835 amorphous carbon, 53,668 hydrogen, and 5,014 sodium). The 50% CVF system contained 113,160 atoms (49,232 crystalline carbon, 26,563 amorphous carbon, 32,353 hydrogen, and 5,012 lithium/sodium). The 10% CVF system contained 689,788 atoms (423,744 crystalline carbon, 131,915 amorphous carbon, 102,814 hydrogen, and 31,278 sodium). The large number of atoms in each system are necessary to capture both the mesoscale structure of LBCC anodes and an accurate CVF with appropriately sized crystallites. These model structures have been extensively compared to synthesized carbon composites.[20] A full table of system details can be found in Table .

For application purposes, knowledge of diffusion rates and ion migration is critical to understanding the performance of an anode material. The self-diffusion coefficient is obtained by using the Einstein relation and calculating a single-particle autocorrelation function, the MSD. The Einstein relation includes the condition that the MSD is linearly proportional to observation time, which occurs in the infinite-time limit. Simulating confined systems that operate with short time scales (1 ns) often do not meet this condition, and thus, the application of the Einstein relation is not valid.[33] A robust solution to this issue is shown by Calvo-Muñoz et al. where the MSD of MD simulations can be extended to reach the infinite-time limit by fitting the MSD of a CRW simulation to the MSD from the MD simulation.[33] The CRW theory uses two physical parameters, cage radius and cage-to-cage hopping probability. These parameters represent the physical system’s dimensions and the activation barrier for diffusion respectively, ensuring an accurate result for the self-diffusion coefficient. This work uses the same CRW simulation code as Calvo-Muñoz et al. to obtain self-diffusion coefficients for lithium and sodium in the LBCC anodes.