Defects in Hard Carbon: Where Are They Located and How Does the Location Affect Alkaline Metal Storage?

50%) in defect formation energy is observed for all defects except the nitrogen substitutional defect. These results confirm the experimentally observed enhanced defect concentration at surfaces. Previous studies have shown that defects have a marked impact on metal storage. This work explores the interplay between position and defect type for lithium, sodium, and potassium adsorption. Regardless of defect location, it is found that the energetic contributions to the metal adsorption energies are princi-pally dictated by the defect type and carbon interlayer distance.


Introduction
The increase in global energy demand has necessitated the development of efficient and sustainable energy storage technologies. Rechargeable batteries play a crucial role in providing energy storage needed to transition from a fossil-fuel-based to a sustainable energy economy. [1,2] Alkali ion batteries including lithium (Li)-ion batteries (LIBs), sodium (Na)-ion batteries (NIBs), and potassium (K)-ion batteries (KIBs) have received increasing attention as rechargeable batteries. [1,3,4] Na and K are among the most abundant elements, making them promising and suitable Pyrolysis of precursors such as sugars, biomass, or phenolic resins releases heteroatoms, and the heteroatom concentration decreases with thermal treatment temperature. [6,7,10,23] Even HCs with similar structure but with different defect concentrations show markedly different NIB cycling behaviors, making the consideration of defects imperative. [23] Furthermore, heteroatom doping, annealing, and defect engineering as a means to improve the anode performance have been reported. [9,15,28,29] However, contrasting reports on the effect of these defects on anode performance exist and further atomic-scale study is required to understand the role of individual defects. [9,30,31] The complex and varying structure of HC makes the modeling of a universal HC material unfeasible. [9] Breaking the HC material down into simpler structural motifs allows the individual structural features found in many HCs to be understood in isolation. Previously, we have studied defects on the basal plane (using graphene as a model for this structural motif), cylindrical pores (using carbon nanotubes as a model for curved motifs), and planar graphitic pores (using bilayer graphite as a model for the graphitic stacks with varying interlayer distances). [7,21,[32][33][34][35] In this paper, we consider more complex HC models, to include the effect of edges, curvature, and strain on the defect formation. These models are based on the reconstructed graphite (10)(11) surface, which we have employed previously to study metal binding at the transition from curved to planar carbon morphologies. The interlayer distances in the planar region (corresponding to stacked graphitic layers) have been chosen to reproduce the experimental mean separation and range. [7,21,33,36] This approach allows for the motifs found in HC materials to be modeled in a systematic way. Furthermore, these models are relevant for other carbon materials based on graphite, carbon nanotubes, and graphene as these all share common structural features. This handling allows us to study the defect formation as a function of lattice position and the effect of defect position on the initial alkali metal (Li/Na/K) incorporation.

Results and Discussion
From experimental evidence, carbon-based anodes have both curved and planar motifs, with surfaces and their defects having a particular influence on electrolyte breakdown and the initial metal adsorption. [7,11,[36][37][38][39] These curved pores have been observed both from pair distribution analysis and transmission electron microscopy imaging. [7,11,[36][37][38] In this study, we investigate the effect of these morphologies on defect formation as a function of lattice position. Based on our previous investigation of defect formation on graphene, [32,34] the following defects are considered: carbon vacancy (V C ), nitrogen substitutional or graphitic nitrogen defect (N C ), nitrogen substitutional defect and carbon vacancy or pyridinic nitrogen defect (N C V C ), oxygen substitutional defect (O C ), double oxygen substitutional defect (2O C ), oxygen substitutional defect and carbon vacancy (O C V C ), triple oxygen substitutional defect (3O C ), nitrogen and oxygen substitutional defect (N C O C ), and double oxygen single nitrogen substitutional defect (N C 2O C ). [7,32,34] The curved carbon models used in this work are based on the reconstructed graphite models previously developed by Lechner et al. and Thinius et al. [20,40] The pristine model consists of 800 carbon atoms, with an average interlayer separation of 3.75 Å (ranging from 4.04 Å at defect position 1 to 3.09 Å at defect position 8) based on previous small angle X-ray scattering/wide angle X-ray scattering characterization of HC anode materials. [7,21,33,36] The defects highlighted above are then considered in this system at various distances from the carbon surface, which will henceforth be referred to as position X (shown in Figure 1). This allows for consideration of defect formation in graphitic stacks (highlighted in Figure 1a as the bulk), in the transition region between intercalation and surface behavior (highlighted in Figure 1a as the near surface), and finally at the curved surface.
The defect formation energy ( ) f defect E can be used as a measurement of defect concentration (the lower the f defect E the higher Figure 1. Schematic representation of a) the different regions (bulk-like, near-surface-like, and surface-like) in the simulation cell as defined in this paper, and b) pristine simulation cell with brown spheres being carbon atoms and green spheres signifying the carbon atoms at c) the general defect positions. Red arrows indicate the x-direction, green the y-direction, and blue the z-direction. as a result of a small expansion in the CC bonds as the surface is approached (+1-2%), reducing the steric repulsion (Figure 2c,d and Table S1 (Supporting Information)). Hence, on an energetic basis, it would be reasonable to expect dramatically more surface or near surface V C defects than bulk V C defects (Figure 2a and Figures S1a and S2 (Supporting Information)). From a structural perspective, removal of a C-atom from the lattice results in three equivalent C-dangling bonds, which after relaxation form a 5-membered ring and a single C-dangling bond (Figure 2c and Figures S6 and S7 (Supporting Information)). The CC separation of the 5-membered ring is 2.0 Å (Figure 2c label 1, and Figure 2d), while the CC separation of the dangling bond is 2.6 Å (Figure 2c labels 2 and 3, and Figure 2d). The most significant geometric change between the bulk and surface V C defect is the reduction in the CC bond length (≈25% contraction) in the 5-membered rings (Figure 2d). This is matched by an increase (≈33%) in CC separation of the dangling bond (Figure 2d and Figure S7 and Table S1 (Supporting Information)). A full list of the CC bond lengths and separations is included in Figure S7 and Table S1 (Supporting Information). In the bulk position, the V C defect has three degenerate Jahn-Teller distortions (which has previously been observed for carbon monovacancies in graphene [51][52][53][54] ), and the 5-membered ring is equally likely to form between any of the C-dangling bonds. The introduction of the curved surface breaks the threefold symmetry of the system and gives rise to two degenerate configurations and one distinct configuration. In defect positions 1-4, there is no energetic difference between these configurations. As the V C defect is moved toward the surface, the differences between configurations becomes significant at positions 8-13 (Figure 2d and Figure S7 and Table S1 (Supporting Information)). It should be noted that after position 8, the initial positional degeneracy breaks further upon relaxation. V C at position 8 is further worthy of more consideration as it is the only configuration that shows crosslinking between adjacent graphene sheets (Figure 2e,f). The degree of cross-linking in HC anodes is from experimental evidence dependent on the initial precursors and the intermediate states in the pyrolysis stages. [1,11] From these density functional theory (DFT) simulations, cross-linking is a surface effect and only observed where the curvature brings the sheet separation to a minimum and the dangling bond is pointing toward the neighboring graphene sheet (Figure 2f). Interlayer distances of <3.7 Å are observed in HC anodes, [7,21] as well as cross-linking between graphitic regions, [1] supporting these computational results. A series of nudged elastic band calculations were performed to identify the transition state and barrier height. These calculations show the process to proceed via a low barrier in the forward direction (CC forming E a-forward = 0.2 eV, and CC breaking E a-reverse = 0.9 eV), although it is important to note that these barriers are highly sensitive to layer separation. Any processes giving rise to these cross-linked motifs are potentially vital in understanding the breakdown of the electrolyte at surface defects in HCs, and more generally for the stability of a given HC and its associated solid electrolyte interphase.

Nitrogen Defects
Nitrogen doping of HCs can increase carbon surface wettability to the electrolyte, improving battery performance and metal ion storage capacity. [55,56] Nitrogen defects in carbon materials are commonly found as pyridinic, graphitic, and pyrrolic. Pyridinic is referred to here as the N C V C defect, where the nitrogen heteroatom is introduced at a carbon lattice site bonded to two carbon atoms in a 6-membered ring together with an adjacent carbon vacancy. Graphitic is referred to as the N C defect, where the nitrogen heteroatom is sitting at a carbon lattice site bonded to three carbons. Pyrrolic is not considered in this work, which is the nitrogen heteroatom bonded to two carbon atoms in a 5-membered ring, opening up the adjacent 6-membered carbon rings. [57,58] Previous work has demonstrated that the N C and N C V C defects are more energetically stable in graphene and carbon nanotubes (which are similarly curved to the surface sites in the model used here) than the pyrrolic defect, [59] and hence the nitrogen defects considered here are limited to the graphitic and pyridinic nitrogen defects. Substituting one carbon atom for a nitrogen forms the N C defect (Figure 3a,c,d and Figure S8 (Supporting Information)). N C is incorporated at the carbon lattice site with negligible distortion and has been widely studied as a dopant in sp 2 carbon lattices. [60,61] The N C V C defect is geometrically similar to the V C defect, with the C 6 ring next to the V C replaced by a C 5 N ring (Figure 3b,e-g and Figures S10 and S11 (Supporting Information)). This defect requires additional configurations (as compared to the N C defect) to be considered (  are found in the region directly below the surface, it would be reasonable to expect an, albeit, small reduction in N C concentration in these regions (Figures S1b and S2, Supporting Information). For the N C defect, the structural changes with respect to defect position are less pronounced (Figure 4a). The N atom sits at the 3-coordinated carbon lattice site with a small contraction of each NC bond (1.41 Å) with respect to the CC bond (1.42 Å), although it is important to note that this contraction of bond length is <1% when compared to the pristine carbon lattice (Figure 4a). The presence of the surface has little effect on the defect geometry with the NC bond showing the same surface bond length expansion as the neighboring CC, with the same symmetry dependence (Figure 4a and Table S2 (Supporting Information)).
For the N C V C defect on the other hand (Figure 3b), there is a clear trend with respect to defect position. As opposed to N C and V C , the N C V C defect shows several cross trends.  for the N C V C defect are still significant, so even the lowest energy configuration has an f defect E of 2.71 eV, substantially higher than the N C f defect E but lower than the V C .
f defect E However, with a significant reduction (≈50%) in f defect E as the N C V C defect position approaches the surface (as compared to the bulk defect positions), a higher concentration of these defects in the surface and near surface regimes would be anticipated (Figures S1c and S3, Supporting Information).
In the N C V C defect, the two C-dangling bonds relax inward and form a 5-membered ring, with a CC separation of ≈1.7 Å in the bulk regime (Figure 3e,f). The CC separation decreases as the surface is approached to a minimum CC separation of 1.5 Å. The C 5 N ring shows only marginal geometry changes as the surface is approached. In the bulk regime, the CN bond length is 1.33 Å and the CC bond length 1.41 Å (Figure 4b and Figure S12 and Table S3 (Supporting Information)), with no deviation observed between the perpendicular and tilted variants. In the near surface region, the previously described trends for the perpendicular and tilted variants are observed. A symmetric extension of the bonds in the direction of the surface is observed in the perpendicular case. The CC extension in the C 5 N rings is very small ≈0.01 Å, whereas the C 5 ring shows a more pronounced extension (0.04 Å), as the CC bonds attempt to equalize the strain around the C 5 ring (Figure 4b).
The tilted configuration shows an extension of the bonds in the direction of the surface, which has a negligible effect on the CN bonds, with an extension of <0.01 Å, the C 5 N CC bonds show a similarly small extension of ≈0.02 Å (Figure 4b and Figure S12 and Table S3 (Supporting Information)). The impact of the surface is observed in the C 5 ring with an extension of ≈0.05 Å. It should be noted that the presence of the surface is only partly responsible, as the main driver for this extension is to equilibrate the CC bonds around the C 5 ring, the presence of the surface provides the structural flexibility to achieve this.

Oxygen Defects
Oxygen-containing defects are prevalent in HC anodes, and can play an important role in increasing the surface wettability, which in turn improves the battery stability and performance. [1,7,15,36,39,[62][63][64][65] Tuning the oxygen defect concentration and position (surface vs bulk) allows for physicochemical properties to be optimized. [1,62,64] Experimental XPS measurements commonly identify an oxygen content of ≈4% in HC, soft carbon, and composite carbon anodes for LIBs, NIBs, and KIBs. [66,67] Here, we consider four oxygen defects: single substitutional oxygen defect (O C ), oxygen substitutional defect and carbon vacancy (O C V C ), double oxygen substitutional defect (2O C ), and triple oxygen substitutional defect (3O C ). The O C and 2O C defects are purely substitutional in nature, whereas the O C V C defect includes a carbon vacancy. The 3O C defect can be considered in two ways, either as a V C with the C-dangling bonds substituted by 2-coordinated O-atoms, removing any unsaturated bonds from the systems, and with them any vacancy-like artefacts. Alternatively, the defect can be viewed purely in substitutional terms with three neighboring C-atoms substituted by O removing the central C-site in the process. Out of the two options, which is the most appropriate depends upon the mechanism of formation/oxidation, and the stage of growth at which the defects are introduced. Direct consideration of which is beyond the scope of this work, but it is important to note the varying ways of viewing the same defect, as in V C terms this represents a passivated defect center, whereas in substitutional terms it would not. The calculated f defect E and optimized structures for the lowest energy configurations are presented in Figure 5a-d for each oxygen defect. For all the oxygen defects, distinct regimes are present; one bulk-like where there is little or no influence from the presence of the surface (blue in Figure 5a (Figure 6a,b). This trend continues until defect position 10, when the elongated CO bond breaks and produces the 2-coordinated O and the C-dangling bond configuration (as observed at defect site 12 in Figure 5e). The accommodation of this is facilitated by the disruption of the π-system in the near surface regime and the favorable incorporation of the 2-coordinated O-center. These factors energetically compensate for the formation of the C-dangling bond and reduce the number of O-electrons in π * -states. In common with the V C defect, the dramatic decrease in f defect E (−1.9 eV) would be suggestive of an increased defect concentration in the surface and near surface regions (with the caveats highlighted above) with respect to the bulk ( Figures S1d and S4, Supporting Information). While the magnitude of the drop in f is 0.02 eV, which is much lower than 3.24 eV for the V C at the same position. This suggests that at these chemical potentials, the HC surface is unstable to oxidation, in good agreement with experiment, [7,36] and one of the vectors of surface oxidation could be expected to occur via the formation of the O C .
The as the defect moves from the bulk to the surface (≈3.0 eV) is observed. Figure 5b shows the same cross trends as seen for N C V C that can be deconvoluted by defect configuration, as shown in the labeled subplots in Figure S14 (Supporting Information). What is interesting for the O C V C defect in contrast to the other considered here, is that the tilted defect configuration at the position directly below the surface is the lowest energy configuration (position 11, as shown in Figure 5f). In this defect position, both portions of the defect are accommodated with minimal strain, with the 5-membered ring showing negligible distortion, resulting in five equivalent CC bonds (1.46 ± 0.03 Å) (Figure 6c,d). The C 5 O ring likewise is incorporated with little distortion (<0.01 Å deviation). In general, a greater O C V C concentration would be expected in the near surface and surface regimes (Figures S1e and S4, Supporting Information). The exception to the general trend is position 11, as described above, the formation energy is f defect E = −0.33 eV at these chemical potentials, which in common with the other oxygen defects suggest that HC is unstable to oxygen.
For the 2O C defect (Figure 5c), the bulk-like region can be seen to show a slight decrease in energy with a reduction of 0.25 eV between defect positions 1-7. In the near surface regime, the decrease in f defect E is far more pronounced with ≈−2 eV reduction between positions 8 and 12. In the bulk-like regime, the for the 3O C defect ( Figure 5h) show a smaller gradual downtrend, 0.17 eV from positions 1 to 7, but is energetically more stable than all the other oxygen defects. Approaching the surface, the energy decrease from positions 8 to 13 is 2.5 eV. The negative f defect E confirms that these oxygen-containing defects are energetically probable, in agreement with the experimental observations [7,15,39] and can be understood with reference to the substitutional oxidation of HC. Initially, a single O replaces a C-atom, producing the O C defect and a C-dangling bond, this dangling bond can then be oxidized to form CO/ CO 2 [23] introducing a second C during an intermediate step. This leaves O in its preferred 2-coordinated configuration with no dangling bonds, and confirms that at these chemical potentials,  Figure S8 and Table S2 (Supporting Information) for the N C defect, Figures S10-S12 and Table S3 (Supporting Information) for the N C V C defect.
the HC surface is prone to both oxygen heteroatoms and oxygen vacancy defects. From experimental evidence, oxygen heteroatoms are beneficial for improving the wettability of HC, and oxygen defects can enhance metal adsorption. [1,7,39,[62][63][64][65] In the bulk-like regime, the O C defects sit at a three coordinated C-lattice site (Figures 5a and 6a, positions 1-9, and Figure S13 (Supporting Information)). As the surface is approached, the CO bond perpendicular to the surface breaks  Figure 1. e) Optimized O C defect at defect site 12 structure, f) optimized O C V C defect at defect site 11a structure, g) optimized 2O C defect at defect site 12 structure, and h) optimized 3O C defect at defect site 13 structure. Carbon atoms are brown spheres, and red oxygen. (Figure 6a), leaving a C-dangling bond. In the bulk-like region, there is a more marked distortion in the lattice than was observed for the N C defect at the same positions ( Figure 6b and Table S4 (Supporting Information)). The CO bond lengths are  Table S4 (Supporting Information) for the O C defect, Figure S17 and Table S5 (Supporting Information) for the O C V C defect, Figure S19 and Table S6 (Supporting Information) for the 2O C defect, and Figure S21 and  Figure 6b). As the defect approaches the surface (Figures 5c  and 6a,b positions 7-9), the threefold symmetry with equal CO bond length is distorted as the CO perpendicular to the surface begins to extend. This extension results in the breaking of the perpendicular CO bond, positions 10-13, resulting in a 2-coordinated O and a C-dangling bond (Figure 6a). The 2-coordinated CO bond is significantly shorter (1.39-1.40 Å) with the broken CO bond relaxing outward to a separation of ≈2.35 Å (Figure 6b). The surface positions result in a rupture of the bond linking the two graphene sheets, resulting in a C-dangling bond projecting from one of the sheets, and an O-capping of the other.
The structural properties of the O C V C defect are similar to the V C defect as it relaxes to give a 5-membered C ring, the C-dangling bond of the V C is replaced by the 2-coordinated O-capped C 5 O ring, leaving no unsaturated bonds in the defect (Figure 6c and Figures S10 and S11 (Supporting Information)). Geometrically, the 5-membered ring is more easily accommodated in the O C V C than the V C with a CC separation of ≈1.8 Å (for V C ≈ 2.0 Å) in the bulk, which decreases to ≈1.5 Å at the surface, as observed for the V C defect (Figure 6d and Figure S17 and Table S5 (Supporting Information)). As there are no unsaturated bonds in the system, the cross-linking observed for the V C at strained edges is not seen in the O C V C defect case. In common with the N C V C defect, there are both perpendicular and tilted variants (Figure 6c), with the same a or b type configurations where the O is positioned at the top or bottom site, respectively. For the a and b type perpendicular configurations, there are no significant distortions in the bulk regime (Figure 6d), whereas a small stretching of the bonds in the direction of the surface in the near surface regime is observed ( Figure S17 and Table S5, Supporting Information).
The 2O C defect is made up of two O C defects sitting at neighboring sites, with an OO separation of ≈2.3 Å in the bulk-like regime, which increases as the surface is approached to give a maximum OO separation of ≈2.5 Å when the defect is directly at the surface (Figure 6e,f and Figures S18 and S19 and Table S6 (Supporting Information)). For the 2O C defect, two configurations are possible with the defect sitting perpendicular to the surface (Figure 6e even), or at a 60° angle with respect to the surface (Figure 6e uneven). In the bulk context, the defect has twofold symmetry with each C 5 O ring equivalent. In the perpendicular case, all CO and CC bonds form equivalent pairs (no distortion within the C 5 O rings observed due to position of a given atom with respect to the surface). The CO bonds show no distortion and a negligible extension (0.01 Å) as the defect is moved from the bulk to the surface (Figure 6f and Table S6 (Supporting Information)). The CC bonds perpendicular to the surface show a similar small extension (0.01 Å) as the surface is approached. A greater relaxation is seen for the final pair of CC bonds opposite the CO, as the surface is approached (Figure 6f) and the OO separation increases the COC flattening with the O sitting of opposite sides of the graphene sheet (defect position 12, Figure S18, Supporting Information), which necessitates the extension of the CC bonds to accommodate (≈2%). This further reduces the f defect E for the 2O C defect at the surface. For the tilted configuration, the picture is more complex as the bond pairs within a C 5 O ring become inequivalent as the surface is approached (Figure 6f). In the bulk context uneven numbers (positions 1-7 Figures 5c and 6f), there is no meaningful impact from the presence of the surface and the defect is symmetric, with twofold symmetry as described above. In the near surface and surface regimes, distortions break the symmetry of the defect as the bonds (CO and CC) closer to the surface extend by 1-3% depending upon position.
The 3O C defect has threefold symmetry with three equivalent O-atoms, the OO separation is 2.58 Å, with CO bond lengths of 1.38 Å (Figure 6g,h). In common with the other defects considered here, as the 3O C is moved toward to the surface, there is an elongation of the bonds perpendicular to the surface (Figure 6h and Figure S20 (Supporting Information)). The OO separation extends to 2.78 Å (2.68 Å for the in-plane O), and the CO bonds become inequivalent with the bond perpendicular to the surface showing a small extension (≈0.02 Å) to 1.41 Å, which is compensated for by an equivalent contraction in the neighboring CC bonds of the C 5 O ring ( Figure S21 and Table S7, Supporting Information).

Nitrogen and Oxygen Defects
From the analysis above, both oxygen and nitrogen defects (or dopants) are expected to be present from their generally low f defect E . From experimental evidence, it has also been observed that extrinsic oxygen defects can arise as a consequence of nitrogen doping. [68] In this section, we investigate two nitrogenand oxygen-containing defects, the N C O C and N C 2O C defects (Figure 7). These defects are geometrically similar to the 2O C and 3O C defects, respectively. The number of possible configurations for the N C O C defect case increases significantly as compared to the previously discussed defects. In common with the 2O C defect, the N C O C defect can sit either perpendicular or tilted with respect to the surface. Additionally, either N (type a) or O (type b) can sit closest to the surface (Figure 7c). This in effect breaks the twofold symmetry that was present for the 2O C defect and requires each configuration to be considered explicitly. The N C 2O C defect represents a variation of the 3O C where one of the O (which was originally a C atom) atoms has been replaced by N. In the same way as the 3O C , this can be viewed as being a purely substitutional defect or based upon a pre-existing V C . The introduction of the N-atom breaks the threefold bulk symmetry of the 3O C (in a similar way to N C O C vs 2O C ) while shifting the principal symmetry axis from perpendicular to in the plane of the graphene sheet (Figure 7d).
The f defect E of the N C O C defect at different positions and configurations are shown in Figure 7a. While the bulk and surface regimes seen for the 2O C defect are present, the distribution of energies is broader with a number of cross trends. To understand these trends, it is instructive to deconvolute the f defect E for the different defect configurations outlined above and shown in Figure 7c and Figure S22 (Supporting Information). The energetic differences of these become more pronounced as the surface is approached, with the symmetric even geometry www.small-journal.com  Figure 1. c) The inequivalent N C O C defect configurations where i) is the type a at uneven defect position configuration, ii) type b at uneven defect position configuration, iii) type a at even defect position configuration, and iv) type b at even defect position configuration in (a). d) The inequivalent N C 2O C defect configurations where i) is the type a at uneven defect position configuration, ii) type b at uneven defect position configuration, iii) type a at even defect position configuration, and iv) type b at even defect position configuration in (b). The optimized structure of the defects at position 11a are shown in e) for N C O C and f) for N C 2O C as an example. All optimized structures of these defects are included in Figures S23, S24, S27, and S28 (Supporting Information). Carbon atoms are brown spheres, gray nitrogen, and red oxygen. lower in energy than in its tilted counterpart at a given depth (Figure 7c and Figure S22 (Supporting Information)). The f defect E for the surface defects give support to the experimental observations that edges are susceptible to nitridation in a similar manner to oxidation, with sites at defect position 10 and above giving negative f defect E at these chemical potentials. In contrast to the 2O C defect, there is a strong location dependence, with the negative f defect E only found in the surface and near surface regimes (Figures S1h and S5, Supporting Information).
The distinct f defect E profile is also clearly visible for N C 2O C (Figure 7b associated with the surface. The bulk energies show a shallow decrease of only 0.1 eV, with 2.9 eV reduction in the surface regime. The same cross trends previously described and due to the presence of the defect configurations shown in Figure 7d give rise to this and can be deconvoluted by separating the configurations ( Figure S26, Supporting Information). There is a hint of a secondary trend along the lines of the uneven-even effect previously described, that can be ascribed to the direction of the defect (Figure 7d), although it is important to note that the energy changes here are very small (∆E f = 0.04 eV) and not as significant as those observed for N C O C .
For the optimized N C O C structures, the N and O atoms deviate from the graphene plane ( Figures S23 and S24, Supporting Information), with the lowest energy configuration placing them on opposite sides of the C-sheet (Figure 7e). The bulk-like configuration has a NO separation of ≈2.3 Å in both arrangements (Figure 8a). The NC bond is ≈4% shorter than the NO, incorporated into the lattice with no appreciable distortion, and is less susceptible to any deviations as the surface is approached (Figure 8a). For the C 5 N and C 5 O units in the perpendicular variant, the bonds are symmetric, forming equivalent pairs as with the 2O C . The tilted variant shows small extensions of the bonding in the direction of the surface (after position 8 in Figure 8a and Figures S23 and S24 (Supporting Information)). In the bulk regime, there is no geometric distinction between the type a and b configurations (Figure 8a and Figure S25 and Table S8 (Supporting Information)), with the absolute location and perpendicular versus tilted configuration being better predictors of the geometric distortions.
Examining the geometric structure of the N C 2O C defect in the bulk-like regime, the different defect configurations are geometrically equivalent (Figure 8b), with an OO separation of 2.60 Å and a NO separation of 2.58 Å. The CN bond is 1.34 Å in common with the other 2-coordinated N-centered defects described here. Likewise, the CO bond is 1.38 Å, again in agreement with the other 2-coordinated O-defects ( Figure S29 and Table S9, Supporting Information). In common with the other defects presented here, as the surface is approached, the bonds perpendicular to the surface extend (Figure 8b), with the OO separation increasing to 2.9 Å at the surface, while the NO increases to 2.65 Å at the same position (12). The changes in the NO and CO bonds are far less pronounced at 0.01 and 0.03 Å, respectively. These extensions are accommodated by matching contractions in the CC bonds of the C 5 O/C 5 N rings (Figure 8b). as a function of position due to it being easily accommodated at the C-lattice site with negligible deviation. These results accord well with the experimental observations, describing an increased defect concentration at the edge sites, along with the observed instability to oxidation that is well-known for HC materials. [1,7,39,62,69] The results presented here allow an understanding of the relative defect concentrations between the bulk and the surface to be understood, along with absolute defect concentrations and defect probability distributions in the case of samples that have reached thermodynamic equilibrium (or that it is a reasonable assumption to make). Hence, it is of importance to not only consider defects on planar motifs (such as graphene and the basal plane as we have considered previously) but also on curved morphologies and surfaces to optimize material performance.

Adsorption of Li, Na, and K on V C , N C , and O C Defects
During electrochemical cycling of HC anodes, the metal incorporation can follow different mechanisms depending on carbon morphology, alkali metal cation, and site saturation. A Interatomic distances for the optimized structures of these defects are included in Figure S25 and Table S8 (Supporting Information) for the N C O C defect, and Figure S29 and Table S9 (Supporting Information) for the N C 2O C defect.
lot of debate remains in the literature about the HC charge/discharge mechanism. Initially, the so-called "falling cards model" was introduced by Stevens and Dahn [4] for Na + storage in HC. Following this model, the initial sodiation (sloping region) was attributed to Na + intercalation in the graphitic stacks, whereas the low voltage plateau region was attributed to Na + storage in nanopores. In recent years, this initial model has been extended and refined, also taking into account the effect of heteroatoms and the fact that not all HCs have the same atomic structure. [9,68,[70][71][72][73] Subsequent studies suggested that the sloping region could also be due to metal ion adsorption, and bonding to carbon vacancies. [7,33,39,65] In general, intercalation, carbon vacancy adsorption, heteroatom site adsorption, and pore filling have all been proposed to contribute to the charge/ discharge mechanisms for the HC anode materials. [7,33,39,65] In situ TEM studies showed that the initial sodiation of a HC anode occurs through surface active site adsorption, followed by intercalation into the graphitic layers. [50] The surface adsorption can also be followed by studying the anode expansion, due to the lower volume expansion associated with adsorption when compared to intercalation. [36,50] The volumetric change of a HC NIB anode was studied by Wang et al. showing that the initial sodiation process took place mainly through surface adsorption (due to the stable volume), whereas intercalation occurred after 300 s where there was a sharp rise in volume. [50] In situ electrochemical dilatometry further confirmed that the anode thickness changes as a function of sodiation time, and is further dependent on the microstructure, pyrolysis temperature, and heteroatom concentration of the HC material. [36] The surface of the anode has been shown to be important, with surface defects (especially when there is a significant oxygen concentration identified from XPS measurements) responsible for providing additional initial metal ion storage sites, and irreversible metal trapping. [7,38,39,74] The effect of surface, interlayer distance, pore size, carbon morphology, defects, and heteroatoms will be dependent on pyrolysis/carbonization temperature (with less surface area, fewer defects, and more graphitic character observed with increasing pyrolysis temperature) and the initial reagents (where oxygen and nitrogen dopants can be introduced to tune material properties). [21,27,36,39,75] Additionally, the main contribution to Li storage in HC comes from intercalation, with initial contribution from surface adsorption as confirmed by in situ Raman analysis of Li in HCs. [33,65,[76][77][78][79] For KIBs, the graphitic stacks' interlayer distances accessible to Na + and Li + storage can be inaccessible, with greater interlayer distances required. [33,61,66,80] Previously, we have studied, in isolation, basal plane defect adsorption, and intercalation in planar graphitic pores with varying interlayer distances (c) as guided by experimental HC characterization. [33,34] From these studies, we showed that metal adsorption is greatly enhanced at defect sites, [34] and that especially the sodium and potassium intercalation is heavily dependent on c, with potassium showing energetically favorable intercalation energies (i.e., negative binding energies) first at c > 3.85 Å, and sodium at c > 3.49 Å. [33] The latter were further confirmed by muon spin rotation spectroscopy. [21] Employing the simulation model described herein, the effect of curvature, graphitic stack interlayer, and defect location on metal adsorption will be explored. This gives important insight adding to the previous knowledge of metal adsorption on defective planar basal plane surfaces, metal intercalation in nondefective graphitic stacks, and metal adsorption on curved motifs. [33,34] Inspecting the binding energy of Li, Na, and K in between planar graphitic layers as a function of c as presented in our previous work [33] at c between 3.3 and 4.0 Å, the metal binding energy as a function of c decreases near linearly. As c increases, Li, Na, and K shift from intercalation to surface adsorption and converge to the graphene metal adsorption energies. [33,34] Performing a linear regression ( Figure S30, Supporting Information) shows that the increase in energy for each Å (in this range of c) is 6.5 eV for K, 2.7 eV for Na, and much lower for Li at 0.59 eV, as Li from the simulated binding energies shows energetically favorable intercalation at all investigated c. Here, we extend our previous studies to investigate how these trends are impacted by the addition of simple point defects (carbon vacancies, oxygen, and nitrogen heteroatoms), to gain insight into the surface adsorption and intercalation mechanisms. The metal adsorption energy (E ads ) at different defect positions was calculated according to Equation (2 Here, E Metal-Defective is the total energy of the defective system with an added metal atom, and μ Metal is the chemical potential of the metal species. The chemical potential of the metal is here taken as the total energy of a single metal ion in a vacuum cell (20 Å × 20 Å × 20 Å). [7,33,34] Hence, the more negative the E ads , the stronger the adsorption of the metal to the carbon system. The resulting adsorption energies are plotted in Figure 9 as a function of defect position (Figure 1), with optimized structures included in the Supporting Information. From Figure 9, it is clear that in agreement with the previously observed trends for both intercalation and surface adsorption, Li, Na, and K behave differently. Hence, we will first discuss the metals on the defect systems individually, and then summarize the intermetal trends separately.
Li shows the strongest interaction with the V C containing carbon lattice (as seen from the most negative adsorption energies in Figure 9a). For this system, Li shows similar adsorption strengths at the bulk-like vacancies (defect positions 1-4) and on the surface vacancies (defect positions 9-12). Comparing the energies at defect positions 1-8, i.e., the intercalation-like region, the Li intercalation energies are much improved as compared to Li intercalation in pristine planar graphitic layers. The Li adsorption energies for V C at defect positions 5-8 are weaker, signifying a transition from surface adsorption to intercalation, where it can be expected that the mechanisms compete. The optimized structures of Li and the V C defects are included in Figure S31 (Supporting Information). The transition (in terms of Li adsorption energies) between the intercalation and surface adsorption behavior is not as clear for the N C and O C defect containing systems. For the N C defect ( Figure S32 For Na, three clear trends are observed for each defect in Figure 9b. In defect positions 1-5, Na is intercalated, and at positions 9-12 adsorbed on the surface. Only for Na at V C defect ( Figure S34, Supporting Information) sites is the sites (and hence transition) between these regimes energetically favorable, although at a very weak adsorption energy (−0.02 eV). This furthers the analysis seen above for Li adsorption, that carbon vacancy defects are highly favorable for metal intercalation. For Na at defect sites 6-8, the Na adsorption energy is positive for both the N C and O C systems. These are the narrowest sites, agreeing with our previous results of Na intercalation in graphitic stacks, with Na intercalation becoming energetically favorable at c = 3.49 Å. The V C defect induces distortions to the carbon lattice, expanding the interlayer distance locally to 4 Å, allowing Na to intercalate. The same structural distortion is observed for Na at N C and O C defect site 8 ( Figures S35 and S36, Supporting Information), but the energetic cost for this distortion is higher, resulting in positive adsorption energies. Combining the knowledge gained above that surface defects have lower defect formation energies, and hence expected to be present in higher concentrations under equilibrium conditions, with the strong Na adsorption energies at the surface defect sites, these results suggest that for sodiation, surface adsorption at curved surface lattice sites forms an important part of the sodiation mechanism. As a general rule, the simulated adsorption (or binding) energies can be directly compared to cell voltage (V) assuming V = (−E ads ). The cutoff voltage for NIB HC anodes is typically 2 V versus Na + /Na, which means that any E ads stronger than −2 eV could lead to metal trapping, leading to irreversible capacity loss. [7,15] To this end, the V C and O C surface defects would be expected from these simulations to lead to irreversible capacity loss. Experimental studies have tested to what extent the cell voltage cutoff affects the cycling performance of Na/HC half cells. [15] Charging the half-cell (made of HC anode containing significant oxygen concentrations and defects) to 3 V for a duration of 2 h led to the release of an additional 100 mAh g −1 specific capacity, suggesting that for HC anodes, careful tuning and employment of cutoff voltage may be necessary when they contain these defects. [15] For K adsorption (Figure 9c), the picture is more complicated. For the V C and O C defect containing systems, K adsorption energies are negative, with V C showing the strongest K adsorption. As opposed to Na and Li, K remains in the bulk (as described in Figure 1) in the systems with these defects at all defect sites except 9-12 ( Figures S37-S39, Supporting Information). Hence, discussing the K adsorption as a function of interlayer distance is for K disingenuous as K remains at the wider interlayer distances. The weak K adsorption energies are further due to the high structural distortions (especially in the N Cand O C -containing systems) induced by the metal. For defect sites 9-11 (which show the strongest K adsorption energies in Figure 9c), K adsorbs above the interlayer space. Inspecting the optimized structures for K in the systems with V C and O C at defect site 10 (with both the strongest K adsorption), K in both these sites interacts with a dangling carbon bond at the strained lattice site ( Figures S37 and S39, Supporting Information). Hence, it can be deduced that these defect sites are particularly important for potassiation, as well as sodiation. This dangling bond is not opened up by the N C defects ( Figure S38, Supporting Information), hence not giving the same enhancement to the K adsorption.
From the above analysis, it is clear that both O C and N C defects enhance the metal surface adsorption, but that the V C defect leads to the strongest intercalation and surface interactions with all the metals. Generally, as discussed previously, [33,34] the energetic ordering of Li, Na, and K is different in the intercalation regime (E ads, Li < E ads, Na < E ads, K ) than in the adsorption regime (E ads, Li < E ads, K < E ads, Na ). This ordering, with Na having stronger interaction with the carbon lattice when in the intercalated regime than K, is due to the larger ionic size of K. [33] K is also the only metal that does not break the sp 3 -hybridized carbon interconnect formed when adding a V C defect to lattice site 8. The Na sits within an interlayer distance of 4 Å, and Li 3.78 Å, both at the defect site, whereas K remains in the bulk regime 7.89 Å away from the sp 3 interlink and the defect site, remaining close to site 1, as defined in Figure 1. Hence, it would be expected that in NIBs, and LIBs, these kinds of sp 3 interlinks would be broken during lithiation and sodiation. For the metals at defect sites 9 and 10, the metals sit above the gap between the carbon sheets. This corresponds to the initial intercalation step. At defect sites 11-13, the metals are adsorbed at the surface. Both Li and K sit closest to the carbon surface (1.03 Å for Li, and 1.12 Å for K) at defect position 13, whereas Na is the closest when V C is at position 11 (1.25 Å). Na adsorbed at V C at site 10 is also the strongest adsorption site for Na, and binds stronger to the surface than K at this site. For Li, no clear preference can be seen for intercalation or surface adsorption, indicating that both these will be equally important in the lithiation mechanism. This is in agreement with previous experimental studies of the lithiation mechanism of HC and other carbon anodes such as graphite. Cell tests comparing the same HC materials for LIBs, NIBs, and KIBs were conducted in a previous publication. [33] The study utilized a series of HC materials synthesized by hydrothermal carbonization from glucose (containing oxygen heteroatom defects, but no traceable amount of nitrogen from XPS analysis) and pyrolysis at different carbonization temperatures. From the charge-discharge curves of the HC materials with the most graphitic character (narrowest interlayer distances at 3.5-3.6 Å, respectively, and the lowest oxygen concentration), it could be observed that both the Li and Na cells exhibited higher specific capacities than the K cells due to the hindrance of the K + ions to intercalate in the smaller interlayer distances (as observed in Figure 9c). Comparing the Li and Na capacity showed that Li had higher specific capacities than Na, agreeing with the trend presented in this work (Figure 9a,b). For the HC materials carbonized at lower temperature with wider interlayer distances (>3.7 Å), the lithiation remained more favorable with higher capacities than both sodiation and potassiation (with the sodiation and potassiation showing similar behavior). This agrees with our analysis that the interlayer distance and defect location have larger impact on the metal incorporation than the type of defect does. The HC materials carbonized at lower temperature also showed more sloping region capacity, which is typically coupled to metal adsorption at defect sites as well as intercalation. This indicates that the presence of defects does lead to higher metal storage (by opening up additional surface adsorption sites), in agreement with our DFT simulations. For optimum performance, there needs to be a balance between accessible bulk intercalation sites, surface defect sites for the initial metal adsorption, and that the defects should not bind the metals too strongly as to hinder the further metal incorporation mechanisms.
From these results, it is clear that Li is the least sensitive to the layer separation and shows little deviation between the intercalation and adsorption regimes. Intercalation becomes less favorable in the near surface region as a result of the reduced layer separation. The energy penalty accords well with the previous defect free study, [33] with the absolute binding energy dictated by the defect type. In the case of Na, the picture is more complex as a result of its larger ionic radius, and weaker binding, which combine to give a much greater sensitivity to layer separation. However, the same trends described for Li can be observed for Na, with the absolute magnitude of the adsorption energy dictated by the defect type. The greater sensitivity to layer separation leads to a much larger energy penalty in the near interface regime, in accordance with the defect free adsorption energy penalty as with Li. In the case of K, the picture becomes less clear as c for the HC characterized here is too narrow to favorably accommodate K. As a result, for K, no stable minima could be identified. This makes it impossible to comment on the energetic penalty due to layer separation for K in these models, although it is still possible to see the shift from intercalation (albeit lacking the c distribution) to adsorption.

Conclusion
In this paper, a systematic DFT study of the effect of defect location on defect formation in nanoporous carbon anode materials was conducted. The defects considered were based on carbon monovacancy, nitrogen substitutional heteroatom defects, oxygen substitutional heteroatom defects, and combinations of these. The defects were studied at different locations in three general lattice regions: bulk-like, near surface, and surface. From this analysis, a strong dependence of defect position was identified, with the defect formation energy dramatically increasing as the surface is approached (as compared to the bulk) due to the increased curvature and strain. To understand what contributions the defect location in respect to the carbon matrix interlayer distance, and surface have on the initial lithiation, sodiation, and potassiation mechanisms, we further calculated Li, Na, and K adsorption at the different single center defect sites. From these simulations, it was shown that the defect (at each defect site) that leads to the strongest metal adsorption is the carbon monovacancy. Furthermore, it was shown that surface adsorption is highly energetically favorable when the O C defect at a strained lattice site leads to a dangling bond, and that careful consideration of the cell voltage would be an important engineering tool to minimize capacity loss due to metal trapping at high adsorption energy defect sites. However, the location of the defect, in terms of intercalation in varying interlayer distances and on surface, was of more importance than the defect. From these simulations, these surfaces have high defect concentrations and disorders, something that needs to be considered when studying both the lithiation/sodiation/potassiation mechanisms, and electrolyte anode interfaces.
Finally, defects can be both beneficial and detrimental to alkali metal ion battery performance. As shown from this and previous work, [7,34,36,81,82] defects can increase the adsorption of Li, Na, and K in these anode materials, both at the surface and in the bulk through intercalation, and thus would enhance the initial lithiation, sodiation, and potassiation, respectively. Heteroatom doping with O and N should be explored either by choosing reagents with high oxygen or nitrogen content when synthesizing these HC compounds, or by heat treatment in an oxidative environment. However, these strong metal adsorption energies at the defect sites could lead to metal trapping, limiting the cycling performance. Hence, from a material design perspective, defect engineering should be considered and explored as a route to