Simulating key properties of lithium-ion batteries with a fault-tolerant quantum computer

There is a pressing need to develop new rechargeable battery technologies that can offer higher energy storage, faster charging, and lower costs. Despite the success of existing methods for the simulation of battery materials, they can sometimes fall short of delivering accurate and reliable results. Quantum computing has been discussed as an avenue to overcome these issues, but only limited work has been done to outline how it may impact battery simulations. In this work, we provide a detailed answer to the following question: how can a quantum computer be used to simulate key properties of a lithium-ion battery? Based on recently introduced first-quantization techniques, we lay out an end-to-end quantum algorithm for calculating equilibrium cell voltages, ionic mobility, and thermal stability. These can be obtained from ground-state energies of materials, which are the core calculations executed by the quantum computer using qubitization-based quantum phase estimation. The algorithm includes explicit methods for preparing approximate ground states of periodic materials in first quantization. We bring these insights together to estimate the resources required to implement a quantum algorithm for simulating a realistic cathode material, dilithium iron silicate.

Figure 1

Quantum computing for battery simulations. (a) Sketches depicting three key properties of lithium-ion batteries that can be obtained from calculations of the ground-state energies of cathode materials and isolated molecules (Sec. 2). (b) Summary of the main steps of the first-quantized quantum algorithm implemented in this work. The ground-state energy $E$ of a given material is obtained by running a qubitization-based quantum phase estimation (QPE) algorithm on a quantum computer (Sec. 3c). The initial state for the QPE method is obtained by calculating Hartree-Fock orbitals and using the quantum computer to prepare the corresponding antisymmetric Hartree-Fock state (Sec. pp4-s1). (c) Examples of measurable quantities that can be derived: the cell voltage is given by the difference between the chemical potentials ($\mu$) of the electrodes computed from the energy variation ($\mathrm{\Delta }E$) of the cathode material; the activation energy ($E_a$), which is used to predict the ionic mobility; and the temperature profile that helps to define the battery thermal stability.

Figure 2

Schematic of a typical lithium-ion battery. The negative electrode is usually a graphitic carbon that holds lithium ions within its layers, whereas the positive electrode is a source of lithium ions. During discharge, lithium ions move within the battery from the negative to the positive electrode. The process is reversed during charging. The electrolyte transports lithium ions between the electrodes and the separator functions as a physical barrier keeping cathode and anode apart. Negative and positive collectors receive electrons from the external circuit during charging and discharging, respectively.

Figure 3

(a) Conventional unit cell for the dilithium iron silicate ${\mathrm{Li}}_2\mathrm{FeSi}{\mathrm{O}}_4$ cathode material [90]. (b) Crystal structure of the ${\beta }_{II}$ polymorph where all tetrahedra point in the same direction. Along the ${\mathit{a}}_3$ direction, chains of ${\mathrm{LiO}}_4$ (green) are parallel to rows of alternating ${\mathrm{SiO}}_4$ (blue) and ${\mathrm{FeO}}_4$ (brown) tetrahedra representing the fourfold coordination of the lithium, silicon, and iron atoms by the oxygen atoms located at the vertices. This figure was produced using the VESTA package [92].

Figure 4

Antisymmetrization circuit. Example of an antisymmetrization circuit for three electrons. The operation $C$ represents a comparison test controlled on the two registers that are being compared. The seed register is measured in order to postselect on the collision-free subspace. The $Z$ gates perform the phase flip when swapping two registers. At the end of the circuit, the auxiliary record qubits (bottom register) can be discarded as they are disentangled [41]. This circuit can be extended to an arbitrary number of electrons $\eta$ by increasing the size of the sorting network and adding additional auxiliary qubits for each required comparison and swap.

Figure 5

Circuit diagram of an example controlled rotation $R_Y\left({\theta }_{pq}\right)$. The rotation is performed on the subspace spanned by $|p\rangle =|0101\rangle$ and $|q\rangle =|0010\rangle$. The following procedure is applied to the bits where they differ, namely, the last three qubits. First, we apply $X$ gates such that $|p\rangle \to |0000\rangle$ and $|q\rangle \to |0111\rangle$. Then cnot gates map these states to $|0000\rangle$ and $|0100\rangle$ respectively. This allows us to perform a rotation on the second qubit controlled on the auxiliary qubit ${|\sum b_{i}mod2\rangle }_a$. Finally, the cnots and $X$ gates are uncomputed, yielding the desired controlled rotation on the subspace $\text{span}\{|p\rangle ,|q\rangle \}$.

Figure 6

Non-Clifford gate cost for initial state preparation and quantum phase estimation. (a) The non-Clifford gate cost due to Givens rotations used in the circuit for initial state preparation. (b) Toffoli gate cost of the quantum phase estimation algorithm. All calculations are done for the unit cell of ${\mathrm{Li}}_2{\mathrm{FeSiO}}_4$ with 156 electrons. The total number of qubits is 2375 for $n_p=4$ and 6652 for $n_p=9$. In the right figure we depict only Toffoli gate count, as the number of T gates is much smaller ($<3\times {10}^5$). The total error $\varepsilon$ includes contributions from different approximations throughout the algorithm, but it does not take into account the error derived from a finite basis set. The range of values for the error $\varepsilon$ was set based on the required accuracy for simulating the battery properties described in Sec. 2a. While an error in the energy of the order of 0.1 eV is tolerable for predicting experimental voltages, simulating the ionic mobility and the thermal stability requires a higher precision of $\varepsilon \sim 0.1\text{mHa}=0.0027\text{eV}$. The slope of the Toffoli gate cost for fixed target precision is a consequence of the leading cost term in (34), $12\eta n_p⌈\left(\pi \lambda \right)/\left(2{\varepsilon }_{\mathrm{QPE}}\right)⌉$, where $n_p=\lceil log(N^{1/3}+1)\rceil$. The indicated $eV$ energy precision directly translates into the precision of the calculated voltage $V$ [see Eq. (13)], as the denominator in that equation represents a multiple of the charge of an electron, in our case $1e$. Note also that in our battery example we would need to compute the energy of both the lithiated and delithiated phase, the latter being cheaper because of a smaller $\eta$. Consequently, the error in both estimations will be added when computing any of the battery properties discussed. Thus a chemical precision error of $\varepsilon =0.043$ eV is just under the target voltage error of 0.1 V. These calculations were performed with the TFermion library [45].

Figure 7

Estimation of the time required to run the algorithm. This figure illustrates total runtime for synthesizing all the Toffoli gates indicated in Fig. 6 for $\varepsilon =0.043$ eV. All calculations are done for the unit cell of ${\mathrm{Li}}_2{\mathrm{FeSiO}}_4$ with 156 electrons, and we assume that the number of plane waves used in the state preparation and quantum phase estimation are the same. The total number of qubits is 2,375 for $n_p=4$ and 6652 for $n_p=9$. We compute the distillation time as the product of the number of Toffoli gates, the surface code distance $d$, and the clock frequency, all divided by a small $n_p$ factor originating from the techniques in [124] that parallelize the CSWAPs and arithmetic computations. We compute $d$ in this figure as in the moderate error case of Ref. [47]. We emphasize that these are rough estimates whose main purpose is to provide a method to interpret the gate cost.

Figure 8

Conceptual circuit diagram for the quantum phase estimation algorithm. The system register is initialized in the eigenstate $|\psi \rangle$ of the target unitary. The auxiliary register consists of $t$ qubits, initialized in the basis state $|0\rangle$. After applying a Hadamard gate (H) on all auxiliary qubits, increasing powers of the target unitary are applied to the system register, controlled on the state of each auxiliary qubit. Concluding with an inverse quantum Fourier transform (${\mathrm{QFT}}^{-1}$) and measuring the output qubits gives a binary string that can be processed to estimate the desired phase.

Figure 9

High-level representation of PREP subroutine. The first two rotations can be attributed to the preparation of registers $a$ and $m$ in Eq. (E26). The joint circuit for ${\mathrm{PREP}}_T$ and the extra step of ${\mathrm{PREP}}_V$ is shown in Fig. 12. Finally, we have the preparation of the momentum state and the QROM routine.

Figure 10

Preparation of a state proportional to ${\sum }_r{2}^{r/2}|r\rangle$, where $r$ is encoded in unary. We prepare exponentially decreasing amplitudes by using a sequence of Hadamard and controlled-Hadamard gates. By flipping all qubits, the desired state with exponentially increasing amplitudes is obtained. The bottom qubit represents a flag qubit that can be measured to project out the invalid all-zero state.

Figure 11

Circuit representation of ${\mathrm{PREP}}_T$ and ${\mathrm{PREP}}_V$. The uniform superpositions can be implemented with a series of single control rotations or with Hadamard gates ($H$) and an inequality test. The creation of exponential superpositions of the form ${\sum }_r{2}^{r/2}|r\rangle$ are implemented as in Fig. 11. The preparation of the uniform superposition over $|j\rangle$ and the equality test with $|i\rangle$ is only required for ${\mathrm{PREP}}_V$, but we include it here due to its conceptual similarity to other preparations in ${\mathrm{PREP}}_T$.

Figure 12

Main technique used in the SEL operator. The strategy used in SEL consists of (i) swapping the $|{\mathit{p}}_j\rangle$ register into an auxiliary register, controlled on the value of $|j\rangle$; (ii) performing the target uncontrolled operation $O$, where some additional register $|a\rangle$ such as $|\mathit{\nu }\rangle$ might intervene; and (iii) reversing the swaps. Each $|{\mathit{p}}_j\rangle$ contains three $\omega$ coordinates $x, y$ and $z$, each with $n_p$ qubits. We use $O$ to represent different potential operations applied during ${\mathrm{SEL}}_H$. In ${\mathrm{SEL}}_T$, the operator $O$ represents the application of a phase ${(-1)}^{b({\mathit{p}}_{\omega ,r}{\mathit{p}}_{\omega ,r}\oplus 1)}$. For ${\mathrm{SEL}}_U$ and ${\mathrm{SEL}}_V$, it may similarly refer to controlled phases or to arithmetic sums for computing $|\mathit{q}-\mathit{\nu }\rangle ,|\mathit{p}+\mathit{\nu }\rangle$.

Figure 13

Quantum circuit for momentum state preparation. The circuit for implementing a state with exponential amplitudes is the same as in Fig. 11. The controlled Hadamard gates correspond to the preparation of ${\mathcal{C}}_{\mu }$ in Fig. 14. There is also a register for the uniform superposition over $|m\rangle$, as well as the three tests. The later one checking the condition ${\left(2^{\mu -2}\right)}^{2}M>m{\parallel \mathit{\nu }\parallel }^2$ constitutes the key step in this procedure.

Figure 14

Preparation of the superposition corresponding to ${\mathcal{C}}_{\mu }$. The first register in $|{\nu }_i\rangle$ is the sign qubit, using controlled Hadamard gates. This procedure has to be repeated for $i\in \{x,y,z\}$. The last multicontrolled not can be understood as part of the detection of $\mathit{\nu }$ having a $-0$ value in one of the components.