Toffoli complexity and qubit cost when constructing the qubitization walk
operator. Unfortunately, $\lambda$ is increased in these cases. The increase
is attributed to the lower variational freedom in constructing non-orthogonal
bases when representing the two-electron integral tensor in factorized form
compared with the non-symmetry adapted setting. For the THC case, no
asymptotic speedup is formally possible. This stems from the linear cost of
unary iteration over all basis states. Nevertheless, due to competing
prefactors between unary iteration and state preparation, we do observe a
$\sqrt{N_{k}}$ scaling improvement in the Toffoli per step and logical qubit
cost for the range of systems studied. This is likely a finite size effect,
but may be a practically important when considering which algorithm to chose
in the future. Thus, improving the $\lambda$ value of THC through more
sophisticated and affordable means is worth further investigation.

Reaching the thermodynamic and complete basis set limit is very challenging,
even for classical wavefunction methods like CCSD and ph-AFQMC. Previous ph-
AFQMC results for simple insulating solids with two-atom unit cells suggest
that at least a $3\times 3\times 3$ and $4\times 4\times 4$ sampling of the
Brillouin zone is required to extrapolate correlation energies to the
thermodynamic limit [108]. Similarly, it has been found that quadruple-zeta
quality basis sets are required to converge the cohesive energy to less than
0.1 eV / atom, while a triple-zeta quality basis is likely sufficient for
quantities such as the lattice constant and bulk modulus [109]. Similar system
sizes and basis sets were found to be required for CCSD simulations of
metallic systems [18]. Although the theory of finite size corrections [110,
111, 112, 113] is still an area of active research [114, 115], the simulation
of bulk systems even with these corrections typically requires on the order of
50 atoms, which in turn corresponds to hundreds of electrons and thousands of
orbitals. For excited state properties, particularly those concerning charged
excitations, even larger system sizes may be required without the use of
sophisticated finite size correction schemes [116]. Thus, we suspect that
simulating large system sizes will continue to be necessary in order to obtain
high accuracy for condensed phase simulations. It is important to note that
high accuracy classical wavefunction methods are often considered too
expensive for practical materials simulation, and DFT is still the workhorse
of the field. Appendix F shows that simulating even simple solids with coarse
$k$-meshes can take on the order of hours, which would otherwise take seconds
for a modern DFT code. From the quantum resource estimates it is clear that
several orders of magnitude of improvement are necessary before practical
materials simulation is possible. Despite this, the fairly low scaling of
phase estimation as a function of system size serves as encouragement to
pursue quantum simulation for materials further.

The aforementioned convergence difficulties are demonstrated in our classical
calculations on the LNO system when attempting to resolve the discrepancy
between band-theory predictions and experimental observations of the ground
state geometry. Furthermore, the variance in energy between CCSD, MP2, ph-
AFQMC, and DMET (and their expenses) make it difficult to select an efficient
method for determining Hamiltonian parameter cutoffs to use in quantum
resource estimation. If anything, this highlights the need for high accuracy
classical computation when performing quantum resource estimates and
ultimately picking an algorithm for quantum simulation. The quantum resource
estimates for LNO simulations are exorbitantly expensive even at small
$k$-mesh; estimated to run in $\mathcal{O}(10^{2})-\mathcal{O}(10^{3})$ days
using the DF LCU. Just as resource estimates for chemistry fell drastically
with algorithmic developments clearly further algorithmic improvements are
needed to make a LNO sized problem feasible on a quantum computer.

Qubitization is a general tool for Hamiltonian simulation and there may be
other simulation scenarios when the improved walk operators yield faster
simulations. There are also areas to further improve the quantum algorithms by
taking advantage of space group symmetry along with translational symmetry. In
classical calculations this can lead to substantial computational savings even
at the mean-field level. Just as in the case of quantum algorithms for
molecular simulations we expect the quantum resource costs to fall with
further exploration of algorithmic improvements.

## Acknowledgements

The authors thank Yuan Su for helpful conversations. FDM thanks Miguel Morales
for helpful discussions on the form of the $k$-point THC factorization. DWB
worked on this project under a sponsored research agreement with Google
Quantum AI. DWB is also supported by Australian Research Council Discovery
Projects DP190102633 and DP210101367.

## References

  * Kassal  _et al._ [2008] I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. Aspuru-Guzik, Proceedings of the National Academy of Sciences 105, 18681 (2008).
  * Babbush  _et al._ [2018a] R. Babbush, C. Gidney, D. Berry, N. Wiebe, J. McClean, A. Paler, A. Fowler, and H. Neven, Physical Review X 8, 041015 (2018a).
  * Su  _et al._ [2021] Y. Su, D. W. Berry, N. Wiebe, N. Rubin, and R. Babbush, PRX Quantum 2, 040332 (2021).
  * Babbush  _et al._ [2018b] R. Babbush, N. Wiebe, J. McClean, J. McClain, H. Neven, and G. K.-L. Chan, Physical Review X 8, 011044 (2018b).
  * Grüneis  _et al._ [2013] A. Grüneis, J. J. Shepherd, A. Alavi, D. P. Tew, and G. H. Booth, The Journal of Chemical Physics 139, 84112 (2013).
  * Kato [1957] T. Kato, Communications on Pure and Applied Mathematics 10, 151 (1957).
  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Reviews of Modern Physics 70, 1039 (1998).
  * Dagotto [1994] E. Dagotto, Reviews of Modern Physics 66, 763 (1994).
  * Sachdev [2003] S. Sachdev, Reviews of Modern Physics 75, 913 (2003).
  * Pisani [2003] C. Pisani, Journal of Molecular Structure (Theochem) 621, 141 (2003).
  * Pisani  _et al._ [2005] C. Pisani, M. Busso, G. Capecchi, S. Casassa, R. Dovesi, L. Maschio, C. Zicovich-Wilson, and M. Schütz, Journal of Chemical Physics 122 (2005), 10.1063/1.1857479.
  * Grüneis  _et al._ [2011] A. Grüneis, G. H. Booth, M. Marsman, J. Spencer, A. Alavi, and G. Kresse, Journal of Chemical Theory and Computation 7, 2780 (2011).
  * Ben  _et al._ [2012] M. D. Ben, J. Hutter, and J. Vandevondele, Journal of Chemical Theory and Computation 8, 4177 (2012).
  * Ben  _et al._ [2013] M. D. Ben, J. Hutter, and J. Vandevondele, Journal of Chemical Theory and Computation 9, 2654 (2013).
  * Booth  _et al._ [2013] G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Nature 493, 365 (2013).
  * Booth  _et al._ [2016] G. H. Booth, T. Tsatsoulis, G. K. L. Chan, and A. Grüneis, Journal of Chemical Physics 145 (2016), 10.1063/1.4961301.
  * McClain  _et al._ [2017] J. McClain, Q. Sun, G. K. L. Chan, and T. C. Berkelbach, Journal of Chemical Theory and Computation 13, 1209 (2017).
  * Neufeld  _et al._ [2022] V. A. Neufeld, H.-Z. Ye, and T. C. Berkelbach, The Journal of Physical Chemistry Letters 13, 7497 (2022).
  * Cui  _et al._ [2020] Z. H. Cui, T. Zhu, and G. K. L. Chan, Journal of Chemical Theory and Computation 16, 119–129 (2020).
  * Zhu  _et al._ [2020] T. Zhu, Z. H. Cui, and G. K. L. Chan, Journal of Chemical Theory and Computation 16, 141 (2020).
  * Cui  _et al._ [2022] Z.-H. Cui, H. Zhai, X. Zhang, and G. K.-L. Chan, Science 377, 1192 (2022).
  * Motta  _et al._ [2019] M. Motta, S. Zhang, and G. K.-L. Chan, Physical Review B 100, 045127 (2019).
  * Cotton [1991] F. A. Cotton,  _Chemical applications of group theory_ (John Wiley & Sons, 1991).
  * Crawford and Di Remigio [2019] T. D. Crawford and R. Di Remigio, in  _Annual Reports in Computational Chemistry_, Vol. 15 (Elsevier, 2019) pp. 79–101.
  * Georges  _et al._ [1996a] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996a).
  * Cortona [1991] P. Cortona, Phys. Rev. B 44, 8454 (1991).
  * Inglesfield [1981] J. Inglesfield, Journal of Physics C: Solid State Physics 14, 3795 (1981).
  * Knizia and Chan [2012a] G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012a).
  * Cui  _et al._ [2019] Z.-H. Cui, T. Zhu, and G. K.-L. Chan, Journal of Chemical Theory and Computation 16, 119 (2019).
  * Pham  _et al._ [2019] H. Q. Pham, M. R. Hermes, and L. Gagliardi, Journal of Chemical Theory and Computation 16, 130 (2019).
  * Zheng  _et al._ [2018] H. Zheng, H. J. Changlani, K. T. Williams, B. Busemeyer, and L. K. Wagner, Frontiers in Physics 6, 43 (2018).
  * Lee  _et al._ [2021] J. Lee, D. W. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush, PRX Quantum 2, 030305 (2021).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Ivanov  _et al._ [2022] A. V. Ivanov, C. Sünderhauf, N. Holzmann, T. Ellaby, R. N. Kerber, G. Jones, and J. Camps, arXiv:2210.02403 (2022).
  * Reiher  _et al._ [2017] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Proceedings of the National Academy of Sciences 114, 7555 (2017).
  * von Burg  _et al._ [2021] V. von Burg, G. H. Low, T. Häner, D. S. Steiger, M. Reiher, M. Roetteler, and M. Troyer, Physical Review Research 3, 033055 (2021).
  * Goings  _et al._ [2022] J. J. Goings, A. White, J. Lee, C. S. Tautermann, M. Degroote, C. Gidney, T. Shiozaki, R. Babbush, and N. C. Rubin, Proceedings of the National Academy of Sciences 119, e2203533119 (2022).
  * Bianchini  _et al._ [2019] M. Bianchini, M. Roca-Ayats, P. Hartmann, T. Brezesinski, and J. Janek, Angewandte Chemie Int. Ed. 58, 2 (2019).
  * Sicolo  _et al._ [2020] S. Sicolo, M. Mock, M. Bianchini, and K. Albe, Chemistry of Materials 32, 10096 (2020).
  * Chen  _et al._ [2011] H. Chen, C. L. Freeman, and J. H. Harding, Physical Review B 84, 085108 (2011).
  * Vandevondele  _et al._ [2005] J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Computer Physics Communications 167, 103 (2005).
  * Kühne  _et al._ [2020] T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, Journal of Chemical Physics 152, 194103 (2020).
  * Dovesi  _et al._ [2020] R. Dovesi, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. Bush, P. D’Arco, Y. Noel, M. Rera, P. Carbonniere, M. Causa, S. Salustro, V. Lacivita, B. Kirtman, A. M. Ferrari, F. S. Gentile, J. Baima, M. Ferrero, R. Demichelis, and M. D. L. Pierre, Journal of Chemical Physics 152, 204111 (2020); (2020).
  * Guidon  _et al._ [2009] M. Guidon, J. Hutter, and J. VandeVondele, Journal of Chemical Theory and Computation 5, 3010 (2009).
  * Guidon  _et al._ [2010] M. Guidon, J. Hutter, and J. VandeVondele, Journal of Chemical Theory and Computation 6, 2348–2364 (2010).
  * Pisani  _et al._ [1988] C. Pisani, R. Dovesi, and C. Roetti, “Hartree-fock ab initio treatment of crystalline systems,” (1988).
  * Blum  _et al._ [2009] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Computer Physics Communications 180, 2175 (2009).
  * Goedecker [1999] S. Goedecker, Reviews of Modern Physics 71, 1085 (1999).
  * Bowler and Miyazaki [2012] D. R. Bowler and T. Miyazaki, Reports on Progress in Physics 75, 036503 (2012).
  * Wu  _et al._ [2009] X. Wu, A. Selloni, and R. Car, Physical Review B - Condensed Matter and Materials Physics 79, 085102 (2009).
  * Lippert  _et al._ [1997] G. Lippert, J. Hutter, and M. Parinello, Molecular Physics 92, 477 (1997).
  * Whitten [1973] J. L. Whitten, The Journal of Chemical Physics 58, 4496 (1973).
  * Mintmire and Dunlap [1982] J. W. Mintmire and B. I. Dunlap, Physical Review A 25, 88 (1982).
  * Dunlap [2000] B. I. Dunlap, Journal of Molecular Structure: THEOCHEM 529, 37 (2000).
  * Weigend [2002] F. Weigend, Physical Chemistry Chemical Physics 4, 4285 (2002).
  * Varga [2005] S. Varga, Physical Review B 71, 073103 (2005).
  * Maschio and Usvyat [2008a] L. Maschio and D. Usvyat, Physical Review B1 78, 073102 (2008a).
  * Burow  _et al._ [2009] A. M. Burow, M. Sierka, and F. Mohamed, Journal of Chemical Physics 131 (2009), 10.1063/1.3267858.
  * Wang  _et al._ [2020] X. Wang, C. A. Lewis, and E. F. Valeev, The Journal of Chemical Physics 153, 124116 (2020).
  * Ye and Berkelbach [2021] H. Z. Ye and T. C. Berkelbach, Journal of Chemical Physics 154 (2021), 10.1063/5.0046617.
  * Hohenstein  _et al._ [2012a] E. G. Hohenstein, R. M. Parrish, and T. J. Martínez, The Journal of Chemical Physics 137, 1085 (2012a).
  * Parrish  _et al._ [2012] R. M. Parrish, E. G. Hohenstein, T. J. Martínez, and C. D. Sherrill, The Journal of Chemical Physics 137, 224106 (2012).
  * Hohenstein  _et al._ [2012b] E. G. Hohenstein, R. M. Parrish, C. D. Sherrill, and T. J. Martínez,  _The Journal of Chemical Physics_ , The Journal of Chemical Physics 137, 221101 (2012b).
  * Ye and Berkelbach [2022] H.-Z. Ye and T. C. Berkelbach, Journal of Chemical Theory and Computation 18, 1595 (2022).
  * Hartwigsen  _et al._ [1998] C. Hartwigsen, S. Goedecker, and J. Hutter, Physical Review B 58, 3641 (1998).
  * Heyd  _et al._ [2005] J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, The Journal of Chemical Physics 123, 174101 (2005).
  * Grüneis  _et al._ [2010] A. Grüneis, M. Marsman, and G. Kresse, The Journal of Chemical Physics 133, 074107 (2010).
  * Nadler and Kempier [1959] M. R. Nadler and C. Kempier, Analytical Chemistry 31, 2109 (1959).
  * Tang  _et al._ [2009] F. L. Tang, X. X. Che, W. J. Lu, G. B. Chen, Y. Xie, and W. Y. Yu, Physica B: Condensed Matter 404, 2489 (2009).
  * Babbush  _et al._ [2018c] R. Babbush, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, A. Paler, A. Fowler, and H. Neven, Physical Review X 8, 041015 (2018c).
  * Berry  _et al._ [2019] D. W. Berry, C. Gidney, M. Motta, J. McClean, and R. Babbush, Quantum 3, 208 (2019).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * Szegedy [2004] M. Szegedy, in  _45th Annual IEEE Symposium on Foundations of Computer Science_ (IEEE, 2004) pp. 32–41.
  * Childs and Wiebe [2012] A. M. Childs and N. Wiebe, Quantum Information and Computation 12, 901–924 (2012).
  * Baldereschi [1973] A. Baldereschi, Physical Review B 7, 5212 (1973).
  * Low  _et al._ [2018] G. H. Low, V. Kliuchnikov, and L. Schaeffer, arXiv: 1812.00954 (2018).
  * Werner  _et al._ [2003] H.-J. Werner, F. R. Manby, and P. J. Knowles, The Journal of chemical physics 118, 8149 (2003).
  * Maschio and Usvyat [2008b] L. Maschio and D. Usvyat, Physical Review B 78, 073102 (2008b).
  * Oumarou  _et al._ [2022] O. Oumarou, M. Scheurer, R. M. Parrish, E. G. Hohenstein, and C. Gogolin, arXiv:2212.07957 (2022).
  * Lu and Ying [2016] J. Lu and L. Ying, Annals of Mathematical Sciences and Applications 1, 321 (2016).
  * Wu  _et al._ [2021] K. Wu, X. Qin, W. Hu, and J. Yang, Journal of Chemical Theory and Computation 18, 206 (2021).
  * Hu  _et al._ [2017] W. Hu, L. Lin, and C. Yang, Journal of Chemical Theory and Computation 13, 5420 (2017).
  * Dong  _et al._ [2018] K. Dong, W. Hu, and L. Lin, Journal of Chemical Theory and Computation 14, 1311 (2018).
  * [84] “Project by BMW, BASF, Samsung SDI and Samsung Electronics to enhance sustainable cobalt mining,” https://www.basf.com/global/en/who-we-are/sustainability/responsible-partnering/cobalt-initiative.html (accessed October 26, 2020).
  * Olivetti  _et al._ [2017] E. A. Olivetti, G. Ceder, G. G. Gaustad, and X. Fu, Joule 1, 229 (2017).
  * Dahn  _et al._ [1990] J. R. Dahn, U. Vonsacken, and C. A. Michal, Solid State Ionics 44, 87 (1990).
  * Ohzuku  _et al._ [1993] T. Ohzuku, A. Ueda, and M. Nagayama, Journal of The Electrochemical Society 140, 1862 (1993).
  * Radin  _et al._ [2020] M. D. Radin, J. C. Thomas, and A. Van der Ven, Physical Review Materials 4, 043601 (2020).
  * Møller and Plesset [1934] C. Møller and M. S. Plesset, Physical Review 46, 618 (1934).
  * Cremer [2011] D. Cremer, Ltd. WIREs Comput Mol Sci 1, 509 (2011).
  * Purvis and Bartlett [1982] G. D. Purvis and R. J. Bartlett, The Journal of Chemical Physics 76, 1910 (1982).
  * Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, Reviews of Modern Physics 79, 291 (2007).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K. L. Chan, Wiley Interdisciplinary Reviews: Computational Molecular Science 8, e1340 (2018).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z. H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. Smith, C. Sun, S. N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K. L. Chan, The Journal of chemical physics 153, 024109 (2020).
  * Kim  _et al._ [2018] J. Kim, A. D. Baczewski, T. D. Beaudet, A. Benali, M. C. Bennett, M. A. Berrill, N. S. Blunt, E. J. L. Borda, M. Casula, D. M. Ceperley, _et al._ , Journal of Physics: Condensed Matter 30, 195901 (2018).
  * Kent  _et al._ [2020] P. R. C. Kent, A. Annaberdiyev, A. Benali, M. C. Bennett, E. J. Landinez Borda, P. Doak, H. Hao, K. D. Jordan, J. T. Krogel, I. Kylänpää, J. Lee, Y. Luo, F. D. Malone, C. A. Melton, L. Mitas, M. A. Morales, E. Neuscamman, F. A. Reboredo, B. Rubenstein, K. Saritas, S. Upadhyay, G. Wang, S. Zhang, and L. Zhao, The Journal of Chemical Physics 152, 174105 (2020).
  * Goedecker  _et al._ [1996] S. Goedecker, M. Teter, and J. Hutter, Physical Review B 54, 1703 (1996).
  * [98] J. Hutter, “New optimization of gth pseudopotentials for pbe, scan, pbe0 functionals. gth pseudopotentials for hartree-fock. nlcc pseudopotentials for pbe,” .
  * VandeVondele and Hutter [2007] J. VandeVondele and J. Hutter, Journal of Chemical Physics 127, 114105 (2007).
  * Georges  _et al._ [1996b] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Reviews of Modern Physics 68, 13 (1996b).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Reviews of Modern Physics 78, 865 (2006).
  * Held [2007] K. Held, Advances in Physics 56, 829 (2007).
  * Vollhardt [2020] D. Vollhardt (2020) p. 11001.
  * Knizia and Chan [2012b] G. Knizia and G. K.-L. Chan, Physical Review Letters 109, 186404 (2012b).
  * Knizia and Chan [2013] G. Knizia and G. K. L. Chan, Journal of Chemical Theory and Computation 9, 1428 (2013).
  * Liu  _et al._ [2023] Y. Liu, O. R. Meitei, Z. E. Chin, A. Dutt, M. Tao, T. Van Voorhis, and I. L. Chuang, arXiv:2301.01457 (2023).
  * Rubin  _et al._ [2022] N. C. Rubin, J. Lee, and R. Babbush, Journal of Chemical Theory and Computation 18, 1480 (2022).
  * Malone  _et al._ [2020] F. D. Malone, S. Zhang, and M. A. Morales, Journal of Chemical Theory and Computation 16, 4286 (2020).
  * Morales and Malone [2020] M. A. Morales and F. D. Malone, The Journal of Chemical Physics 153, 194111 (2020).
  * Chiesa  _et al._ [2006] S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holzmann, Physical Review Letters 97, 076404 (2006).
  * Drummond  _et al._ [2008] N. Drummond, R. Needs, A. Sorouri, and W. Foulkes, Physical Review B 78, 125106 (2008).
  * Azadi and Foulkes [2015] S. Azadi and W. Foulkes, The Journal of Chemical Physics 143, 102807 (2015).
  * Holzmann  _et al._ [2016] M. Holzmann, R. C. Clay III, M. A. Morales, N. M. Tubman, D. M. Ceperley, and C. Pierleoni, Physical Review B 94, 035126 (2016).
  * Dagrada  _et al._ [2016] M. Dagrada, S. Karakuzu, V. L. Vildosola, M. Casula, and S. Sorella, Physical Review B 94, 245108 (2016).
  * Mihm  _et al._ [2021] T. N. Mihm, T. Schäfer, S. K. Ramadugu, L. Weiler, A. Grüneis, and J. J. Shepherd, Nature Computational Science 1, 801 (2021).
  * Yang  _et al._ [2020] Y. Yang, V. Gorelov, C. Pierleoni, D. M. Ceperley, M. Holzmann, _et al._ , Physical Review B 101, 085115 (2020).
  * Sanders  _et al._ [2020] Y. R. Sanders, D. W. Berry, P. C. S. Costa, L. W. Tessler, N. Wiebe, C. Gidney, H. Neven, and R. Babbush, PRX Quantum 1, 020312 (2020).
  * DeYonker  _et al._ [2007] N. J. DeYonker, K. A. Peterson, G. Steyl, A. K. Wilson, and T. R. Cundari, Journal of Physical Chemistry A 111, 11269 (2007).
  * Lee and Taylor [1989] T. J. Lee and P. R. Taylor, International Journal of Quantum Chemistry 36, 199 (1989).
  * Nielsen and Janssen [1999] I. M. Nielsen and C. L. Janssen, Chemical Physics Letters 310, 568 (1999).
  * Lee [2003] T. J. Lee, Chemical Physics Letters 372, 362 (2003).
  * Lee  _et al._ [2022] J. Lee, H. Q. Pham, and D. R. Reichman, Journal of Chemical Theory and Computation 18, 7024 (2022).

##  Appendix A Sparse representation derivations

###  A.1 The Pauli operator representation of the one-body term

Here we derive the Pauli operator form of the one-body operator amenable to
implementation as a Majorana select operation. The one-body operator is
rewritten as

$\displaystyle\sum_{p,q=1}^{N/2}h_{p\mathbf{k},q\mathbf{k}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}\mapsto\frac{1}{4}\sum_{p,q=1}^{N/2}h_{p\mathbf{k},q\mathbf{k}}[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}+iY_{q\mathbf{k}\sigma})]$
$\displaystyle=\frac{1}{8}\sum_{p,q=1}^{N/2}h_{p\mathbf{k},q\mathbf{k}}[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}+iY_{q\mathbf{k}\sigma})]+\frac{1}{8}\sum_{p,q=1}^{N/2}h_{q\mathbf{k},p\mathbf{k}}[\vec{Z}(X_{q\mathbf{k}\sigma}-iY_{q\mathbf{k}\sigma})][\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})]$
$\displaystyle=\frac{1}{8}\sum_{p\neq
q=1}^{N/2}h_{p\mathbf{k},q\mathbf{k}}[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}+iY_{q\mathbf{k}\sigma})]-\frac{1}{8}\sum_{p\neq
q=1}^{N/2}h^{*}_{p\mathbf{k},q\mathbf{k}}[\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}-iY_{q\mathbf{k}\sigma})]$
$\displaystyle\quad+\frac{1}{4}\sum_{p=1}^{N/2}h_{p\mathbf{k},p\mathbf{k}}[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})]$
$\displaystyle=\frac{1}{8}\sum_{p\neq q=1}^{N/2}{\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})\left\\{[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}+iY_{q\mathbf{k}\sigma})]-[\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}-iY_{q\mathbf{k}\sigma})]\right\\}$
$\displaystyle\quad+\frac{i}{8}\sum_{p,q=1}^{N/2}{\rm
Im}(h_{p\mathbf{k},q\mathbf{k}})\left\\{[\vec{Z}(X_{p\mathbf{k}\sigma}-iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}+iY_{q\mathbf{k}\sigma})]+[\vec{Z}(X_{p\mathbf{k}\sigma}+iY_{p\mathbf{k}\sigma})][\vec{Z}(X_{q\mathbf{k}\sigma}-iY_{q\mathbf{k}\sigma})]\right\\}$
$\displaystyle\quad+\frac{1}{2}\sum_{p=1}^{N/2}h_{p\mathbf{k},p\mathbf{k}}(\openone_{p\mathbf{k}\sigma}-Z_{p\mathbf{k}\sigma})$
$\displaystyle=\frac{1}{4}\sum_{p\neq q=1}^{N/2}{\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})\left\\{-i\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+i\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}$
$\displaystyle\quad+\frac{i}{4}\sum_{p,q=1}^{N/2}{\rm
Im}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}+\frac{1}{2}\sum_{p=1}^{N/2}h_{pp}(\mathbf{k})(\openone_{p\mathbf{k}\sigma}-Z_{p\mathbf{k}\sigma})$
$\displaystyle=\frac{i}{4}\sum_{p\neq q=1}^{N/2}{\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}X_{q\mathbf{k}\sigma}\vec{Z}Y_{p\mathbf{k}\sigma}+\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}$
$\displaystyle\quad+\frac{i}{4}\sum_{p,q=1}^{N/2}{\rm
Im}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}+\frac{1}{2}\sum_{p=1}^{N/2}h_{pp}(\mathbf{k})(\openone_{p\mathbf{k}\sigma}-Z_{p\mathbf{k}\sigma})$
$\displaystyle=\frac{i}{2}\sum_{p,q=1}^{N/2}{\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}+\frac{i}{4}\sum_{p,q=1}^{N/2}{\rm
Im}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}+\frac{1}{2}\sum_{p=1}^{N/2}h_{p\mathbf{k},p\mathbf{k}}\openone.$
(74)

In the last line we have used the symmetry of ${\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})$ to combine
$\vec{Z}X_{q\mathbf{k}\sigma}\vec{Z}Y_{p\mathbf{k}\sigma}$ and
$\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}$, then used the fact
that $iXY=-Z$ to combine the sum with $p\neq q$ with that for $p$. The
complete expression for the Hamiltonian has the sum over $\sigma$ and
$\mathbf{k}$, which we have left out for simplicity here. Including those
gives the expression in Eq. (III.1).

###  A.2 One-body correction for sparse case

Next we derive the effective one-body term from the two-electron part of the
Hamiltonian. In the case $p=q$ and $\mathbf{Q}=0$, the second term in square
brackets in Eq. (20) can be written as

$\displaystyle-V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle=V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle\quad-V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}(a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}+a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma})a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle=V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle\quad-V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}.$
(75)

In the last line we have used the fact that for $p=q$ and $\mathbf{Q}=0$,
$a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}+a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}$
is just the identity, so this becomes a one-body operator.

Similarly, if $r=s$ and $\mathbf{Q}=0$ (but $p\neq q$), the second term in
square brackets in Eq. (20) can be written as

$\displaystyle-V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}$
$\displaystyle=V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle\quad-V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}(a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}+a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau})$
$\displaystyle=V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle\quad-V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}.$
(76)

Thus we see that in either case ($p=q$ or $r=s$), we have the same expression
as in Eq. (21), plus a one-body operator. Moreover, because of the symmetry of
$V$ (in swapping the $pq$ pair with the $rs$ pair), these corrections are
equal. Note also that we can relabel swapping $p$ with $q$ and $r$ with $s$ to
replace
$V^{*}_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q\mathbf{k}\sigma}^{\dagger}$
with (now explicitly taking $\mathbf{Q}=0$)

$V^{*}_{q\mathbf{k},p\mathbf{k},s\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}a_{q\mathbf{k}\sigma}a_{p\mathbf{k}\sigma}^{\dagger}=-V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}.$
(77)

This means that the contribution of these corrections is

$\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p,q=1}^{N/2}\left(\sum_{r=1}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}.$
(78)

In this expression the constant factor is determined as follows. There is a
factor of $1/4$ in Eq. (20). Next, there is a factor of 2 because we have the
contribution from $p=q$ as well as that from $r=s$. Last, there is the factor
of 2 from the summation over the spin $\tau$. As a result, these factors
cancel to give 1 above. Therefore, for $p\neq q$, we can combine this one-body
term with $h_{pq}$ as

$h^{\prime}_{pq}=h_{pq}+\sum_{r=1}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}.$
(79)

Next we consider the case where $p=q$, $r=s$, and $\mathbf{Q}=0$. Then the
second term in square brackets in Eq. (20) can be written as

$\displaystyle
V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}$
$\displaystyle=V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle\quad+V^{*}_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}(a_{p\mathbf{k}\sigma}^{\dagger}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}^{\dagger}a_{s\mathbf{k}^{\prime}\tau}-a_{p\mathbf{k}\sigma}a_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}^{\dagger}a_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q})\tau}a_{s\mathbf{k}^{\prime}\tau}^{\dagger}).$
(80)

The operators in brackets in the final line can be written as, taking $p=q$,
$r=s$, and $\mathbf{Q}=0$,

$\displaystyle
a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}a_{r\mathbf{k}^{\prime}\tau}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}+a_{p\mathbf{k}\sigma}a_{p\mathbf{k}\sigma}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}-a_{p\mathbf{k}\sigma}a_{p\mathbf{k}\sigma}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}-a_{p\mathbf{k}\sigma}a_{p\mathbf{k}\sigma}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}a_{r\mathbf{k}^{\prime}\tau}^{\dagger}$
$\displaystyle=a_{r\mathbf{k}^{\prime}\tau}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}-a_{p\mathbf{k}\sigma}a_{p\mathbf{k}\sigma}^{\dagger}$
$\displaystyle=a_{r\mathbf{k}^{\prime}\tau}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}+a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}-\openone.$
(81)

By symmetry of swapping $p$ and $q$, and swapping $r$ and $s$, we must be able
to simplify the final line of (80) to

$V_{p\mathbf{k},p\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}(a_{r\mathbf{k}^{\prime}\tau}^{\dagger}a_{r\mathbf{k}^{\prime}\tau}+a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}-\openone).$
(82)

That is, this value of $V$ is real. We can also relabel $r$ and $p$ and use
symmetry to show the contribution from the first term in Eq. (82) is
equivalent to

$V_{r\mathbf{k}^{\prime},r\mathbf{k}^{\prime},p\mathbf{k},p\mathbf{k}}a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}=V_{p\mathbf{k},p\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}.$
(83)

Hence the contribution of these corrections is

$\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p=1}^{N/2}\left(\sum_{r=1}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}V_{p\mathbf{k},p\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right)(a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}-\openone/2).$
(84)

In this case, the constant factor comes from $1/4$ in Eq. (20), and a factor
of 2 from the sum over $\tau$. As a result the expression in Eq. (82) is
divided by 2 here, and apart from the identity we have the same expression as
that accounting for only one of the pairs $p,q$ and $r,s$ being equal. The
operator proportional to the identity can be ignored in the implementation of
the Hamiltonian because it just gives a global shift in the eigenvalues.

###  A.3 Complexity for sparse implementation

The fundamental operator we are aiming to implement is in the form of Eq.
(III.1) for the one-body term and Eq. (III.1) for the two-body term. In both
we have a real part and an imaginary part; for the one-body term this is
$h_{pq}$, and for the two-body term this is
$V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime}}$.
We need to perform a state preparation that provides _real_ amplitudes for the
real and imaginary parts of $h$ and $V$ on separate basis states (not just
real and imaginary parts of an amplitude on each basis state). This means the
number of items of data to output is doubled in order to give the real and
imaginary parts. The state preparation is otherwise essentially unchanged from
that in [71], as described in Eq. (48) of that work and the accompanying
explanation.

Recall that in the sparse state preparation procedure, we use a register
indexing the nonzero entries (see Eq. (43) of [71]). That is used to output
“ind”, “alt”, and “keep” values via QROM (see Eq. (44) of [71]). The “ind”
values are values of $p,q,r,s$, as well as the sign needed, and a qubit
distinguishing between the one- and two-body terms. The “alt” values are
alternate values of these quantities, and “keep” governs the probability of
swapping these registers for the state preparation via coherent alias
sampling. Since we need a bit to flag whether the amplitude being produced is
for the real or imaginary part, that would indicate we need two extra bits
output, one for the “ind” value and one for the “alt” value. However, we can
use one bit in the register indexing the nonzero entries to flag between real
and imaginary parts. It is most convenient to make this register the least
significant bit. Then we just need to produce “alt” values of this register,
so the output size is only increased by 1 bit instead of 2. A requirement for
this approach is that the non-zero entries of $V$ that are retained are the
same for the real and imaginary parts.

A further increase in the size of the output register is because we need to
output values of $\mathbf{k}$, $\mathbf{k}^{\prime}$, and $\mathbf{Q}$. The
number of bits needed to store $\mathbf{k}$ is not simply $\lceil\log
N_{k}\rceil$ because $\mathbf{k}$ is a vector. The number of bits will be
denoted $n_{k}$. If we assume that the number of values is given by the
product of numbers in the three dimensions $N_{k}=N_{x}N_{y}N_{z}$, then

$n_{k}=\lceil\log N_{x}\rceil+\lceil\log N_{y}\rceil+\lceil\log N_{z}\rceil.$
(85)

Therefore $\mathbf{k}$, $\mathbf{k}^{\prime}$, and $\mathbf{Q}$ increase the
size of both the ind and alt registers by $3n_{k}$, for a total of $6n_{k}$.
The size of the output is given in Eq. (A13) of [32] as
$m=\aleph+8\lceil\log(N/2)\rceil+4$, and would here be increased to

$\aleph+8\lceil\log(N/2)\rceil+6n_{k}+5,$ (86)

where we have also increased the size of the output by 1 to account for
selecting between real and imaginary parts, as discussed above. The quantity
$\aleph$ is the number of bits for the “keep” register.

The remaining consideration for the sparse state preparation is the symmetry.
In prior work there were three symmetries, with swap of $p,q$ with $r,s$ as
well as swaps within the $p,q$ and $r,s$ pairs. The method to take advantage
of this was described from about Eq. (49) on in [71]. There you only perform
the preparation for a restricted range of $p,q,r,s$, then use three qubits to
control swaps to generate the symmetries.

Here we have the symmetry with swap of $p,q$ with $r,s$, but we can only swap
the $p,q$ and $r,s$ pairs simultaneously. We also need to take the complex
conjugate when performing that swap. In order to implement the symmetries
here, we will have two control qubits. One qubit will be in a
$\mathinner{|{+}\rangle}$ state and control swap of the $p,q$ with $r,s$ as
before, except we now have the registers containing
$\mathbf{k},\mathbf{k}^{\prime},\mathbf{k}{\ominus}\mathbf{Q},\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$
to swap. That qubit is only set to $\mathinner{|{+}\rangle}$ for the two-body
term, since that symmetry does not make sense for the one-body term. The
second qubit is used to simultaneously swap the $p,q$ and $r,s$ pairs, as well
as the registers containing $\mathbf{k}$, etc. It will also be used as a
control for a $Z$ phase gate on a qubit flagging imaginary components. That is
a Clifford gate and is not included in the Toffoli count.

The net result is that the cost of the swaps to produce these symmetries is
unchanged from that in [71], except in that we are counting the qubits needed
to store $\mathbf{k}$, etc, as well as $p,q,r,s$. Since a controlled swap of
two qubits can be performed with a single Toffoli (and Clifford gates), the
Toffoli cost of the two controlled swaps of registers is the total number of
qubits used to store $p,q,r,s$ as well as
$\mathbf{k},\mathbf{k}^{\prime},\mathbf{k}{\ominus}\mathbf{Q},\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$,
which is $4\lceil\log(N/2)\rceil+4n_{k}$. Note that in the state preparation
we will be producing $\mathbf{k},\mathbf{k}^{\prime},\mathbf{Q}$, and need to
compute $\mathbf{k}{\ominus}\mathbf{Q}$ and
$\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$ before performing the swaps for these
symmetries.

Assuming for the moment that $N_{x},N_{y},N_{z}$ are all powers of two, then
the number of Toffolis needed for the modular subtractions of the three
components will be $n_{k}-3$, unless one or more of $N_{x},N_{y},N_{z}$ is
equal to 1. It is simpler to give the cost as $n_{k}$ Toffolis, to avoid
needing to address special cases. A further complication is when one or more
of $N_{x},N_{y},N_{z}$ are not powers of two. In this case, the subtraction
can be performed in the usual way for two’s complement binary. Then you can
check if the result for any component is negative, and if it is then add the
appropriate $N_{x},N_{y},N_{z}$ to make it non-negative. The controlled
addition of a classically given number has complexity $n_{k}$, so this at
worst doubles the complexity to $2n_{k}$ for the modular subtraction.

The other major feature that we need to account for is the modified select
operation needed. The basic circuit primitive was given in Figure 13 of [32],
in order to apply $\vec{Z}Y_{p,\sigma}$ followed by $\vec{Z}X_{q,\sigma}$. A
more complicated circuit primitive was given in Figure 1 of [71], which
included testing $p=q$ which is not needed in the approach of [32]. Here the
scheme is more complicated, because instead of having a fixed sequence where
we need to apply $Y$ followed by $X$ we have every combination. This can be
achieved by simply performing each twice; once with a controlled $\vec{Z}Y$
and once with a controlled $\vec{Z}X$, with a doubling of the Toffoli
complexity. That can be seen easily from the diagram in Figure 9 of [2]. There
a control qubit is used, so that can be used to control application of this
circuit with $Y$, then to control application of this circuit with $X$.

To understand how $X$ versus $Y$ is selected, note that there are effectively
five bits controlling here. Let us call the bit selecting between the one- and
two-body terms $b_{0}$; this is created in the sparse state preparation. Let
us call the bit selecting real versus imaginary parts $b_{1}$; this is again
created in the state preparation. There also needs to be a bit $b_{2}$ for
selecting between the two lines for real and the two lines for imaginary in
the expression in Eq. (III.1). Then we have $b_{3}$ to select between the two
terms in the first set of square brackets in each line of Eq. (III.1), and a
bit $b_{4}$ selecting between the two terms in the second set of square
brackets.

Now, considering the operators indexed by $r,s$ first, these are applied for
the two-body terms but not the one-body term. This control of the operations
adds only one Toffoli to the cost. For the first operation,
$\vec{Z}X_{s\mathbf{k}^{\prime}\tau}$ or
$\vec{Z}Y_{s\mathbf{k}^{\prime}\tau}$, we can see that the selection between
$X$ and $Y$ depends only on bit $b_{4}$. For the second operation, the
selection is independent of whether we have the real or imaginary part. We
select $X$ if we have $b_{4}=0$ (the first term) and $b_{2}=0$ (the first
line), or if we have $b_{4}=b_{2}=1$. To create a bit selecting between $X$
and $Y$ we can simply perform a CNOT between these bits, with no Toffoli cost.

Next, consider the operators indexed by $p,q$. For simplicity we will first
consider just the two-body terms. Again the first operation can select between
$X$ and $Y$ just by using the bit $b_{3}$ selecting between the terms. Then
for the second operation, we select $X$ if we have $b_{1},b_{2},b_{3}$ equal
to $0,0,0$, or $0,1,1$, or $1,0,1$, or $1,1,0$. It is easily seen that if we
apply CNOTs with $b_{1}$ then $b_{2}$ as control and $b_{3}$ as target, then
we should apply $X$ if we have $b_{3}=0$. This selection can be performed
without Toffolis again.

Now to take account of how the one-body terms are applied, it is convenient to
rewrite the first line of Eq. (III.1) as

$-\frac{i}{4}\sum_{p,q=1}^{N/2}{\rm
Re}(h_{p\mathbf{k},q\mathbf{k}})\left\\{\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}-\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right\\}.$
(87)

Then the selection between the operations is identical to that for $b_{2}=1$
(second lines) for the two-body part. Therefore, for the above analysis of the
two-body implementation, we can replace $b_{2}$ with a bit that is 1 if
$b_{2}=1$ OR $b_{0}=0$. This operation requires one more Toffoli.

Another modification we need to make is to compute
$\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$ and $\mathbf{k}{\ominus}\mathbf{Q}$
to use in the selection for the two-body operations. As explained above, these
modular subtractions have complexity at worst $2n_{k}$. The calculation
$\mathbf{k}{\ominus}\mathbf{Q}$ needs to be controlled on the bit $b_{0}$
selecting between the one- and two-body terms, which increases its complexity
by $n_{k}$. Therefore the complexity of this arithmetic is $3\lceil\log
N_{k}\rceil$ Toffolis. We can keep the working qubits in order to uncompute
this arithmetic with Clifford gates.

Finally, we should account for the phase factors needed in the implementation.
The phase factors needed are as follows.

  1. 1.

We should apply an $i$ phase factor on the one-body term. That can be
implemented with an $S$ gate which is Clifford.

  2. 2.

If we have the one-body term (flagged by $b_{0}=0$) we should flip the sign of
the real part (flagged by $b_{1}=0$). This can be done with a controlled
phase, which is again Clifford.

  3. 3.

For the two-body term ($b_{0}=1$), real ($b_{1}=0$), and second line
($b_{2}=1$) we should flip the sign. This doubly controlled phase has a cost
of one Toffoli.

  4. 4.

We should flip the sign with $b_{3}=1$ if we have the two-body term
($b_{0}=1$) and the second line for real ($b_{1}=0,b_{2}=1$) or the first line
for imaginary ($b_{1}=1,b_{2}=0$). We should also flip the sign with $b_{3}=1$
if we have the one-body term ($b_{0}=0$) and real ($b_{1}=0$). To achieve this
we can first perform a CNOT with $b_{1}$ as control and $b_{2}$ as target.
Then, if $b_{0}=0,b_{1}=0$ OR $b_{0}=1,b_{2}=1$ we should apply a $Z$ gate to
the qubit containing $b_{3}$. This can be achieved with two double controlled
phase gates, so has Toffoli cost 2.

  5. 5.

We should flip the sign for $b_{4}=1$ if we have the second line $b_{2}=1$.
That is just a controlled phase with no non-Clifford cost.

As a result, the total complexity of implementing these phase factors is 3
Toffoli gates.

The total additional complexity is therefore $2$ Toffolis for the selection of
$X$ versus $Y$ when we account for needing to perform the one-or two-body
term, $3\lceil\log N_{k}\rceil$ Toffolis for subtractions, 3 Toffolis for
phase factors, and doubling the selection cost to select between $X$ and $Y$.
The two Toffolis to account for the one-body term were one for selecting
performing the operators indexed by $r,s$, and another Toffoli to perform an
OR between $b_{0}$ and $b_{2}$ for the operators indexed by $p,q$.

A further complication arises where the $h$ and $V$ are dependent on the spins
$\sigma$ and $\tau$. This is easily accounted for by outputting the values of
$\sigma,\tau$ as part of the state preparation. This means that the size of
both the “ind” and “alt” outputs are increased by 2, making the total size of
the output increase by 4 to be

$\aleph+8\lceil\log(N/2)\rceil+6n_{k}+9.$ (88)

Often there is the symmetry that for $V$ the value with
$\sigma=\uparrow,\tau=\downarrow$ are the same as for
$\sigma=\downarrow,\tau=\uparrow$. This means that we can omit the case
$\sigma=\downarrow,\tau=\uparrow$, and use a swap of these two qubits
controlled by an ancilla qubit in the usual way for obtaining symmetries. In
the detailed costing below, we give results for the case where $h$ and $V$ are
not dependent on spin for simplicity.

The QROM output size is

$m=\aleph+8n_{N}+6n_{k}+5,$ (89)

where $n_{N}=\lceil\log(N/2)\rceil$. This output size is increased above that
analysed in [32]. Then, using that output size, the formula for the cost of
the preparation with $d$ unique nonzero entries is

$\lceil d/k_{1}\rceil+m(k_{1}-1)$ (90)

and of the inverse preparation is

$\lceil d/k_{2}\rceil+k_{2}.$ (91)

Here $k_{1}$ and $k_{2}$ must be chosen as powers of 2. This formula is the
same as in [32], but with the modified value of $m$.

To begin the state preparation, we need to prepare an equal superposition
state over $d$ basis states. The analysis is described in [32], which gives
the costing $3\lceil\log d\rceil-3\eta+2b_{r}-9$ Toffoli gates. Here $\eta$ is
a number such that $2^{\eta}$ is a factor of $d$, and $b_{r}$ is a number of
bits used for rotation of an ancilla qubit to improve the amplitude of
success. This is a cost needed both for the preparation and inverse
preparation.

Other minor Toffoli costs are as follows. We use extra ancillas to save cost,
because a large number of ancillas were used for the QROM, and can be reused
here without increasing the maximum number of ancillas needed. In the
following we use the notation $n_{N}=\lceil\log(N/2)\rceil$.

  1. 1.

Perform select as shown in Figure 13 of [32] twice, but controlling between
$X$ and $Y$. This complexity is $4NN_{k}-6$, since we have $8$ times a
complexity of $NN_{k}/2-1$ for each of the selected operations, plus 2
Toffolis to generate the qubits we need for the control. There were two
Toffolis needed to account for selecting between one- and two- body terms, and
otherwise the selection to account for the various terms can be performed
using Clifford gates. In addition to this, we need to perform swaps controlled
by spin qubits twice for each of the two spin qubits, with a complexity
$2NN_{k}$. That then gives a total complexity of this step $6NN_{k}-6$.

  2. 2.

The state preparation needs an inequality test on $\aleph$ qubits, as well as
controlled swaps. The controlled swaps are on $n_{N}+3n_{k}+2$ qubits. Here
$4\lceil\log(N/2)\rceil$ are for the values of $p$, $q$, $r$, and $s$, the
$3n_{k}$ is for the $\mathbf{k}$, $\mathbf{k}^{\prime}$, and $\mathbf{Q}$
values, a $+1$ is for the qubit which distinguishes between the one- and two-
electron terms, and a further $+1$ comes from the qubit for selecting between
the real and imaginary parts. There are also ind and alt values of the sign,
but the correct phase can be applied with Clifford gates, so this does not add
to the Toffoli cost. The cost of the inequality test on $\aleph$ qubits is
$\aleph$. As in [32], we can eliminate the non-Clifford cost of the inverse
preparation using ancillas and measurements, so the Toffoli cost is
$\aleph+4n_{N}+3n_{k}+2$.

  3. 3.

The controlled swaps used to generate the symmetries have a cost of
$4n_{N}+4n_{k}$. This is increased by $4n_{k}$ over that in [32], since we
need to swap $\mathbf{k}$ registers as well. Although there are only two
controlled swaps rather than three in [32], two of the controlled swaps in
[32] together act on as many qubits as one controlled swap here, so the factor
of 4 is the same as in [32]. A further $4n_{k}$ cost is for computing
$\mathbf{k}{\ominus}\mathbf{Q}$ and $\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$
(or $2n_{k}$ if $N_{x},N_{y},N_{z}$ are powers of 2), and an extra $n_{k}$ is
needed to make the computation of $\mathbf{k}{\ominus}\mathbf{Q}$ controlled.
Again these controlled swaps can be inverted for the inverse preparation with
measurements and Clifford gates. Thus the total Toffoli cost here is
$4n_{N}+9n_{k}$.

  4. 4.

For the qubitization construction a reflection on the ancilla is needed as
well. The qubits that need to be reflected on are

    1. (a)

the $\lceil\log d\rceil$ qubits for preparing the state,

    2. (b)

$\aleph$ qubits for the equal superposition state in coherent alias sampling,

    3. (c)

two qubits that are used for controlled swaps to generate the symmetries of
the state,

    4. (d)

the two spin qubits,

    5. (e)

the ancilla qubit that is rotated to produce the equal superposition state,

    6. (f)

and the qubits storing $b_{2},b_{3},b_{4}$, which are also used in the linear
combination of unitaries.

There is no non-Clifford Toffoli cost for the preparation on
$b_{2},b_{3},b_{4}$, since an equal superposition may be prepared with a
Hadamard. They are control qubits that need to be reflected upon for the
qubitization, so add a cost of 3 Toffolis to the reflection giving a total
cost $\lceil\log d\rceil+\aleph+6$.

  5. 5.

As before, the control for the phase estimation uses unary iteration on the
control registers, with one more Toffoli for each step. The control by these
registers is implemented simply by controlling the reflection, which needs
just one Toffoli per step.

  6. 6.

An extra three Toffolis are needed for the phase factors.

Adding all these minor costs together gives, in the spin-independent case

$\displaystyle 2(3\lceil\log
d\rceil-3\eta+2b_{r}-9)+(6NN_{k}-6)+(\aleph+4n_{N}+3n_{k}+2)+4n_{N}+9n_{k}+\lceil\log
d\rceil+\aleph+6+2+3$ $\displaystyle=6NN_{k}+8n_{N}+10\lceil\log
N_{k}\rceil+2\aleph+7\lceil\log d\rceil-6\eta+4b_{r}-8.$ (92)

The total cost for a single step is then

$\left\lceil\frac{d}{k_{1}}\right\rceil+m(k_{1}-1)+\left\lceil\frac{d}{k_{2}}\right\rceil+k_{2}+6NN_{k}+8n_{N}+12n_{k}+2\aleph+7\lceil\log
d\rceil-6\eta+4b_{r}-8,$ (93)

with $m=\aleph+8n_{N}+6\lceil\log N_{k}\rceil+5$,
$n_{N}=\lceil\log(N/2)\rceil$, $\eta$ an integer such that $2^{\eta}$ is a
factor of $d$, and $b_{r}$ the number of bits used for rotation of an ancilla
qubit.

We may count the qubit costs by considering the maximum used during the QROM,
as the advanced QROM has a high qubit usage that will not be exceeded in other
parts of the algorithm. The qubit costs are therefore as follows.

  1. 1.

The control register for the phase estimation uses
$\lceil\log(\mathcal{I}+1)\rceil$ qubits, and there are
$\lceil\log(\mathcal{I}+1)\rceil-1$ qubits for the unary iteration.

  2. 2.

The system uses $NN_{k}$ qubits.

  3. 3.

The $\lceil\log d\rceil+\aleph+8$ qubits that need to be reflected upon listed
above.

  4. 4.

A qubit is needed to flag success of the equal superposition state
preparation.

  5. 5.

The phase gradient state uses $b_{r}$ qubits.

  6. 6.

The QROM uses qubits (including the output)
$mk_{1}+\lceil\log(d/k_{1})\rceil$.

This gives a total number of logical qubits

$2\lceil\log(\mathcal{I}+1)\rceil+NN_{k}+\lceil\log
d\rceil+b_{r}+\aleph+mk_{1}+\lceil\log(d/k_{1})\rceil+8,$ (94)

with $m=\aleph+8n_{N}+6\lceil\log N_{k}\rceil+5$.

##  Appendix B Single-factorization derivations

###  B.1 One-body correction for single factorization

For the single factorized form of the Hamiltonian, we may use the same
expressions for $\hat{A}$ and $\hat{B}$ for the case $\mathbf{Q}=0$ as for
$\mathbf{Q}\neq 0$, with an additional correction proportional to the
identity. This yields a one-body correction in the case of $\hat{A}$ but not
$\hat{B}$. For $\hat{A}_{n}(\mathbf{Q}=0)$ we obtain a term proportional to
the identity, as follows

$\displaystyle\hat{A}_{n}(\mathbf{Q}=0)$
$\displaystyle=\frac{1}{2}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}\sum_{p\neq
q}\left(L_{p\mathbf{k},q\mathbf{k},n}a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}+L_{p\mathbf{k},q\mathbf{k},n}^{*}a_{q\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}\right)+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}\sum_{p}L_{p\mathbf{k}p\mathbf{k},n}a_{p\mathbf{k}\sigma}^{\dagger}a_{p\mathbf{k}\sigma}$
(95)
$\displaystyle=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p\neq
q}^{N/2}\left(\frac{i{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right.$
$\displaystyle\quad+\left.\frac{i{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right)+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p}^{N/2}\frac{L_{p\mathbf{k},p\mathbf{k},n}}{2}(\openone-Z)$
$\displaystyle=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{pq}^{N/2}\left(\frac{i{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right.$
$\displaystyle\quad+\left.\frac{i{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]}{4}\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q(\mathbf{k}{\ominus}\mathbf{Q})\sigma}\right)\right)+\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p}^{N/2}\frac{L_{p\mathbf{k}p\mathbf{k},n}}{2}\openone.$
(96)

Here we have used the symmetry of $L$, so $L_{p\mathbf{k}p\mathbf{k},n}$ is
real. This derivation is similar to that for the one-body term in Appendix
A.1.

Because $\hat{A}_{n}(\mathbf{Q}=0)$ is squared, the identity term gives rise
to a one-body correction

$\displaystyle\frac{i}{4}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{n}^{M}\sum_{\mathbf{k}}^{N_{k}}\sum_{p,q}^{N/2}\left({\rm
Re}[L_{p\mathbf{k}q\mathbf{k},n}]\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}\right)\right.$
$\displaystyle\quad+\left.{\rm
Im}[L_{p\mathbf{k}q\mathbf{k},n}]\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right)\right)\sum_{\mathbf{k}^{\prime}}^{N_{k}}\sum_{r=1}^{N/2}L_{r\mathbf{k}^{\prime}r\mathbf{k}^{\prime},n}$
$\displaystyle=\frac{i}{4}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\mathbf{k}}^{N_{k}}\sum_{p,q}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}\sum_{r=1}^{N/2}\left({\rm
Re}[V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}]\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}-\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}\right)\right.$
$\displaystyle\quad+\left.{\rm
Im}[V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}]\left(\vec{Z}X_{p\mathbf{k}\sigma}\vec{Z}X_{q\mathbf{k}\sigma}+\vec{Z}Y_{p\mathbf{k}\sigma}\vec{Z}Y_{q\mathbf{k}\sigma}\right)\right).$
(97)

Here there was a factor of $1/2$ on the square of $\hat{A}_{n}(\mathbf{Q}=0)$,
a factor of 2 from the cross term in the square, a factor of 2 from the sum
over the spin on the identity, and so a factor of $1/2$ has been cancelled.
The form of this correction is identical to that for the one-body term, except
$h_{p\mathbf{k},q\mathbf{k}}$ is replaced with

$\sum_{r=1}^{N/2}\sum_{\mathbf{k}^{\prime}}^{N_{k}}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}.$
(98)

For $\hat{B}_{n}(\mathbf{Q}=0)$, it is easily seen that the symmetry
$L_{p\mathbf{k}q\mathbf{k},n}=L^{*}_{q\mathbf{k}p\mathbf{k},n}$ implies that
$\hat{B}_{n}(\mathbf{Q}=0)=0$. If we use the form for
$\hat{B}_{n}(\mathbf{Q}=0)$ in terms of Pauli operators given in Eq. (III.2),
then it will be proportional to the identity due to the case $p=q$. Squaring
then just gives a correction proportional to the identity (which can be
ignored in the implementation because it is just an energy shift), and it
gives no one-body correction. As a result we add the expression in Eq. (98) to
$h_{pq}$ to obtain the complete one-body Hamiltonian given in Eq. (38).

###  B.2 Complexity for single-factorized representation

To see the changes we need to make to the algorithm for the single-factorized
representation, recall that the two-body term was of the form [32]

$W^{\prime}=\frac{1}{8}\sum_{\ell=1}^{L}\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{p,q=1}^{N/2}W^{(\ell)}_{pq}Q_{pq\sigma}\right)^{2},$
(99)

where $Q_{pq\sigma}$ was an individual Pauli string. So the changes in the
representation are

  * •

The sum over $\ell$ up to $L$ has been replaced with a sum over $\mathbf{Q}$
and $n$, as well as a sum over the squares of $A$ and $B$.

  * •

Inside the square, the sum over just $\sigma,p,q$ now also has a sum over
$\mathbf{k}$.

  * •

Inside the sum, instead of just having a single Pauli string, we have a sum
over $4$, with real and imaginary parts of
$L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}$.

The amendments we will make to the original algorithm (according to the
description in [32]) to implement the block encoding are as follows.

  * •

For the sum over $\mathbf{Q}$ and $n$ we can combine them into $\ell$, and use
the same state preparation method as before. The value of $\mathbf{Q}$ will
need to be used in the select operation, so needs to be output as part of that
state preparation.

  * •

In the preparation for the block encoding of $A$ and $B$, the index
$\mathbf{k}$ will be needed as well as $p$ and $q$.

  * •

We no longer take advantage of $p,q$ symmetry.

  * •

We need to perform arithmetic to compute $\mathbf{k}{\ominus}\mathbf{Q}$ and
$\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$, with a cost of $4n_{k}$ (or $2n_{k}$
if $N_{x},N_{y},N_{z}$ are powers of 2).

  * •

A number of qubits can be used for selecting between the parts of the linear
combination of unitaries, similar to the sparse case. We have $b_{0}$ to
select between the one- and two-body terms, $b_{1}$ for selecting between the
real and imaginary parts, and $b_{3}$ selecting between the two terms in one
application of $A$ or $B$. The qubit $b_{2}$ can be used for selecting between
$A$ and $B$, which is a change from the sparse case, where it was used for
selecting between lines. We do not need $b_{4}$ because we are implementing
$A$ or $B$ twice (and creating the bit $b_{3}$ both times).

  * •

There needs to be a doubling of the selection cost to select between $X$ and
$Y$ as in the sparse case.

  * •

The creation of the qubits for controlling between $X$ and $Y$ can be
performed with one additional Toffoli. Note first that the terms in $A$ are
equivalent to the one-body part, and the terms in $B$ are the same except with
the real and imaginary lines swapped around. This means that we can use
$b_{0}$ and $b_{2}$ as a control to flip $b_{1}$, which effectively swaps the
real and imaginary parts for $B$ so it can be implemented in the same way.
Now, for the first selection of $X$ versus $Y$, we can apply a CNOT with
$b_{1}$ as control and $b_{3}$ as target, and use that as control For the
second selection we can simply use $b_{3}$ as control.

  * •

For the phase factors, we just need a sign flip if $b_{1}=0$ and $b_{3}=1$,
which is a Clifford controlled phase.

To explain the modifications needed for the costings, here we give the
sequence of steps with the same numbering as in [32], explaining the
differences.

  1. 1.

We first prepare a state as

$\frac{1}{\sqrt{\lambda}}\left(\mathinner{|{0,0,0,0}\rangle}\sqrt{\sum_{p,q}\left(|{\rm
Re}(h^{\prime}_{pq})|+|{\rm
Im}(h^{\prime}_{pq})|\right)}+\frac{1}{\sqrt{2}}\sum_{\mathbf{Q},n}\mathinner{|{\ell,\mathbf{Q},n,1}\rangle}\sum_{\mathbf{k},pq}(|{\rm
Re}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]|+|{\rm
Im}[L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n}]|)\right),$ (100)

where $\mathinner{|{\ell,\mathbf{Q},n}\rangle}$ indicates $\ell$ which starts
from 1 indexing values of $\mathbf{Q},n$, but $\mathbf{Q}$ and $n$ are also
output in registers. That is, we will be preparing $\ell$ while outputting
values of $\mathbf{Q},n$. We are assuming the more difficult case where the
number of values of $\mathbf{Q}$ or $n$ are not powers of 2, but if they are
then further simplifications are possible. This has complexity as follows.

    1. (a)

Preparing an equal superposition on $MN_{k}+1$ basis states has complexity
$3n_{MN}+2b_{r}-9$, where $b_{r}$ is the number of bits used for the rotation
on the ancilla,

$n_{MN}=\lceil\log(MN_{k}+1)\rceil.$ (101)

    2. (b)

A QROM is applied with output size

$b_{MN}=\aleph_{1}+n_{MN}+2n_{k}+2,$ (102)

with $\aleph_{1}$ being the number of bits used for the keep values (which
govern the precision of the state preparation via the inequality test). Here
$n_{MN}$ and $2n_{k}$ are for $\ell$ and $\mathbf{Q}$, with the factor of 2
accounting for ind and alt values of $\mathbf{Q}$. The extra 2 qubits are for
outputting a qubit showing if $\ell=0$ (for selecting between the one- and
two-body parts). The complexity is

$\left\lceil\frac{MN_{k}+1}{k_{MN}}\right\rceil+b_{MN}(k_{MN}-1).$ (103)

    3. (c)

An inequality test is performed with complexity $\aleph_{1}$.

    4. (d)

A controlled swap is performed with complexity $n_{k}+\lceil\log M\rceil+1$.

  2. 2.

Next, we prepare a state on the second register as

$\displaystyle\frac{1}{\sqrt{\lambda}}\left(\mathinner{|{0,0,0,0}\rangle}\sum_{p,q}\left[\sqrt{2{|{\rm
Re}(h^{\prime}_{pq})|}}\mathinner{|{\theta_{pq0}^{(0)}}\rangle}\mathinner{|{0,p,q,0}\rangle}+\sqrt{2{|{\rm
Im}(h^{\prime}_{pq})|}}\mathinner{|{\theta_{pq1}^{(0)}}\rangle}\mathinner{|{0,p,q,1}\rangle}\right]\right.$
$\displaystyle+\frac{1}{\sqrt{2}}\sum_{\mathbf{Q},n}\mathinner{|{\ell,\mathbf{Q},n,1}\rangle}\sqrt{\sum_{\mathbf{k},rs}(|{\rm
Re}[L_{r\mathbf{k}s(\mathbf{k}{\ominus}\mathbf{Q}),n}]|+|{\rm
Im}[L_{r\mathbf{k}s(\mathbf{k}{\ominus}\mathbf{Q}),n}]|)}$
$\displaystyle\times\left.\sum_{\mathbf{k},p,q}\left[\sqrt{|{\rm
Re}(L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n})|}\mathinner{|{\theta_{\mathbf{k}pq0}^{(\ell)}}\rangle}\mathinner{|{\mathbf{k},p,q,0}\rangle}+\sqrt{|{\rm
Im}(L_{p\mathbf{k}q(\mathbf{k}{\ominus}\mathbf{Q}),n})|}\mathinner{|{\theta_{\mathbf{k}pq1}^{(\ell)}}\rangle}\mathinner{|{\mathbf{k},p,q,1}\rangle}\right]\right)\mathinner{|{+}\rangle}\mathinner{|{+}\rangle},$
(104)

where $\theta_{\mathbf{k}pq0}^{(\ell)}$, $\theta_{\mathbf{k}pq1}^{(\ell)}$ are
used to obtain the correct signs on the terms, and the
$\mathinner{|{+}\rangle}$ states at the end are used to select the spin and
control the swap between the $p$ and $q$ registers.

Now we have a distinction from [32] in that we have separate real and
imaginary parts, and a separate prepared qubit to flag between the real and
imaginary parts. Because of the large number of variables, we will again use a
single variable for iteration, and use it to output $\mathbf{k},p,q$. The
complexity of this state preparation is then as follows.

    1. (a)

First, prepare an equal superposition over the variable for iteration. There
are $P=N_{k}N^{2}/2$ values to take, which includes a factor of $2$ for the
real and imaginary parts, $N_{k}$ for $\mathbf{k}$, and $N^{2}/4$ for the
values of $p,q$. Then the complexity of preparing the equal superposition is
$3n_{P}-3\eta+2b_{r}-9$, where $n_{P}=\lceil\log P\rceil$, with $\eta$ being
the largest number such that $2^{\eta}$ is a factor of $P$.

    2. (b)

The size of the QROM output is

$b_{p}=2n_{k}+4n_{N}+\aleph_{2}+3,$ (105)

where the first term is for the three components of $\mathbf{k}$, the second
is for $p$ and $q$. The third is for ind and alt values of the qubit to store
the correct sign, as well as an alt value of the extra qubit for selecting
between the real and imaginary parts. We do not include an ind value for that
qubit, because it is part of the register we are iterating over. The
complexity of this QROM will be

$\left\lceil\frac{MN_{k}+1}{k_{p1}}\right\rceil\left\lceil\frac{P}{k_{p2}}\right\rceil+b_{p}(k_{p1}k_{p2}-1),$
(106)

where we are accounting for the cost to select based on both the index from
the factorization and the index for $\mathbf{k},p,q$, and using the result for
the complexity of QROM on two registers from Appendix G of [32].

    3. (c)

Perform the inequality test with cost $\aleph_{2}$, which is the bits of
precision for this state preparation.

    4. (d)

Perform the controlled swap with the alt values with cost $n_{k}+2n_{N}+1$.
Here we are swapping the ind and alt values of $\mathbf{k},p,q$, as well as
the qubit selecting between real and imaginary parts. The sign required for
the sign qubits can be implemented with Cliffords as in [32], so does not add
to this Toffoli cost.

  3. 3.

We no longer perform swaps of $p$ and $q$ for symmetry, but we do need to
perform arithmetic to compute $\mathbf{k}{\ominus}\mathbf{Q}$ and
$\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$, with a cost of $4n_{k}$.

  4. 4.

Perform select by performing the sequence of four controlled
$\vec{Z}X_{p,\sigma}$ or $\vec{Z}Y_{p,\sigma}$ operations. The cost is
$4(NN_{k}/2-1)$ Toffolis since it must be controlled, and there is a cost of
one more Toffoli to create the qubits to control on. In order to select the
spin we also perform a swap controlled by the spin selection qubit before and
after, with a cost of $NN_{k}$ Toffolis.

  5. 5.

Reverse steps 2 and 3, where the complexities are the same except the QROM
complexity which is changed to

$\left\lceil\frac{MN_{k}+1}{k^{\prime}_{p1}}\right\rceil\left\lceil\frac{P}{k^{\prime}_{p2}}\right\rceil+k^{\prime}_{p1}k^{\prime}_{p2}.$
(107)

  6. 6.

Reflect on the qubits that were prepared in step 2. The qubits we need to
reflect on are as follows.

    1. (a)

The $n_{P}$ qubits for the variable of iteration.

    2. (b)

We need to reflect on the $\aleph_{2}$ registers that are used for the equal
superposition state for the state preparation.

    3. (c)

One that is rotated for the preparation of the equal superposition state.

    4. (d)

One for the spin.

    5. (e)

One for controlling the swap between the $p$ and $q$ registers.

    6. (f)

One for selecting between the real and imaginary part.

    7. (g)

One for selecting between $A$ and $B$.

That gives a total of $n_{P}+\aleph_{2}+5$ qubits. The reflection needs to be
controlled on the success of the preparation on the $\ell$ register, and
$\ell\neq 0$, making the total cost $n_{P}+\aleph_{2}+5$ Toffolis.

  7. 7.

Perform steps 2 to 5 again, but this time $MN_{k}+1$ is replaced with $MN_{k}$
in Eq. (106) and Eq. (107). Also, the select operation needs to be controlled
on $\ell\neq 0$, which flags the one-body term. That requires another 4
Toffolis.

  8. 8.

Invert the state preparation on the $\ell$ register, where the complexity of
the QROM is reduced to

$\left\lceil\frac{MN_{k}+1}{k^{\prime}_{P}}\right\rceil+k^{\prime}_{P}.$ (108)

  9. 9.

To complete the step of the quantum walk, perform a reflection on the ancillas
used for the state preparation. There are
$n_{MN}+n_{P}+\aleph_{1}+\aleph_{2}+5$, where the qubits we need to reflect on
are as follows.

    1. (a)

The $n_{MN}$ qubits for the $\ell$ register.

    2. (b)

The $n_{P}$ qubits for the registers in the state preparation for $A$ and $B$.

    3. (c)

The $\aleph_{1}$ qubits for the equal superposition state used for preparing
the state on the $\ell$ register using the coherent alias sampling.

    4. (d)

The $\aleph_{2}$ qubits for the equal superposition state for preparing the
state for $A$ and $B$.

    5. (e)

Two qubits rotated for the boosting the success probability for the equal
superposition states.

    6. (f)

One qubit for the spin.

    7. (g)

One qubit for controlling the swap of the $p$ and $q$ registers.

    8. (h)

One for selecting between the real and imaginary part.

    9. (i)

One for selecting between $A$ and $B$.

This reflection has cost $n_{MN}+n_{P}+\aleph_{1}+\aleph_{2}+4$.

  10. 10.

The steps of the walk are made controlled by using unary iteration on an
ancilla used for the phase estimation. Each step requires another two Toffolis
for the unary iteration and making the reflection controlled.

In this list of steps we have not explicitly included the part for applying
the phase factors, but that has no non-Clifford cost.

Next we consider the total number of logical qubits needed for the simulation
via this method.

  1. 1.

The control register for the phase estimation, and the ancillas for the unary
iteration, together need $2\lceil\log\mathcal{I}\rceil-1$ qubits.

  2. 2.

There are $NN_{k}$ qubits for the target system.

  3. 3.

There are $n_{MN}+2$ qubits for the $\ell$ register, the qubit rotated in
preparing the equal superposition, and the qubit flagging success of preparing
the equal superposition.

  4. 4.

The state preparation on the $\ell$ register uses $b_{MN}=2n_{k}+2\lceil\log
M\rceil+2\aleph_{1}+2$ qubits. Here $2n_{k}+2\lceil\log M\rceil$ is for the
ind and alt values of $\mathbf{Q}$ and $n$, $\aleph_{1}$ are for keep values,
$\aleph_{1}$ are for the equal superposition state, 1 is for the output of the
inequality test, and 2 are for the qubit flagging $\ell\neq 0$ and its
alternate value.

  5. 5.

There are $n_{P}+2$ qubits needed for the register preparing $p,q,\mathbf{k}$
values, a qubit that is rotated for the equal superposition, and a qubit
flagging success of preparing the equal superposition.

  6. 6.

The equal superposition state used for the second preparation uses
$\aleph_{2}$ qubits.

  7. 7.

The phase gradient register uses $b_{r}$ qubits.

  8. 8.

The qubits for the spin, controlling the swap of $p$ and $q$, selection
between the real and imaginary parts, and selection between $A$ and $B$ for a
total of 4.

  9. 9.

The QROM needs a number of qubits
$b_{p}k_{p1}k_{p2}+\lceil\log[(MN_{k}+1)/k_{p1}]\rceil+\lceil\log[L/k_{p2}]\rceil$.

The QROM for the state preparation on the second register uses a large number
of temporary ancillas, which can be reused by later parts of the algorithm, so
those later parts of the algorithm do not need the number of qubits counted.
The total number of qubits used is then

$2\lceil\log\mathcal{I}\rceil+NN_{k}+n_{MN}+n_{P}+2n_{k}+2\lceil\log
M\rceil+2\aleph_{1}+\aleph_{2}+b_{r}+9+b_{p}k_{p1}k_{p2}+\lceil\log[(MN_{k}+1)/k_{p1}]\rceil+\lceil\log[L/k_{p2}]\rceil$
(109)

with $b_{p}=2n_{k}+2n_{N}+\aleph_{2}+3$, $n_{N}=\lceil\log(N/2)\rceil$,
$n_{P}=\lceil\log P\rceil$, $L=N_{k}N(N+2)/4$. This completes the costing of
the low rank factorization method.

##  Appendix C Double-factorization derivations

###  C.1 One-body correction

Here we derive the correction for the one-body Hamiltonian as given in Eq.
(51). The lambda value for the Hamiltonian can be calculated by determining
the total L1-norm using the second factorization

$\hat{H}^{\prime}_{2}=\frac{1}{2}\sum_{\mathbf{Q}}^{N_{k}}\sum_{n}^{M}\left(\hat{A}^{2}_{n}(\mathbf{Q})+\hat{B}^{2}_{n}(\mathbf{Q})\right)$
(110)

with

$\displaystyle 2\hat{A}_{n}(\mathbf{Q})$
$\displaystyle=\sum_{\mathbf{k}}\left[U^{A}_{n}(\mathbf{Q},\mathbf{k})\left(\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}f^{A}_{p}(\mathbf{Q},n,\mathbf{k})(\mathbb{1}-Z_{p\mathbf{k}\sigma})\right)U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\right]$
$\displaystyle=\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{\mathbb{1}}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}-\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{Z}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}$
(111)

where
$\hat{\mathbb{1}}^{A}_{\mathbf{k}}=\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}f^{A}_{p}(\mathbf{Q},n,\mathbf{k})\mathbb{1}$
and
$\hat{Z}^{A}_{\mathbf{k}}=\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}f^{A}_{p}(\mathbf{Q},n,\mathbf{k})Z_{p\mathbf{k}\sigma}$,
and

$\displaystyle 2\hat{B}_{n}(\mathbf{Q})$
$\displaystyle=\sum_{\mathbf{k}}\left[U^{B}_{n}(\mathbf{Q},\mathbf{k})\left(\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},B}}f^{B}_{p}(\mathbf{Q},n,\mathbf{k})(\mathbb{1}-Z_{p\mathbf{k}\sigma})\right)U^{B}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\right]$
$\displaystyle=\sum_{\mathbf{k}}U^{B}_{n}(\mathbf{Q},\mathbf{k})\hat{\mathbb{1}}^{B}_{\mathbf{k}}U^{B}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}-\sum_{\mathbf{k}}U^{B}_{n}(\mathbf{Q},\mathbf{k})\hat{Z}^{B}_{\mathbf{k}}U^{B}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}$
(112)

where
$\hat{\mathbb{1}}^{B}_{\mathbf{k}}=\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},B}}f^{B}_{p}(\mathbf{Q},n,\mathbf{k})\mathbb{1}$
and
$\hat{Z}^{B}_{\mathbf{k}}=\sum_{\sigma}\sum_{p}^{\Xi_{\mathbf{Q},n,\mathbf{k},B}}f^{B}_{p}(\mathbf{Q},n,\mathbf{k})Z_{p\mathbf{k}\sigma}$.
The factor of 1/2 from the Jordan-Wigner transform is squared to 1/4, which is
moved outside each term and combined with the prefactor 1/2 to produce a
prefactor of 1/8. We note that $\hat{A}_{n}(\mathbf{Q})^{2}$ can be written as

$\displaystyle 4\hat{A}_{n}(\mathbf{Q})^{2}$
$\displaystyle=\left(\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{\mathbb{1}}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}-\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{Z}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\right)$
$\displaystyle\quad\times\left(\sum_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})\hat{\mathbb{1}}^{A}_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})^{\dagger}-\sum_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})\hat{Z}^{A}_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})^{\dagger}\right)$
$\displaystyle=2\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{\mathbb{1}}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\hat{A}_{n}(\mathbf{Q})+2\hat{A}_{n}(\mathbf{Q})\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{\mathbb{1}}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}$
$\displaystyle\quad+\sum_{\mathbf{k},\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{Z}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})\hat{Z}^{A}_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})^{\dagger}$
$\displaystyle\quad-\sum_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{\mathbb{1}}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\sum_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})\hat{\mathbb{1}}^{A}_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})^{\dagger}.$
(113)

The last term in the above equation is proportional to the identity and is
ignored. A similar expression can be derived for $\hat{B}_{n}(\mathbf{Q})^{2}$
and thus the component of the two-body term involving two Pauli $Z$ operators
is written as

$\displaystyle V$
$\displaystyle=\frac{1}{8}\sum_{\mathbf{Q},n,\mathbf{k},\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k})\hat{Z}^{A}_{\mathbf{k}}U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})\hat{Z}^{A}_{\mathbf{k}^{\prime}}U^{A}_{n}(\mathbf{Q},\mathbf{k}^{\prime})^{\dagger}$
$\displaystyle\quad+\frac{1}{8}\sum_{\mathbf{Q},n,\mathbf{k},\mathbf{k}^{\prime}}U^{B}_{n}(\mathbf{Q},\mathbf{k})\hat{Z}^{B}_{\mathbf{k}}U^{B}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}U^{B}_{n}(\mathbf{Q},\mathbf{k}^{\prime})\hat{Z}^{B}_{\mathbf{k}^{\prime}}U^{B}_{n}(\mathbf{Q},\mathbf{k}^{\prime})^{\dagger}$
(114)

which implies the two-body L1-norm, $\lambda_{\mathrm{DF},2}$, is

$\displaystyle\lambda_{\mathrm{DF},2}=\frac{1}{4}\sum_{\mathbf{Q},n}\left[\left(\sum_{\mathbf{k},p}^{N_{k}\Xi_{\mathbf{Q},n,\mathbf{k},A}}|f^{A}_{n}(p,\mathbf{Q},\mathbf{k})|\right)^{2}+\left(\sum_{\mathbf{k},p}^{N_{k}\Xi_{\mathbf{Q},n,\mathbf{k},B}}|f^{B}_{n}(p,\mathbf{Q},\mathbf{k})|\right)^{2}\right]$
(115)

where the factor of $1/8$ becomes a factor of $1/2$ accounting for spin. This
factor of $1/2$ is further divided by two because we perform oblivious
amplitude amplification–i.e. the inner step of qubitization evolving by
$2\hat{A}_{n}(\mathbf{Q})^{2}-\mathbb{1}$ and
$2\hat{B}_{n}(\mathbf{Q})^{2}-\mathbb{1}$.

Next, the one-body terms in the third line of Eq. (113) can be rewritten as

$2\sum_{\mathbf{k}}\hat{\mathbb{1}}^{A}_{\mathbf{k}}\hat{A}_{n}(\mathbf{Q})+2\hat{A}_{n}(\mathbf{Q})\sum_{\mathbf{k}}\hat{\mathbb{1}}^{A}_{\mathbf{k}}\,.$
(116)

This expression needs to be divided by 8 to give the contribution to the
Hamiltonian, and there is a similar contribution from
$\hat{B}_{n}(\mathbf{Q})^{2}$ to give the overall contribution to the one-body
Hamiltonian

$\displaystyle\frac{1}{2}\sum_{n,\mathbf{Q}}\left(\sum_{\mathbf{k}}\hat{\mathbb{1}}^{A}_{\mathbf{k}}\hat{A}_{n}(\mathbf{Q})+\sum_{\mathbf{k}}\hat{\mathbb{1}}^{B}_{\mathbf{k}}\hat{B}_{n}(\mathbf{Q})\right).$
(117)

Taking the trace of Eq. (C.1) and (C.1) then implies

$\displaystyle\sum_{\mathbf{k}}\hat{\mathbb{1}}^{A}_{\mathbf{k}}$
$\displaystyle=\mathbb{1}\operatorname{Tr}(\hat{A}_{n}(\mathbf{Q}))=\mathbb{1}\,\frac{1}{2}\left[\operatorname{Tr}(\hat{\rho}_{n}(\mathbf{Q}))+\operatorname{Tr}(\hat{\rho}^{\dagger}_{n}(\mathbf{Q}))\right],$
(118) $\displaystyle\sum_{\mathbf{k}}\hat{\mathbb{1}}^{B}_{\mathbf{k}}$
$\displaystyle=\mathbb{1}\operatorname{Tr}(\hat{B}_{n}(\mathbf{Q}))=\mathbb{1}\,\frac{i}{2}\left[\operatorname{Tr}(\hat{\rho}_{n}(\mathbf{Q}))-\operatorname{Tr}(\hat{\rho}^{\dagger}_{n}(\mathbf{Q}))\right].$
(119)

The trace of $\hat{\rho}_{n}(\mathbf{Q})$ is non-zero only for $\mathbf{Q}=0$.
In that case

$\operatorname{Tr}(\hat{\rho}_{n}(0))=2\sum_{\mathbf{k}}\left(\sum_{r}^{N/2}L_{r\mathbf{k}r\mathbf{k},n}\right)$
(120)

which is real. Moreover, it is easily seen that $\hat{\rho}(0)$ is Hermitian
using the symmetry
$L_{p\mathbf{k}q\mathbf{k},n}=L_{q\mathbf{k}p\mathbf{k},n}^{*}$, so
$\hat{A}_{n}(0)=\hat{\rho}_{n}(0)$ and $\hat{B}_{n}(0)=0$. Therefore

$\displaystyle\sum_{\mathbf{k}}\hat{\mathbb{1}}^{A}_{\mathbf{k}}\hat{A}_{n}(0)+\sum_{\mathbf{k}}\hat{\mathbb{1}}^{B}_{\mathbf{k}}\hat{B}_{n}(0)=2\sum_{\mathbf{k},p,q,\sigma}\left(\sum_{\mathbf{k}^{\prime},r}L_{p\mathbf{k}q\mathbf{k},n}L_{r\mathbf{k}^{\prime}r\mathbf{k}^{\prime},n}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}.$
(121)

Therefore the contribution to the one-body Hamiltonian becomes

$\displaystyle\frac{1}{2}\sum_{n,\mathbf{Q}}\left(\sum_{\mathbf{k}}\hat{\mathbb{1}}^{A}_{\mathbf{k}}\hat{A}_{n}(\mathbf{Q})+\sum_{\mathbf{k}}\hat{\mathbb{1}}^{B}_{\mathbf{k}}\hat{B}_{n}(\mathbf{Q})\right)=\sum_{\mathbf{k},p,q,\sigma}\left(\sum_{\mathbf{k}^{\prime},r}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}.$
(122)

As a result, the complete one-body Hamiltonian is

$\displaystyle
H_{1}^{\prime}=\sum_{\mathbf{k},p,q,\sigma}\left(h_{p\mathbf{k},q\mathbf{k}}+\sum_{\mathbf{k}^{\prime},r}V_{p\mathbf{k},q\mathbf{k},r\mathbf{k}^{\prime},r\mathbf{k}^{\prime}}\right)a_{p\mathbf{k}\sigma}^{\dagger}a_{q\mathbf{k}\sigma}.$
(123)

This is identical to the result that was obtained in the single-factorization
case as in Eq. (98). Thus the L1-norm of $H_{1}^{\prime}$ is the sum

$\displaystyle\lambda_{\mathrm{DF},1}=\sum_{\mathbf{k}}\sum_{p}|\lambda_{\mathbf{k},p}|$
(124)

where $\lambda_{\mathbf{k},p}$ is an eigenvalue of the matrix representing
$H_{1}^{\prime}(\mathbf{k})$ which are the coefficients in the parenthesis of
Eq. (123).

###  C.2 Complexity of the double-factorized representation

Our form of the two-body part of the Hamiltonian is

$\hat{H}^{\prime}_{2}=\frac{1}{2}\sum_{\mathbf{Q}}^{N_{k}}\sum_{n}^{M}\left(\hat{A}^{2}_{n}(\mathbf{Q})+\hat{B}^{2}_{n}(\mathbf{Q})\right).$
(125)

with

$\displaystyle\hat{A}_{n}(\mathbf{Q})=\sum_{\mathbf{k}}\left[U^{A}_{n}(\mathbf{Q},\mathbf{k})\left(\sum_{\sigma}\sum_{r}^{\Xi_{\mathbf{Q},n,\mathbf{k},A}}f^{A}_{r}(\mathbf{Q},n,\mathbf{k})n_{r,\mathbf{k},\sigma}\right)U^{A}_{n}(\mathbf{Q},\mathbf{k})^{\dagger}\right]$
(126)

and similarly for $\hat{B}_{n}(\mathbf{Q})$. In comparison, the double-
factorized Hamiltonian from [36, 32] is

$F^{\prime}=\frac{1}{8}\sum_{\ell=1}^{L}U_{\ell}\left(\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{p=1}^{\Xi^{(\ell)}}f_{p}^{(\ell)}Z_{p,\sigma}\right)^{2}U_{\ell}^{\dagger}.$
(127)

So, in contrast to the decomposition before, instead of a sum over $\ell$, we
have a sum over $\mathbf{Q},n$, and a qubit indexing over $\hat{A},\hat{B}$.
This difference can be accounted for easily in the method as presented in
[32]. That method may be summarized as follows.

  1. 1.

Perform a state preparation over $\ell$ for the first factorisation.

  2. 2.

Use a QROM on $\ell$ to output some parameters needed for the state
preparation for the second factorisation (the operator that is squared).

  3. 3.

Perform the inner state preparation over $p$.

  4. 4.

Apply a QROM to output the sequence of rotations dependent on $\ell$ and $p$.

  5. 5.

Apply the Givens rotations.

  6. 6.

Apply a controlled $Z$.

  7. 7.

Invert the Givens rotations, QROM, and state preparation over $p$.

  8. 8.

Perform a reflection on the ancilla qubits used for the state preparation over
$p$.

  9. 9.

Perform steps 3 to 7 again.

  10. 10.

Invert the QROM from step 2.

  11. 11.

Invert the state preparation from step 1.

Note that this is distinct from the procedure in [36] which combined the
$\ell$ and $p$ preparations.

To account for the changes here, the index $\ell$ can be used to iterate
through all possible values of $\mathbf{Q},n$, and the qubit indexing over
$\hat{A},\hat{B}$. Most of the steps can be performed ignoring these values,
but we will need to know $\mathbf{Q}$ before performing the Givens rotations.
It is convenient to output this value in the QROM used in step 2, which
slightly increases the output size of this QROM. We will also need to output
$\mathbf{k}$ values, and these will be given in the second state preparation
used in step 3. But, that preparation will produce a joint index of $p$ and
$\mathbf{k}$ without giving $\mathbf{k}$ explicitly (similar to our
preparation over $\ell$ not giving $\mathbf{Q}$ explicitly. This can be output
by the QROM in step 4.

In order to apply the Givens rotations, we will need to perform controlled
swaps of system registers $\mathbf{k},\mathbf{k}{\ominus}\mathbf{Q}$ into
working registers, then apply the Givens rotations on those working registers.
Since $\mathbf{k}{\ominus}\mathbf{Q}$ is not given directly by the state
preparation, it needs to be computed with cost $2n_{k}$ (or $n_{k}$ if
$N_{x},N_{y},N_{z}$ are powers of 2). The controlled swaps have a Toffoli cost
of $2n_{k}$ for the unary iteration, and $NN_{k}$ for the controlled swaps.
The cost of $NN_{k}$ is because we need to run through $NN_{k}/2$ system
qubits twice. These controlled swaps are performed 4 times, because they need
to be performed before and after each application of the Givens rotations.
That gives a total cost from this part

$4NN_{k}+12n_{k}.$ (128)

Then for the QROM outputting the Givens rotations, the number of items of data
can be given as

$\sum_{\mathbf{Q},n,\mathbf{k}}(\Xi_{\mathbf{Q},n,\mathbf{k},A}+\Xi_{\mathbf{Q},n,\mathbf{k},B}),$
(129)

where $\Xi_{\mathbf{Q},n,\mathbf{k},A}$ and $\Xi_{\mathbf{Q},n,\mathbf{k},B}$
are the cutoffs in the sums for $\hat{A}$ and $\hat{B}$. As per Eq. (47), we
define $\Xi$ as this quantity divided by $L=2N_{k}M$, so we can write the
number of items of data as $L\Xi$.

The Givens rotations need to be on $2N$ orbitals, so there are $2N$ Givens
rotations. For each of these rotations two angles need to be specified, in
contrast to one in [36, 32]. The size of the data output for the QROM for the
Givens rotations is increased to $2N\beth$, because there are 2 registers of
size $N/2$, and there are 2 rotations of $\beth$ of precision for each Givens
rotation. The total complexity of applying the Givens rotations is increased
to $16N(\beth-2)$. This is an increase of a factor of 4 over that in [32],
with a factor of 2 from using 2 working registers, and a factor of 2 because
there are two rotations for each Givens rotation.

The other changes in the cost are relatively trivial. There is a swap on the
system registers controlled on the spin register. Since this is now on
$NN_{k}$ qubits instead of $N$, the cost is multiplied by $N_{k}$.

So, to summarize the complexity using the same numbering of steps as in [32],
we have the following.

  1. 1.

The cost of the state preparation over $\ell$ is

$(3n_{L}-3\eta+2b_{r}-9)+\left\lceil\frac{L+1}{k_{p1}}\right\rceil+b_{p1}(k_{p1}-1)+\aleph_{1}+n_{L},$
(130)

where $L$ is now $2N_{k}M$ and as before $b_{p1}=n_{L}+\aleph_{1}$,
$n_{L}=\lceil\log L\rceil$.

  2. 2.

The complexity of the QROM on $\ell$ is now

$\left\lceil\frac{L+1}{k_{o}}\right\rceil+b_{o}(k_{o}-1),$ (131)

with

$b_{o}=n_{k}+n_{\Xi}+n_{L,\Xi}+b_{r}+1,$ (132)

with the extra $n_{k}$ being to output $\mathbf{Q}$. Here $n_{\Xi}$ is the
number of bits needed for $\Xi_{\mathbf{k},p}$ values of $p$, and $n_{L,\Xi}$

$n_{L,\Xi}=\lceil\log(L\Xi+N_{k}N/2)\rceil$ (133)

is the number of bits needed for the offset.

  3. 3.

The cost of the second stage of state preparation is

$\displaystyle
4(7n_{\Xi}+2b_{r}-6)+4(n_{L,\Xi}-1)+\left(\left\lceil\frac{L\Xi+NN_{k}/2}{k_{p2}}\right\rceil+\left\lceil\frac{L\Xi}{k_{p2}}\right\rceil+2b_{p2}(k_{p2}-1)\right)+4(\aleph_{2}+n_{\Xi}),$
(134)

where the brackets are used to indicate the cost of parts (a) to (d) of step
3. As well as using our modified definition of $\Xi$, the only change over the
costing in [32] is replacing $N/2$ with $NN_{k}/2$ for the range of values for
the one-body term. In this cost we are including the second use of the
preparation in part 7.

  4. 4.

The cost of the number operators via QROM is

$\left\lceil\frac{L\Xi+NN_{k}/2}{k_{r}}\right\rceil+\left\lceil\frac{L\Xi}{k_{r}}\right\rceil+(4N\beth+n_{k})(k_{r}-1)+\left\lceil\frac{L\Xi+NN_{k}/2}{k^{\prime}_{r}}\right\rceil+\left\lceil\frac{L\Xi}{k^{\prime}_{r}}\right\rceil+2k^{\prime}_{r}+4(n_{L,\Xi}-1)+16N(\beth-2)+2NN_{k}+2.$
(135)

Here the term $4N\beth(k_{r}-1)$ has been increased by a factor of 4 over that
in [32], because we have 2 times as many qubits that the Givens rotations need
to act on, and there are twice as many rotations needed for each Givens
rotation. (This term is corresponding to the output size for the QROM.) We
have also added $n_{k}$ for the output size so we can output the value of
$\mathbf{k}$ needed to select the register. Again $N/2$ is replaced with
$NN_{k}/2$ for the one-body term. The quantity $16N(\beth-2)$ is for the cost
of the Givens rotations, and is also multiplied by a factor of 4 over that in
[32]. The $2NN_{k}$ for the controlled swaps for spin, and is increased over
$2N$ in [32] because we now have $\mathbf{k}$.

  5. 5.

The inversion of the state preparation has cost

$\displaystyle
2(7n_{\Xi}+2b_{r}-6)+2(n_{L,\Xi}-1)+\left(\left\lceil\frac{L\Xi+NN_{k}/2}{k^{\prime}_{p2}}\right\rceil+\left\lceil\frac{L\Xi}{k^{\prime}_{p2}}\right\rceil+2k^{\prime}_{p2}\right)+2(\aleph_{2}+n_{\Xi}).$
(136)

This cost is the same as in part 3, except the cost of erasing the QROM is
reduced. We are again including both uses (with the second described in step
7).

  6. 6.

The reflection for the oblivious amplitude amplification has an unchanged cost

$n_{\Xi}+\aleph_{2}+2.$ (137)

  7. 7.

The cost of the second use of the block encoding to give the square are
already accounted for above.

  8. 8.

The cost of inverting step 1 is

$(3n_{L}-3\eta+2b_{r}-9)+\left\lceil\frac{L+1}{k^{\prime}_{p1}}\right\rceil+k^{\prime}_{p1}+\aleph_{1}+n_{L},$
(138)

and for inverting step 2 is

$\left\lceil\frac{L+1}{k^{\prime}_{o}}\right\rceil+k^{\prime}_{o},$ (139)

where we are using the improved cost for erasing QROM.

  9. 9.

The reflection cost is unchanged at

$n_{L}+n_{\Xi}+\aleph_{1}+\aleph_{2}+1.$ (140)

  10. 10.

The extra cost of unary iteration on the control register and of controlling
the reflection on that register is 2 Toffolis.

  11. 11.

The new costs of performing controlled swaps into working registers and
arithmetic to compute $\mathbf{k}{\ominus}\mathbf{Q}$ and
$\mathbf{k}^{\prime}{\ominus}\mathbf{Q}$ are

$4NN_{k}+12n_{k}.$ (141)

Adding all these costs together gives the total cost for block encoding the
Hamiltonian.

The cost in terms of logical qubits is very similar to that for the original
double-factorized approach. The differences are as follows.

  1. 1.

There are registers needed to store
$\mathbf{k},\mathbf{Q},\mathbf{k}{\ominus}\mathbf{Q}$. Because
$\mathbf{k}{\ominus}\mathbf{Q}$ can be computed in place in the $\mathbf{Q}$
register, we only need storage for 2. Moreover, because $\mathbf{k}$ is given
in the QROM output in part 4 above, it does not need to be added to that qubit
costing.

  2. 2.

There are $N$ qubits used for the working registers (2 of size $N/2$).

  3. 3.

A number of parameters are changed, in particular the number of system qubits
is now $NN_{k}$, and $L$ is computed from the number of values of $\mathbf{Q}$
and $n$.

  4. 4.

The size of the output for the Givens rotations is multiplied by a factor of
4.

##  Appendix D Tensor hypercontraction derivations

###  D.1 THC symmetries

In this section we derive the symmetry relationships for the central tensor
based on the four-fold symmetry of the two-electron integral tensor as used in
Eq. (62). Recall an element of the two-electron integral tensor can be
represented in THC form as

$\displaystyle
V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r\mathbf{k}^{\prime}{\ominus}\mathbf{Q},s\mathbf{k}^{\prime}}=\sum_{\mu,\nu}\chi_{p\mathbf{k},\mu}^{*}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\mu}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\nu}^{*}\chi_{s\mathbf{k}^{\prime},\nu}$
(142)

where $\mathbf{G}_{1}$ is shorthand for
$\mathbf{G}_{\mathbf{k},\mathbf{k}-\mathbf{Q}}$ and $\mathbf{G}_{2}$ is
shorthand for $G_{\mathbf{k}^{\prime},\mathbf{k}^{\prime}-\mathbf{Q}}$. The
four fold symmetry of the complex valued two-electron integral tensor is
reflected in the central tensor $\zeta$. We recover the symmetry by first
noting the four equivalent two-electron integrals

$\displaystyle
V_{p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q}),r\mathbf{k}^{\prime}-\mathbf{Q},s\mathbf{k}^{\prime}}=$
$\displaystyle\sum_{\mu,\nu}\chi_{p\mathbf{k},\mu}^{*}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\mu}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\nu}^{*}\chi_{s\mathbf{k}^{\prime},\nu}$
$\displaystyle
V_{q(\mathbf{k}{\ominus}\mathbf{Q}),p\mathbf{k},s\mathbf{k}^{\prime},r\mathbf{k}^{\prime}-\mathbf{Q}}^{*}=$
$\displaystyle\left(\sum_{\mu,\nu}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\mu}^{*}\chi_{p\mathbf{k},\mu}\zeta_{\mu\nu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{1},!\mathbf{G}_{2}}\chi_{s\mathbf{k}^{\prime},\nu}^{*}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\nu}\right)^{*}$
$\displaystyle
V_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),s\mathbf{k}^{\prime},p\mathbf{k},q(\mathbf{k}{\ominus}\mathbf{Q})}=$
$\displaystyle\sum_{\mu,\nu}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\mu}^{*}\chi_{s\mathbf{k}^{\prime},\mu}\zeta_{\mu\nu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{2},!\mathbf{G}_{1}}\chi_{p\mathbf{k},\nu}^{*}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\nu}$
$\displaystyle\
V_{s\mathbf{k}^{\prime},r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),q(\mathbf{k}{\ominus}\mathbf{Q}),p\mathbf{k}}=$
$\displaystyle\left(\sum_{\mu,\nu}\chi_{s\mathbf{k}^{\prime},\mu}^{*}\chi_{r(\mathbf{k}^{\prime}{\ominus}\mathbf{Q}),\mu}\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{2},\mathbf{G}_{1}}\chi_{q(\mathbf{k}{\ominus}\mathbf{Q}),\nu}^{*}\chi_{p\mathbf{k},\nu}\right)^{*}$

which, implies

$\displaystyle\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}=\left(\zeta_{\mu\nu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{1},!\mathbf{G}_{2}}\right)^{*}=\zeta_{\nu\mu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{2},!\mathbf{G}_{1}}=\left(\zeta_{\nu\mu}^{\mathbf{Q},\mathbf{G}_{2},\mathbf{G}_{1}}\right)^{*}.$
(143)

Here $({\ominus}\mathbf{Q})$ is used to indicate a modular negative of
$\mathbf{Q}$, similar to modular subtraction. In the above expression the
complement of $\mathbf{G}_{1}$, $!\mathbf{G}_{1}$, is defined through

$\displaystyle\mathbf{k}_{p}-\mathbf{k}_{q}=\mathbf{Q}+\mathbf{G}_{1}$
$\displaystyle\mathbf{k}_{q}-\mathbf{k}_{p}=({\ominus}\mathbf{Q})+!\mathbf{G}_{1}$
$\displaystyle!\mathbf{G}_{1}=-\left(\mathbf{Q}+\mathbf{G}_{1}+({\ominus}\mathbf{Q})\right),$
(144)

and it is important to note that $({\ominus}\mathbf{Q})$ is defined to be in
the original set of $k$-points and it is useful as we only build
$\zeta^{\mathbf{Q}}$. A similar expression can be derived for
$!\mathbf{G}_{2}$. It is helpful to consider some concrete examples, which are
given in Table 7.

$k$-mesh | $\mathbf{k}_{p}$ | $\mathbf{k}_{q}$ | $\mathbf{k}_{p}-\mathbf{k}_{q}$ | $\mathbf{Q}$ | $\mathbf{G}$ | $\mathbf{k}_{q}-\mathbf{k}_{p}$ | $({\ominus}\mathbf{Q})$ | $!\mathbf{G}$
---|---|---|---|---|---|---|---|---
$[1,1,4]$ | (0, 0, 3) | (0, 0, 1) | (0, 0, 2) | (0, 0, 2) | (0,0,0) | $(0,0,-2)$ | (0, 0, 2) | $(0,0,-4)$
$[1,4,4]$ | (0, 2, 1) | (0, 3, 1) | $(0,-1,0)$ | (0, 3, 0) | $(0,-4,0)$ | $(0,1,0)$ | (0, 1, 0) | $(0,0,0)$
$[1,4,4]$ | (0, 2, 1) | (0, 3, 3) | $(0,-1,-2)$ | (0, 3, 2) | $(0,-4,-4)$ | $(0,1,2)$ | (0, 1, 2) | $(0,0,0)$
$[1,4,4]$ | (0, 1, 2) | (0, 1, 3) | $(0,0,-1)$ | (0, 0, 3) | $(0,0,-4)$ | $(0,0,1)$ | (0, 0, 2) | $(0,0,0)$
$[1,4,4]$ | (0, 1, 3) | (0, 1, 2) | (0, 0, 1) | (0, 0, 1) | (0,0,0) | $(0,0,-1)$ | (0, 0, 3) | $(0,0,-4)$
$[4,4,4]$ | (2, 1, 3) | (3, 1, 2) | $(-1,0,1)$ | (3, 0, 1) | $(-4,0,0)$ | $(1,0,-1)$ | (0, 0, 3) | $(0,0,-4)$
$[4,4,4]$ | (2, 1, 2) | (3, 3, 3) | $(-1,-2,-1)$ | (3, 2, 3) | $(-4,-4,-4)$ | $(1,2,1)$ | (1, 2, 1) | $(0,0,0)$
Table 7: Some examples of the values that the different momentum labels can
take in Eq. 144. We restrict $\mathbf{k},\mathbf{Q},({\ominus}\mathbf{Q})$ to
be in the original $k$-point set.

###  D.2 Complexity of the tensor hypercontraction representation

The following is a detailed costing for the qubitization oracles using the THC
LCU. In the initial state preparation, we need to prepare a superposition over
$\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2},\mu,\nu$ with weights
$\sqrt{|\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}|}$. The
state can be prepared via the coherent alias sampling procedure, starting with
QROM to output keep and alt values. One option here is to produce an equal
superposition over $\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2},\mu,\nu$, then
calculate a contiguous register from these values to use for the QROM. That
procedure is fairly complicated, because it requires preparing equal
superpositions over three components of $\mathbf{Q}$ as well as
$\mathbf{G}_{1},\mathbf{G}_{2},\mu$ and $\nu$, then arithmetic for the
contiguous register. To simplify the procedure we give the complexity for
giving ind values like for sparse state preparation. That is, we prepare the
contiguous register, and use the QROM to output both ind (index) and alt
values of $\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2},\mu,\nu$.

There is the symmetry
$\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}=(\zeta_{\nu,\mu}^{\mathbf{Q},\mathbf{G}_{2},\mathbf{G}_{1}})^{*}$,
which indicates only half the range of $\mu,\nu,\mathbf{G}_{1},\mathbf{G}_{2}$
values need be prepared. It is convenient to prepare the full range, but use
part of the range for real and part for imaginary components. If we only were
considering $\mu,\nu$, we could use $\mu\leq\nu$ for real components and
$\mu>\nu$ for imaginary components. To account for
$\mathbf{G}_{1},\mathbf{G}_{2}$ as well, we can combine them with $\mu,\nu$ as
least-significant bits for combined integers to use in inequality tests. This
inequality test between $\mu,\mathbf{G}_{1}$ and $\nu,\mathbf{G}_{2}$ is used
to give a qubit flagging that the component should be imaginary. A further
qubit in a $\mathinner{|{+}\rangle}$ state is used to control a swap of
$\mu,\mathbf{G}_{1}$ with $\nu,\mathbf{G}_{2}$ registers, and a controlled $Z$
gate on the qubit flagging the imaginary component gives the desired complex
conjugate. As a result the range for $\mu,\nu$ is $M^{2}$ taking account of
giving real and imaginary components.

For $\mathbf{Q}$ there is also the symmetry where
$\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}=(\zeta_{\mu\nu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{1},!\mathbf{G}_{2}})^{*}$,
so it is only necessary to produce approximately half as many values of
$\mathbf{Q}$. This is complicated by the cases where
$\mathbf{Q}={\ominus}\mathbf{Q}$. If $N_{x},N_{y},N_{z}$ are odd, then the
only case where this can be true is that $\mathbf{Q}=0$, so the number of
values of $\mathbf{Q}$ that need be considered is $(N_{x}N_{y}N_{z}+1)/2$. If
one of $N_{x},N_{y},N_{z}$ is even and the other two are odd, then for the one
that is even there will be a second values of that component of $\mathbf{Q}$
that is equal to its negative. That means there are two value of $\mathbf{Q}$
overall satisfying $\mathbf{Q}={\ominus}\mathbf{Q}$, and the number of unique
values is $N_{x}N_{y}N_{z}/2+1$. Similarly, if there are two even values of
$N_{x},N_{y},N_{z}$, then there are four values of $\mathbf{Q}$ satisfying
$\mathbf{Q}={\ominus}\mathbf{Q}$, and the number of unique values is
$N_{x}N_{y}N_{z}/2+2$. For all three of $N_{x},N_{y},N_{z}$ even the number of
unique values is $N_{x}N_{y}N_{z}/2+4$. We also need $NN_{k}/2$ values for the
one-body term. The number of values is then

$d=32[N_{x}N_{y}N_{z}+2^{v}]M^{2}+NN_{k}/2,$ (145)

where $v$ is the number of even values of $N_{x},N_{y},N_{z}$.

The size of the output is then

$m=2(2\lceil\log M\rceil+n_{k}+8)+\aleph.$ (146)

where $\aleph$ is the number of bits for the keep register. There is a factor
of 2 at the front to account for ind and alt values, then $\lceil\log M\rceil$
for each of $\mu$ and $\nu$, and $n_{k}$ for the components of $\mathbf{Q}$.
There is a further qubit distinguishing between the one and two-electron
terms, a qubit giving the sign of the real or imaginary component of $\zeta$,
and 6 qubits for $\mathbf{G}_{1},\mathbf{G}_{2}$, for a total $+8$.

  1. 1.

There is a cost of $N_{k}N/2$ for controlled swaps for the spin. In principle
this is performed four times, because it is performed before and after the two
$c^{\dagger}c$ operators. The middle pair can be combined, with the single
controlled swap being controlled by the parity of the two spin qubits, for a
total cost of $3N_{k}N/2$.

  2. 2.

Before the state preparation, we need to prepare an equal superposition over
$d$ basis states, with costing $3\lceil\log d\rceil-3\eta+2b_{r}-9$ Toffoli
gates. As before, $\eta$ is a number such that $2^{\eta}$ is a factor of $d$,
and $b_{r}$ is a number of bits used for rotation of an ancilla qubit to
improve the amplitude of success. This cost is incurred twice, once for the
preparation and once for the inverse preparation.

  3. 3.

The complexity of the QROM being used for the state preparation is

$\left\lceil\frac{d}{k_{p}}\right\rceil+m(k_{p}-1),$ (147)

with $k_{p}$ being a power of 2. The inverse preparation then has a cost

$\left\lceil\frac{d}{k^{\prime}_{p}}\right\rceil+k^{\prime}_{p}.$ (148)

  4. 4.

We perform an inequality test with cost $\aleph$. Accounting for the inverse
of the preparation gives a total cost $2\aleph$.

  5. 5.

The controlled swap based on the result of the inequality test is on

$2\lceil\log M\rceil+n_{k}+7$ (149)

pairs of qubits, so has this Toffoli cost. Note that we have $+7$ here rather
than $+8$. This is because we do not need to swap the sign qubits; the sign
can be applied with $Z$ gates controlled on the result of the inequality test,
not adding to the Toffoli cost (as usual). This cost is incurred again in the
inverse preparation for a total of $4\lceil\log M\rceil+2n_{k}+14$.

  6. 6.

As described above, we perform an inequality test between $\mu,\mathbf{G}_{1}$
and $\nu,\mathbf{G}_{2}$ to give the qubit flagging whether we have a real or
imaginary component. Then we perform a controlled swap of $\mu,\mathbf{G}_{1}$
with $\nu,\mathbf{G}_{2}$ to generate one symmetry for the state preparation,
with the complex conjugate applied using a Clifford gate. This part therefore
has Toffoli cost $2\lceil\log M\rceil+6$. This cost is incurred again in the
inverse preparation giving a total cost $4\lceil\log M\rceil+12$. In addition
to this controlled swap, we perform a controlled swap in the middle, but it is
not controlled so does not add to the Toffoli complexity.

  7. 7.

For the symmetry where
$\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}=(\zeta_{\mu\nu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{1},!\mathbf{G}_{2}})^{*}$,
we can use a second control qubit to flip the sign on $\mathbf{Q}$, negate
$\mathbf{G}_{1}$ and $\mathbf{G}_{2}$, and apply the complex conjugate. The
complex conjugate again can be applied with a Clifford gate, and so can the
controlled-NOT gates on $\mathbf{G}_{1},\mathbf{G}_{2}$. A controlled sign
flip of $\mathbf{Q}$ can be performed with $2n_{k}$ Toffolis, simply by
flipping the sign in usual two’s complement binary, then controlling addition
of $N_{x},N_{y},N_{z}$ in each component.

  8. 8.

Next we need to prepare a superposition over allowed values of $\mathbf{k}$,
because $\mathbf{k}-\mathbf{Q}-\mathbf{G}$ needs to be in the allowed range of
$\mathbf{k}$ values (using $\mathbf{G}$ to indicate $\mathbf{G}_{1}$ or
$\mathbf{G}_{2}$ depending on which part we are performing). In particular,
for the $x$-component we have an allowed range for $k_{x}$ from $Q_{x}$ to
$N_{x}-1$ when $G_{x}=0$, or $0$ to $Q_{x}-1$ when $G_{x}\neq 0$. It is
similar for the other two components. We can therefore prepare a superposition
over the appropriate range then add $Q_{x}$ if $G_{x}=0$.

Creating an equal superposition requires Hadamards on the appropriate subset
of qubits, as well as a $Q_{x},G_{x}$-dependent rotation to give a high
success probability for the amplitude amplification. This information can be
output with Toffoli cost $2N_{x}-2$ on the qubits representing $Q_{x},G_{x}$.
The complexity of the controlled Hadamards is then $\lceil\log N_{x}\rceil$
Toffolis, assuming we use a catalytic T state as in [32]. The complexity of
preparing the equal superposition is then $6\lceil\log N_{x}\rceil+2b_{r}-6$,
including $3\lceil\log N_{x}\rceil$ for three rounds of $\lceil\log
N_{x}\rceil$ controlled Hadamards. The reason why there is $-6$ rather than
$-9$ is the inequality test is with a value in a quantum register (in each of
three tests), which requires one more Toffoli than an inequality test with a
classically given value.

The controlled addition of $Q_{x}$ has complexity $2\lceil\log N_{x}\rceil$.
The total complexity of the preparation of the superposition for the three
components of $\mathbf{k}$ is therefore

$N_{x}+N_{y}+N_{z}+8n_{k}+6b_{r}-24.$ (150)

This cost is incurred four times for the preparation and inverse preparation
of $\mathbf{k}$ and $\mathbf{k}^{\prime}$.

  9. 9.

In order to account for the one-body term, we note than the one-body term has
a single $\mu$ and $\mathbf{k}$ rather than $\mathbf{Q}$. We also do not want
the operations we perform in the two-body part for the symmetry to affect the
one-body part. We can therefore output $\mu=\nu$ for the one-body part in the
QROM, so the swap of $\mu$ and $\nu$ has no effect. The value of $\mathbf{k}$
for the one-body part can be stored in the same register as used for
$\mathbf{Q}$ for the two-body part. To prevent the operations used to generate
the symmetry
$\zeta_{\mu\nu}^{\mathbf{Q},\mathbf{G}_{1},\mathbf{G}_{2}}=(\zeta_{\mu\nu}^{({\ominus}\mathbf{Q}),!\mathbf{G}_{1},!\mathbf{G}_{2}})^{*}$
being applied for the one-body part, we can simply apply a Toffoli to produce
a new control qubit. The remaining part above is the preparation of the
superposition over the $\mathbf{k}$ values controlled on $\mathbf{Q}$; this
does not need to be amended to account for the one-body part because there we
will not be using this value in the extra register.

  10. 10.

Now that we have prepared the register that is in an equal superposition over
the appropriate range of $\mathbf{k}$, we need to use that in combination with
$\mathbf{Q}$ and $\mu$ to prepare a superposition with the correct weights. To
do this, we will use coherent alias sampling in the usual way, but will need
to construct an appropriate register to iterate over from registers
$\mathbf{k},\mathbf{Q},\mathbf{G},\mu$. First we compute
$\mathbf{k}-\mathbf{Q}-\mathbf{G}$ in an ancilla register. These two
subtractions have cost $2n_{k}$. Since it needs to be computed and uncomputed
for each of the two factors in the Hamiltonian, the total cost is $8n_{k}$.

Now, because $\mathbf{k}$ and $\mathbf{k}{\ominus}\mathbf{Q}$ uniquely specify
$\mathbf{Q},\mathbf{G}$, these two registers can be used for the iteration
instead of $\mathbf{Q}$, with the additional advantage that they are both over
the full range of the Brillouin zone. Now we need to compute a contiguous
register

$(((((\mathbf{k}_{x}N_{y}+\mathbf{k}_{y})N_{z}+\mathbf{k}_{z})N_{x}+\mathbf{k}^{\prime}_{x})N_{y}+\mathbf{k}^{\prime}_{y})N_{z}+\mathbf{k}^{\prime}_{z})M+\mu,$
(151)

where we are using $\mathbf{k}^{\prime}$ for $\mathbf{k}{\ominus}\mathbf{Q}$.
This contiguous register includes many multiplications by classically chosen
constants, which has complexity depending on how many ones are in these
constants. The worst case is where these numbers are all ones, so we will give
the cost for that case even though it is rare.

As discussed in [117] the cost of multiplying two integers when one is given
classically is no more than the product of the numbers of bits. For the
additions, the cost is no more than the number of bits on the larger number.
So, we have a cost as follows.

    1. (a)

For multiplying $\mathbf{k}_{x}N_{y}$ a cost of $\lceil\log
N_{x}\rceil\lceil\log N_{y}\rceil$. Here $N_{y}$ would have more bits if it
were a power of 2, but then the multiplication cost would be zero.

    2. (b)

For adding $+\mathbf{k}_{y}$ a cost of $\lceil\log N_{x}N_{y}\rceil$.

    3. (c)

For multiplying $\times N_{z}$ the cost is $\lceil\log
N_{x}N_{y}\rceil\lceil\log N_{z}\rceil$.

    4. (d)

For adding $+\mathbf{k}_{z}$ the cost is $\lceil\log N_{k}\rceil$.

    5. (e)

For multiplying $\times N_{x}$ the cost is $\lceil\log N_{k}\rceil\lceil\log
N_{x}\rceil$.

    6. (f)

For adding $+\mathbf{k}^{\prime}_{x}$ the cost is $\lceil\log
N_{x}N_{k}\rceil$.

    7. (g)

For multiplying $\times N_{y}$ the cost is $\lceil\log
N_{x}N_{k}\rceil\lceil\log N_{y}\rceil$.

    8. (h)

For adding $+\mathbf{k}^{\prime}_{y}$ the cost is $\lceil\log
N_{x}N_{y}N_{k}\rceil$.

    9. (i)

For multiplying $\times N_{z}$ the cost is $\lceil\log
N_{x}N_{y}N_{k}\rceil\lceil\log N_{z}\rceil$.

    10. (j)

For adding $+\mathbf{k}^{\prime}_{z}$ the cost is $\lceil\log
N_{k}^{2}\rceil$.

    11. (k)

For multiplying $\times M$ the cost is $\lceil\log N_{k}^{2}\rceil\lceil\log
M\rceil$.

    12. (l)

For finally adding $+\mu$ the cost is $\lceil\log N_{k}^{2}M\rceil$.

We need to add all these items together to give the total cost, and it needs
to be multiplied by 4 because we compute and uncompute for each of the two
factors in the Hamiltonian.

Next we have a QROM on this contiguous register with cost

$\left\lceil\frac{N_{k}^{2}M}{k_{\rm nrm}}\right\rceil+(k_{\rm
nrm}-1)(n_{k}+\aleph).$ (152)

with $k_{\rm nrm}$ a power of 2. This is because there are $N_{k}^{2}M$ items
to iterate over, and we need to output $n_{k}$ bits for the alternate value of
$k$ and $\aleph$ for the keep value. We have twice this cost because of the
two factors in the Hamiltonian, but the erasure cost for each factor is

$\left\lceil\frac{N_{k}^{2}M}{k_{\rm era}}\right\rceil+k_{\rm era}.$ (153)

The last two steps of the coherent alias sampling are an inequality test with
cost $\aleph$ and a controlled swap with cost $n_{k}$. These costs are
incurred 4 times, once for preparation and once for inverse preparation for
each of the two factors for the Hamiltonian.

  11. 11.

We will need to prepare a register that is $\mathbf{k}-\mathbf{Q}-\mathbf{G}$
again. We previously computed this, but we need to compute it again because we
have performed a state preparation on $\mathbf{k}$. This has a cost of
$2n_{k}$ again, and needs to be done 4 times for a total cost of $8n_{k}$.

A further subtlety is that we are storing the value of $\mathbf{k}$ to use in
the $\mathbf{Q}$ register in the one-body case. We can perform a controlled
swap into the working register for $\mathbf{k}$ or
$\mathbf{k}{\ominus}\mathbf{Q}$, which has a total cost of $4n_{k}$. Combined
with the arithmetic cost this is $12n_{k}$.

  12. 12.

To use the register with $\mathbf{k}$ or $\mathbf{k}{\ominus}\mathbf{Q}$ to
control the swap of system registers into working registers, we can use each
qubit to control swaps of the system registers in a similar way as is used for
advanced QROM. The cost for selecting each qubit out of $N_{k}$ is $N_{k}-1$,
similar to the use in advanced QROM, despite the use of multiple components.
In particular, we can perform swaps of system registers based on the
$x$-component of $\mathbf{k}$ with cost $(N_{x}-1)N_{y}N_{z}$. Then swapping
the registers based on the $y$-component out of the subset of $N_{y}N_{z}$ has
cost $(N_{y}-1)N_{z}$. Then the cost of swapping based on the $z$-component
has cost $N_{z}-1$. Adding these three costs together gives
$N_{x}N_{y}N_{z}-1=N_{k}-1$. This is performed for each of the $N/2$ qubits we
need, for cost $N(N_{k}-1)/2$. We need to swap and inverse swap 8 times on
$NN_{k}/2$ system registers, for a total cost of $4N(N_{k}-1)$.

  13. 13.

Next, we consider the output of the rotations for the $c$ modes. These will be
controlled by the registers with $\mu$ and $\mathbf{k}$ (or
$\mathbf{k}{\ominus}\mathbf{Q}$), as well as the qubit selecting between the
one- and two-body terms. A difficulty is that we would need a contiguous
register in order to be able to effectively apply the advanced QROM. A method
around this is to use a QROM on the selection qubit and $\mathbf{k}$ to output
an offset, then add $\mu$. The complexity of that QROM is $2N_{k}$, then the
complexity of the addition is $\lceil\log N_{k}(M+N/2)\rceil$.

That is in the case where we need to apply the one-body term as part of the
implementation. Recall that we only need to do that once when we are applying
$c^{\dagger}c$ twice for the two-body term. In the part where we are not
applying the one-body term we instead have complexity $N_{k}+\lceil\log
N_{k}M\rceil$.

The size of the QROM output for the rotations is then $N\beth$. That is again
because we need Givens rotations on $N/2$ qubits, and need two angles for each
Givens rotation with $\beth$ each. The complexity is then, in the case with
the one-body term,

$\left\lceil\frac{N_{k}(M+N/2)}{k_{r}}\right\rceil+N\beth(k_{r}-1),$ (154)

and for the case where we do not need the one-body term

$\left\lceil\frac{N_{k}M}{k_{r}}\right\rceil+N\beth(k_{r}-1),$ (155)

with $k_{r}$ being a power of 2. The cost of the Givens rotations is
$2N(\beth-2)$, because we have $N/2$ qubits and 2 angles for each Givens
rotation.

The cost of erasing the QROM is for the two cases is

$\displaystyle\left\lceil\frac{N_{k}(M+N/2)}{k^{\prime}_{r}}\right\rceil+k^{\prime}_{r},$
(156)
$\displaystyle\left\lceil\frac{N_{k}M}{k^{\prime}_{r}}\right\rceil+k^{\prime}_{r}.$
(157)

Lastly we note that we need to apply the sequence of Givens rotations 8 times
to account for the four $c^{\dagger}$ and $c$ operators, and similarly we need
to apply the QROM and invert it four times. That gives a total complexity for
the QROM-based basis rotations

$\displaystyle
2\left\lceil\frac{N_{k}(M+N/2)}{k_{r}}\right\rceil+2N\beth(k_{r}-1)+2\left\lceil\frac{N_{k}M}{k_{r}}\right\rceil+2N\beth(k_{r}-1)+2\left\lceil\frac{N_{k}(M+N/2)}{k^{\prime}_{r}}\right\rceil+2k^{\prime}_{r}$
$\displaystyle+2\left\lceil\frac{N_{k}M)}{k^{\prime}_{r}}\right\rceil+2k^{\prime}_{r}+16N(\beth-2)+12N_{k}+4\lceil\log
N_{k}(M+N/2)\rceil+4\lceil\log N_{k}M\rceil,$ (158)

where $16N(\beth-2)$ is for the Givens rotations themselves, and the terms
from $12N_{k}$ on are for creating and erasing contiguous registers.

  14. 14.

The next part of the complexity that needs to be accounted for is the
selection of $\vec{Z}X$ and $\vec{Z}Y$ for the implementation of the
$c^{\dagger}$ and $c$ operators. We only need to select between $\vec{Z}X$ and
$\vec{Z}Y$, but not select the location the $X$ or $Y$ is performed since
these are applied in the working registers. This selection can therefore be
performed entirely with Clifford gates. However, we do need additional control
to avoid performing these operators for one of the $c^{\dagger}c$ for the one-
body component. That is a cost of just one Toffoli for each of $c^{\dagger}$
and $c$.

A complication is that we need to perform $Z$ gates on remaining system
registers (that have not been swapped into the working registers). We perform
unary iteration on the register containing $\mathbf{k}$ (or
$\mathbf{k}{\ominus}\mathbf{Q}$), and use that to apply the appropriate $Z$
operators with a Toffoli cost $N_{k}-1$. The Toffoli complexity is independent
of $N$, because we only have a Toffoli cost for the unary iteration, not for
the controlled-$Z$ gates.

  15. 15.

The last part to the complexity when constructing the qubitised operator is
the reflection on the control ancillas. The qubits we need to reflect on are
as follows.

    1. (a)

The $\lceil\log d\rceil$ qubits used for the equal superposition state for the
initial preparation, and another qubit rotated for this preparation.

    2. (b)

The $\aleph$ used for the equal superposition state for the inequality test in
state preparation.

    3. (c)

The two qubits for the two spins in the sum.

    4. (d)

There are $n_{k}$ qubits that an equal superposition of $\mathbf{k}$ values is
prepared in, and this is done twice for $2n_{k}$ qubits. To save on qubit use
we can also reuse these qubits and flag that they are equal to zero in between
preparing $\mathbf{k}$ and $\mathbf{k}^{\prime}$, but that does not affect the
Toffoli count.

    5. (e)

There are $\aleph$ used for the equal superposition state for two rounds of
state preparation for $\mathbf{k}$. Again these can be reused but that does
not affect the qubit count.

    6. (f)

There are 4 qubits to select between $X$ and $Y$ for each of the
$c^{\dagger},c$.

    7. (g)

There are two qubits used to generate the symmetries.

This is a total of

$\lceil\log d\rceil+3\aleph+2n_{k}+9$ (159)

qubits for the control, and the reflection cost is 2 less than this. We do
require an additional Toffoli for the control by the phase estimation
registers, and another for unary iteration on the phase estimation registers.
Therefore this expression can be used for the additional Toffoli cost.

Next we account for the qubit costs.

  1. 1.

There are $NN_{k}$ system qubits.

  2. 2.

The control register for the phase estimation, and the ancillas for the unary
iteration, together need $2\lceil\log\mathcal{I}\rceil-1$ qubits.

  3. 3.

There are all the qubits needed for controls as described in the last item in
the Toffoli costing above. Taking the inversion of the equal superposition
over $\mathbf{k}$ to be flagged by a single qubit, the $2n_{k}$ can be
replaced with $n_{k}+1$. Similarly flag qubits can be used on the qubits used
for the equal superposition state for the inequality test to replace $3\aleph$
with $\aleph+2$ qubits. That gives a number of qubits

$\lceil\log d\rceil+\aleph+n_{k}+12.$ (160)

  4. 4.

A phase gradient state of size $\beth$ is used for rotations.

  5. 5.

A single T state is used for the controlled Hadamards.

  6. 6.

The QROM used for the first state preparation outputs $m$ qubits.

  7. 7.

This first state preparation also uses
$m(k_{p}-1)+\lceil\log(d/k_{p})\rceil-1$ temporary qubits.

  8. 8.

A single qubit is used for the result of the inequality test for the coherent
alias sampling.

  9. 9.

A single qubit is used for the result of the inequality test between
$\mu,\mathbf{G}_{1}$ and $\nu,\mathbf{G}_{2}$.

  10. 10.

There are $n_{k}+3b_{r}$ qubits used as the output of the QROM used to give
the information needed for the preparation of the equal superposition over
$\mathbf{k}$, as well as $b_{r}$ for $\mathbf{k}$.

  11. 11.

In the preparation of the state for $\mathbf{k}$ we compute
$\mathbf{k}{\ominus}\mathbf{Q}$ in an ancilla needing $n_{k}$ qubits, and
compute a contiguous register that needs $\lceil\log(N_{k}^{2}M)\rceil$
qubits.

  12. 12.

The state preparation for $\mathbf{k}$ uses $n_{k}+\beth$ output qubits.

  13. 13.

The state preparation uses

$(k_{\rm nrm}-1)(n_{k}+\beth)+\left\lceil\log\left(\frac{N_{k}^{2}M}{k_{\rm
nrm}}\right)\right\rceil-1$ (161)

temporary qubits.

  14. 14.

We also use $\aleph$ in this state preparation for a superposition state, and
another qubit for the result of the inequality test.

  15. 15.

There are $\lceil\log N_{k}(M+N/2)\rceil$ qubits used for the contiguous
register for the QROM for the qubit rotations.

  16. 16.

There are $N\beth k_{r}$ used for the QROM for the rotations, with another

$\left\lceil\log\left(\frac{N_{k}(M+N/2)}{k_{r}}\right)\right\rceil-1$ (162)

temporary qubits.

Accounting for the maximum total involves determining the maximum number used
which depends on the use of temporary ancillas. We first need to determine
whether the temporary qubits described in part 13 or the total qubits in parts
14 to 16 are larger. We take the maximum of these, and add it to the qubits
used in parts 8 to 12. Then we compare that to the number of temporary qubits
in part 7. The maximum of that is added to the qubits in parts 1 to 6.

Figure 15: This shows how to construct a self-inverse procedure for block
encoding two-electron terms, as in Fig. 6 of [32]. The left side is the
manifestly self-inverse form, and the right side is a more intuitive form
where the $\mathinner{|{+}\rangle}$ state is used to generate the symmetry
between $\mu$ and $\nu$, and the two $V$ operations are controlled by $\mu$
and $\nu$ in succession.

Next we give a little more detail on how the construction is made self-
inverse. As explained in Fig. 6 of [32] (shown above as Figure 15) the THC
construction may be made self-inverse by using the qubit in the
$\mathinner{|{+}\rangle}$ state which controls the swap of the
$\mathinner{|{\mu}\rangle}$ and $\mathinner{|{\nu}\rangle}$ registers. As can
be seen in Figure 6, we are using a similar procedure, where the
$\mathinner{|{+}\rangle}$ state controls the swap of $\mu$ and $\nu$ at the
top left and at the lower right. The $X$ gate and swaps on the lower left
corresponds to the $X$ gate and swap in the middle of Figure 15.

We have currently just drawn the quantum circuit as having $c$ and
$c^{\dagger}$, but these would be implemented using ancilla qubits to control
the election between $X$ and $Y$ in $X\pm iY$. For the implementation to be
self-inverse, we would want the qubit used to control the first $c$ to be the
same as that for the final $c^{\dagger}$. This can be achieved by taking the
four qubits for control of each of the $c$ and $c^{\dagger}$ so the top one
can be used as the control each time. In particular, after the first $c$, swap
the first two qubits, then after the $c^{\dagger}$ swap the first pair with
the second pair, then after the next $c$ swap the first two qubits again. As a
result, the first qubit can be used as the control each time. Moreover, this
arrangement of swaps is obviously self-inverse.

The application of the controlled $X$ and $iY$ gates is also automatically
self-inverse. The reason is that $iY$ squared is $-\openone$. In the block
encoding we perform $Z$ gates for the control qubits for $c$ but not
$c^{\dagger}$. Then performing the unitaries for the block encoding twice, we
have first that the final controlled $c^{\dagger}$ in the first block encoding
is matched with the first $c$ for the second block encoding. The same control
qubit is used, so the two controlled $X$ operations cancel, and the two
controlled $iY$ operations give $-\openone$. This cancels with the $Z$ gate on
that control qubit. In this way all the operations can be cancelled, and
because each time we have controlled $c$ and $c^{\dagger}$ matched, which
cancel the $Z$ gate on the control qubit. There is no additional Toffoli cost
for these swaps and phase gates, because they are all Clifford gates.

##  Appendix E Correlation diagnostics for LNO

structure | max($|t_{1}|$) | max($|t_{2}|$) | T1 | D1 | UHF $S^{2}$
---|---|---|---|---|---
C2/m | 0.2538 | 0.0330 | 0.0482 | 0.1912 | 0.7783
P21/c | 0.2313 | 0.0322 | 0.0472 | 0.2178 | 0.8027
P2/c | 0.1688 | 0.0571 | 0.0371 | 0.2089 | 1.0447
Table 8:  Some common diagnostics of strong correlation from the ROHF-CCSD
calculations for each of the distorted LiNiO2 structures. The UHF $S^{2}$
values in the final column are given per formula unit and should be compared
with the exact doublet ($\langle S^{2}\rangle=0.75$). The basis set (GTH-DZVP)
and other details of the calculations are described in Section V.2

For insulators like the distorted structures of LiNiO2, there are several
common diagnostics that are used in molecular calculations to identify cases
of “strong correlation.” Here we examine the maximum elements of the CCSD T1
and T2 tensors (max($|t_{1}|$) and max($|t_{2}|$)), which are commonly used as
a measure of correlation[118], as well as the T1[119] and D1[120, 121]
diagnostics computed from the ROHF $k$-point CCSD calculations. We also show
the expectation value of the $S^{2}$ operator for the UHF solution because
spin contamination in the UHF calculation can be a signature of strong
correlation. The results are shown in Table 8. As expected, these results do
not suggest particularly strong correlation. Only the max($|t_{1}|$) values
for the C2/m and P21/c structures and the UHF $S^{2}$ for the P2/c structure
are larger than might be expected. The larger max($|t_{1}|$) is likely an
indication that the ROHF orbitals are far from optimal, and the symmetry
breaking in the UHF calculations does not mean that CCSD cannot provide
reliable energies.

##  Appendix F Classical timing benchmarks

In order to compare the quantum algorithm run time to state of the art
classical algorithms we measured the cost to compute the CCSD and ph-AFQMC
total energy for the benchmark systems listed in Table 1 in double- and
triple-zeta basis sets. The results of these timings are presented in Fig. 16.

For CCSD we used pyscf [93, 94] and timed the data on a node with 30 3.1 GHz
Intel Xeon CPUs (30 OpenMP threads).

For ph-AFQMC, which is considerably more expensive than CCSD for small system
sizes, we estimated the run time by performing a short ph-AFQMC calculation
for each system using 8 Nvidia V100 GPUs. From this data we can then estimate
the run time to achieve an statistical error bar (per atom) of $1\times
10^{-4}$ Ha through the assumption that the statistical error of ph-AFQMC
decays like $N_{s}^{-1/2}$, where $N_{s}$ is the number of Monte Carlo
samples. Formally, ph-AFQMC should asymptotically scale like
$\mathcal{O}(N_{k}^{3})$ [22] assuming the number of samples required to reach
the desired precision does not scale with the system size. Interestingly, we
found that the statistical error bar per atom for a fixed number of samples
actually decreased with $N_{k}$, which implies the variance of ph-AFQMC is
increasing sub-linearly with the system size. For smaller system sizes it is
important to note that practically one can saturate the GPU with walkers with
nearly no loss speed [108], thus reducing the error bar given a fixed wall
time. As a result the small $N_{k}$ AFQMC numbers represent a large
overestimation in runtime one would practically need to obtain the desired
statistical error. Another confounding factor which may affect the scaling of
ph-AFQMC is the time step error (we fixed the time step at 0.005 Ha-1 for all
systems). Recent results suggest that ph-AFQMC suffers from a size extensivity
error [122], which is practically remedied through time step extrapolation.