Interactive Multiscale Modeling to Bridge Atomic Properties and Electrochemical Performance in Li-CO₂ Battery Design

We performed first-principles-based Density Functional Theory (DFT) calculations to optimize the crystal structures of Sb_0.67Bi_1.33Te₃ and Li₂CO₃ compounds using the Vienna Ab$-$initio Simulation Package (VASP) [53, 54]. The Sb_0.67Bi_1.33Te₃ structure is constructed based on the Sb₂Te₃ structure from the Materials Project database [55]. Multiple possible structures for Sb_0.67Bi_1.33Te₃ were optimized, and the most energetically favorable configuration (i.e., with the lowest optimization energy) is presented in Figure 1 (B). The Li₂CO₃ structure presented in 1B is obtained from the Materials Project database [55]. The projector augmented wave (PAW) formalism described the valence electrons of Bi, Sb, Te, Li, C, and O atoms using plane wave-based wave functions was employed [56]. The structure optimization, involving the minimization of ground state energy, utilized the generalized gradient approximation of Perdew and Wang method that brings the correct exchange-correlation part of electronic energy [57, 58]. A kinetic cutoff energy of 500 eV was set to enhance calculation accuracy. The unit cell of both Sb_0.67Bi_1.33Te₃ and Li₂CO₃ compounds (Figure 1 (B)) extended in all three directions such that the lattice constants become more than 20Å to sample $1\times 1\times 1$ k-meshes of the Brillouin zones, respectively. Ionic and electronic optimizations were alternately performed until the forces on each ion reached less than $\pm 10$ meV/Å.

For the Sb_0.67Bi_1.33Te₃ surface interaction (Figure 1 (D)), the Perdew-Burke-Ernzerhof (PBE) functional, within the generalized gradient approximation (GGA), was employed to describe the exchange-correlation energy. A plane-wave basis set was used with an appropriate energy cutoff to ensure energy convergence, and the interaction between valence and core electrons was treated using the projector augmented wave (PAW) method. A cut-off energy was set at 500 eV and the Brillouin zone sampled at gamma point with a k-mesh of 4 x 4 x 1. The convergence criterion of the total energy was set to be within 1x10^-5 eV within the K-point integration and optimizations were performed until the forces on each atom were below 0.001 eV/Å. Finally, to account for the van der Waals interactions, the van der Waals D3 corrections were implemented [59]. Zero-point energy, heat capacity, and entropy corrections were computed using statical mechanics within the harmonic approximation.

Using optimized structures of Sb_0.67Bi_1.33Te₃ and Li₂CO₃ compounds from our DFT calculation (Figure 1 (B)), we carried out the Ab initio molecular dynamics (AIMD) simulations at 300K up to 1ps with NVT ensemble and 1fs time step to select 10 different frames in the 1ps trajectory data. Subsequently, we generated wavefunctions ($\Psi_{i,k}$ where $i$ and $k$ are band and k-point indices, respectively) data for the selected frames using DFT calculations as explained by Knyazev et. al. [60, 48] to compute the electronic conductivity based on Kubo-Greenwood (KG) formula. During the calculation of wavefunction data from VASP, we included 100% more unoccupied bands than occupied bands for both Sb_0.67Bi_1.33Te₃ and Li₂CO₃ compounds to account for excited state wavefunctions and their corresponding energies. The dynamic electrical conductivity $\kappa(\omega)=\kappa_{1}(\omega)+i\kappa_{2}(\omega)$ is represented as a coefficient of the current density $\mathbf{J}(\omega)$ and electric field $\mathbf{E}(\omega)$:

$$\mathbf{J}(\omega)=[\kappa_{1}(\omega)+i\kappa_{2}(\omega)]\mathbf{E}(\omega)$$

(1)

where $\omega$ is the frequency of the external electric field $\mathbf{E}(\omega)$. The important component $\kappa_{1}(\omega)$ relates to the energy absorbed by the electrons and is computed by the Kubo-Greenwood formula as [48, 49, 50]

$$\kappa_{1}(\omega)=\frac{2\pi e^{2}\hbar^{2}}{3m\omega\Omega}\sum_{i,j,\alpha,k}W(k)|\langle\Psi_{i,k}|\nabla_{\alpha}|\Psi_{j,k}\rangle|^{2}[f(\varepsilon_{i,k})-f(\varepsilon_{j,k})]\delta(\varepsilon_{j,k}-\varepsilon_{i,k}-\hbar\omega)$$

(2)

Here, the sum is over all $\mathbf{k}$ points in the Brillouin zone, band indices $i$ and $j$ and three spatial dimensions $\alpha$. $\kappa_{1}(\omega)$ is the frequency-dependent electrical conductivity. $e$ is the elementary charge. $\hbar$ is the reduced Planck constant. $m$ is the electron mass. $\omega$ is the angular frequency. $\Omega$ is the volume of the system. $f_{i}$ and $f_{j}$ are the Fermi-Dirac distribution functions for states $i$ and $j$. $|\psi_{i}\rangle$ and $|\psi_{j}\rangle$ are the wavefunctions of states $i$ and $j$. $\varepsilon_{i,k}$ and $\varepsilon_{i,k}$ i are the energies of states $i$ and $j$. $\delta$ is the Dirac delta function. $W(k)$ denotes a weight function at a particular $\mathbf{k}$ point. The wave function $\Psi_{i,\mathbf{k}}$ is generated from our AIMD simulations at room temperature as explained above. $f(\varepsilon_{i,\mathbf{k}})$ is the Fermi function. In practical calculations, $\delta$-function is written in terms of the Gaussian function:

$$\delta(\varepsilon_{j,k}-\varepsilon_{i,k}-\hbar\omega)\rightarrow\frac{1}{\sqrt{2\pi}\Delta E}\exp\left(-\frac{(\varepsilon_{j,\mathbf{k}}-\varepsilon_{i,\mathbf{k}})}{2(\Delta E)^{2}}\right)$$

(3)

where $\Delta E$ is broadening of the $\delta$-function. This technical parameter plays an essential role in controlling the algorithm in the Kubo-Greenwood formula, particularly at low frequencies.

Classical Molecular Dynamics (MD) simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS [61]. A cut-off distance of $r_{cut}=14$ Å is used for both Lennard Jones and Coulomb interactions, and a time-step of $\delta t=1$ fs is utilized. Using a particle-particle particle-mesh solver (PPPM), long-range electrostatic interactions are resolved with an accuracy of 10^-4 by mapping atomic charges to a 3D mesh, solving Poisson’s equation on the mesh via 3D Fast Fourier Transforms (FFTs), and interpolating electric fields from mesh points back to atoms. [62, 63]. The electrolyte is prepared using an DMSO:(EMIM⁺, BF${}_{4}^{-}$) ratio of 3:2, 1M (Li⁺, TFSI^-), and includes 6000 DMSO molecules, 1857 (EMIM⁺, BF${}_{4}^{-}$) pairs , 566 (Li⁺, TFSI^-), 14 ZnI₂ and 14 CO₂ molecules, giving a total number of atoms of 113142. Each molecule geometry was optimized using the Generalized Amber Force Field (GAFF) [64] and the Universal Force Field force fields (UFF) [65], successively. Then the Automated force field Topology Builder (ATB) is used to obtain the final GROMOS FF parameters for TFSI^-, DMSO, EMIM⁺, BF${}_{4}^{-}$, and CO₂ [66, 67, 68]. We note that the ATB parameter generation considers the following: Lennard-Jones parameters refined against experimental solvation and pure liquid properties, ATB parametrization and validation using a single-range 14 Å cutoff for both Lennard-Jones and Coulomb interactions, atomic charges fitted to quantum mechanical (QM) electrostatic potentials and bonded parameters are assiged using force constants estimated at B3LYP/6-31G* level of theory, resulting in robust symmetry routines for assigning identical parameters to chemically equivalent atoms, bonds, and angles. For a comprehensive understanding of the ATB parameter generation process, refer to the detailed explanations provided in the original research papers [66, 67, 68]. The UFF interaction potential is used for ZnI₂. For the Li⁺, the interaction potential parameters were taken from the Canongia Lopes & Padua (CL&P) force field for ionic liquids [69, 70]. The initial distribution was generated using the Moltemplate package [71]. In this process, the simulated system constituents are placed randomly into an oversized cubic box. Then, the system was energetically minimized to ensure optimal particle distribution. The energy minimization algorithm consists of adjusting particle coordinates iteratively [72]. A finite, unitless stopping tolerance for energy (Estop=10^-7) that represents the energy variation between consecutive iterations divided by the energy magnitude, and a stopping tolerance for force (Fstop=10^-8 kcal/Å), are chosen. After that, the system is agitated using the Langevin thermostat [73] at T=900 K for t=1 ns, followed by the same process using the Nosé-Hoover thermostat [61, 73]. Then, to reach the desired ambient condition of temperature and pressure (T=300 K, P=1 bar), the simulated systems were equilibrated in the NPT statistical ensemble with isotropic constant pressure control using Berendsen barostat [61, 74]. Firstly, from T=900K to T=300K at a constant pressure of P=500 bar for t=3 ns. Secondly, the pressure decreases from P=500 bar to P=1 bar at a constant temperature of T=300 K over a time period of t=3 ns. In the final NPT equilibration stage, the systems were equilibrated at room temperature and ambient pressure (T=300 K, P=1 bar) for t=10 ns. To investigate the transport properties under conditions similar to experimental settings, we performed simulations in the NVT (constant number of particles, volume, and temperature) statistical ensemble. We employed the Nosé-Hoover thermostat to maintain a constant temperature of 300 K while keeping the volume fixed.

A mathematical framework has been developed for modeling the behavior of prismatic Li-CO₂ batteries, drawing inspiration from previous studies on Li-O₂ and Li-ion batteries [75, 76, 77]. The model integrates key principles such as mass and current conservation, species transport phenomena, and reaction kinetics within the cathode and separator regions, with the aim of providing deeper insights into the complex mechanisms that occur within the cell during operation. The Li-CO₂ cell under investigation consists of several components: a lithium metal sheet serving as the negative electrode, a separator, a porous carbon cathode coated with Sb_0.67Bi_1.33Te₃, and an organic electrolyte that permeates the porous structures, as illustrated schematically in Figure 1 (A). The used electrolyte composition is described in section 2.3. To facilitate electron transfer, current collectors are located at the rear of both electrodes. This configuration enables a thorough examination of the electrochemical and transport processes that govern battery cell performance and limitations. The main objective of this mathematical framework is to provide a comprehensive understanding of the intricate interactions among various components and processes within the Li-CO₂ battery cell system.

Calculated dynamic electrical conductivities for both Sb_0.67Bi_1.33Te₃ and Li₂CO₃ are shown in Figure 2, for decremented values of $\Delta E$ from 0.9 to 0.54 with a step of 0.04. The electrical conductivities calculated for different values of $\Delta E$ show the conductive behavior of Sb_0.67Bi_1.33Te₃ (a conductivity range of the order of $10^{5}$ S/m) and insulative behavior of Li₂CO₃ (a conductivity range of the order of $10^{-3}$ S/m). Quantitatively, the electrical conductivities converge with decreasing $\Delta E$, i.e., the dependence of conductivity on $\Delta E$ is smaller. The static electrical conductivities at zero frequency are calculated by interpolating small frequency data shown in Figure 2 (a2) and (b2). Using $\Delta E$ values where conductivity is converged, below 0.66 for Sb_0.67Bi_1.33Te₃ and below 0.86 for Li₂CO₃, our computational analysis reveals a static electrical conductivity of $(1.43\pm 0.3)\times 10^{5}$ S/m for Sb_0.67Bi_1.33Te₃, indicating its strong electrical conducting properties. This finding aligns with the well-established characteristics of bismuth telluride (Bi₂Te₃) and its alloys, which are extensively employed in thermoelectric devices due to their good electrical conductivity [83, 84]. The high conductivity can be ascribed to the semiconductor nature of Bi₂Te₃ and the antimony (Sb) doping effect, which augments the carrier concentration [85]. Experimental investigations on Bi₂Te₃ have reported electrical conductivity values as high as $1.22\times 10^{5}$ S/m [84], further corroborating the computational results. Our calculations reveal that lithium carbonate (Li₂CO₃) exhibits a static electrical conductivity range of $0.1-2\times 10^{-3}$ S/m, indicating its poor electrical conductivity and insulating behavior [86]. This calculated value aligns closely with experimentally reported values, $2\times 10^{-3}$ S/m [87] and $3.16\times 10^{-4}$ S/m [88]. The insulating nature of Li₂CO₃ can be attributed to its ionic nature and the strong ionic bonds that hinder electron mobility [89].

Calculated dynamic electrical conductivities for both Sb_0.67Bi_1.33Te₃ and Li₂CO₃ are shown in Figure 2 (a1) and (b1) for a large frequency range and in Figures 2 (a2) and (b2) for a small frequency range, for decremented values of $\Delta E$ from 0.9 to 0.54 with a step of 0.04. Figures 2 (a2) and (b2), the electrical conductivities converge, for small frequencies, with decreasing $\Delta E$, i.e., the dependence of conductivity on $\Delta E$ becomes smaller. The static electrical conductivities at zero frequency are calculated by interpolating small frequency data shown in Figures 2 (a2) and (b2), using $\Delta E$ values where conductivity has converged, below 0.66 for Sb_0.67Bi_1.33Te₃ and below 0.86 for Li₂CO₃. The static electrical conductivities calculated for different values of $\Delta E$ show the conductive behavior of Sb_0.67Bi_1.33Te₃ (a conductivity range on the order of $10^{5}$ S/m) and insulative behavior of Li₂CO₃ (a conductivity range on the order of $10^{-4}$-$10^{-3}$ S/m). Our computational analysis reveals a static electrical conductivity of $(1.43\pm 0.3)\times 10^{5}$ S/m for Sb_0.67Bi_1.33Te₃, indicating its strong electrical conducting properties. This finding aligns with the well-established characteristics of bismuth telluride (Bi₂Te₃) and its alloys, which are extensively employed in thermoelectric devices due to their good electrical conductivity [83, 84]. The high conductivity can be attributed to the semiconductor nature of Bi₂Te₃ and the antimony (Sb) doping effect, which augments the carrier concentration [85]. Experimental investigations on Bi₂Te₃ have reported electrical conductivity values as high as $1.22\times 10^{5}$ S/m [84], further corroborating the computational results. Our calculations reveal that lithium carbonate (Li₂CO₃) exhibits a static electrical conductivity range of $0.1-2\times 10^{-3}$ S/m, indicating its poor electrical conductivity and insulating behavior [86]. This calculated value aligns closely with experimentally reported values of $2\times 10^{-3}$ S/m [87] and $3.16\times 10^{-4}$ S/m [88]. The insulating nature of Li₂CO₃ can be attributed to its ionic nature and the strong ionic bonds that hinder electron mobility [89]. In Li₂CO₃, Li⁺ ions act as the primary charge carriers rather than electrons, a characteristic typical of ionic compounds. The low conductivity of Li₂CO₃ presents challenges to battery performance, including high internal resistance and passivation of the electrode surface [90, 91]. To address these challenges, potential solutions include enhancing Li₂CO₃ conductivity, promoting its decomposition through conductive additives or catalysts, designing porous electrodes for increased surface area, and exploring alternative Li₂CO₃ structures with improved conductivity or decomposition properties.

The DFT calculations of the reaction path for Li₂CO₃/C formation on the Sb_0.67Bi_1.33Te₃ as Li-CO₂ battery cathode catalyst have revealed intriguing insights into the adsorption behavior of reaction products and the step-wise mechanism of the discharge process. To comprehensively understand the formation of Li₂CO₃/C and C, we calculated the Gibbs free energy ($\Delta$G) changes for various elementary reaction steps on the Sb_0.67Bi_1.33Te₃ catalyst surface, revealing a sequence of intermediate stages that elucidate the complex reaction pathway.

The Gibbs free energy $G$ of the reaction was calculated using the following equation:

$$G=E+\text{ZPE}+\int_{0}^{T}C_{p}\,dT-TS$$

(38)

where ZPE is the zero-point energy, $C_{p}$ is the heat capacity, $S$ is the entropy, $T$ is the temperature, and $E$ is the electronic energy.
We considered ZPE and entropy changes in Step-1 and assumed that these changes for other elementary steps are small and therefore negligible. The potential-dependent reaction energies were calculated based on the computational hydrogen electrode (CHE) model with reference to the Li/Li⁺ electrode [92, 93]. The potential-dependent reaction energies profile is shown in Figure 2 (c), and the respective reaction energies of each step are shown in Table 1. To level up the electrochemical step, a potential was applied to the formation of LiCO₂. As can be seen in Figure 2 (c), this potential occurs at 2.3 V versus Li/Li⁺ electrode.

The reaction begins with the adsorption of CO₂^∗ on the Bi site, with a Bi-O distance of 3.5 Å. This initial step is marginally endergonic, with a $\Delta$G of 0.002 eV, indicating a slight energy barrier to CO₂ adsorption. Following this step, the coupled addition of Li⁺ and an electron to the adsorbed CO₂^∗ forms the LiCO₂^∗ intermediate. This step is highly exergonic, with a $\Delta$G of -2.3 eV, suggesting a thermodynamically favorable process that likely drives the overall reaction forward. The reaction pathway continues with the interaction between two LiCO₂^∗ intermediates, producing Li₂CO₂^∗ and CO₂^∗. This step has a $\Delta$G of -1.2 eV, indicating another energetically favorable process. Subsequently, a pivotal step occurs: the conversion of Li₂CO₂^∗ and CO₂^∗ into Li₂CO₃ and CO^∗, with a $\Delta$G of $-$1.7 eV. This step is crucial as it marks the formation of one of the main discharge products, Li₂CO₃. The final step in the pathway involves the generation of carbon, facilitated through the reaction of CO^∗ and CO₂^∗. Interestingly, this step is energetically uphill, with a $\Delta$G of 2.8 eV, and also produces CO₃^∗ on the surface. This energy barrier for carbon formation could explain the observed preferential adsorption of carbon on the catalyst surface, as the system may need to overcome this energy barrier to complete the carbon formation process.
The complete 5-step reaction pathway can be summarized as follows:

	${}^{*}+CO_{2}(g)\rightarrow CO_{2}{}^{*}$		(Step-1)
	$\displaystyle O_{2}{}^{*}+Li^{+}+e^{-}\rightarrow LiCO_{2}{}^{*}$		(Step-2)
	$\displaystyle 2LiCO_{2}{}^{*}\rightarrow CO_{2}{}^{*}+Li_{2}CO_{2}{}^{*}$		(Step-3)
	$\displaystyle CO_{2}{}^{*}+Li_{2}CO_{2}{}^{*}\rightarrow Li_{2}CO_{3}+CO{}^{*}{}+{}^{*}$		(Step-4)
	$\displaystyle CO_{2}{}^{*}+CO{}^{*}\rightarrow C+CO_{3}{}^{*}{}+{}^{*}$		(Step-5)

These results highlight the essential role of coupled electron-cation transfer processes in the growth of discharge products. The strong exergonic nature of steps 2-4 suggests that these processes occur readily, while the endergonic nature of the final carbon formation step aligns with the observed preferential adsorption of carbon on the Sb_0.67Bi_1.33Te₃ catalyst surface. The lack of Li₂CO₃ adsorption on the catalyst surface, as indicated by the simulation, can be explained by the reaction pathway. Li₂CO₃ is formed in step 4, but this step also releases CO^∗ and a free surface site (^∗). This suggests that Li₂CO₃ formation occurs near, but not directly on, the catalyst surface, reducing the cathode porosity and possibly clogging the cathode pores.
These insights into the reaction mechanism have significant implications for catalyst design and overall battery performance. The strong affinity of carbon for the catalyst surface, coupled with the energetically unfavorable step of its formation, suggests that catalyst degradation through carbon accumulation could be a concern over long-term battery operation. This highlights the need for strategies to mitigate carbon build-up, such as developing self-cleaning catalysts or implementing periodic regeneration protocols. Furthermore, the spatial separation of Li₂CO₃ formation from the immediate catalyst surface underscores the importance of electrolyte engineering in Li-CO₂ battery design. Since Li₂CO₃ formation appears to occur away from the catalyst, the electrolyte’s ability to facilitate Li₂CO₃ formation, dissolution, and transport becomes crucial. This suggests that future research should focus on developing electrolytes or additives that can enhance Li₂CO₃ solubility or promote its formation and decomposition.

This section presents the results of MD simulations that elucidate the atomic-level dynamic and solvation properties of the Li-CO₂ battery electrolyte depicted in Figure 3 (a). These simulations provide microscopic insights into the electrolyte’s behavior, bridging the gap between atomic and macroscopic transport properties in our multiscale approach.

The simulated voltage-capacity curves of the Li-CO₂ battery cell during the charge-discharge phase is presented in Figure 4 (a). This curve was derived from finite element analysis (FEA) using the atomistic simulation results shown in Figure 2 and Table 2, along with other parameters detailed in Tables S1, S2, and S3. The FEA results were compared with the experimental data obtained from the set-up described in the Supporting Information (SI), at a current density of 1 mA/cm². The data illustrated in Figure 4 (a) indicate a strong correlation between the simulated results and the experimental voltage. During discharge, the cell potential exhibited a rapid decline from 3.1 V to a plateau of approximately 2.65 V, followed by a gradual reduction to the cutoff voltage of 2.2 V. During charging, the cell potential increased rapidly from 3.7 V to a plateau of 4.2 V, followed by a gradual rise to the cutoff voltage of 4.5 V. Furthermore, the specific capacity, calculated on the basis of the catalyst mass, was about 6200 mAh/g (2.48 mAh/cm²).

Figure 4 (b) illustrates the impact of applied current density of 1 mA/cm², 0.5 mA/cm², 0.25 mA/cm², and 0.1 mA/cm² on the simulated discharge curves of Li-CO₂ battery cell and Figure S1 (a) shows the same curves for current densities of 0.2 mA/cm², 0.3 mA/cm², 0.4 mA/cm², 0.6 mA/cm², 0.7 mA/cm², 0.8 mA/cm² and 0.9 mA/cm²). The discharge capacity exhibits a substantial reduction, decreasing from 81,570 mAh/g at a low current density of 0.1 mA/cm² to 6,200 mAh/g at a higher current density of 1 mA/cm². This observation aligns with prior experimental results, indicating that current density has a significant influence on cell capacity [98]. Furthermore, the discharge voltage plateau decreases as the current density increases (Figure 4 (b)). The capacity decrease at elevated discharge rates can be attributed to limitations in CO₂ transport through the cathode, which becomes saturated with electrolyte and unable to sustain the electrochemical reaction. Consequently, as the discharge rate escalates, CO₂ reduction is confined to a localized area near the cathode-current collector interface. Furthermore, the rapid reduction in porosity due to Li₂CO₃ and C deposition on the active area’s surface further impedes CO₂ transport into the cell, preventing complete utilization of the electrode’s porosity. Figures 5, S3 and S4 illustrate and analyze these phenomena in detail. Figure S1 (b) illustrates the overpotential behavior of the studied Li-CO₂ battery cell during discharge, revealing a transition from a positive overpotential to a negative overpotential. The initial positive overpotential can be attributed to kinetic and mass transport limitations associated with the CO₂ reduction reaction (ORR) at the cathode. This phenomenon is consistent with observations in other battery systems, such as vanadium-manganese flow batteries, where kinetic barriers result in a positive overpotential during discharge [99]. As the discharge progresses, the overpotential becomes negative, indicating a shift in the reaction dynamics. This transition is due to morphological changes in the cathode surface or discharge products (Figure 5 (a1) and (a2)) that could create a more favorable electrochemical environment [22]. Furthermore, as the discharge continues, the distribution of reactants and products becomes more balanced, alleviating initial transport limitations. This phenomenon is related to the morphological changes observed in the cathode surface during discharge (Figure 5 (a1) and (a2)). We note that Li et al. reported that the formation of $\text{Li}_{2}\text{C}_{2}\text{O}_{4}$ (lithium oxalate) particles on the cathode surface created a more favorable electrochemical environment by providing additional active sites, improving conductivity and increasing the effective surface area [100]. These changes resulted in higher discharge capacity and improved cycling stability, underscoring the complexity of the discharge process and the importance of optimizing materials and design to enhance overall battery cell performance.

The observed behavior where the overpotential in the Li-CO₂ battery cell during discharge decreases to 0 V as current density increases from 0.1 to 1 mA/cm² can be explained by changes in the morphology and formation mechanism of Li₂CO₃/C (Figure 5 (a1) and (a2)). It has been found for the Li-air battery that at lower current densities, the deposited species (Li₂O₂) tends to form larger, toroidal structures, which can lead to higher overpotentials due to less efficient electron transport [101, 102, 103]. Similarly to what was observed for Li-O₂ battery, as the current density increases in Li-CO₂ battery, the deposited species, Li₂CO₃ and carbon form a thin film (see Figure 5 (a1) and (a2)), which maintains closer contact with the surface of the conductive catalyst, thus reducing the overpotential (Figure S1 (b)) because Li₂CO₃ has low electrical conductivity (Figure 2). Figure 4 (c) shows that the capacity of the studied Li-CO₂ battery cell decreases exponentially with increasing current density, as can be described by the following equation [102]:

$$C=C_{0}+A\exp\left(-\frac{i}{i_{0}}\right)$$

(47)

where: $C_{0}$ = 7500 mAh/g is the baseline capacity, $A$ = 110000 mAh/g is a scaling factor, $i_{0}$ = 0.157 $mA/cm^{2}$ is a characteristic current density. At low current density (i<<0.157 $mA/cm^{2}$), the capacity $C$ tends towards its maximum value of 117500 mAh/g. This is the sum of the baseline capacity ($C_{0}$) and the additional capacity represented by $A$. As the current density $i$ increases, there is a rapid decrease in the additional capacity provided by $A$, resulting in a decrease in the total capacity $C$. The parameter $i_{0}$ determines the rate at which the capacity decreases with increasing current density. A smaller $i_{0}$ would mean a faster decrease in capacity for a given increase in current density, while a larger $i_{0}$ would indicate a slower decrease. As the current density increases further, the capacity approaches the baseline value $C_{0}$ of 7500 mAh/g. This suggests that at high current densities, the battery’s capacity is primarily limited to its baseline capacity, as the contribution from the exponential term becomes negligible. The exponential decrease in capacity with increasing current density highlights a trade-off between power output and energy capacity [104]. At higher current densities (i>>0.157 $mA/cm^{2}$), the battery can deliver more power, but at the cost of reduced capacity. This is a critical consideration in applications where both high-power and high-capacity are required. The connection between cell voltage ($V$) and current density ($i$) during the discharge of the studied Li-CO₂ battery cell in Figure 4 (d) is represented by the formula [102, 105]:

$$V=A\cdot\log\left(\frac{i}{\text{mA/cm${}^{2}$}}\right)+B$$

(48)

where $A=-0.256$ V and $B=2.696$ V. The parameter $A$ represents the slope of the voltage change with the logarithm of the current density/(mAcm^-2), while $B$ represents the voltage at a reference current density of 1 mA/cm². This formula indicates that as the current density increases, the cell voltage decreases logarithmically, a common characteristic in electrochemical cells due to increased polarization and decreased voltage at higher current densities [102, 105]. The performance of Li-CO₂ batteries is governed by a complex interplay of structural and operational parameters, including the initial cathode porosity, CO₂ diffusion in the battery electrolyte, initial CO₂ concentration in the cathode, and the kinetics of the cathodic reaction. These factors collectively influence key performance metrics, such as capacity and voltage plateau presented in Figure S2 and Figure 4 (e) and (f). Our analysis reveals a strong positive correlation between porosity ($\varepsilon_{l,0}$) and capacity ($C$), as demonstrated by the best linear fit equation:

$$C/\text{(mAh/g)}=59479.29\,\,\varepsilon_{l,0}-25199.46$$

(49)

This trend is attributed to the increased surface area and improved CO₂ accessibility in higher porosity structures, which facilitate more efficient electrochemical reactions. Similarly, the CO₂ diffusion coefficient ($D_{\text{CO}_{2}}$) and initial CO₂ concentration ($c_{\text{CO}_{2,0}}$) exhibit significant positive correlations with capacity, as shown by the best linear fit equations:

$$C/\text{(mAh/g)}=1.38\times 10^{13}\,(D_{\text{CO}_{2}}/\text{m}^{2}\text{s}^{-1})+1087.06$$

(50)

$$C/\text{(mAh/g)}=3840.00\,(c_{\text{CO}_{2,0}}/(mol/m^{3}))-603.20$$

(51)

These results highlight the critical role of CO₂ transport and availability in determining battery capacity. Enhanced diffusion and higher CO₂ concentrations improve the accessibility of CO₂ to active sites, thereby promoting more efficient electrochemical reactions. This is consistent with recent findings by Chen et al., who demonstrated that a Li-O₂/CO₂ battery can achieve a full discharge capacity of 6628 mAh/g and a long life of 715 cycles, surpassing the performance of pure Li-O₂ batteries [106]. However, the cathodic reaction rate constant ($k_{c}$) shows a logarithmic relationship with both capacity and voltage plateau, as evidenced by:

$$C/\text{(mAh/g)}=763.86\,\,log(K_{c})+28999.81$$

(52)

$$V/(V)=0.06\,\,log(K_{c})+3.64$$

(53)

This logarithmic relationship underscores the significant impact of reaction kinetics on battery performance, emphasizing the need for efficient catalysts to optimize reaction rates. Furthermore, both CO₂ diffusion and concentration demonstrate positive correlations with the voltage plateau, as shown by:

$$V/(V)=2.14\times 10^{7}\,(D_{\text{CO}_{2}}/\text{m}^{2}\text{s}^{-1})+2.80$$

(54)

$$V/(V)=0.01\,(c_{\text{CO}_{2,0}}/(mol/m^{3}))+2.75$$

(55)

These findings indicate that improving CO₂ transport and availability not only enhances capacity but also stabilizes the voltage plateau, which is critical for achieving stable and efficient battery operation. A recent study demonstrated a high-performance Li-CO₂ battery using a Cu(II) solid redox mediator, achieving a higher discharge voltage of 2.8 V and a lower charge potential of 3.7 V [25]. Therefore, optimizing the cathode porosity, CO₂ diffusion and concentration, while improving reaction kinetics through advanced catalysts, are essential strategies for enhancing the performance of Li-CO₂ batteries. Future research should focus on developing innovative electrode materials, electrolyte and catalytic systems to leverage these relationships and advance the practical application of Li-CO₂ batteries as a sustainable energy storage technology. For instance, a recent study reported a high-rate Li-CO₂ battery enabled by a 2D medium-entropy catalyst, which operates at a high current density of 5000 mA/g and capacity of 5000 mAh/g for up to 125 cycles, far exceeding previously reported values [32]. Figure 5 (a1) and (a2) illustrate how porosity ($\varepsilon_{l}$) and volume fraction of Li₂CO₃/C deposited species ($\varepsilon_{s}$) change over space and time with current densities of i=1 mA/cm² and i=0.1 mA/cm², respectively. Figure S3 presents a comparison of these parameters at i=0.1, i=0.25, i=0.5, and i=1 mA/cm². For high current density (i=1 mA/cm²), porosity decreases gradually, primarily near the CO₂ feed side of the cathode (X/L${pos}$=1), from its initial value of $\varepsilon_{l}=0.73$ to 0.48 at the end of discharge. This decrease is due to the formation of a thin film of Li₂CO₃/C, represented by the gradual increase of $\varepsilon_{s}$ from 0 to 0.25 near the cathode feed side. This film formation leads to electrode clogging, impeding efficient CO₂ diffusion into the cathode. For low current density (i=0.1 mA/cm²), Li₂CO₃/C forms larger toroidal structures [101, 102, 103], evidenced by the gradual increase of $\varepsilon_{s}$ across a wider region of the cathode from 0 to 0.64. This results in a distributed porosity decrease from 0.73 to 0.09, initially allowing more CO₂ transport into the cathode. However, as discharge progresses, the substantial porosity decrease hinders CO₂ diffusion even at low current densities. Figure 5 (b1) and (b2) present the CO₂ concentration profiles within the cell during discharge at i=1 mA/cm² and i=0.1 mA/cm², respectively. Figure S4 compares concentration profiles at various current densities. At i=1 mA/cm², the CO₂ concentration ($c_{CO_{2}}$) decreases sharply from its initial value of 4.8 mol/m³ to 0 mol/m³ in the porous cathode, indicating strongly impeded CO₂ diffusion. At i=0.1 mA/cm², the decrease is more gradual. The impediment to CO₂ transport, which is inversely proportional to the current density, results from the accumulation of Li₂CO₃ and carbon on the active cathode surface. This accumulation reduces pore availability for the electrolyte, increases the amount of insulating species, and obstructs pores, collectively restricting CO₂ transport and limiting cell capacity (Figure 4 (c)). To address these limitations, an ideal cathode should feature: large surface area for reaction sites and product deposition, large pores to facilitate active species transport and prevent clogging, high porosity to accommodate discharge product growth, and high electrical conductivity and stability. Such a cathode design could mitigate the observed limitations and potentially enhance the overall performance of Li-CO₂ batteries [107].