A First-Principles Study on the Multilayer Graphene Nanosheets Anode Performance for Boron-Ion Battery

Advanced battery materials are urgently desirable to meet the rapidly growing demand for portable electronics and power. The development of a high-energy-density anode is essential for the practical application of B3+ batteries as an alternative to Li-ion batteries. Herein, we have investigated the performance of B3+ on monolayer (MG), bilayer (BG), trilayer (TG), and tetralayer (TTG) graphene sheets using first-principles calculations. The findings reveal significant stabilization of the HOMO and the LUMO frontier orbitals of the graphene sheets upon adsorption of B3+ by shifting the energies from −5.085 and −2.242 eV in MG to −20.08 and −19.84 eV in 2B3+@TTG. Similarly, increasing the layers to tetralayer graphitic carbon B3+@TTG_asym and B3+@TTG_sym produced the most favorable and deeper van der Waals interactions. The cell voltages obtained were considerably enhanced, and B3+/B@TTG showed the highest cell voltage of 16.5 V. Our results suggest a novel avenue to engineer graphene anode performance by increasing the number of graphene layers.


Introduction
The past two decades have witnessed impressive improvements in lithium-ion battery (LIB) technologies for portable consumer electronics, electrical devices, and energy storage [1,2]. LIBs have demonstrated outstanding performance, including superior energy density, operating voltage, life cycle, and minimal rate of self-discharge, as well as low volume [3]. Generally, LIBs have exhibited exceptional performance over other known rechargeable ion battery systems [4][5][6]. More significantly, the 2019 Nobel Prize award in Chemistry was received by researchers in the field of lithium-ion batteries [7,8]. Despite the appealing applicability of LIBs, limitations such as narrow lifetime, inadequate performance at low temperatures, and most severely, the swift exhaustion of the lithium mineral reserves, may represent a setback for LIB technology. For example, around one-quarter of the global Li precursor materials are utilized to produce LIBs, thereby delivering a strong hike in the price of lithium carbonate [9][10][11]. On this premise, a wide range of research has been dedicated to alternative elements for ion battery technology [2,8,[12][13][14]. Among many candidates, sodium-ion batteries (SIBs) are expected to substitute the LIBs due to their low cost, non-toxicity, and nearly limitless sodium mineral reserve [15,16]. Moreover, sodium possesses similar chemical properties as lithium since they belong to the same alkali metal Nanomaterials 2022, 12, 1280 3 of 16 superior electrochemical performance for ion battery anodic applications. Therefore, we investigated the electrochemical performance of single-layer to four-layer graphene sheets. With this investigation, we aim to pave the way for the successful design of extremely effective materials for energy storage.

Theoretical Methods
All structural optimizations and electronic properties calculations were performed employing Density Functional Theory (DFT) as implemented on Gaussian 09 suite. The Perdew-Burke-Ernzerhof (PBE) functional belonging to the generalized gradient approximation (GGA) functional was used to account for the exchange-correlation energy because it provides reasonable accuracy without prolonged computational times, while the 6-311 G(d,p) basis set was adopted [8,[59][60][61]. The PBE is an effective standardized functional because, by design, each component adheres to some exact conditions [62]. It follows the spin-scaling relationship exactly and reclaims the linear response of the LDA for a small-scale gradient [63]. The PBE functional is a full ab initio functional which relies on µ, β, K, and γ parametric values fixed from theoretical factors. Similarly, the 6-311 G(d,p) basis set was adopted because it provided a reliable assessment of the energies of solvation when employing an implicit solvent standard/prototype when compared to other highly revered basis sets [64]. The influence of dispersion-corrected functionals on the geometry of the layers was similarly examined using Grimme's D3 approach [65][66][67][68]. Typically, all atoms were set free in the single-and double-layer slab. The resulting layers slab models are shown in Figure 1. However, to stabilize the computed surface energies for the three and four layers, the atoms were frozen. Studies have shown that freezing layers does not affect surface energies [59]. The adsorption energies (E ad ) were calculated according to the following equation [69]: where E complex , E adsorbate , and E substrate represent the energies of the absorbate and boron ion(s), graphitic layer(s), and the substrate, respectively. Likewise, the final cell voltages (V cell ) were calculated utilizing the Nernst equation: where zF represents the charge on the metal ions and Faraday constant, and ∆G cell stands for the Gibbs free energy of the overall cell reactions as where ∆G corresponds to the change in the internal energies of the cell, because the influences of entropy and the volume effects are insignificant (<0.01 V) to the cell voltage (V cell ) [1]. We investigated the effect of the number of layers and boron ions on the cell voltage. Our model of the systems with pristine layered graphene sheets is illustrated in Figure 1. The models of the boron-intercalated layered graphene sheets are shown in Figure 2. The GaussSum program was used to depict the density of states (DOS) plot [70].

Structural Properties
Pristine graphene nanosheets: The optimized electronic structures of the graphene sheet and the different multilayer orientations of the sheets are shown in Figure 1. In the multilayer orientations, the graphene sheets were differently orientated: AB and AA for bilayer; ABA and AAB for trilayer; and ABAB, AABB, and AABA for tetralayer before optimization. However, after optimization, it was observed that the layers favored AA orientation in all the arrangements, as similarly observed by other reports [8].

Structural Properties
Pristine graphene nanosheets: The optimized electronic structures of the graphene sheet and the different multilayer orientations of the sheets are shown in Figure 1. In the multilayer orientations, the graphene sheets were differently orientated: AB and AA for bilayer; ABA and AAB for trilayer; and ABAB, AABB, and AABA for tetralayer before optimization. However, after optimization, it was observed that the layers favored AA orientation in all the arrangements, as similarly observed by other reports [8].

Structural Properties
Pristine graphene nanosheets: The optimized electronic structures of the graphene sheet and the different multilayer orientations of the sheets are shown in Figure 1. In the multilayer orientations, the graphene sheets were differently orientated: AB and AA for bilayer; ABA and AAB for trilayer; and ABAB, AABB, and AABA for tetralayer before optimization. However, after optimization, it was observed that the layers favored AA orientation in all the arrangements, as similarly observed by other reports [8].
Boron-intercalated and boron-free graphene nanosheets: Figure 2 shows the optimized structures of the boron atom and its corresponding ion (B 3+ ) adsorbed on the monolayer graphene sheets (MG) (Figure 2a,b,i,j). The boron species preferred a central position on the MG after optimization even after being positioned at different points on the sheets before optimization. In the multilayered graphene orientations, the intercalated B or B 3+ preferred the edges of the multilayer sheets after optimization. The reduced density gradient (RDG) maps for the complexes imply weak van der Waals interactions predominantly among the boron species and the graphene sheets. It may offer explanations for the preferred arrangements of the complexes. As shown in Figure S1 in the Supplementary Materials, the TG and TTG have two boron atoms and/or ions (B 3+ ) alternatively intercalated within two layers of the graphene sheets.
We examined the distances between the layers of the graphene (α-layer and β-layer) and the boron species (B/B 3+ ). As shown in Figure 3, the red bars represent the distance or degree of compactness between the β-layer and the boron species. Similarly, the green bars represent the degree of compactness between the α-layer and the boron species (B/B 3+ ). The monolayer graphene arrangement results from the "α-layer"-B/B 3+ arrangement with the subsequent multilayer arrangements resulting from the successive addition of two or more graphene layers in the αand/or β-layered direction(s). The "α-layer"-B/B 3+ distances (the green bars) are generally shorter than the "β-layer"-B/B 3+ distances (the red bars), probably due to stronger interaction between the α-layer and the boron species before the subsequent addition of the β-layer(s) to generate the multilayers. The β-layer distances from the boron atoms (B) are also longer than their corresponding distances from the ionic boron species (B 3+ ), which further amplifies the ionic interaction effect of the boron with the graphene layer(s). Boron-intercalated and boron-free graphene nanosheets: Figure 2 shows the optimized structures of the boron atom and its corresponding ion (B 3+ ) adsorbed on the monolayer graphene sheets (MG) (Figure 2a,b,i,j). The boron species preferred a central position on the MG after optimization even after being positioned at different points on the sheets before optimization. In the multilayered graphene orientations, the intercalated B or B 3+ preferred the edges of the multilayer sheets after optimization. The reduced density gradient (RDG) maps for the complexes imply weak van der Waals interactions predominantly among the boron species and the graphene sheets. It may offer explanations for the preferred arrangements of the complexes. As shown in Figure S1 in the Supplementary Materials, the TG and TTG have two boron atoms and/or ions (B 3+ ) alternatively intercalated within two layers of the graphene sheets.
We examined the distances between the layers of the graphene (α-layer and β-layer) and the boron species (B/B 3+ ). As shown in Figure 3, the red bars represent the distance or degree of compactness between the β-layer and the boron species. Similarly, the green bars represent the degree of compactness between the α-layer and the boron species (B/B 3+ ). The monolayer graphene arrangement results from the "α-layer"-B/B 3+ arrangement with the subsequent multilayer arrangements resulting from the successive addition of two or more graphene layers in the α-and/or β-layered direction(s). The "α-layer"-B/B 3+ distances (the green bars) are generally shorter than the "β-layer"-B/B 3+ distances (the red bars), probably due to stronger interaction between the α-layer and the boron species before the subsequent addition of the β-layer(s) to generate the multilayers. The β-layer distances from the boron atoms (B) are also longer than their corresponding distances from the ionic boron species (B 3+ ), which further amplifies the ionic interaction effect of the boron with the graphene layer(s). The distance between the graphene layers (α-layer and β-layer) and the boron species (B and/or B 3+ ). The green bars represent the "α-layer"-B/B 3+ distance while the red bars represent the "β-layer"-B/B 3+ distance.
In the cases of arrangements having more than one boron atom/ion, similarly, distance trends and effects are pertinent. Generally, the presence of the boron atoms/ions intercalation influences the degree of compactness of the graphene layer arrangement; however, this may likely be insignificant to the progressive drifting of the boron atoms and/ions during the charging and discharging process of the ion battery, as illustrated in the RDG surface analysis. The distance between the graphene layers (α-layer and β-layer) and the boron species (B and/or B 3+ ). The green bars represent the "α-layer"-B/B 3+ distance while the red bars represent the "β-layer"-B/B 3+ distance.
In the cases of arrangements having more than one boron atom/ion, similarly, distance trends and effects are pertinent. Generally, the presence of the boron atoms/ions intercalation influences the degree of compactness of the graphene layer arrangement; however, this may likely be insignificant to the progressive drifting of the boron atoms and/ions during the charging and discharging process of the ion battery, as illustrated in the RDG surface analysis.

Electronic Properties of Graphene Sheets and Absorbed B 3+ on the Graphene Sheets
The energies of HOMO, LUMO, and HOMO-LUMO bandgap of the graphene sheets and the intercalated B 3+ ions are calculated in Figure 4a-g and Figures S2-S5 in the Supplementary Materials. The presence of higher HOMO energy is a characteristic associated with donating tendencies, while low LUMO energy is considered as accepting ability.
with donating tendencies, while low LUMO energy is considered as accepting ability.
The HOMO energies (Table S1 in Supplementary Materials) of the graphitic layers relatively decrease from monolayer to tetralayer, whereas the LUMO energies recorded no significant changes. This suggests that the HOMO and LUMO levels of graphene sheets were significantly stabilized by the adsorption of the B 3+ ion. The HOMO and LUMO levels shift from −5.085 and −2.242 eV in MG to −20.08 and −19.84 eV in 2B 3+ @TTG. Although there are no significant differences between the energies of the HOMO and LUMO of the singly adsorbed B 3+ ions, relative variations are seen with doubly adsorbed B 3+ ions: −20.08 and −19.84 eV for 2B 3+ @TG and −18.33 and −17.66 eV, corresponding to 2B 3+ @TTG, which can be attributed to the addition of B 3+ into the graphene sheets. Accordingly, upon adsorption of B 3+ , the global electronic energy bandgap (Eg) of the graphene sheets is reduced significantly to about 99.89%, and the same trend is observed with adsorbed neutral boron (Table S1 in Supplementary Materials). The observed order with adsorbed neutral boron could be due to the existence of unpaired electrons in the valence shell of the neutral boron; the HOMO level of the graphitic layer is largely impacted by altering to higher energies, suggesting a considerable destabilization [71]. Similarly, the Eg of the doubly adsorbed B 3+ increases in comparison with the singly adsorbed B 3+ , which could be attributed to the decreased LUMO level and increased HOMO level. In electrostatic potential maps (ESPs), colors are used as indicators for different electrostatic potential values: blue demonstrates high positive (electron-deficient) regions of the species, while green displays the region of zero potential [72]. As illustrated in Figures S6 and S7, the ESPs of the B 3+ form of B 3+ @MG, B 3+ @BG, B 3+ @TG, B 3+ @TTG_asym, B 3+ @TTG_sym, 2B 3+ @TG, and 2B 3+ @TTG reveal predominantly more electropositive regions (blue color). This is more severe with B 3+ @MG, B 3+ @BG, 2B 3+ @TG, and 2B 3+ @TTG configurations. However, with the addition of layers, a slightly neutral region with B 3+ @TG, B 3+ @TTG_asym, and B 3+ @TTG_sym was observed. The HOMO energies (Table S1 in Supplementary Materials) of the graphitic layers relatively decrease from monolayer to tetralayer, whereas the LUMO energies recorded no significant changes. This suggests that the HOMO and LUMO levels of graphene sheets were significantly stabilized by the adsorption of the B 3+ ion. The HOMO and LUMO levels shift from −5.085 and −2.242 eV in MG to −20.08 and −19.84 eV in 2B 3+ @TTG. Although there are no significant differences between the energies of the HOMO and LUMO of the singly adsorbed B 3+ ions, relative variations are seen with doubly adsorbed B 3+ ions: −20.08 and −19.84 eV for 2B 3+ @TG and −18.33 and −17.66 eV, corresponding to 2B 3+ @TTG, which can be attributed to the addition of B 3+ into the graphene sheets. Accordingly, upon adsorption of B 3+ , the global electronic energy bandgap (E g ) of the graphene sheets is reduced significantly to about 99.89%, and the same trend is observed with adsorbed neutral boron (Table S1 in Supplementary Materials).
The observed order with adsorbed neutral boron could be due to the existence of unpaired electrons in the valence shell of the neutral boron; the HOMO level of the graphitic layer is largely impacted by altering to higher energies, suggesting a considerable destabilization [71]. Similarly, the E g of the doubly adsorbed B 3+ increases in comparison with the singly adsorbed B 3+ , which could be attributed to the decreased LUMO level and increased HOMO level. In electrostatic potential maps (ESPs), colors are used as indicators for different electrostatic potential values: blue demonstrates high positive (electron-deficient) regions of the species, while green displays the region of zero potential [72]. As illustrated in Figures S6 and S7, the ESPs of the B 3+ form of B 3+ @MG, B 3+ @BG, B 3+ @TG, B 3+ @TTG_asym, B 3+ @TTG_sym, 2B 3+ @TG, and 2B 3+ @TTG reveal predominantly more electropositive regions (blue color). This is more severe with B 3+ @MG, B 3+ @BG, 2B 3+ @TG, and 2B 3+ @TTG configurations. However, with the addition of layers, a slightly neutral region with B 3+ @TG, B 3+ @TTG_asym, and B 3+ @TTG_sym was observed.

Adsorption of B 3+ Ions on the Wall of Graphene Sheets
The adsorption of B 3+ ions on the wall of the solid host is considered an alternative B 3+ -storage mechanism. The model sheets studied with adsorbed B 3+ ions are depicted in Figure 2. The role of the numbers of graphene sheets and the number of B 3+ ions are examined. We studied the adsorption characteristics of graphitic carbon sheets by intercalating B 3+ ions between the bilayer, trilayer, and tetralayer ( Figure 2). The positions of the intercalated B 3+ ions between the graphene sheets were selected thoroughly. All the B 3+ ions were placed in the central plane defined by the graphene sheet and allowed to relax in all directions.

Adsorption of B 3+ Ions on the Wall of Graphene Sheets
The adsorption of B 3+ ions on the wall of the solid host is considered an alternative B 3+ -storage mechanism. The model sheets studied with adsorbed B 3+ ions are depicted in Figure 2. The role of the numbers of graphene sheets and the number of B 3+ ions are examined. We studied the adsorption characteristics of graphitic carbon sheets by intercalating B 3+ ions between the bilayer, trilayer, and tetralayer ( Figure 2). The positions of the intercalated B 3+ ions between the graphene sheets were selected thoroughly. All the B 3+ ions were placed in the central plane defined by the graphene sheet and allowed to relax in all directions.
The adsorption energies (Ead) for the singly intercalated B 3+ ions were calculated as 3.987, −320.1, −477.0, −638.4, and −637.6 kcal/mol for B 3+ @MG, B 3+ @BG, B 3+ @TG, 2B 3+ @TTG_asym, and B 3+ @TTG_sym, respectively ( Figure 5). Furthermore, the tetralayer graphitic carbon B 3+ @TTG_asym and B 3+ @TTG_sym exhibited the most favorable and stable adsorption configuration of B 3+ intercalated within different layers of graphene sheets. Similarly, the Ead corresponding to the formation of 2B 3+ @TG and 2B 3+ @TTG configurations having double B 3+ ion intercalations are −467.9 and −628.7 kcal/mol, respectively. Thus, the addition of an extra sheet to form 2B 3+ @TTG was more energetically favored than the formation of the 2B 3+ @TG with fewer layers. Meanwhile, a very weak interaction was observed on B 3+ @MG with an Ead of 3.987 kcal/mol, in agreement with some research findings on Li [73]. The reduced density gradient (RDG) isosurfaces of the interactions are shown in Figures 6 and S8. Bader's theory of atoms in molecules categorized van der Waals interactions, strong steric effects, and hydrogen bond interactions to exhibit low ρ and relatively larger ρ values, respectively [74]. RDG isosurfaces plots reveal the binding interaction region and the modes [75]. The RDG isosurfaces analysis results of BC6 and BC2 mode of binding interactions manifest the weak van der Waals interaction between B 3+ and graphene sheets, as displayed in Figures 6 and S8, suggesting the physical nature of the reaction and high diffusion affinity of the B 3+ ions [76]. Except for B 3+ @MG, in which the intercalated B 3+ is equidistant with six carbons of the graphene ring to form BC6 (Figure 6a,b), The reduced density gradient (RDG) isosurfaces of the interactions are shown in Figure 6 and Figure S8. Bader's theory of atoms in molecules categorized van der Waals interactions, strong steric effects, and hydrogen bond interactions to exhibit low ρ and relatively larger ρ values, respectively [74]. RDG isosurfaces plots reveal the binding interaction region and the modes [75]. The RDG isosurfaces analysis results of BC 6 and BC 2 mode of binding interactions manifest the weak van der Waals interaction between B 3+ and graphene sheets, as displayed in Figure 6 and Figure S8, suggesting the physical nature of the reaction and high diffusion affinity of the B 3+ ions [76]. Except for B 3+ @MG, in which the intercalated B 3+ is equidistant with six carbons of the graphene ring to form BC 6 (Figure 6a,b), all other configurations revealed the BC 2 geometry, where the B 3+ ions are adjacent to two carbon atoms, one above and one below the graphene layer (Figure 6c-f and Figure S8). all other configurations revealed the BC2 geometry, where the B 3+ ions are adjacent to two carbon atoms, one above and one below the graphene layer (Figures 6c-f and S8). To gain in-depth knowledge of the interactions between the intercalated B 3+ /B atom and the numbers of graphene sheets, we calculated the density of states (DOS) and the projected density of states (PDOS) of the B 3+ /B atom adsorbed graphene sheets as displayed in Figures 7, S9-S12 . The projected density of states (PDOS) plots revealed the reduction in Eg from 2.843 eV in MG to 2.472 eV with TTG, which is likely due to an increase in the number of graphene sheets (Figure 7). Similarly, as seen in the PDOS and the frontier molecular orbital, the number of states at the HOMO is influenced by the carbon atoms, while the contribution at the LUMO increases with an increase in the graphene sheet. In general, the valence and conduction bands of MG, BG, TG, and TTG are dominated by the p orbitals of both types of carbon atoms (either bonded with C or H). Unlike the boron atoms, the adsorption of B 3+ significantly reduces the levels of HOMO and LUMO in the graphene sheets and limits their Eg, which would simplify the diffusion of electrons or charges. Analysis of the PDOS and frontier molecular orbital implies that the LUMO is localized on the B 3+ while the HOMO is confined on the graphitic sheets but improved with the addition of the sheets. Furthermore, the adsorption of two B 3+ into TG and TTG followed the same trend. However, upon adsorption of boron, a new electronic peak (indicated by the red arrow) is generated, which mainly comes from the contribution of 2p unpaired electrons of the B atom ( Figures S9a-c, S10c,d and S12c,d ). The presence of unpaired electrons at the HOMO generates unstable singly occupied molecular orbitals, resulting in the lowering of the HOMO from −5.085 eV in B@MG to 0.040 eV in 2B@TG (Table S1). In general, the PDOS and the frontier molecular orbital analysis illustrate that the contribution of the number of states at the HOMO is mainly influenced by the boron atoms, while the LUMO is controlled by the atoms of the graphene sheet. To gain in-depth knowledge of the interactions between the intercalated B 3+ /B atom and the numbers of graphene sheets, we calculated the density of states (DOS) and the projected density of states (PDOS) of the B 3+ /B atom adsorbed graphene sheets as displayed in Figure 7 and Figures S9-S12. The projected density of states (PDOS) plots revealed the reduction in E g from 2.843 eV in MG to 2.472 eV with TTG, which is likely due to an increase in the number of graphene sheets (Figure 7). Similarly, as seen in the PDOS and the frontier molecular orbital, the number of states at the HOMO is influenced by the carbon atoms, while the contribution at the LUMO increases with an increase in the graphene sheet. In general, the valence and conduction bands of MG, BG, TG, and TTG are dominated by the p orbitals of both types of carbon atoms (either bonded with C or H). Unlike the boron atoms, the adsorption of B 3+ significantly reduces the levels of HOMO and LUMO in the graphene sheets and limits their E g , which would simplify the diffusion of electrons or charges. Analysis of the PDOS and frontier molecular orbital implies that the LUMO is localized on the B 3+ while the HOMO is confined on the graphitic sheets but improved with the addition of the sheets. Furthermore, the adsorption of two B 3+ into TG and TTG followed the same trend. However, upon adsorption of boron, a new electronic peak (indicated by the red arrow) is generated, which mainly comes from the contribution of 2p unpaired electrons of the B atom ( Figures S9a-c, S10c,d and S12c,d). The presence of unpaired electrons at the HOMO generates unstable singly occupied molecular orbitals, resulting in the lowering of the HOMO from −5.085 eV in B@MG to 0.040 eV in 2B@TG (Table S1). In general, the PDOS and the frontier molecular orbital analysis illustrate that the contribution of the number of states at the HOMO is mainly influenced by the boron atoms, while the LUMO is controlled by the atoms of the graphene sheet. Nanomaterials 2022, 12, x FOR PEER REVIEW 9 of 17 Natural bond orbitals (NBO) population analysis was conducted to estimate the charge transfer in the complexes and is presented in Figures S13 and S14. As shown in Table 1, the complexes exhibit σ donation between the graphene sheets and the B 3+ ion, and π back-donation from the B 3+ ion to the vacant orbitals of the graphene atoms. The NBO results likewise indicate that the charge transfer is more sensitive to the sum of adsorbed B 3+ ions in 2B 3+ @TG and 2B 3+ @TTG than the other complexes, suggesting greater stability in the B 3+ @MG, B 3+ @BG, B 3+ @TG, B 3+ @TTG_asym, and B 3+ @TTG_sym configurations in contrast to the 2B 3+ @TG and 2B 3+ @TTG complexes.  Natural bond orbitals (NBO) population analysis was conducted to estimate the charge transfer in the complexes and is presented in Figures S13 and S14. As shown in Table 1, the complexes exhibit σ donation between the graphene sheets and the B 3+ ion, and π back-donation from the B 3+ ion to the vacant orbitals of the graphene atoms. The NBO results likewise indicate that the charge transfer is more sensitive to the sum of adsorbed B 3+ ions in 2B 3+ @TG and 2B 3+ @TTG than the other complexes, suggesting greater stability in the B 3+ @MG, B 3+ @BG, B 3+ @TG, B 3+ @TTG_asym, and B 3+ @TTG_sym configurations in contrast to the 2B 3+ @TG and 2B 3+ @TTG complexes. C1 and C2 represent the shortest distant carbon of the graphene sheet above or below, respectively, and B 3+ is the adsorbed boron ion, π* represent the anti-bonding π-orbitals.

Ion Battery Applications
The graphene sheets (single and multilayers) were considered as anodic electrodes for the boron-ion batteries with the half-cell reactions taking place at the cathode and anode shown in Equations (4)- (7). Equation (8) also shows the overall cell reaction.
Anode: B@sheet → B 3+ @sheet + 3e − Charging Process: B 3+ @sheet + 3e − → B@sheet Discharging Process: B@sheet → B 3+ @sheet + 3e − Overall cell reaction: B 3+ + B@sheet → B 3+ @sheet + B + ∆G cell (8) where B/B 3+ is the boron atom/ion and "sheet" is the monolayer graphene sheet (MG), bilayer graphene sheet (BG), trilayer graphene sheet (TG), and tetralayer graphene sheet (TTG). The cell voltage (V cell ) and the cell energy are computed according to Equations (2) and (3). As depicted in Figure 5, the adsorption energies of the boron atom(s) and/or ion(s) with the graphene sheets increased negatively as the number of the layers increased, with the boron ion (B 3+ ) adsorption being more negative than the atomic counterparts. The monolayer graphene sheet (B 3+ /B@MG) has the least negative adsorption energy, while the tetralayer graphene sheet (B 3+ /B@TTG) showed the highest negative adsorption energy. The complexes with more negative adsorption energies express more favorable and deeper van der Waals interactions [65,77]. Table 2 shows the overall cell energies (∆E cell ) and the corresponding voltages of the cells. The cells' energy increases with an increasing number of graphene sheets. However, for the intercalation of two atoms/ions of boron, the overall cell adsorption energy is higher. This can be attributed to the increased weak interactions between the graphene sheets and the two boron ions, which still permit free movement of the boron ions within the graphene layers during the charging and discharging processes. Moreover, compared to sodium and lithium-ion batteries (SIBs and LIBs) [1] and other battery materials reported (Table 3), the voltages obtained here are significantly improved. The B 3+ /B@TTG showed the highest cell voltage of 16.5 V, which indicates that in addition to the complexes being efficient alternatives for the anodic electrode of boron-ion batteries, increasing the number of the graphene layers also improved the electrochemical performance of the storage system. Generally, the boron-graphene sheets have shown excellent voltage outputs, hinting at their suitability as anodic material(s) for enhancing the energy storage performance of boron-ion batteries (BIBs) if extensively studied.

Conclusions
We have studied the high-energy-density anodic compartment for boron-ion batteries as an alternative to lithium-ion and sodium-ion batteries using the first-principles calculations within the framework of Density Functional Theory. We investigated the electrochemical performance of the boron ion(s) on the monolayer, bilayer, trilayer, and tetralayer graphene sheet electrodes. Significant decreases in the HOMO-LUMO energy gap from −5.085 eV to −2.242 eV for B 3+ @MG and from −20.08 eV to −19.84 eV for 2B 3+ @TTG were recorded. The predominant interaction force for the graphene layers was van der Waals, as predicted by the reduced density gradient isosurface analyses, with the B 3+ @TTG_asym and B 3+ @TTG_sym configurations showing the most favorable interactions.
Furthermore, the electrochemical cell voltages obtained with the single-layer (B 3+ /B@MG) and/or multilayer (B 3+ /B@BG, B 3+ /B@TG, and B 3+ /B@TTG) graphene sheets were significantly improved, with B 3+ /B@MG (13.7 V) showing the least voltage. On the contrary, B 3+ /B@TTG (16.5 V) showed the highest voltage. The results revealed that increasing the number of graphene layers improves the electrochemical performance of the anodic electrode of the boron-ion battery. Therefore, the number of layers or the thickness of the graphene nanosheets is another effective parameter to tune the anode performance for boron-ion batteries. This might suggest an additional dimension to further engineer the graphene anode performance by increasing the number of graphene layers. Our theoretical investigations demonstrate the suitability of the graphene-based anodic electrodes for boron-ion batteries when extensively investigated for large-scale application as a highly efficient energy storage substitute for lithium-and sodium-based ion batteries, which could be beneficial to the material design of ion batteries.

Data Availability Statement:
The raw/processed data required to reproduce these findings cannot be shared at this time due to time limitations and will be made available by the corresponding author upon request.
Acknowledgments: For computing resources, we would like to acknowledge the ITC unit of the King Fahd University of Petroleum and Minerals (KFUPM) for providing access to the HPC computer resources.

Conflicts of Interest:
The authors declare no conflict of interest. The graphene layer directly beneath the intercalated boron atom/ion β-layer

RDG
The graphene layer directly above the intercalated boron atom/ion AA Similar graphene layer stacking for bilayer AB Different graphene layer stacking for bilayer ABA Alternated graphene layer stacking for trilayer AAB Random graphene layer stacking for trilayer ABAB Alternated graphene layer stacking for tetralayer AABA Random graphene layer stacking for tetralayer AABB Regular graphene layer stacking for tetralayer