Molecular Understanding of Energy Storage in Supercapacitors with Porous Graphyne Electrodes: From Semiconductor to Conductor

Two-dimensional (2D) porous materials with high specific surface area and ordered 5 morphology exhibit great potential as supercapacitor electrodes. The fundamental 6 understanding of the charge storage and charging dynamics of 2D porous materials can 7 help the optimal design of supercapacitors. Herein, we investigated the energy storage, 8 including the double layer and quantum capacitances, of supercapacitors with typical 2D 9 porous graphynes in the ionic liquid electrolyte by combining molecular dynamics 10 simulation and density functional theory. Simulations revealed that supercapacitors with 11 porous graphyne electrodes could obtain excellent double-layer capacitances, but their total 12 capacitances are limited by the low quantum capacitances. We further predicted boron- 13 /nitrogen-doped graphynes and found that the new porous graphynes turn into good 14 conductors after doping and could achieve a quite high quantum capacitance. The charging 15 dynamics in nanoscale and capacitive performance in macroscale based on the predicted 16 graphyne electrodes were evaluated by combining molecular simulation and transmission 17 line model. Results demonstrate that both outstanding gravimetric and volumetric energy 18 and power densities could be obtained in doped porous graphyne supercapacitors. These 19 findings pave the way for understanding energy storage mechanisms and designing high- 20 performance supercapacitors.

Supercapacitors can provide fast charging and high power density but moderate energy density 22 compared to batteries 1,2 . Enhancing their energy density without compromising the power 23 density can make their more widespread applications. Porous carbons have been widely used as 24 electrode materials in supercapacitors to enhance their energy storage 3-5 . The amorphous porous 25 carbons (e.g., activated carbons 6 and carbide-derived carbons 7 ) are the most popular due to their 26 large specific surface area and high electronic conductivity 5 . However, the random morphology 27 structure would lead to long ion motion paths and severe ion collisions, thereby slowing down 28 the charging dynamics and reducing power density 5, 8, 9 . 29 Conductive 2D porous materials, with regular porosities and superior chemical stability, 30 have attracted growing research interest as electrodes for supercapacitors [10][11][12][13][14] . Sheberla et al. 31 revealed that the supercapacitor with electrodes made entirely of the conductive metal-organic 32 framework (MOF) could obtain high areal capacitance and low cell resistance in an organic 33 electrolyte 10 . Supercapacitors with hexaaminobenzene-derived MOF electrodes and aqueous 34 electrolytes were measured to have both high volumetric and high areal capacitances 12 . 35

Combining molecular dynamics (MD) simulations and experimental measurements, conductive 36
MOFs in an ionic liquid (IL) were uncovered to achieve unprecedentedly high volumetric energy 37 and power densities 13 . Nevertheless, those conductive MOFs only have a moderate specific 38 surface area, and therefore, more 2D porous electrodes need to be developed. Porous graphynes, 39 such as hydrogen substituted graphyne (HsGY) 15 , hydrogen substituted graphdiyne (HsGDY) 40 with AA stacking 16 and AB stacking 17 structures, and fluoride substituted graphdiyne 18 , with 41 extended π-conjugated carbon skeleton of porous structure and high specific surface area are 42 synthesized and have been applied to lithium and sodium ion storage [17][18][19] and supercapacitor 20 . 43 For example, supercapacitors based on HsGDY electrodes and 6M KOH aqueous electrolyte 44 could achieve a high specific capacitance (230 F/g) 20 . 45 However, the band gaps have been found in both HsGY 15 and HsGDY 21 , implying that they 46 are not good conductors. Previous studies showed that low-dimensional materials with low 47 density of states (DOS) near Fermi level, such as graphene 22

64
In-pore double-layer structure 65 We begin our analysis by examining the microstructure of electrolytes inside pores. anions become more heterogenous: counterions (cations at negative electrode and anions at 70 positive electrode) are adsorbed closer to the pore surface, while co-ions are repulsed from the 71 pore surface to the center. The radial distribution exhibits more detailed layering, as shown in Fig.  72 1b. Co-ions form a two-layer structure inside the pore under weaker polarization: one layer near 73 the pore surface and the other near the center ( Fig. 1b and Fig. S2a), but this distribution will 74 transit to a one-layer structure as the polarization gets stronger. As for axial distribution shown in 75  Moreover, cations distribute more homogenously in axial direction under negative polarization 79 than at PZC or under positive polarization due to their more random orientation (Fig. S2c).

88
For AB stacking HsGDY, the pore is divided into two parts ( Fig. S3a): the larger part can be 89 wetted by the ions at the PZC and thus defined as a channel; while the smaller part is named as a 90 cage because it shows ionophobicity and only counterions could be adsorbed into this part under 91 polarization. Interestingly, the in-plane ion distributions in the smaller pore (HsGY, Fig. S4) and 92 larger pore (HsGTY, Fig. S5) exhibit a petal pattern and a round pattern, respectively. The 93 distributions of cations and anions in HsGY are homogeneous and almost immune to polarization 94 (Fig. S4c); while in HsGTY, more layers of ions along the radial direction are observed, 95 especially under negative polarization, and the orientation of cation is distributed more randomly 96 ( Fig. S5b vs. Fig. S2b and Fig. S4b)

106
∆ is the interaction energy difference between ions in the pores and ions in the bulk.

107
We further scrutinized the EDL capacitance of porous graphyne supercapacitors. The differential 108 EDL capacitance, , is calculated by d /d , where is surface charge density of 109 electrode and is the electrostatic potential of the electrode. The AA stacking HsGDY shows a 110 camel shape of the capacitance-potential curve with two maximal values of 8.36 F/cm 2 at -0.71 V 111 and 7.95 F/cm 2 at 1.05 V (Fig. 2a). Such shapes of capacitance-potential curves were also 112 observed for HsGY and HsGTY; however, their specific capacitances are smaller than AA 113 stacking HsGDY (Fig. S6). But for AB stacking HsGDY, the capacitance-potential curve turns 114 into a bell shape with a maximum value of 13.3 F/cm 2 at 0 V. The area-specific EDL capacitance 115 of those porous graphyne supercapacitors is comparable to that of general porous carbon 116 materials. While, porous graphynes could bring excellent gravimetric EDL capacitance and 117 volumetric EDL capacitance (e.g., 220 F/g and 130 F/cm 3 for AA stacking HsGDY, 297 F/g and 118 192 F/cm 3 for AB stacking HsGDY) ( Fig. 2b and Fig. S7), contributed by their high specific 119 surface area and specific volume. 120 The capacitance-potential curves could be understood by the potential-dependent ion 121 number density inside pores (Fig. 2c). Under weaker polarization (between -0.75 V and 1 V), the 122 numbers of cations and anions change faster as polarization is intensified. Under stronger 123 polarization (< -0.75 V and > 1 V), the number changes become less pronounced with increasing 124 polarization. Especially when the potential is lower than -2 V, the anions are almost excluded 125 from the pores, and their number is nearly immune to the potential variation. That accounts for 126 the dual-peak capacitance-potential curve, and smaller EDL capacitance below -2 V. For AB 127 stacking HsGDY, the number density of ions changes the fastest near PZC (Fig. 2d), 128 corresponding to the maximum EDL capacitance obtained at 0 V. Surprisingly, counterions in the 129 channel of AB stacking HsGDY would abnormally decrease because some counterions enter into 130 the cage. The smaller anions in the positive electrode are more likely to enter the cage than the 131 bigger cations in the negative electrode. Therefore, the capacitance of EDL near the positive 132 electrode is higher than that near the negative electrode (Fig. 2a). 133 The analyses of ion packing inside pores reveal that ions in AA stacking pores gather more 134 closely to each other than those at bulk (Fig. 2e), leading to stronger interactions (Fig. 2f). While 135 ion distribution and interaction energy in AB stacking pores are closer to those at bulk due to the 136 screening effect of the electrode. Weaker interactions between ions in AB stacking HsGDY make 137 for a lower steric effect, resulting in more convenient ion transport. 26 That explains why the 138 capacitance of EDL near AB stacking HsGDY is higher than that near AA stacking HsGDY. 139

Quantum capacitance 140
We then investigated the quantum capacitance of porous graphyne supercapacitors. The applied 141 voltage to the device equals the difference of electrochemical potentials between two electrodes, 142 and the electrochemical potential is composed of chemical potential and electrostatic potential 27 . 143 The influence of chemical potential is negligible for well-conductive electrodes (e.g., metal) but 144 is significant for semiconductors (Fig. 3a). For the latter, the electrostatic potential difference 145 equals electrochemical potential difference minus chemical potential difference. As a result, the 146 actual electrostatic potential drop of the EDL is smaller than the applied voltage, and it is 147

154
We developed quantum capacitance theory in supercapacitors and clarified the related 155 terminologies (details see Methods and Supplementary part 4). In short, the chemical potential 156 of the electrode is a function of its excess charge, and the proportionality coefficient is quantum 157 capacitance, which can be calculated from the DOS of the electrode. The quantum capacitance is 158 calculated as: 159 (1) 160 where is the DOS, is the thermal broadening function and ℎ is the chemical 161 potential. is the relative energy with respect to Fermi level, and is the elementary charge. 162 Notably, the sign of ℎ must be positive. Therefore, the total capacitance, , can be 163 regarded as the sum of quantum capacitance and EDL capacitance in series (Fig. S8): 164 .
(2) 165 As the DOS of graphynes shown in Fig. 3b, band gaps are observed in both HsGDYs (1.41 166 eV for AA stacking one and 2.07 eV for AB stacking one), evidencing that they are typical 167 semiconductors with non-negligible quantum capacitance. Similar band gaps are also observed in 168 HsGY (1.72 eV) and HsGTY (1.21 eV) (Fig. S9). Figure 3c depicts the quantum capacitances of 169 both HsGDYs. One can find that, different from the graphene 22 , the quantum capacitance of 170 porous graphynes exhibits obviously asymmetric feature. The quantum capacitances of porous 171 graphynes are quite large under positive polarization, even 6000 F/g, much higher than that of 172 graphene 22 , and thus have little effect on total capacitances. But under negative polarization, the 173 quantum capacitances are so small due to the band gap that they would dominate the total 174 capacitance. 175 As predicted, the total capacitances under the positive polarization approximate the EDL 176 capacitances, while those for the negative electrode are lowered quite a lot (e.g., 25 F/g at -1.25 177 V for AA stacking HsGDY, Fig. 3d), which would eventually limit the practical application of 178 porous graphynes as electrodes. 179

Doping enhances capacitance 180
The relatively wide band gaps of HsGDYs will lead to low quantum capacitance and 181 conductivity 17 , and thus reduce the energy and power densities of the device. It has been revealed 182 that the quantum capacitance 28, 29 and conductivity 30 of graphene could be enhanced by boron (B) 183 or nitrogen (N) doping. We predicted the structures of AA stacking HsGDY with B or N-doped 184 by first-principles calculations to improve their capacitive performance (Fig. S10a). The results 185 of DOS reveal that Fermi levels of all doped HsGDYs shift to the conduction band or valence 186 band (Fig. 4a, Fig. S10b and S11a), meaning that HsGDYs transit from semiconductor to 187 conductor after doping. Specifically, the electronegativity of B (N) is lower (higher) than carbon, 188 leading to electron deficiency (electron injection) and down-shift (up-shift) of the Fermi level. 189 The higher doping concentration of B (N) leads to more electron deficiency (electron injection), 190 resulting in further down-shift (up-shift) of the Fermi level ( Fig. 4a and Fig. S11a).

196
As exhibited in Fig. 4b, Fig S10c and S11b, the quantum capacitance of doped HsGDY is 197 much higher than that of the pristine one and could reach 1000 F/g at 0 V even with lower 198 doping concentrations (i.e., B1-HsGDY and N1-HsGDY). The quantum capacitance of B6-199 HsGDY has a rather high quantum capacitance, reaching 6246 F/g at 0.31 V. Therefore, the 200 doping of B/N could noticeably improve the quantum capacitance and, thereby, reduce the 201 influence on the total capacitance. 202 We also examined how doping influences EDL structure and capacitance. It can be noted 203 that, despite its high doping concentration, the structure and capacitance of double layers in N6-204 HsGDY supercapacitors are almost the same as those in the pristine one ( Fig. S12a-b) capacitance of N-doped graphene supercapacitor is similar to that of pristine graphene) 31 . 207 Therefore, the total capacitance of the doped HsGDY can be simply evaluated from the quantum 208 capacitance of the doped HsGDY (e.g., Fig. 4b) and EDL capacitance of the pristine HsGDY 209 supercapacitor (e.g., Fig. 2b) using Eq. (2). Unlike pristine HsGDY, the total capacitance of 210 doped HsGDYs is dominated by EDL capacitance rather than quantum capacitance ( Fig. 4c and  211   Fig. S11c). Consequently, the capacitance-potential curves display double-humped shapes. (3) 224 Meanwhile, the EDL energy density, , is defined as: 225 .
(4) 226 Therefore, the total energy density, , is regarded as the sum of the band-structure energy 227 density and EDL energy density: 228 As seen in Fig. 4d, the energy density of the pristine HsGDY is very small under negative 230 polarization. After doping, the energy density is significantly enhanced, even for low doping 231 concentrations. The energy density of B1-HsGDY is 57.93 Wh/kg (28.60 Wh/kg for the pristine 232 one) at 3 V, and an outstanding energy density with 98.89 Wh/kg could be achieved at 3.82 V of 233

B6-HsGDY. 234
It is worth noting that we above focused on quaternary doping. In practice, the N-doped has 235 different N bonding states, including quaternary doping, pyridinic doping, and pyrrole doping. 29 236 Pyridinic N-doped HsGDY has been synthesized and showed an excellent catalytic 237 performance 33 . However, as shown in Fig. S13, there is still a large band gap in pyridinic N-238 where ∞ is the charge when the pore is fully charged, is the pore volume divided by its surface 248 area, is the full length of the pore ( = 0.2, =6.27 nm for B6-HsGDY), is the intrinsic 249 relaxation time and can be calculated by fitting MD data. The ionic conductivity inside the pore, 250 , is further calculated by: 251 .
(7) 252 Therefore, the charging curve obtained from MD simulations can be fitted with TLM to quantify 253 the charging dynamics (Fig. 5a). The intrinsic relaxation time decreases as the voltage increases 254 (Fig. 5b), which can be explained by in-pore ion packing (Fig. 1b). Under lower polarization, the 255 radial distributions of cations and anions overlap as both counterions and co-ions reside near the 256 wall, causing the conflicting motion paths of adsorbed counterions and desorbed co-ions, which 257 blocks ion transport and thus slowing down the charging dynamics (Fig. 5c). Under higher 258 polarization, cations and anions reside in different regions because co-ions will be excluded from 259 the pore center, so that their motion paths are separate, resulting in faster charging (Fig. 5d). B6-HsGDY-based device are 6-16 s (Fig. S14b), much smaller than that of porous carbon (about 275 tens of minutes) 35 . The ionic transport resistance in the electrode can be evaluated by 13 : 276 where denotes the porosity (here is 0.33), and the resistances are found to be quite small as 278 1.47~4.24 Ω cm (Fig. S14b). 279 The energy and power densities of the device are further calculated based on the 280 capacitance and resistance 36,37 . As shown in the Ragone plot (Fig. 5e) This work built a multi-scale simulation method from the atomic scale to macroscale for 2D 304 porous electrodes. Although we took porous graphynes as an example, the method may be 305 applied to other 2D porous materials, such as MOFs 41 and COFs 42 . 306

Structure optimization and property calculation of graphynes 308
The geometry optimization of each graphyne was derived from DFT calculations as implemented 309 in the Vienna ab-initio simulation package (VASP) 43 (Supplementary Ta ble 1 ). The Perdew-310 DFT calculations. The projector augmented wave (PAW) method 46 with a cutoff energy of 550 313 eV was used to describe the interaction between nuclei and electrons. The periodic images of 314 atoms were used in all three directions. A k-point mesh of 3 × 3× 11 and a convergence criterion 315 for the electronic self-consistent loop of 10 -5 eV were adopted in both structural optimization and 316 DOS calculation. All the atomic positions are fully optimized until all components of the residual 317 forces are smaller than 0.03 eV Å -1 . 318 The structure properties (i.e., pore size, specific surface area, and pore volume) for porous 319 graphynes were calculated by Monte Carlo methods using Zeo++ software 47 with a nitrogen 320 probe radius of 0.186 nm (Supplementary Table 2). 321

MD simulation 322
As shown in In our simulations, the applied potential between two porous graphynes electrodes is 333 maintained by the constant potential method (CPM), which allows the fluctuations of changes on 334 electrode atoms during simulation 52-54 . Apparently, this approach is more realistic than the often 335 used constant-charge method, as it considers the electronic polarizability of the electrode and the 336 image force effects 55 . The porous graphynes typically have a 3D rough surface, so that the local 337 redistribution of charges on the electrode due to interaction with electrolyte could not be 338 neglected 37, 52 . Besides, the constant charge simulation approach leads to an unphysical 339 enhancement of the charging dynamics 56, 57 , so it is essential to use the CPM to study the 340 charging dynamics. The CPM method was implemented in a methodology proposed by 341 Siepmann et al. 58 and refined by Reed et al. 59 . The system was first held for 20 ns for 342 equilibration without polarization. To obtain EDL structure and capacitance inside the pore, MD 343 simulations were performed under the cell voltage ranging from 0 to 6 V. In this approach, the 344 charge on each electrode was updated at each time step (2 fs). 345

Calculation of quantum capacitance 346
This theory of quantum capacitance is based on fixed-band approximation (i.e., the electronic 347 structure remains unchanged during charging 60 ). We assume that the electrochemical potential of 348 electrons is the summation of chemical potential and electrostatic potential, and that chemical 349 potential only depends on the number of electrons 27 . 350 We first discuss the condition at the temperature of 0 K. The numbers of electrons inside the 351 electrode at the neutral state and the polarized state are respectively described as: where D(E) is the density of states, is the Fermi level, and are the chemical potential at 355 the neutral state and the polarized state, respectively, ̅ is the electrochemical potential, and is 356 the electrostatic potential of the electrode. 357 The net charge of the electrode is 358 .
(11) 359 Here we introduce a new concept ℎ called chemical potential in volt as δE N /δQ, 360 which means the change in free energy as the change of the charge. The aim is to introduce a 361 physical quantity that has the same dimension as electrostatic potential so that we can use them 362 in an equivalent circuit 363 ℎ .
(13) 366 The main difference between this quantum capacitance from that of the same name in the field of 367 transistors is the negative sign. Only in this way can we set up the equivalent circuit model as 368 quantum capacitance and EDL capacitance in series. 369 To consider the effects of temperature, the thermal broadening function, 370 sech /2 /4 , needs to be included as shown in Eq. (1). More information could be found 371 in Supplementary part 4. 372