ADVANCES IN SIMULATING FERMIONS ON A QUANTUM COMPUTER By Gabriel Given A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics—Doctor of Philosophy 2024 ABSTRACT One of the principal challenges in simulating fermions on a quantum computer is that qubits lack the anti-symmetry of fermions. The simplest solution, the Jordan-Wigner transformation, converts local interactions into non-local ones. I will describe a method based on Majorana fermions that preserves locality, and propose some improvements to it that reduce the CNOT gate cost and make the algorithm more suited to simulating nuclear matter. I will also suggest how a perturbation theory-based approach can be useful for studies in nuclear physics. Finally, I will discuss contributions I have made involving time fractals and quantum algorithms such as the rodeo algorithm, an eigenvalue estimation algorithm that can obtain precise results even on noisy quantum computers. ACKNOWLEDGEMENTS Thank you to everyone at FRIB for being part of a wonderful community to work in. Thank you also to my mom, Hannah, and Chessie for your support over the past several years. iii LIST OF ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Mathematical Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Quantum Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Variational Principle 1.4 Many-Body Wavefunctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Second Quantization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Lattice Effective Field Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Quantum Circuits . . . v 1 1 1 3 4 4 6 8 CHAPTER 2 . 15 FERMION TO QUBIT MAPPINGS . . . . . . . . . . . . . . . . . . . Jordan-Wigner Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Parity Basis . . . . . . . . . . . . . . . . . . . . . . 17 2.3 General Linear Transformation Mappings 2.4 Auxiliary Qubit Mappings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Majorana Fermion Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 . 21 2.6 Equivalence of Fermion to Qubit Mappings . . . . . . . . . . . . . . . . . . . . . CHAPTER 3 . . THE 𝜂 𝜒 METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1 The 𝜂 𝜒 Method . 3.2 Exponentiation of 𝑋 𝑋 + 𝑌𝑌 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3 Single-particle creation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 The 𝜂 𝜒 Method in 3 Dimensions . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.5 Using Perturbation Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.6 Perturbative Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Implementation of the 𝜂 𝜒 Method . . . . . . . . . . . . . . . . . . . . . . . . 38 3.7 CHAPTER 4 CONTRIBUTIONS TO VARIOUS PAPERS . . . . . . . . . . . . . . . 40 4.1 Scale Invariance and Time Fractals . . . . . . . . . . . . . . . . . . . . . . . . 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2 Projected Cooling . . 4.3 Rodeo Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.4 Quantum Monte Carlo and Superfluidity . . . . . . . . . . . . . . . . . . . . . 46 . . CHAPTER 5 CONCLUSIONS AND OUTLOOK . . . . . . . . . . . . . . . . . . . . 51 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 iv LIST OF ABBREVIATIONS CNOT Controlled Not Gate CZ NISQ Controlled 𝑍 Gate Noisy Intermediate-Scale Quantum (Computers) v CHAPTER 1 INTRODUCTION In this chapter, I will introduce nuclear lattice effective field theory, and Hamiltonians used in it. I will also provide some background about circuits used on quantum computers. In the second chapter, I will discuss various algorithms for converting Hamiltonians into quantum circuits. In the third chapter, I will discuss one particular algorithm in depth, including some improvements I have made to it and future applications of it. In the final chapter, I will discuss contributions that I have made to other papers authored by members of Dean Lee’s research group. 1.1 Mathematical Preliminaries Here I will define some mathematical notation that I will use. The commutator of two operators is defined [ 𝐎, 𝐵] = 𝐎𝐵 − 𝐵𝐎. The commutator is zero if the two operators commute. The anti-commutator of two operators is defined {𝐎, 𝐵} = 𝐎𝐵 + 𝐵𝐎. The anti-commutator is zero if the two operators anti-commute, i.e. 𝐎𝐵 = −𝐵𝐎. The Kronecker delta function is defined 𝛿𝑖 𝑗 = 1, 0, 𝑖 = 𝑗 𝑖 ≠ 𝑗 .    (1.1) (1.2) (1.3) An operator 𝑈 acting on a Hilbert space 𝐻 is unitary if 𝑈†𝑈 = 𝑈𝑈† = 𝐌, where 𝑈† is the adjoint, i.e. complex conjugate transpose, and 𝐌 is the identity matrix. Unitary operators have the property that they preserve inner products. 1.2 Quantum Mechanics In classical mechanics, the time evolution of an object is given by Newton’s 2nd law of motion, F = 𝑚a. 1 (1.4) If we can write down the forces acting on an object, we can solve this differential equation for the position of the object at any future time. Alternatively, we can write down kinetic energy 𝑇 and potential energy 𝑉 as functions of some generalized coordinates 𝑞 𝑗 and their time derivatives (generalized velocities) (cid:164)𝑞 𝑗 . Then we can solve the (mathematically equivalent to Newton’s 2nd law) Euler-Lagrange equations, d d𝑡 𝜕L 𝜕 (cid:164)𝑞 𝑗 = 𝜕L 𝜕𝑞 𝑗 , (1.5) where L = 𝑇 −𝑉. Solving either of these differential equations is not always practical. One example where they don’t work well is chaotic systems, systems that have exponential dependence on the initial conditions. This means that the initial conditions must be known exactly to get accurate predictions of the future behavior of the system. Additionally, there are situations where classical mechanics does not provide accurate predictions, such as in double-slit interference patterns of electrons. This is because, to the best of our knowledge, the world obeys quantum mechanics instead of classical mechanics. Classical mechanics remains useful as a limiting case of quantum mechanics for large systems. While classical mechanics is deterministic, quantum mechanics is probabilistic. Instead of solving for position as a function of time, we solve for a complex-valued wave function that represents a probability amplitude, Κ(𝑥, 𝑡). From the wave function we can calculate a probability density, which leads to the natural normalization of the wave function 𝜌(𝑥) = |Κ(𝑥, 𝑡)|2, ∫ ∞ −∞ |Κ(𝑥, 𝑡)|2 = 1. (1.6) (1.7) The state of a system need not be written as a wave function in the position basis; we can write the general state as |Κ(𝑡)⟩. The dynamics of the state are given by the time-dependent Schrödinger equation, 𝑖 𝑑 𝑑𝑡 |Κ(𝑡)⟩ = ˆ𝐻 |Κ(𝑡)⟩ , (1.8) 2 where ˆ𝐻 = ˆ𝑇 + ˆ𝑉 is the sum of kinetic and potential energy operators. If the dependence of |Κ⟩ on position and time are independent, we can use separation of variables to obtain the time-independent Schrödinger equation, and obtain the time evolution operator, ˆ𝐻 |Κ⟩ = 𝐞 |Κ⟩ , 𝑈 (𝑡) = 𝑒−𝑖 ˆ𝐻𝑡 . (1.9) (1.10) 1.3 Variational Principle Suppose that we have some Hamiltonian 𝐻 with eigenstates |𝑛⟩. Then we can define the eigenvalues (energies) of the eigenstates as 𝐻 |𝑛⟩ = 𝐞𝑛 |𝑛⟩ . (1.11) If we don’t know the ground state, we can still get an upper bound on the ground state energy 𝐞0 using the variational principle. We can write down a state that depends on some parameters 𝜃, |𝜓(𝜃)⟩. Since the eigenstates form a complete basis, we can decompose |𝜓(𝜃)⟩ into a linear combination over the eigenstates, |𝜓(𝜃)⟩ = ∑ 𝑛 |𝑛⟩ ⟚𝑛|𝜓(𝜃)⟩ = 𝑐𝑛 |𝑛⟩ . ∑ 𝑛 (1.12) Then the expectation value of the Hamiltonian will give an upper bound on the ground state energy: ⟚𝜓(𝜃)| 𝐻 |𝜓(𝜃)⟩ = = = 𝑐∗ 𝑚𝑐𝑛 ⟚𝑚| 𝐻 |𝑛⟩ 𝑐∗ 𝑚𝑐𝑛𝐞𝑛𝛿𝑛𝑚 |𝑐𝑛|2𝐞𝑛 ∑ 𝑛𝑚 ∑ 𝑛𝑚 ∑ 𝑛 ≥ 𝐞0. (1.13) (1.14) (1.15) (1.16) So minimizing the expectation value of the Hamiltonian over 𝜃 gives an estimate of the ground state energy. If there exist 𝜃 such that |𝜓(𝜃)⟩ is the ground state, then the ground state energy will be exact. 3 1.4 Many-Body Wavefunctions A single particle has a wave function 𝜓𝑖 (𝑥𝑖) that depends on the potential that the particle is in. If there are multiple particles, the product 𝜓𝑖 (𝑥𝑖) (cid:214) 𝑖 (1.17) approximates a solution to the Schrödinger equation [41]. For bosons, this is satisfactory. However, fermions obey the Pauli exclusion principle through wave functions that are antisymmetric under particle exchange. This requires, for example, that 𝜓1(𝑥1)𝜓2(𝑥2) = −𝜓1(𝑥2)𝜓2(𝑥1). (1.18) Slater notes that any permutation of the 𝑥𝑖 still results in an approximate solution to the Schrödinger equation. It is possible to construct a unique antisymmetric linear combination of these permuta- tions. This is most conveniently written in terms of a Slater determinant, 1 √ 𝑁! (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 𝜓1(𝑥1) 𝜓2(𝑥1) . . . 𝜓𝑁 (𝑥1) 𝜓1(𝑥2) 𝜓2(𝑥2) . . . 𝜓𝑁 (𝑥2) ... ... ... 𝜓1(𝑥𝑁 ) 𝜓2(𝑥𝑁 ) . . . 𝜓𝑁 (𝑥𝑁 ) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) . (1.19) Permuting two 𝑥𝑖 exchanges two columns of the matrix. By the properties of the determinant, this flips its sign. Thus the Slater determinant wave function has the antisymmetry we desire. It is also easy to verify that the Slater determinant satisfies the Pauli exclusion principle: if two particles are in the same state, two rows of the matrix will be identical, causing the determinant to be zero. It is possible to construct more complex many-body wave functions. For example, the Pfaffian wave function is an antisymmetric product that contains two-particle pairing functions as well as single-particle states [1]. 1.5 Second Quantization Before I introduce nuclear lattice effective field theory, I will introduce the language that its Hamiltonians are written in, second quantization. Second quantization describes many-body 4 systems by the number of particles in each single-particle state. This is very convenient for defining simulations on a lattice. Suppose that we have some (possibliy infinite) set of orthonormal single-particle states. These states could be, for example, lattice sites or electron orbitals. Then the many-body state can be written as |𝑛1, 𝑛2, 𝑛3, . . . ⟩ , (1.20) where 𝑛𝑖 is the number of particles in state 𝑖. For bosons (particles with integer spin), 𝑛𝑖 can be any nonnegative integer, but for fermions (particles with half-integer spin) 𝑛𝑖 is restricted to be 0 or 1 by the Pauli exclusion principle. The total number of particles in the system can be found by summing the occupation numbers of each state, 𝑁 = (cid:205)𝑖 𝑛𝑖. We introduce creation and annihilation operators that add or remove one particle from a many- body state. For bosons, we denote the creation operator as 𝑏† and the annihilation operator as 𝑏, and they are defined such that 𝑖 | . . . , 𝑛𝑖−1, 𝑛𝑖, 𝑛𝑖+1, . . . ⟩ = √𝑛𝑖 + 1| . . . , 𝑛𝑖−1, 𝑛𝑖 + 1, 𝑛𝑖+1, . . . ⟩ 𝑏† 𝑏𝑖 | . . . , 𝑛𝑖−1, 𝑛𝑖, 𝑛𝑖+1, . . . ⟩ = √ 𝑛𝑖 | . . . , 𝑛𝑖−1, 𝑛𝑖 − 1, 𝑛𝑖+1, . . . ⟩. (1.21) (1.22) Note that if we try to annihilate a particle that doesn’t exist, the state is replaced by 0. Due to the symmetry of bosonic wavefunctions under particle exchange, the creation and annihilation operators obey the following commutation relations: [𝑏† 𝑖 , 𝑏† 𝑗 ] = 0, [𝑏𝑖, 𝑏 𝑗 ] = 0, [𝑏𝑖, 𝑏† 𝑗 ] = 𝛿𝑖 𝑗 . (1.23) (1.24) (1.25) Nuclear matter is made of protons and neutrons, both particles with spin 1 2 and thus fermions. So we will primarily deal with fermionic creation and annihilation operators, denoted 𝑎† and 𝑎 or 𝑓 † and 𝑓 , respectively. Fermion wavefunctions are antisymmetric under particle exchange, so 5 these fermionic operators obey the anti-commutation relations { 𝑓 † 𝑖 , 𝑓 † 𝑗 } = 0, { 𝑓𝑖, 𝑓 𝑗 } = 0, { 𝑓𝑖, 𝑓 † 𝑗 } = 𝛿𝑖 𝑗 . (1.26) (1.27) (1.28) The fermionic operators act on many-particle states by 𝑓 † 𝑖 |𝑛1, . . . , 𝑛𝑖−1, 𝑛𝑖, 𝑛𝑖+1, . . . ⟩ = (−1) 𝑓𝑖 |𝑛1 . . . , 𝑛𝑖−1, 𝑛𝑖, 𝑛𝑖+1, . . . ⟩ = (−1) (cid:205)𝑖−1 𝑗=1 𝑛 𝑗 (1 − 𝑛𝑖))|𝑛1, . . . , 𝑛𝑖−1, 𝑛𝑖 + 1, 𝑛𝑖+1, . . . ⟩ (cid:205)𝑖−1 𝑗=1 𝑛 𝑗 𝑛𝑖 |𝑛1, . . . , 𝑛𝑖−1, 𝑛𝑖 − 1, 𝑛𝑖+1, . . . ⟩. (1.29) (1.30) Note that unphysical application of the creation and annihilation operators, trying to create a second particle in a state in violation of the Pauli exclusion principle or annihilating a particle that doesn’t exist, sends the state to 0. Starting from the vacuum state |0⟩ with no particles, any many body state can be constructed by repeated application of creation operators: (𝑎† 𝑝𝑖 )𝑛𝑖 = |𝑛1, 𝑛2, . . . ⟩ . (cid:214) 𝑖 (1.31) One useful operator constructed from creation and annihilation operators is the number operator, 𝑁𝑖 = 𝑓 † 𝑖 𝑓𝑖, which gives the number of particles in state 𝑖. 1.6 Lattice Effective Field Theory Nuclear lattice effective field theory is an established method for simulating nuclear matter [24, 25]. Lattice quantum chromodynamics, our best theory explaining matter, is too computationally expensive for systems of more than a few nucleons. So effective field theories are introduced, which describe interactions between protons and neutrons instead of between quarks. The Hamiltonian for leading order pionless effective field theory is 𝐻LO = 𝐻free + 𝑉 + 𝑉𝐌2 + 𝑉 (3𝑁). The free part of the Hamiltonian is 𝐻free = 1 2𝑚 ∫ ∑ 𝜎,𝜏 d3r ∇𝑎† 𝜎𝜏 (r) · ∇𝑎𝜎𝜏 (r), 6 (1.32) (1.33) which when discretized for a lattice simulation becomes 𝐻free = 3 𝑚 ∑ n,𝜎,𝜏 𝑎† 𝜎𝜏 (n)𝑎𝜎𝜏 (n)− 1 2𝑚 ∑ ∑ n,𝜎,𝜏 𝑙=𝑥,𝑊,𝑧 (cid:2)𝑎† 𝜎𝜏 (n)𝑎𝜎𝜏 (n + e𝑙) + 𝑎† 𝜎𝜏 (n)𝑎𝜎𝜏 (n − e𝑙)(cid:3) . (1.34) Under periodic boundary conditions, we can rearrange the terms in the second sum to obtain 𝐻free = 3 𝑚 ∑ n,𝜎,𝜏 𝑎† 𝜎𝜏 (n)𝑎𝜎𝜏 (n) − 1 2𝑚 ∑ ∑ n,𝜎,𝜏 𝑙=𝑥,𝑊,𝑧 (cid:2)𝑎† 𝜎𝜏 (n)𝑎𝜎𝜏 (n + e𝑙) + 𝑎† 𝜎𝜏 (n + e𝑙)𝑎𝜎𝜏 (n)(cid:3) , (1.35) which makes clear that the second term is a "hopping term" that moves particles between adjacent lattice sites. The potential terms are written in terms of nucleon density and local isospin density, 𝜌(r) = ∑ 𝜎,𝜏 𝑎† 𝜎𝜏 (r)𝑎𝜎𝜏 (r), 𝜌𝐌 (r) = ∑ 𝜎,𝜏,𝜏′ 𝑎† 𝜎𝜏 (r) [𝜏𝐌]𝜏𝜏′𝑎𝜎𝜏′ (r), (1.36) (1.37) where 𝜏𝐌 are Pauli matrices in isospin space. Using :: to indicate normal ordering, i.e. annihilation operators are to the right of creation operators, we can write 𝑉 = 𝑉𝐌2 = 𝑉 (3𝑁) = ∫ 𝐶 2 𝐶𝐌2 2 𝐷 6 ∫ d3r : [𝜌(r)]2 :, ∫ d3r : [𝜌𝐌 (r)]2 :, d3r : [𝜌(r)]3 : . (1.38) (1.39) (1.40) Recently, smeared creation and annihilation operators have been introduced to enforce conser- vation of orbital angular momentum [25]. The smeared operators also reduce lattice artifacts. The result is that short-range interactions now have a non-local term. The local portion of the short-range interactions is 𝑉𝐿 = 𝑐𝐿 2 + ∑ n 𝑐𝐌,𝐿 2 ∑ n,𝐌 : 𝜌𝐿 (n) 𝜌𝐿 (n) : + 𝑐𝑆,𝐿 2 : 𝜌𝐌,𝐿 (n) 𝜌𝐌,𝐿 (n) : + n,𝑆 𝑐𝑆,𝐌,𝐿 2 ∑ : 𝜌𝑆,𝐿 (n) 𝜌𝑆,𝐿 (n) : (1.41) : 𝜌𝑆,𝐌,𝐿 (n) 𝜌𝑆,𝐌,𝐿 (n) : , (1.42) ∑ n,𝑆,𝐌 and the non-local portion is 𝑉𝑁 𝐿 = 𝑐𝑁 𝐿 2 ∑ n : 𝜌𝑁 𝐿 (n) 𝜌𝑁 𝐿 (n) : + 𝑐𝐌,𝑁 𝐿 2 ∑ n,𝐌 : 𝜌𝐌,𝑁 𝐿 (n) 𝜌𝐌,𝑁 𝐿 (n) : , (1.43) 7 where : : indicates normal ordering. The interactions are written using smeared density operators defined as follows: 𝜌𝑁 𝐿 (n) = 𝑎† 𝑁 𝐿 (n)𝑎𝑁 𝐿 (n), 𝜌𝐌,𝑁 𝐿 (n) = 𝑎† 𝑁 𝐿 (n) [𝜏𝐌]𝑎𝑁 𝐿 (n), ∑ 𝜌𝐿 (n) = 𝑎†(n)𝑎(n) + 𝑠𝐿 𝑎†(n′)𝑎(n′), 𝜌𝑆,𝐿 (n) = 𝑎†(n) [𝜎𝑆]𝑎(n) + 𝑠𝐿 ∑ 𝑎†(n′) [𝜎𝑆]𝑎(n′), ⟹n′n⟩ 𝜌𝐌,𝐿 (n) = 𝑎†(n) [𝜏𝐌]𝑎(n) + 𝑠𝐿 ⟹n′n⟩ ∑ 𝑎†(n′) [𝜏𝐌]𝑎(n′), ⟹n′n⟩ 𝜌𝑆,𝐌,𝐿 (n) = 𝑎†(n) [𝜎𝑆 ⊗ 𝜏𝐌]𝑎(n) + 𝑠𝐿 𝑎†(n′) [𝜎𝑆 ⊗ 𝜏𝐌]𝑎(n′), ∑ ⟹n′n⟩ (1.44) (1.45) (1.46) (1.47) (1.48) (1.49) where the sum over ⟹n′n⟩ is over nearest neighbors on the lattice. We use 𝜎𝑆 to represent the spin Pauli matrices and 𝜏𝐌 to represent the isospin Pauli matrices. The smeared creation and annihilation operators are 𝑎𝑁 𝐿 (n) = 𝑎(n) + 𝑠𝑁 𝐿 ∑ 𝑎(n′), 𝑎† 𝑁 𝐿 (n) = 𝑎†(n) + 𝑠𝑁 𝐿 ⟹n′n⟩ ∑ 𝑎†(n′). ⟹n′n⟩ (1.50) (1.51) 1.7 Quantum Circuits In 1982, Richard Feynman proposed that suitable devices could by created to simulate any quantum system: quantum computers [19]. There are certain problems where it is believed that quantum computers can outperform classical computers. One famous example is the factorization of large numbers. RSA encryption relies on the inability of classical computers to find the prime factors of large numbers. However, in 1994, Shor proposed an algorithm that would allow a sufficiently large quantum computer to find the prime factorization of a large number [40]. Fortunately for the security of the internet, quantum computers are not yet capable of running Shor’s algorithm on large numbers, as we are still in the Noisy Intermediate-Scale Quantum (NISQ) era [34]. Two- qubit gates typically have error rates of order 0.1%. There is also an error rate when taking a 8 measurement of a quantum circuit, of order 1% for superconducting qubits. Intermediate-Scale refers to that most quantum computers have no more than a few hundred qubits. IBM has recently introduced a quantum computer with over 1000 qubits [8], but error rates continue to be one of the primary challenges in quantum computing. I wish to determine how to use this rapidly-improving technology for simulations of nuclear lattice effective field theory. Qubits store data in a quantum two-level system, with states |0⟩ = (cid:169) (cid:173) (cid:173) (cid:171) of a qubit can be written as a superposition of those two states, 1 0 (cid:170) (cid:174) (cid:174) (cid:172) 0 and |1⟩ = (cid:169) (cid:173) (cid:173) 1 (cid:171) (cid:170) (cid:174) (cid:174) (cid:172) . The state 𝑎|0⟩ + 𝑏|1⟩, (1.52) where 𝑎 and 𝑏 are complex numbers with |𝑎|2 + |𝑏|2 = 1. Qubits can store much more information than classical bits: it takes 2𝑁 classical bits to store the same amount of information as 𝑁 qubits. Analogously to logic gates on classical computers, there are one- and two-qubit quantum gates. The three Pauli matrices, 0 1 1 0 (cid:170) 𝑋 = (cid:169) (cid:174) (cid:173) (cid:174) (cid:173) (cid:171) (cid:172) 0 −𝑖 𝑌 = (cid:169) (cid:173) (cid:173) (cid:171) 𝑖 0 , , (cid:170) (cid:174) (cid:174) (cid:172) and 𝑍 = (cid:169) (cid:173) (cid:173) (cid:171) 1 0 0 −1 , (cid:170) (cid:174) (cid:174) (cid:172) (1.53) (1.54) (1.55) are used as single-qubit gates. The Hadamard gate creates a superposition of the |0⟩ and |1⟩ states, 𝐻 = 1 √ 2 1 1 1 −1 (cid:169) (cid:173) (cid:173) (cid:171) . (cid:170) (cid:174) (cid:174) (cid:172) The phase shift gate leaves the |0⟩ state unchanged and applies a phase 𝑒𝑖𝜑 to the |1⟩ state, 𝑃(𝜑) = (cid:169) (cid:173) (cid:173) (cid:171) 1 0 0 𝑒𝑖𝜑 . (cid:170) (cid:174) (cid:174) (cid:172) 9 (1.56) (1.57) Two important special cases of the phase gate are the 𝑆 gate, and the 𝑇 gate, 𝑆 = 𝑃(𝜋/2) = (cid:169) (cid:173) (cid:173) (cid:171) 1 0 0 𝑖 , (cid:170) (cid:174) (cid:174) (cid:172) I will also use a rotation about the 𝑧-axis, 𝑇 = 𝑃(𝜋/4) = (cid:169) (cid:173) (cid:173) (cid:171) 1 0 0 𝑒𝑖𝜋/4 . (cid:170) (cid:174) (cid:174) (cid:172) 𝑅𝑧 (𝜃) = 𝑒−𝑖𝑍𝜃/2 = (cid:169) (cid:173) (cid:173) (cid:171) 𝑒−𝑖𝜃/2 0 0 𝑒𝑖𝜃/2 . (cid:170) (cid:174) (cid:174) (cid:172) Gates for rotations about the 𝑥 and 𝑊 axes are defined similarly, 𝑅𝑥 (𝜃) = (cid:169) (cid:173) (cid:173) (cid:171) cos(𝜃/2) −𝑖 sin(𝜃/2) −𝑖 sin(𝜃/2) cos(𝜃/2) (cid:170) (cid:174) (cid:174) (cid:172) 𝑅𝑊 (𝜃) = (cid:169) (cid:173) (cid:173) (cid:171) cos(𝜃/2) − sin(𝜃/2) sin(𝜃/2) cos(𝜃/2) (cid:170) (cid:174) (cid:174) (cid:172) (1.58) (1.59) (1.60) (1.61) (1.62) In classical computation, the NAND gate is universal, i.e. any circuit can be written using only NAND gates. For quantum computation, universality requires a set of gates. One such universal set is {𝐻, 𝑇, 𝐶𝑁𝑂𝑇 } [4]. An important class of two-qubit gates is controlled 1-qubit gates. If the control qubit (indicated by a dot) is in the |1⟩ state, a single qubit gate is applied to the target qubit. Otherwise, the identity is applied. As the Pauli 𝑋 gate acts as a not gate (flipping |1⟩ to |0⟩ and vice versa), the controlled 𝑋 gate is also called a controlled not gate, CNOT. The controlled 𝑍 gate, CZ, has the special characteristic that it is symmetric in the control and target qubits. The circuit symbols for the CNOT and CZ gates are shown in Fig. 1.1. It is often useful (or necessary) to decompose two qubit gates into circuits composed of single qubit gates and a certain subset of two-qubit gates, such as only CNOT gates. This can be done by the transpilation tools in Qiskit [36]. A controlled 𝑅𝑧 gate is equivalent to the circuit shown in Fig. 1.2, composed of 𝑅𝑧 gates and CNOT gates. 10 Figure 1.1 From left to right, a CNOT gate, CZ gate, and alternate method of displaying a CZ gate. The upper qubit is the control qubit. 𝑍 Figure 1.2 Decomposition of controlled 𝑅𝑧 (2) gate into 𝑅𝑧 and CNOT gates. The lower qubit is the control qubit. This circuit is an example of the utility of controlled reversal gates [3]. Let 𝑅 be a product of single-qubit gates that anti-commutes with a Hamiltonian 𝐻. Then 𝑅𝐻𝑅† = −𝐻 follows by right-multiplication by 𝑅†. Thus a controlled 𝑅 gate, denoted 𝐶𝑅, flips the sign of time evolution (exponentiation) of 𝐻. In the circuit in Fig. 1.2, 𝑅 = 𝑋 and 𝐻 = 𝑍. From its definition, the 𝑅𝑧 (𝜃) gate is the time evolution of 𝑍 over time 𝜃/2. When the control qubit is on, both 𝑅𝑧 gates contribute a rotation over time 1; when the control qubit is off, the two 𝑅𝑧 gates cancel out, making the circuit into the identity. Thus the circuit creates the desired controlled 𝑅𝑧 (2). One other useful circuit is to exponentiate a Pauli string [38]. Fig. 1.3 shows the exponentiation of 𝑌 𝑋 𝑍. Single-qubit change of basis gates are applied to each qubit so that a single 𝑅𝑧 can rotate all of the qubits at once. We then return each qubit to the computational basis. This circuit can be controlled by simply controlling the 𝑅𝑧 gate: if the 𝑅𝑧 gate is turned off, the left and right sides of the circuit cancel each other. 1.7.1 Quantum Algorithms One old algorithm for using quantum computers to find eigenvalues is Quantum Phase Estima- tion [15]. Suppose that there is some unitary operator 𝑈 with eigenvalue 𝑒𝑖2𝜋𝜙, and that we know how to implement controlled 𝑈2𝑘 for 𝑘 = 0, 1, 2, . . . . The 𝑗th auxiliary qubit is used to control 11 Figure 1.3 Circuit for 𝑒−𝑖𝑌 (0) 𝑋 (1) 𝑍 (2) /2. This circuit can be controlled by controlling the 𝑅𝑧 gate. 𝑈2 𝑗 , so it ends up in the state |0⟩ + 𝑒𝑖2𝜋(2 𝑗 𝜙) |1⟩ . (1.63) Then an inverse Quantum Fourier Transform can extract the best 𝑚-bit approximation of 𝜙 with a high probability, where 𝑚 is the number of auxiliary qubits. This algorithm is not fault tolerant, so it is not practical to implement on NISQ-era devices. In 2014, a new eigensolver was introduced based on the variational principle: the Variational Quantum Eigensolver [32]. This is a hybrid quantum-classical algorithm. First, the Hamiltonian is written as a sum of terms such that the expectation value of each term can be calculated efficiently on a quantum computer (e.g. by writing the Hamiltonian as a sum of Pauli strings). Then, by linearity, the expectation value of the Hamiltonian is the sum of expectation values for each term. An initial state with variational parameters is chosen, and the expectation value of the Hamiltonian is calculated. The ground state will minimize the expectation value. The classical computer is then used to perform the minimization and choose the next set of parameters for the initial state. Compared to Quantum Phase Estimation, the Variational Quantum Eigensolver requires a larger number of smaller circuits. This reduces the time it takes to run each circuit, making it more practical given the short decoherence times of NISQ era hardware. It has been proposed that, for large systems that can be decomposed into subsystems and interaction terms, it is possible to use variational quantum eigensolvers to construct ground states of each subsystem, and then a second variational quantum eigensolver to find the ground state of the entire system [20]. 12 In 2012, a method was proposed to implement a linear combination of unitary operators on a quantum computer [12]. The group of unitary operators is not closed under addition, so a linear combination of unitaries cannot in general be implemented exactly on a quantum computer. However, it is possible to implement the linear combination 𝜅𝑈1 + 𝑈2 (up to normaliation) with an error rate at most 4𝜅/(𝜅 + 1)2. Iterating this method allows any general linear combination of unitaries to be implemented. The quantum approximation optimization algorithm is another method for finding eigenvalues [18, 23]. It involves two Hamiltonians: 𝐻𝑃, the Hamiltonian that we desire to find the eigenvalues of, and 𝐻𝑀 that is easy to solve. The algorithm takes advantage of the adiabatic theorem to start with an eigenvector of 𝐻𝑀 and end up with a good approximation of an eigenvector 𝐻𝑃. Time evolutions of the two Hamiltonians are applied alternately, according to a set of parameters that are optimized classically. 1.7.2 Quantum Hardware The algorithms I work with are hardware-independent. However it is still useful to understand some of the technologies used to build quantum computers. IBM is developing quantum computers based on superconducting qubits [43]. An 𝐿𝐶 circuit, formed of an inductor and a capacitor, creates a harmonic oscillator potential. This has evenly- spaced energy levels. An 𝐿𝐶 circuit therefore doesn’t work for a qubit, as there is no way to distinguish transitions between the ground state and the first excited state from higher-level transitions. So an anharmonicity is added to the circuit in the form of a Josephson junction, a very thin (insulator width of order 2 nm) superconductor-insulator-superconductor circuit element. This makes the transition energies different, so that the two-state system of the ground state and first excited state can be isolated. Low-loss cavity resonators are used to mediate interactions between qubits. The Josephson junction circuits and cavity resonators can be manufactured in large quantities, allowing large quantum computers to be constructed. One of the leading alternatives to superconducting quantum computers is trapped ion quantum computers, such as the Quantinuum H2 series [16]. The H2 computers use hyperfine states of 13 ytterbium. The ions, paired with barium coolant ions, are stored in a race track-shaped surface- electrode trap. It is possible to rearrange the order of the ions on the trap, and two-qubit gates can be applied to adjacent pairs of ions. This allows two-qubit gates to be applied to any pair of qubits, unlike for superconducting qubits, where two-qubit gates can only be applied to qubits connected by a resonator. 14 CHAPTER 2 FERMION TO QUBIT MAPPINGS In lattice effective field theory, we have a lattice of sites where the fermions can be located. We also have a Hamiltonian written in terms of fermionic creation and annihilation operators 𝑓 † and 𝑓 . Our goal now is to implement time evolution of this Hamiltonian on a quantum computer. We first rewrite the Hamiltonian as a sum of Pauli strings acting on qubits, 𝐻 = (cid:205) 𝑗 ℎ 𝑗 where ℎ𝑖 are Pauli strings. Then we can use the first order Lie-Trotter formula [13], 𝑒−𝑖𝐻𝑡 ≈ (cid:214) 𝑒−𝑖ℎ 𝑗 𝑡, 𝑗 (2.1) along with the circuit from Fig. 1.3 to exponentiate the Pauli strings to create a quantum circuit to represent our Hamiltonian. The first order Lie-Trotter formula has error 𝑂 (𝑡2). Higher order Lie-Trotter formulas have errors that scale with higher powers of 𝑡. I will primarily consider lattice effective theory in two dimensions. 2.1 Jordan-Wigner Transformation The simplest method to rewrite Hamiltonians in terms of Pauli strings is the Jordan-Wigner transformation [49, 38]. We use a lattice of the same size as the physical lattice, and use a one-to- one mapping between the fermion sites and qubits. We let the qubit states represent the occupation of the corresponding lattice sites: |1⟩𝑖 means that there is a fermion on site 𝑖 and |0⟩𝑖 means that there is no fermion on site 𝑖. The fermionic creation operator satisfies 𝑗 | 𝑓1 . . . 𝑓 𝑗−10 𝑓 𝑗+1 . . . 𝑓𝑛⟩ = (−1)(cid:205) 𝑗 −1 𝑓 † 𝑖=0 𝑓𝑖 | 𝑓1 . . . 𝑓 𝑗−11 𝑓 𝑗+1 . . . 𝑓𝑛⟩ and 𝑓 † 𝑗 | 𝑓1 . . . 𝑓 𝑗−11 𝑓 𝑗+1 . . . 𝑓𝑛⟩ = 0, (2.2) (2.3) where 𝑓𝑖 is the number of fermions on lattice site 𝑖 [5]. Note that Eq. 2.2 has a sign determined by parity; this is due to the anti-symmetry of fermion wave functions. Equation 2.3 enforces the Pauli exclusion principle by forbidding a second fermion from being added to a particular lattice site. 15 We want to find an equivalent operator written using Pauli matrices in order to construct a quantum circuit. A reasonable guess is ˆ𝑄+ 𝑗 = |1⟩⟚0| = ˆ𝑄− 𝑗 = |0⟩⟚1| = 1 2 1 2 0 0 (𝑋 𝑗 − 𝑖𝑌𝑗 ) = (cid:169) (cid:173) (cid:173) 1 0 (cid:171) 0 1 (𝑋 𝑗 + 𝑖𝑌𝑗 ) = (cid:169) (cid:173) (cid:173) (cid:171) 0 0 (cid:170) (cid:174) (cid:174) (cid:172) 𝑗 (cid:170) (cid:174) (cid:174) (cid:172) 𝑗 , . (2.4) (2.5) However these operators do not reproduce the parity term of Eq. 2.2 required due to the anti- commutation of fermions. To make acceptable qubit creation and annihilation operators, we must restore the anti-commutation relation of fermionic creation and annihilation operators, (cid:8)𝑎 𝑗 , 𝑎𝑘 (cid:9) = 0. (2.6) We note that (cid:110) ˆ𝑄± 𝑗 , 𝑍 𝑗 (cid:111) = 0, so appending a string of 𝑍 operators will restore anti-commutation, 𝑓 † 𝑗 → 𝑓 𝑗 → 𝑗−1 (cid:204) 𝑘=1 𝑗−1 (cid:204) 𝑘=1 𝑍𝑘 ⊗ ˆ𝑄+ 𝑗 , 𝑍𝑘 ⊗ ˆ𝑄− 𝑗 . (2.7) (2.8) The downside to the Jordan-Wigner transformation (in more than one dimension) is that the creation and annihilation operators are non-local, and are Pauli strings with length that scales with the system size. In two (or more) dimensions, the lattice sites need to be numbered in some order for the operators to be well-defined. For a hopping term on an 𝐿 by 𝐿 lattice, with the qubits numbered left to right and top to bottom, the horizontal hop will be efficient, 𝑂 (1), as the two sites have consecutive indices. However the vertical hop will have a Pauli 𝑍 string of length 𝑂 (𝐿) due to the gap in indices. In three dimensions, hops between layers of qubits would have 𝑍 strings of length 𝑂 (𝐿2). 16 2.2 Parity Basis Since the computational cost of the Jordan-Wigner transformation is due to calculating parity, one might suspect that it would be more efficient to store parity, 𝑝 𝑗 = (cid:205) 𝑗 𝑘=1 𝑓𝑘 , rather than occupation in the qubits [38]. This complicates the qubit creation and annihilation operators, 𝑗 = |0⟩⟚0| 𝑗−1 ⊗ ˆ𝑄± ˆ𝑃± 𝑗 − |1⟩⟚1| 𝑗−1 ⊗ ˆ𝑄∓ 𝑗 1 2 = (𝑍 𝑗−1 ⊗ 𝑋 𝑗 ∓ 𝑖𝑌𝑗 ), (2.9) (2.10) as we need to use ˆ𝑄− 𝑗 to create a particle if the parity of site 𝑗 is already 1. But we will also need to update the parity of all sites with indices greater than 𝑗 with Pauli 𝑋 matrices. So we have not solved the problem of the Pauli 𝑍 string, but merely moved it to a Pauli 𝑋 string on different qubits. For a hopping term on a 𝐿 by 𝐿 lattice, the horizontal hop will be efficient, 𝑂 (1), and the vertical hop will have a long Pauli 𝑍 string, 𝑂 (𝐿), the same as for the Jordan-Wigner transformation. 2.3 General Linear Transformation Mappings In two dimensions, the Jordan-Wigner transformation and parity basis methods are two examples of a general class of linear transformation-based methods, described in Ref. [44]. One particularly interesting case is a compromise between the occupation and parity bases that reduces the length of the Pauli strings: the Bravyi-Kitaev method gives strings of length 𝑂 (log 𝐿) instead of 𝑂 (𝐿), where 𝐿 is the size of the lattice [5, 38]. To construct the general linear transformation mapping, let 𝑛 = 𝐿2 be the total number of fermion sites. Given an invertible 𝑛 × 𝑛 matrix 𝐎 with elements in F2, we define a fermion to qubit transformation 𝝎 = 𝐎 𝒇 , 𝒇 = 𝐎−1𝝎, (2.11) where 𝒇 = ( 𝑓1 𝑓2 . . . 𝑓𝑛) is the vector of fermion mode occupation numbers, 𝝎 ∈ Z⊗𝑛 2 is the vector of qubit states, and all math is done in F2. For the Jordan-Wigner transformation, 𝐎 is the identity 17 matrix. For the parity transformation, 𝐎 is an upper-triangular matrix of 1s, 𝐎 = 1 1 1 1 (cid:170) (cid:174) 0 1 1 1 (cid:174) (cid:174) (cid:174) (cid:174) 0 0 1 1 (cid:174) (cid:174) (cid:174) 0 0 0 1 (cid:172) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) , 𝐎−1 = 1 1 0 0 (cid:170) (cid:174) 0 1 1 0 (cid:174) (cid:174) (cid:174) (cid:174) 0 0 1 1 (cid:174) (cid:174) (cid:174) 0 0 0 1 (cid:172) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) . (2.12) We will define several sets of qubits that will be convenient for writing the creation and annihilation operators. First are the update sets 𝑈 ( 𝑗), which is the set of qubits that must be updated when the occupation of qubit 𝑗 changes. 𝑈 ( 𝑗) is the set of qubits with nonzero entries in the 𝑗th column of 𝐎. Notice that for the parity basis the update set 𝑈 ( 𝑗) is all qubits with indices ≥ 𝑗, and for the Jordan-Wigner transformation 𝑈 ( 𝑗) contains only 𝑗. Similarly to the parity basis method, we will need to apply Pauli 𝑋 operators to qubits in the update set. To satisfy the anti-commutation of fermions, we need to calculate parity, 𝑝 𝑗 = (cid:205) 𝑗−1 𝑖=0 𝑓𝑖. We define the parity set 𝑃( 𝑗) as the set of qubits needed to calculate parity, which requires ∑ 𝜔𝑖 = 𝑖∈𝑃(𝐜) 𝑗−1 ∑ 𝑖=0 𝑓𝑖, (2.13) If we define 𝑅 as an upper-triangular matrix of 1s with zeros on the diagonal, 𝑃( 𝑗) contains the column indices of nonzero entries in row 𝑗 of 𝑅 𝐎−1. We will need to apply Pauli 𝑍 operators to qubits in the parity set. The third useful set is the flip set 𝐹 ( 𝑗), which contains the set of qubits necessary to determine the occupation of qubit 𝑗. This allows us to determine whether 𝑄+ or 𝑄− is necessary. 𝐹 ( 𝑗) contains the nonzero entries in row 𝑗 of 𝐎−1. We can then use all three sets to find qubit representations of the fermionic creation and annihilation operators: 𝑓 † 𝑗 → 𝑓 𝑗 → 1 2 1 2 (cid:204) 𝑘∈𝑈 ( 𝑗) (cid:204) 𝑘∈𝑈 ( 𝑗) (cid:169) (cid:173) (cid:171) (cid:169) (cid:173) (cid:171) 𝑋𝑘 (cid:170) (cid:174) (cid:172) 𝑋𝑘 (cid:170) (cid:174) (cid:172) (cid:169) (cid:173) (cid:171) (cid:169) (cid:173) (cid:171) (cid:204) I + 𝑙∈𝐹 ( 𝑗) (cid:204) 𝑚∈𝑃( 𝑗) 𝑍𝑙(cid:170) (cid:174) (cid:172) (cid:169) (cid:173) (cid:171) 𝑍𝑚(cid:170) (cid:174) (cid:172) (cid:204) I − 𝑙∈𝐹 ( 𝑗) (cid:204) 𝑚∈𝑃( 𝑗) 𝑍𝑙(cid:170) (cid:174) (cid:172) (cid:169) (cid:173) (cid:171) . 𝑍𝑚(cid:170) (cid:174) (cid:172) 18 (2.14) (2.15) 2.3.1 Bravyi-Kitaev Transformation The Bravyi-Kitaev transformation is a sort of compromise between the Jordan-Wigner and Parity bases [5, 38, 45]. The transformation matrix 𝐎𝐵𝐟 has a binary tree structure. It can be 𝐎𝐵𝐟,1 = (1) (2.16) recursively defined by and 𝐎𝐵𝐟,2𝑥+1 = 𝐎𝐵𝐟,2𝑥 0 0 𝐎𝐵𝐟,2𝑥 (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172) (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) , (2.17) ← 1 → where ← 1 → represents a row of ones. Recursive formulae for the update, parity, and flip sets are given in Ref. [45]. The parity and update sets have Pauli weight (string length) 𝑂 (log 𝐿), and hence the hopping terms also have weight 𝑂 (log 𝐿). However, for large systems, it is not possible to embed the Bravyi-Kitaev transformation into a square lattice of qubits [44]. 2.4 Auxiliary Qubit Mappings Auxiliary Qubit Mappings are an improvement on the Jordan-Wigner transformation that use auxiliary qubits to store the parity of sets of lattice qubits. These auxiliary qubits essentially provide shortcuts, avoiding the long Pauli strings of the Jordan-Wigner transformation [44]. It’s not clear to me how much the need to update the auxiliary qubits offsets this advantage. Ref. [44] claims a Pauli string length of 𝑂 (1). 2.5 Majorana Fermion Methods Several methods based on Majorana fermions also give constant scaling of Pauli string length with the lattice size. Mathematically, a Dirac fermion 𝑓 can be split into two Majorana fermions 𝛟 that are their own antiparticles: where 𝛟2 𝑗−1 = 𝑓 𝑗 + 𝑓 † 𝑗 𝛟2 𝑗 = −𝑖 (cid:16) 𝑓 𝑗 − 𝑓 † 𝑗 (cid:17) , 𝛟 𝑗 𝛟 𝑗 = 1 and 𝛟 𝑗 = 𝛟† 𝑗 . 19 (2.18) (2.19) (2.20) The Majorana fermions satisfy the anti-commutation relation (cid:8)𝛟 𝑗 , 𝛟𝑘 (cid:9) = 2𝛿 𝑗 𝑘 . (2.21) The 𝜂 𝜒 method that I use is a Majorana fermion method. First, I will describe a couple others. 2.5.1 Bravyi-Kitaev Superfast Several methods put qubits on the edges of a graph. Bravyi-Kitaev Superfast is one of them [5, 44]. The vertices of the graph represent the fermionic modes, and the edges represent interaction terms in the Hamiltonian. For example, if the term 𝑓 † 𝑖 𝑓 𝑗 appears in the Hamiltonian, then (𝑖, 𝑗) ∈ 𝐞, the set of edges. We use the Majorana fermions to define edge and vertex operators, 𝐎 𝑗 𝑘 = −𝑖𝛟2 𝑗−1𝛟2𝑘−1 (edge) 𝐵 𝑗 = −𝑖𝛟2 𝑗−1𝛟2 𝑗 (vertex). (2.22) (2.23) These operators can generate any parity-preserving fermionic operator. For example, a hopping term can be written as 𝑗 𝑓𝑘 + 𝑓 † 𝑓 † 𝑘 𝑓 𝑗 = 𝑖 2 ( 𝐎 𝑗 𝑘 𝐵𝑘 + 𝐵 𝑗 𝐎 𝑗 𝑘 ). (2.24) We choose an (arbitrary) orientation 𝜖 𝑗 𝑘 = ±1 for every edge. For each vertex k, we order the incident edges ( 𝑗, 𝑘), denoted by the symbol < 𝑘 . This allows us to define the edge and vertex operators: 𝐵𝑘 = (cid:204) 𝑍𝑎𝑘 , 𝑎:(𝑎,𝑘)∈𝐞 (2.25) . (2.26) 𝐎 𝑗 𝑘 = 𝜖 𝑗 𝑘 𝑋 𝑗 𝑘 (cid:204) 𝑍𝑏𝑘 𝑏:(𝑏,𝑘)< 𝑘 ( 𝑗,𝑘) (cid:204) 𝑍 𝑗 𝑐 𝑐:( 𝑗,𝑐)< 𝑗 ( 𝑗,𝑘) (cid:170) (cid:174) (cid:174) (cid:172) (cid:169) (cid:173) (cid:173) (cid:171) (cid:169) (cid:173) (cid:173) (cid:171) (cid:170) (cid:174) (cid:174) (cid:172) We observe that the maximum Pauli weight of a hopping term is a constant, independent of number of particles or lattice size. To prevent the Majoranas from interacting with themselves, it is necessary to apply a loop constraint. For any closed length-𝑙 loop in the graph, ( 𝑗1, . . . , 𝑗𝑙, 𝑗1), we have the constraint 𝐶 𝑗1,..., 𝑗𝑙 = 𝑖𝑙 𝐎 𝑗1 𝑗2 . . . 𝐎 𝑗𝑙−1 𝑗𝑙 𝐎 𝑗𝑙 𝑗1. (2.27) 20 Downsides of this method include that the number of fermions is required to be even and that there is no direct correspondence between qubit occupation and fermionic mode occupation. 2.5.2 Verstraete-Cirac Encoding Another Majorana fermion method is the Verstraete-Cirac encoding [47]. Auxiliary qubits are added for each lattice site, and set up in the ground state of the auxiliary Hamiltonian Haux = 𝑖𝑐𝑖𝑑 𝑗 ∑ ⟚𝑖, 𝑗⟩ (2.28) for Majoranas 𝑐 and 𝑑, and the sum over links in a graph. The operators 𝑖𝑐𝑖𝑑 𝑗 leave the ground state of Haux unchanged, so they can be appended to the physical Hamiltonian without changing the calculated energies. This conversion from quadratic to quartic operators is similar to what I will do with the 𝜂 𝜒 method, which I find easier to understand. Another presentation of the Verstraete-Cirac encoding is in Ref. [48]. They show that in 3 dimensions, the maximum Pauli weight of a hopping term is 12, independent of the lattice size. This is the same as I will find for the 𝜂 𝜒 method. 2.5.3 Other Majorana-Based Methods The compact mapping of Ref. [17] constructs edge and vertex operators in a way similar to Bravyi-Kitaev superfast. It applies only to a square lattice, and adds auxiliary qubits to every other face of the lattice. This results in requiring fewer qubits than the Verstraete-Cirac encoding, as well as having a lower Pauli string length for hopping terms. Bravyi-Kitaev superfast was generalized to the custom fermionic codes of Ref. [11]. It uses local Majoranas in a similar way to the 𝜂 𝜒 method. A discussion of parallelization of custom fermionic codes is in Ref. [6]. 2.6 Equivalence of Fermion to Qubit Mappings A recent paper has shown, using finite-depth generalized local unitary transformations, that two-dimensional fermion to qubit mappings are mathematically equivalent [10]. Starting with an exact bosonization (using Majorana fermions), they obtain the Jordan-Wigner transformation, Verstraete-Cirac encoding, Bravyi-Kitaev superfast, and other encodings. They find a “supercom- pact” mapping that has a qubit to fermion ratio of only 1.25 for spinless fermions. This is better 21 than the ratio of 2 for Verstraete-Cirac, Bravyi-Kitaev superfast, and 𝜂 𝜒, but comes at the expense of a potentially higher (although still constant) Pauli string length for hopping terms. I chose to work with the 𝜂 𝜒 method because I find it relatively easy to understand and implement, especially since it contains a direct correspondence between qubit occupation and fermionic mode occupation. 22 CHAPTER 3 THE 𝜂 𝜒 METHOD 3.1 The 𝜂 𝜒 Method The fermion to qubit mapping that I choose to use is the 𝜂 𝜒 method of Ref. [28]. It requires more qubits than some of the most efficient mappings, but it has more consistently efficient hopping terms and in my opinion it is simpler to understand and implement. A nice presentation of the 𝜂 𝜒 method for spinless fermions is in Appendix B of [33]. This section will be a similar presentation, but for the more complicated case of nuclear matter (protons and neutrons with spin) on a two-dimensional rectangular lattice with periodic boundary conditions. For each site, there can be spin-up and spin-down protons and neutrons, with at most one of each fermion species at a site due to the Pauli exclusion principle. This creates a 16-dimensional Hilbert space. We will use the frequently-used mathematical trick of rewriting the fermion creation and annihilation operators in terms of Majorana fermions, which satisfy {𝛟𝑖 r , 𝛟 𝑗 r′ } = 2𝛿𝑖 𝑗 𝛿r,r′: (𝛟1 r − 𝑖𝛟2 r ) (𝛟3 r − 𝑖𝛟4 r ) (𝛟5 r − 𝑖𝛟6 r ) 𝑓 𝑝 ↑r = 𝑓 𝑝 ↓r = 𝑓 𝑛 ↑r = 𝑓 𝑛 ↓r = 1 2 1 2 1 2 1 (𝛟7 2 r − 𝑖𝛟8 r ). (3.1) (3.2) (3.3) (3.4) We will next introduce two sets of auxiliary Majorana fermions on each lattice site. There will be 8 𝜂𝑖, corresponding to the internal degrees of freedom (e.g. site occupation, spin, isospin), and 4 𝜒𝑖 corresponding to links between adjacent lattice sites. We define bosonic operators, which (since they are products of two Majoranas) commute at different sites. We call these operators bosonic due to this commutation. 23 𝑖 𝑗 𝑖 𝑗 r = 𝑖𝜂𝑖 Θ r = 𝑖𝜂𝑖 r = −𝑖 𝜒𝑖 Λ 𝑖 𝑗 r r 𝜂 𝑗 r r 𝜒 𝑗 r 𝜒 𝑗 r Ί (3.5) (3.6) (3.7) We can think of Θ as a single-site operator, describing interactions between nucleon species on the same site. Similarly, Λ describes interaction between nucleons on different lattice sites. Ί is a “pass-through” operator, and will be used to connect lattice sites that are not nearest neighbors. We can write some useful identities about the Ί operator. Ί𝑖𝑙 = −Λ 𝑗𝑖Θ 𝑗 𝑘 Λ𝑘𝑙 = −Λ 𝑗𝑖𝑖𝜂 𝑗 𝜂𝑘𝑖𝜂𝑘 𝜒𝑙 = −𝑖Λ 𝑗𝑖𝑖𝜂 𝑗 𝜒𝑙 = −𝑖Λ 𝑗𝑖Λ 𝑗𝑙 . (3.8) (3.9) (3.10) (3.11) We define an association between our initial Majoranas and these bosonic operators, so that we can rewrite our Hamiltonian in terms of these bosonic operators: 𝑖𝛟𝑖 r 𝑖 𝑗 𝛟 𝑗 r = Θ r 𝑖𝛟𝑖 r 𝛟 𝑗 𝑖𝜇 r′ = Λ r Λ 𝑗 𝜈 r′ , (3.12) (3.13) where the indices 𝜇, 𝜈 are determined by the orientation of the edge connecting the lattice sites. We follow the convention of Ref. [30] that 4 is the +𝑥 direction, 3 is the −𝑥 direction, 2 is the +𝑊 direction, and 1 is the −𝑊 direction. Note that 𝛟𝑖 ≠ 𝜂𝑖, since multiplication by 𝜒𝜇 is necessary for the terms involving multiple sites. We now rewrite our 𝜂 𝜒 Majoranas as Dirac fermions. It may feel as though we are going in 24 circles, but as I just showed, there are subtle differences between the types of Majoranas. 𝜂1 = 𝑎 + 𝑎†, 𝜂2 = 𝑖(𝑎 − 𝑎†) 𝜂3 = 𝑏 + 𝑏†, 𝜂4 = 𝑖(𝑏 − 𝑏†) 𝜂5 = 𝑐 + 𝑐†, 𝜂6 = 𝑖(𝑐 − 𝑐†) 𝜂7 = 𝑑 + 𝑑†, 𝜂8 = 𝑖(𝑑 − 𝑑†) 𝜒1 = 𝑒 + 𝑒†, 𝜒2 = 𝑖(𝑒 − 𝑒†) 𝜒3 = 𝑔 + 𝑔†, 𝜒4 = 𝑖(𝑔 − 𝑔†) Because there are 6 fermion modes, we will need 6 qubits for each site (which we will soon reduce to 5). We apply the Jordan-Wigner transformation to write each fermionic operator in terms of Pauli matrices applied to qubits: 𝑎† → 𝑄 (1)+ 𝑏† → 𝑍 (1)𝑄 (2)+ 𝑐† → 𝑍 (1) 𝑍 (2)𝑄 (3)+ 𝑑† → 𝑍 (1) 𝑍 (2) 𝑍 (3)𝑄 (4)+ 𝑒† → 𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4)𝑄 (5)+ 𝑔† → 𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5)𝑄 (6)+ , where the indices are qubit indices, an occupied qubit is represented by |1⟩ = (cid:169) (cid:173) (cid:173) (cid:171) , and 0 (cid:170) (cid:174) (cid:174) 1 (cid:172) 𝑄+ = 1 2 (𝑋 − 𝑖𝑌 ). 25 This gives us 𝜂1 → 𝑋 (1), 𝜂3 → 𝑍 (1) 𝑋 (2), 𝜂5 → 𝑍 (1) 𝑍 (2) 𝑋 (3), 𝜂7 → 𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑋 (4), 𝜂2 → −𝑌 (1) 𝜂4 → −𝑍 (1)𝑌 (2) 𝜂6 → −𝑍 (1) 𝑍 (2)𝑌 (3) 𝜂8 → −𝑍 (1) 𝑍 (2) 𝑍 (3)𝑌 (4) 𝜒1 → 𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑋 (5), 𝜒2 → −𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4)𝑌 (5) 𝜒3 → 𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5) 𝑋 (6), 𝜒4 → −𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5)𝑌 (6) (3.14) (3.15) (3.16) (3.17) (3.18) (3.19) The 𝜂 and 𝜒 operators that we have constructed span a 64-dimensional Hilbert space, 4 times as large as the physical Hilbert space we are modeling. Following Ref. [30], we can restrict the on-site parton parity, (−𝑖)6Π8 𝑖=1𝜂𝑖Π4 𝑗=1 𝜒 𝑗 → −1 −𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5) 𝑍 (6) → 1. This shows that qubit 6 is redundant, and that we can eliminate operators acting on it by applying these substitutions: 𝑋 (6) → 1 𝑍 (6) → −𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5) 𝑌 (6) → −𝑖𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5). This changes two of our Pauli string expressions for the Majoranas: 𝜒3 → 𝑍 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑍 (5) 𝜒4 → 𝑖 We now are able to write the bosonic operators in terms of Pauli operators. (3.20) (3.21) (3.22) (3.23) (3.24) 26 Θ 𝑖 𝑗 r = 𝑖 r r (cid:169) (cid:173) −𝑍 (1) (cid:173) (cid:173) r (cid:173) (cid:173) r 𝑋 (2) −𝑌 (1) (cid:173) (cid:173) (cid:173) r 𝑌 (2) 𝑌 (1) (cid:173) (cid:173) (cid:173) −𝑌 (1) r 𝑍 (2) (cid:173) (cid:173) (cid:173) (cid:173) 𝑌 (1) r 𝑍 (2) (cid:173) (cid:173) (cid:173) r 𝑍 (2) −𝑌 (1) (cid:173) (cid:173) (cid:173) 𝑌 (1) r 𝑍 (2) (cid:173) (cid:171) r 𝑋 (3) r 𝑌 (3) r 𝑍 (3) r 𝑍 (3) r r r r r 𝑋 (4) r 𝑌 (4) 𝑍 (1) r r r 𝑖 r 𝑋 (2) −𝑋 (1) 𝑋 (1) r 𝑌 (2) r 𝑍 (2) −𝑋 (1) 𝑋 (1) r 𝑍 (2) −𝑋 (1) r 𝑍 (2) 𝑋 (1) r 𝑍 (2) r 𝑋 (3) r 𝑌 (3) r 𝑍 (3) r 𝑍 (3) r r r r 𝑋 (4) r 𝑌 (4) r r 𝑌 (1) r 𝑋 (2) r 𝑋 (2) 𝑋 (1) 𝑖 r r −𝑌 (1) r 𝑌 (2) r 𝑌 (2) −𝑋 (1) 𝑍 (2) r r r −𝑍 (2) r −𝑌 (2) r 𝑋 (3) 𝑌 (2) r 𝑌 (3) r 𝑍 (3) −𝑌 (2) 𝑌 (2) r 𝑍 (3) r 𝑋 (4) r 𝑌 (4) r r r r 𝑖 r 𝑋 (3) −𝑋 (2) 𝑋 (2) r 𝑌 (3) −𝑋 (2) r 𝑍 (3) r 𝑍 (3) 𝑋 (2) r 𝑋 (4) r 𝑌 (4) r r r r r r 𝑋 (3) r 𝑋 (3) 𝑌 (1) r 𝑍 (2) r 𝑍 (2) 𝑋 (1) 𝑌 (2) r 𝑋 (3) r 𝑋 (3) 𝑋 (2) 𝑖 r r −𝑍 (3) r r 𝑋 (4) −𝑌 (3) 𝑌 (3) r 𝑌 (4) r r −𝑌 (1) −𝑋 (1) r r r 𝑌 (3) r 𝑌 (3) r 𝑍 (2) r 𝑍 (2) r 𝑌 (3) −𝑌 (2) −𝑋 (2) r 𝑌 (3) 𝑍 (3) r r r 𝑖 −𝑋 (3) r 𝑋 (4) r 𝑌 (4) 𝑋 (3) r r r r 𝑋 (4) r 𝑋 (4) r r r 𝑍 (3) r 𝑍 (3) r 𝑋 (4) r 𝑋 (4) 𝑌 (1) r 𝑍 (2) r 𝑍 (2) 𝑋 (1) 𝑌 (2) r 𝑍 (3) 𝑋 (2) r 𝑍 (3) 𝑌 (3) r 𝑋 (4) r 𝑋 (4) 𝑋 (3) 𝑖 r r r −𝑌 (1) −𝑋 (1) r r 𝑌 (4) r 𝑌 (4) r r r 𝑍 (3) r 𝑍 (3) r 𝑌 (4) r 𝑌 (4) r 𝑍 (2) r 𝑍 (2) r 𝑍 (3) −𝑌 (2) −𝑋 (2) r 𝑍 (3) −𝑌 (3) r 𝑌 (4) −𝑋 (3) r 𝑌 (4) 𝑍 (4) r r r r −𝑍 (4) r 𝑖 (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172)𝑖 𝑗 r r 𝑋 (5) r 𝑋 (5) r r r 𝑍 (4) r 𝑍 (4) r 𝑋 (5) r 𝑋 (5) r 𝑍 (3) r 𝑍 (3) r 𝑍 (4) r 𝑍 (4) r 𝑋 (5) r 𝑋 (5) r 𝑍 (2) 𝑌 (1) r 𝑍 (2) 𝑋 (1) r 𝑍 (3) 𝑌 (2) r 𝑍 (3) 𝑋 (2) r 𝑍 (4) 𝑌 (3) r 𝑍 (4) 𝑋 (3) 𝑌 (4) r 𝑋 (5) r 𝑋 (5) 𝑋 (4) r r r r r Λ 𝑖 𝑗 r = (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) −𝑌 (1) −𝑋 (1) r r 𝑌 (5) r 𝑌 (5) r r r 𝑍 (4) r 𝑍 (4) r 𝑌 (5) r 𝑌 (5) r 𝑍 (3) r 𝑍 (3) r 𝑍 (4) r 𝑍 (4) r 𝑌 (5) r 𝑌 (5) r 𝑍 (2) r 𝑍 (2) r 𝑍 (3) −𝑌 (2) r 𝑍 (3) −𝑋 (2) −𝑌 (3) r 𝑍 (4) r 𝑍 (4) −𝑋 (3) r 𝑌 (5) −𝑌 (4) r 𝑌 (5) −𝑋 (4) r r r r r r r 𝑍 (5) r 𝑍 (5) r r r 𝑍 (4) r 𝑍 (4) r 𝑍 (5) r 𝑍 (5) r 𝑍 (3) r 𝑍 (3) r 𝑍 (4) r 𝑍 (4) r 𝑍 (5) r 𝑍 (5) r 𝑍 (2) 𝑌 (1) r 𝑍 (2) 𝑋 (1) r 𝑍 (3) 𝑌 (2) r 𝑍 (3) 𝑋 (2) r 𝑍 (4) 𝑌 (3) r 𝑍 (4) 𝑋 (3) r 𝑍 (5) 𝑌 (4) r 𝑍 (5) 𝑋 (4) r r r r r r −𝑋 (1) r 𝑌 (1) r r 𝑋 (2) −𝑍 (1) 𝑍 (1) r 𝑌 (2) −𝑍 (1) r 𝑍 (2) 𝑍 (1) r 𝑍 (2) −𝑍 (1) r 𝑍 (2) r 𝑍 (2) 𝑍 (1) r 𝑋 (3) r 𝑌 (3) r 𝑍 (3) r 𝑍 (3) r r r r r r 𝑋 (4) r 𝑌 (4) (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172)𝑖 𝑗 27 −𝑖 𝑍 (5) r 𝑌 (5) r r 𝑍 (4) r 𝑍 (3) r 𝑍 (2) 𝑍 (1) r 𝑋 (5) r Ί 𝑖 𝑗 r = (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) −𝑍 (5) r −𝑖 𝑋 (5) r r 𝑍 (3) −𝑍 (1) r 𝑍 (2) r 𝑍 (4) r 𝑌 (5) r −𝑌 (5) r −𝑋 (5) r −𝑖 r 𝑍 (4) r 𝑍 (3) r 𝑍 (5) r r 𝑍 (2) 𝑍 (1) r r 𝑍 (2) 𝑍 (1) r 𝑍 (2) −𝑍 (1) r 𝑍 (2) 𝑍 (1) r 𝑋 (5) r 𝑌 (5) r 𝑍 (5) r 𝑍 (4) r 𝑍 (3) r 𝑍 (4) r 𝑍 (3) r 𝑍 (4) r 𝑍 (3) 𝑖 r r (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172)𝑖 𝑗 We will now consider how to set up an initial state by creating fermion pairs. We first consider a pair of particles adjacent in the x-direction. As an example, 𝑓 𝑝† ↑,r 𝑓 𝑛† ↑,r+x = r + 𝑖𝛟2 r ) (𝛟5 1 2 r+x + Λ14 r Λ53 (−𝑖Λ14 r+x + 𝑖𝛟6 r+x) r Λ63 r+x + Λ24 r Λ53 r+x + 𝑖Λ24 r Λ63 r+x) = = (𝛟1 1 2 1 4 1 4 1 4 1 4 = −𝑄+(1) = = (−𝑖Λ14 r (Λ53 r+x + 𝑖Λ63 r+x) + Λ24 r (Λ53 r+x + 𝑖Λ63 r+x)) (−𝑖Λ14 r + Λ24 r + 𝑌 (1) r (𝑖𝑋 (1) r ) (Λ53 r+x + 𝑖Λ63 r+x) r+x + 𝑖𝑋 (3) 𝑍 (5) 𝑍 (4) r+x r+x ) (𝑌 (3) r+x 𝑍 (4) r+x 𝑍 (5) r+x) r 𝑄+(3) r+x 𝑍 (4) r+x 𝑍 (5) r+x . (3.25) (3.26) (3.27) (3.28) (3.29) (3.30) To obtain the corresponding terms for other nucleon species, the first indices on the Λ operators are changed according to particle type. The 2nd indices are changed from {4,3} to {2,1} for particles adjacent in the y-direction. For two particles of different types on the same lattice site, we obtain 𝑓 𝑝† ↑ 𝑓 𝑛† ↑ = = = = = 1 2 1 4 1 4 1 4 1 4 (𝛟1 + 𝑖𝛟2) 1 2 (𝛟5 + 𝑖𝛟6) (−𝑖Θ15 + Θ16 + Θ25 + 𝑖Θ26) (−𝑖𝑌 (1) 𝑍 (2) 𝑋 (3) − 𝑌 (1) 𝑍 (2)𝑌 (3) + 𝑋 (1) 𝑍 (2) 𝑋 (3) − 𝑖𝑋 (1) 𝑍 (2)𝑌 (3)) (−𝑖𝑌 (1) 𝑍 (2) (𝑋 (3) − 𝑖𝑌 (3)) + 𝑋 (1) 𝑍 (2) (𝑋 (3) − 𝑖𝑌 (3))) (𝑋 (1) − 𝑖𝑌 (1))𝑍 (2) (𝑋 (3) − 𝑖𝑌 (3)) = 𝑄+(1) 𝑍 (2)𝑄+(3). (3.31) (3.32) (3.33) (3.34) (3.35) (3.36) The parity of the site is unchanged, so no operation on the auxiliary qubit is needed. 28 Now that we have constructed the initial state, we will write our Hamiltonian in terms of Pauli operators. The kinetic portion of Hamiltonians realistic to nuclear physics consists primarily of hopping terms. 𝑓 𝑝† ↑,r ↑,r′ + 𝑓 𝑝† 𝑓 𝑝 ↑,r′ 𝑓 𝑝 ↑,r = = = 1 4 1 2 1 2 (𝛟1 r + 𝑖𝛟2 r ) (𝛟1 r′ − 𝑖𝛟2 r′) + (𝑖𝛟2 r r′ − 𝑖𝛟1 𝛟1 r 𝛟2 r′) (Λ 2𝜇 r Λ1𝜈 r′ − Λ 1𝜇 r Λ2𝜈 r′ ) 1 4 (𝛟1 r′ + 𝑖𝛟2 r′) (𝛟1 r − 𝑖𝛟2 r ) (3.37) (3.38) (3.39) For a hop in the +x-direction, 𝜇 = 4 and 𝜈 = 3. This gives us 𝑓 𝑝† ↑,r ↑,r′ + 𝑓 𝑝† 𝑓 𝑝 ↑,r′ 𝑓 𝑝 ↑,r = 1 2 (𝑋 (1) r 𝑋 (1) r′ + 𝑌 (1) r 𝑌 (1) r′ )𝑍 (2) r′ 𝑍 (3) r′ 𝑍 (4) r′ 𝑍 (5) r′ . (3.40) For a hop in the +y-direction, 𝜇 = 2 and 𝜈 = 1. A Hamiltonian with smearing (as in Ref. [25]) would also have next-to-nearest neighbor hopping. As an example, we take a hop in the x-direction, where r′ = r + 2x. 𝑓 𝑝† ↑,r ↑,r′ + 𝑓 𝑝† 𝑓 𝑝 ↑,r′ 𝑓 𝑝 ↑,r = = (𝛟1 r + 𝑖𝛟2 r ) (𝛟1 r′ − 𝑖𝛟2 r′) + 1 4 (𝛟1 r′ + 𝑖𝛟2 r′) (𝛟1 r − 𝑖𝛟2 r ) 1 4 1 2 (𝑖𝛟2 r 1 2 − 𝑖𝛟1 r 1 2 − Λ14 1 2 (Λ24 = − = r′ − 𝑖𝛟1 𝛟1 r 𝛟2 r′) = − (𝑖𝛟2 r 𝛟1 r+x 𝑖𝛟1 r+x 𝛟2 r+x 𝑖𝛟2 r+x 𝛟1 r+2x 𝑖𝛟2 𝑖𝛟1 𝛟1 r+x r Λ13 𝛟2 r+x r+xΛ24 r+x r+xΘ12 r+x r+xΛ13 r+2x 𝛟2 r+2x) (Λ24 r Λ13 r Ί34 r+xΘ12 r+xΛ13 r+xΛ23 r+xΛ24 r+2x − Λ14 r+2x) r Ί34 r+xΛ23 r+2x) This shows the usefulness of the pass-through operators Ί. We can similarly consider next-to-nearest neighbor hopping in the 𝑊-direction, with r′ = r + 2y. 29 𝑓 𝑝† ↑,r ↑,r′ + 𝑓 𝑝† 𝑓 𝑝 ↑,r′ 𝑓 𝑝 ↑,r = = r′ − 𝑖𝛟1 𝛟1 r 𝛟2 r′) (𝛟1 r + 𝑖𝛟2 r ) (𝛟1 r′ − 𝑖𝛟2 r′) + 1 4 (𝛟1 r′ + 𝑖𝛟2 r′) (𝛟1 r − 𝑖𝛟2 r ) 1 4 1 2 (𝑖𝛟2 r 1 2 − 𝑖𝛟1 r 1 2 − Λ12 1 2 (Λ22 = − (𝑖𝛟2 r 𝛟1 r+y 𝑖𝛟1 r+y 𝛟2 r+y 𝑖𝛟2 r+y 𝛟1 r+2y 𝛟1 r+y 𝑖𝛟1 r+y 𝛟2 r+y 𝑖𝛟2 r+y 𝛟2 r+2y) = − (Λ22 r Λ11 r+yΘ12 r+yΛ22 r+yΛ11 r+2y r Λ11 r Ί12 r+yΘ12 r+yΛ11 r+yΛ21 r+yΛ22 r+2y − Λ12 r+2y) r Ί12 r+yΛ21 r+2y) = Typical Hamiltonians in nuclear physics also contain terms proportional to number operators: 𝑓 𝑝† ↑,r 𝑓 𝑝 ↑,r = = = r − 𝑖𝛟2 r ) 𝛟1 r + 1) r + 𝑖𝛟2 r ) (𝛟1 1 2 1 4 1 − Θ12 r 2 (1 − 𝑖𝛟1 r (𝛟1 1 2 r + 𝑖𝛟2 𝛟2 r 1 − 𝑍 (1) r 2 = . So we have seen how to write the Hamiltonian as as sum of Pauli strings. Even after eliminating a redundant qubit, the Hilbert space of the qubits at each lattice site is larger than the physical Hilbert space by a factor of 2. We will resolve this by introducing more constraints [28]. The first type is plaquette constraints that arise from constructing an identity operator around a plaquette in the lattice: r 𝑍 (2) 𝑍 (1) r 𝑍 (3) r 𝑍 (4) r 𝑌 (5) r 𝑋 (5) r+x 𝑌 (5) r+x+y r+xΊ13 r Ί32 Ί24 𝑍 (3) 𝑍 (2) r+y r+y r+x+yΊ41 𝑍 (4) r+y r+y = −1 𝑋 (5) r+y = 1 𝑍 (1) r+y We need to apply this constraint before creating any particles, to initialize the vacuum state. When there are no particles in the lattice, we can impose this plaquette constraint by applying the following operator to each plaquette: 𝐶r = 𝑌 (5) r 𝑋 (5) r+x 𝑌 (5) r+x+y 𝑋 (5) r+y . 30 There are not quite enough plaquette constraints to constrain the Hilbert space to its physical size due to the periodicity of the lattice, so we also have Wilson loop constraints. These arise from constructing identities along horizonal or vertical loops around the toroidal lattice, r Ί12 Ί12 r+y r Ί34 Ί34 r+x . . . Ί12 r−y = −1, . . . Ί34 r−x = −1. (3.41) (3.42) This constraint can be applied to the vacuum state by applying 𝑋 gates to the auxiliary qubits along the diagonal of the lattice [29]. 3.2 Exponentiation of 𝑋 𝑋 + 𝑌𝑌 As we just saw, the hopping terms in the Hamiltonian, once rewritten as Pauli strings, are often of the form (𝑋 𝑋 + 𝑌𝑌 )𝑍 𝑍, (3.43) with varying number of 𝑍 operators. I figured out a method to exponentiate these terms using one circuit instead of two, as would be needed for the basic exponentiation method. We start by noting that in general a two-qubit circuit can be implemented using only 3 CNOT gates [46]. Operators meeting certain conditions can be implemented even more efficiently: for example 𝑋 𝑋 + 𝑍 𝑍 can be exponentiated using only two CNOT gates [42]. This can be converted to 𝑋 𝑋 + 𝑌𝑌 or 𝑋𝑌 + 𝑌 𝑋, that occur in hopping terms, by applying single-qubit gates. We then use our idea of "controlled reversal gates" [3]. We note that the eigenvalues of the 𝑍 operator are 1 and −1, depending on the occupation of the qubit. Hence the effect of the 𝑍 string is to possibly apply a negative sign to the overall exponential. We know that 𝑍 anti-commutes with 𝑋 𝑋 + 𝑌𝑌 , so 𝐶 𝑍 will act as our controlled reversal gate. This gives the circuit shown in Fig. 3.1. Similarly to the controlled version of the canonical Pauli string exponential circuit, if the rotation gates are removed from the circuit in Fig. 3.1, the rest of the circuit cancels out. Using Qiskit [36], I was able to determine that the controlled version of my circuit takes 4 CNOT gates. This is the same number of CNOT gates as would be needed for controlling two copies of the canonical circuit. My controlled circuit is shown in Fig. 3.2. 31 Figure 3.1 Exponentiation of (𝑋 𝑋 + 𝑌𝑌 )𝑍 𝑍. The parity of the qubits that 𝑍 acts on is collected. Then the CZ gates flip the sign of the 𝑋 𝑋 + 𝑌𝑌 circuit if necessary. The H and S gates are needed to perform 𝑋 𝑋 + 𝑌𝑌 instead of 𝑋 𝑋 + 𝑍 𝑍. Figure 3.2 Controlled exponentiation of (𝑋 𝑋 + 𝑌𝑌 )𝑍 𝑍. The top qubit is the control qubit. 3.3 Single-particle creation So far, we have only considered the case where the initial state has an even number of particles. However by changing the direction of one link in the lattice and making the association 𝛟1 0 → Λ11 0 we can access the odd parity sector [33, 28]. At the origin, where the lattice defect is, 𝑓 𝑝† ↑,0 = = = 1 2 1 2 1 2 (𝛟1 0 + 𝑖𝛟2 0) 0 (1 + Θ12 Λ11 0 ) 𝑌 (1) 𝑍 (2) 𝑍 (3) 𝑍 (4) 𝑋 (5) (1 + 𝑍 (1)). (3.44) (3.45) (3.46) (For the rest of this section, all operators are implicitly at the origin.) The other cases are slightly more complicated, since we only have that 𝛟1 → Λ11. But we can insert the identity 𝛟1𝛟1 to add additional Θ terms and eliminate single 𝛟𝑖 for 𝑖 ≠ 1. 32 = − 𝑌 (2) 𝑍 (3) 𝑍 (4) 𝑋 (5) (1 + 𝑍 (2)) 𝑓 𝑝† ↓ = = = 𝑖 𝑓 𝑛† ↑ = = = 𝑖 𝑓 𝑛† ↓ = = = 𝑖 (𝛟3 + 𝑖𝛟4) 1 2 1 𝛟3(1 + Θ34) 2 1 𝛟1Θ13(1 + Θ34) 2 1 2 (𝛟5 + 𝑖𝛟6) 1 2 1 𝛟5(1 + Θ56) 2 1 𝛟1Θ15(1 + Θ56) 2 1 2 (𝛟7 + 𝑖𝛟8) 1 2 1 𝛟7(1 + Θ78) 2 1 𝛟1Θ17(1 + Θ78) 2 1 2 = − 𝑌 (4) 𝑋 (5) (1 + 𝑍 (4)) = − 𝑌 (3) 𝑍 (4) 𝑋 (5) (1 + 𝑍 (3)) Single particles can be created away from the origin by constructing chains of operators to connect to the origin, but this is unnecessary due to translational invariance. Further work needs to be done to determine the impact on the constraints of an odd number of particles. 3.4 The 𝜂 𝜒 Method in 3 Dimensions The 𝜂 𝜒 method can be extended to work on a three-dimensional lattice. The primary change is that there are now two more nearest neighbors to each lattice site. This requires adding two additional 𝜒 Majoranas for each lattice site. Accordingly, one more qubit is required per lattice site, and the lengths of the Pauli strings for hopping terms increase by up to 2. Additionally, due to the increased degree of freedom, there are many more plaquette constraints, since they now occur in three planes. There are also additional Wilson loop constraints. This approximately triples the computational cost of setting up the vacuum state. 33 Figure 3.3 Circuit used to calculate (cid:68) 𝜓 (𝑏) 0 (cid:12) (cid:12) (cid:12) O (cid:12) (cid:12) (cid:12) 𝜓 (𝑎) 0 (cid:69). Credit: Nicholas Cariello. 3.5 Using Perturbation Theory A common method used in physics is perturbation theory, where a complicated, hard to solve, Hamiltonian is decomposed into a part that is easy to solve plus a part that is harder to solve. Then observables, such as energy, can be approximated using the wave functions for the simple Hamiltonian. We want to do perturbation theory on a quantum computer. My colleague Nicholas Cariello constructed a circuit that can calculate (cid:68) (cid:69) is the 𝑎th-order (cid:69), where 𝜓 (𝑏) 0 (cid:12) (cid:12) (cid:12) O (cid:12) (cid:12) (cid:12) 𝜓 (𝑎) 0 (cid:12) (cid:12) (cid:12) 𝜓 (𝑎) 0 correction to the ground state wave function [7]. The circuit for this is shown in Fig. 3.3. The circuit contains single-qubit gates, time evolution of the base Hamiltonian, and controlled application of the perturbation. Monte Carlo methods will be used to handle perturbations that are sums of unitary operators. We consider a Hamiltonian of the form 𝐻0 + 𝑐𝐻1, where 𝐻0 has terms that preserve nucleon species and 𝐻1 has more complicated terms that mix nucleon species. This way we can do a simpler time evolution with only the 𝐻0 terms, and never have to exponentiate 𝐻1. Instead of using the full 𝜂 𝜒 method with 5 qubits per lattice site and the associated long Pauli strings, we could use something like the stacked lattices from Ref. [48] (a separate lattice with local encoding for each nucleon species) to reduce the CNOT cost. However, I’m not sure how we would then handle the perturbative terms, due to the separate local encodings. We can instead use a modification I propose: we can construct a composite lattice as shown in Fig. 3.4 and construct a local encoding using the 𝜂 𝜒 method on the entire composite lattice. Then we can easily define the perturbative terms. 34 𝑝 ↑ 𝑝 ↓ 𝑛 ↑ 𝑛 ↓ Figure 3.4 The layout of the composite lattice I propose. Proton and neutron spin-up and -down each are located in their own quadrant. We will use as the base Hamiltonian a Hamiltonian from leading order pionless effective field theory without a 3-body term [24]. where 𝐻0 = 𝐻free + 𝑉, 𝑉 = ∫ 𝐶 2 d3(cid:174)𝑟 : [𝜌((cid:174)𝑟)]2 :, normal ordering is represented by ::, and 𝜌((cid:174)𝑟) = ∑ 𝜎,𝜏 𝑓 † 𝜎𝜏 ((cid:174)𝑟) 𝑓𝜎𝜏 ((cid:174)𝑟). 𝐻free = 2 𝑚 ∑ ∑ (cid:174)𝑛 𝜎,𝜏 𝑓 † 𝜎𝜏 ( (cid:174)𝑛) 𝑓𝜎𝜏 ( (cid:174)𝑛) − 1 2𝑚 ∑ ∑ ∑ (cid:104) (cid:174)𝑛 𝜎𝜏 𝑙=𝑥,𝑊 𝜎𝜏 ( (cid:174)𝑛) 𝑓𝜎𝜏 ( (cid:174)𝑛 + (cid:174)𝑙) + 𝑓𝜎𝜏 ( (cid:174)𝑛) 𝑓 † 𝑓 † 𝜎𝜏 ( (cid:174)𝑛 + (cid:174)𝑙) (cid:105) . We note that the hopping terms preserve nucleon species, i.e. each hopping term is contained within one quadrant of the composite lattice. When we apply the 𝜂 𝜒 method, we have only one nucleon species per lattice site, as the four nucleon species are spread across lattice sites in the four quadrants of the composite lattice. This means that we will need two 𝜂 Majoranas and four 𝜒 Majoranas for each site. This initially requires three qubits. But, as in the case where we consider full nuclear matter at each lattice site, one qubit is redundant: 𝑖𝑍 (1) 𝑍 (2) 𝑍 (3) → 1. We can replace the Pauli operators acting on the 3rd qubit as follows: 𝑋 (3) → 1 𝑍 (3) → 𝑖𝑍 (1) 𝑍 (2) 𝑌 (3) → −𝑍 (1) 𝑍 (2) 35 We need a total of eight qubits for each physical lattice site: two (including one auxiliary qubit) for each nucleon species. However the bosonic operators have Pauli strings of length at most 2 instead of 5: Θ 𝑖 𝑖 𝑗 r = (cid:169) (cid:173) −𝑍 (1) (cid:173) r (cid:171) −𝑌 (1) −𝑋 (1) r 𝑌 (2) r 𝑌 (2) r r 𝑍 (1) r 𝑖 (cid:170) (cid:174) (cid:174) (cid:172)𝑖 𝑗 r 𝑍 (2) 𝑌 (1) r 𝑍 (2) 𝑋 (1) r r r r 𝑋 (2) 𝑌 (1) r 𝑋 (2) 𝑋 (1) r Λ 𝑖 𝑗 r = (cid:169) (cid:173) (cid:173) (cid:171) Ί 𝑖 𝑗 r = −𝑖 𝑍 (2) r 𝑌 (2) r r 𝑋 (2) −𝑍 (1) r (cid:169) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:173) (cid:171) −𝑍 (2) r −𝑖 𝑋 (2) r r 𝑌 (2) 𝑍 (1) r −𝑌 (2) r −𝑋 (2) r −𝑖 r 𝑍 (2) −𝑍 (1) r A hopping term in the +𝑥 direction gives 𝑋 (1) r −𝑌 (1) r (cid:170) (cid:174) (cid:174) (cid:172)𝑖 𝑗 r 𝑋 (2) −𝑍 (1) r 𝑌 (2) 𝑍 (1) r 𝑍 (2) 𝑖 −𝑍 (1) r r r (cid:170) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:174) (cid:172)𝑖 𝑗 𝑓 𝑝† ↑,r ↑,r′ + 𝑓 𝑝† 𝑓 𝑝 ↑,r′ 𝑓 𝑝 ↑,r = = 1 2 1 2 (Λ24 r Λ13 r 𝑋 (1) r′ − Λ14 r′ + 𝑌 (1) r Λ23 r′ ) r′ )𝑍 (2) r 𝑌 (1) r′ (𝑋 (1) . (3.47) (3.48) This is a shorter Pauli string than without using the composite lattice. At the cost of the extra auxiliary qubits, we now have a maximum Pauli string length for hopping terms of 4, while it was previously 10. As we have seen, number operators can be written as So ↑,r = 𝑓 𝑝† 𝑁 𝑝 ↑,r 𝑓 𝑝 ↑,r = 1 − Θ12 r 2 = 1 − 𝑍 (1) r 2 . 𝑁r𝑁r′ = 1 4 (1 − 𝑍 (1) r ) (1 − 𝑍 (1) r′ ). Exponentiating (𝐶/2)𝑁r𝑁r′ gives us, up to a constant prefactor, 𝑒−𝑖𝑡𝐶 𝑍 (1) r /8𝑒−𝑖𝑡𝐶 𝑍 (1) r′ /8𝑒𝑖𝑡𝐶 𝑍 (1) r 𝑍 (1) r′ /8. 36 We also should consider the special case where both number operators act on the same site: 𝑁r𝑁r = 1 4 (1 − 𝑍 (1) r ) (1 − 𝑍 (1) r ) = 1 4 (2 − 2𝑍 (1) r ) = 𝑁r. This is the expected result, as both eigenvalues of the number operator (0 and 1) are idempotent under multiplication. 3.6 Perturbative Term So far we have considered how to implement the base Hamiltonian, composed of hopping terms and number operator terms. We now consider as a perturbation an isospin term, 𝑉𝐌2 = 𝐶𝐌2 2 ∑ : r (cid:34) ∑ ( 𝜎 𝑓 𝑝† 𝜎,r 𝑓 𝑛 𝜎,r + 𝑓 𝑛† 𝜎,r 𝑓 𝑝 𝜎,r)2 + ( ∑ 𝜎 −𝑖 𝑓 𝑝† 𝜎,r 𝑓 𝑛 𝜎,r + 𝑖 𝑓 𝑛† 𝜎,r 𝑓 𝑝 𝜎,r)2 + ( ∑ 𝜎 𝑓 𝑝† 𝜎,r 𝑓 𝑝 𝜎,r − 𝑓 𝑛† 𝜎,r 𝑓 𝑛 𝜎,r)2 (cid:35) : 𝐶𝐌2 2 = ∑ ∑ r 𝜎≠𝜎′ (4 𝑓 𝑝† 𝜎,r 𝑓 𝑛 𝜎,r 𝜎,r 𝑓 𝑝 𝑓 𝑛† 𝜎,r − 2 𝑓 𝑝† 𝜎,r 𝑓 𝑝 𝜎,r 𝑓 𝑛† 𝜎,r 𝑓 𝑛 𝜎,r + 2 𝑓 𝑝† 𝜎,r 𝑓 𝑛 𝜎,r 𝑓 𝑛† 𝜎′,r 𝑓 𝑝 𝜎′,r + 2 𝑓 𝑛† 𝜎,r 𝑓 𝑝 𝜎,r 𝑓 𝑝† 𝜎′,r + 𝑓 𝑝† 𝜎,r 𝑓 𝑝 𝐶𝐌2 ∑ 2 𝜎,r 𝑓 𝑝† 𝜎′,r ∑ (−6 𝑓 𝑝† 𝑓 𝑝 𝜎,r 𝑓 𝑝 𝜎′,r − 𝑓 𝑝† 𝜎,r 𝑓 𝑝 𝜎,r 𝑓 𝑛† 𝜎,r 𝑓 𝑛† 𝜎′,r 𝜎,r + 4 𝑓 𝑝† 𝜎,r 𝑓 𝑛 𝜎′,r − 𝑓 𝑛† 𝜎,r 𝑓 𝑛 𝑓 𝑛 𝜎,r 𝑓 𝑛† 𝜎,r 𝑓 𝑛 𝜎,r 𝜎′,r = 𝑓 𝑝 𝜎′,r 𝑓 𝑝† 𝜎′,r 𝑓 𝑝 𝜎′,r + 𝑓 𝑛† 𝜎,r 𝑓 𝑛 𝜎,r 𝑓 𝑛† 𝜎′,r 𝑓 𝑛 𝜎′,r) r + 𝑓 𝑝† 𝜎,r 𝑓 𝑝 𝜎≠𝜎′ 𝜎,r 𝑓 𝑝† 𝜎′,r 𝑓 𝑝 𝜎′,r + 𝑓 𝑛† 𝜎,r 𝑓 𝑛 𝜎,r 𝑓 𝑛† 𝜎′,r 𝜎′,r − 2 𝑓 𝑝† 𝑓 𝑛 𝜎,r 𝑓 𝑝 𝜎,r 𝑓 𝑛† 𝜎′,r 𝑓 𝑛 𝜎′,r). (3.49) 𝑓 𝑛 𝜎′,r (3.50) (3.51) (3.52) (3.53) Except for the second term, these are products of two number operators. (We choose to not write the terms in normal ordering to make the number operators apparent.) So they are simple in principle to implement. Applying the 𝜂 𝜒 method to the second term, using a 2x2 composite lattice, we find (for 𝜎 =↑): 37 𝑓 𝑝 𝑓 𝑛† ↓,r ↓,r r+2x + 𝑖𝛟2 + 𝑓 𝑝† 𝑓 𝑛 ↓,r ↓,r r+2x)(𝛟1 (𝛟1 1 𝑓 𝑝 𝑓 𝑛† ↑,r = ↑,r 16 r+2y+2x) (𝛟1 r+2y+2x − 𝑖𝛟2 r + 𝑖𝛟2 (𝛟1 r+2y + 𝑖𝛟2 r+2y) (𝛟1 r − 𝑖𝛟2 r ) r ) (𝛟1 r+2y − 𝑖𝛟2 r+2y) (𝛟1 r+2y+2x + 𝑖𝛟2 r+2y+2x) (𝛟1 r+2x − 𝑖𝛟2 r+2x) = r+2y − 𝑖𝛟1 𝛟1 r r+2y + 𝑖𝛟2 𝛟2 r r+2y + 𝛟2 𝛟1 r 𝛟2 r+2y) r+2y+2x r+2x − 𝑖𝛟1 𝛟1 r+2y+2x − 𝑖𝛟1 𝛟1 r+2y+2x r+2x + 𝑖𝛟2 𝛟2 r+2y+2x + 𝑖𝛟2 𝛟2 r+2x r+2y+2x r+2x + 𝛟2 𝛟1 r+2y+2x + 𝛟2 𝛟1 𝛟2 r+2x) 𝛟2 r+2y+2x) r+2x r+2x r+2x r+2y r − 𝑖𝛟1 𝛟1 r Ί12 r+yΛ11 r+2y r + 𝑖𝛟2 𝛟2 r+2y − Λ12 r+2y r + 𝛟2 𝛟1 r+yΛ21 𝛟2 r ) r+2y + Λ22 r Ί12 r Ί12 r+yΛ11 r+2y − 𝑖Λ22 r Ί12 r+yΛ21 r+2y) r+y+2xΛ11 r+2xΊ12 r+2y+2x + Λ22 r+y+2xΛ11 r+2xΊ12 r+2y+2x − Λ22 r+y+2xΛ11 r+2xΊ12 r+2y+2x − Λ12 r+y+2xΛ11 r+2xΊ12 r+2y+2x + Λ12 r+y+2xΛ21 r+2xΊ12 r+2y+2x + 𝑖Λ22 r+y+2xΛ21 r+2xΊ12 r+2y+2x − 𝑖Λ22 r+y+2xΛ21 r+2xΊ12 r+2y+2x) r+y+2xΛ21 r+2y+2x) r+yΛ11 r + Λ22 r+2y + Λ12 r )(Ί12 r Ί12 r+yΛ21 r+y)(Λ11 r+2y − Λ22 r+2y − 𝑖Λ21 r Ί12 r+yΛ11 r+2y) (−𝑖Λ12 r+2y + 𝑖Λ22 r+2x + Λ22 r Ί12 r+yΛ21 r+2x) (−Ω12 r+2y) r+y+2x) (Λ11 r+2y+2x − 𝑖Λ21 r+2y+2x) r + Λ22 r ) (−Ω12 r+y) (Λ11 r+2y − 𝑖Λ21 r+2y) r+2x + Λ22 r − 𝑋 (1) r+2x)(Ί12 r 𝑌 (2) r+y+2x)(Λ11 r+y)(𝑌 (1) )(−𝑍 (2) r+2y+2x − 𝑖Λ21 𝑋 (2) r+2y r r 𝑌 (2) − 𝑋 (1) 𝑌 (2) r+2x r+2x )(−𝑌 (2) r + 𝑋 (1) r r 𝑍 (2) r+y 𝑋 (2) r+2y )(−𝑍 (2) r+y+2x r+2y+2x) (−𝑖Λ12 𝑋 (2) − 𝑖𝑋 (1) r+2y r+2y − 𝑖𝑋 (1) ) r+2y )(𝑌 (1) 𝑋 (2) r+2y+2x + 𝑋 (1) r+2y ) r+2y+2x )(−𝑖) (𝑖𝑌 (1) r+2y ) (−𝑖) (𝑖𝑌 (1) + 𝑋 (1) r+2x )(𝑌 (2) r+2x 𝑍 (2) r+y+2x r 𝑍 (2) r+y 𝑋 (2) r+2y 𝑄−(1) r+2y )(𝑄+(1) r+2x 𝑍 (2) r+y+2x 𝑋 (2) r+2y+2x 𝑌 (2) r+2x + 𝑋 (1) r+2y+2x ) 𝑄−(1) r+2y+2x ). r+2y+2x 𝑋 (2) r+2y+2x 𝑋 (2) r+2y+2x ) r+2y+2x 𝑓 𝑝† ↑,r 𝑓 𝑛 ↑,r + 1 16 (𝛟1 r 1 16 (𝛟1 r+2y+2x 1 (𝛟1 16 + + = = r+2y (−𝑖Λ12 (𝛟1 1 16 r+2xΊ12 (𝑖Λ12 1 (−𝑖Λ12 16 (𝑖Λ12 r Ί12 1 (−𝑖Λ12 16 1 16 1 8 (𝑖𝑌 (1) r+2x 1 (−𝑖𝑌 (1) 8 (−𝑖𝑌 (1) r+2x =2(𝑄+(1) r 𝑌 (2) (−𝑖Λ12 𝑌 (2) r+2x (𝑖𝑌 (1) = − = + The 𝜂 𝜒 method allowed us to determine which Pauli operators on auxiliary qubits are required. 3.7 Implementation of the 𝜂 𝜒 Method I prepared code to implement the (controlled and non-controlled) time evolution of nuclear lattice effective field theory Hamiltonians using the 𝜂 𝜒 method. I am initially using a 2 by 2 composite lattice, but writing the code so that it also works for larger lattice sizes. The general structure of the code is based on the code from Ref. [29]. I start by constructing the grids of physical and auxiliary qubits, and helper functions for finding 38 the qubit indices corresponding to specific lattice sites. The physical qubits are numbered 0 to 15, and the auxiliary qubits are numbered 16 to 31. (For example, the 2 by 2 lattice for proton spin-up uses physical qubits 0, 1, 4, and 5, and auxiliary qubits 16, 17, 20, and 21.) Qiskit initializes every qubit in the |0⟩ state, which is not the vacuum state for the auxiliary qubits. I initialize the auxiliary qubits into their vacuum state by applying the plaquette and Wilson loop constraints. The Wilson loop constraints are implemented with 𝑋 gates along the diagonal of the auxiliary lattice, in this case qubits 16, 21, 26, and 31. The plaquette constraints are implemented by applying the operator r 𝑋 (2) 𝑌 (2) r+x 𝑌 (2) r+x+y 𝑋 (2) r+y (3.54) to each plaquette. Qubit 2 is the auxiliary qubit for each lattice site. Now that the system is initialized in the vacuum state, I can create pairs of particles for the initial state. The two particle creation operators are products of qubit creation operators on each physical qubit and Pauli operators on the auxiliary qubits. Since we know that we are starting from the vacuum state, we note that the action of 𝑄+ on the vacuum is the same as the action of 𝑋 on the vacuum. I can finally construct the Hamiltonian. For each lattice site, I apply the (controlled if necessary) hopping terms in the +𝑥 and +𝑊 directions, using the circuits from Fig. 3.1 and 3.2. I can then apply the number operator terms using (controlled) 𝑅𝑧 gates and the circuit from Ref. [42]. To assist with constructing the hopping terms, I wrote code to calculate the Pauli strings for each bosonic operator. This code also allowed me to easily construct the bosonic operator matrices earlier in the chapter. 39 CHAPTER 4 CONTRIBUTIONS TO VARIOUS PAPERS 4.1 Scale Invariance and Time Fractals We consider bosons on a one-dimensional lattice, and show that with the correct choice of parameters there is discrete scale invariance and fractal-like time dependence [27]. 4.1.1 Discrete Scale Invariance If we introduce a single-site deep trapping potential to immobilize one boson at the origin, we obtain that for low energies, the Hamiltonian for the second boson is 𝐻 ( 𝑝, 𝑟) = 2𝐜0 sin(𝛌𝜋/2)Γ(1 − 𝛌)| 𝑝|𝛌−1 + 𝑉0 |𝑟 |𝛌−1 , (4.1) where 𝛌 is a parameter which we will set to 2. In the zero energy limit, we obtain wave functions 𝜓+(𝑟) = 𝜓−(𝑟) = 1 2 1 2 (cid:16) |𝑟 |𝑖𝛿+ + |𝑟 |−𝑖 ¯𝛿+ (cid:17) (cid:16) , sgn(𝑟) |𝑟 |𝑖𝛿− + |𝑟 |−𝑖 ¯𝛿− (cid:17) , where sgn(𝑟) is the sign function and 𝛿± satisfy the constraints 𝛿+ = 𝛿− = 𝑉0 𝐜0𝜋 𝑉0 𝐜0𝜋 coth(𝛿+𝜋/2), tanh(𝛿−𝜋/2). The values of 𝛿± are real when 𝑉0/𝐜0 is positive. When 𝑉0/𝐜0 ≫ 𝜋, 𝛿+ ≈ 𝛿− ≈ 𝑉0 𝐜0𝜋 . In the case where 𝛿± are real, we can rewrite the wave functions as 𝜓+(𝑟) = cos[𝛿+ ln(|𝑟 |)], 𝜓−(𝑟) = sgn(𝑟) cos[𝛿− ln(|𝑟 |)]. If we apply discrete scale transformations 𝑟 → 𝜆±𝑟, we obtain 𝜓+(𝑟) → cos[𝛿+ ln(|𝑟 |) + 𝛿+ ln(𝜆+)], 𝜓−(𝑟) → sgn(𝑟) cos[𝛿− ln(|𝑟 |) + 𝛿− ln(𝜆−)]. 40 (4.2) (4.3) (4.4) (4.5) (4.6) (4.7) (4.8) (4.9) (4.10) We obtain discrete scale invariance (up to a minus sign) if 𝜆+ = exp(𝜋/𝛿+), 𝜆− = exp(𝜋/𝛿−). (4.11) The energies scale as 𝐞± → 𝜆−1 ± 𝐞±. 4.1.2 Time Fractals We consider a state that is a linear combination of the first 𝑁 states from Eq. 4.2, We calculate the amplitude |𝑆⟩ = 𝑁−1 ∑ 𝑛=0 |𝜓 (𝑛) + ⟩. 𝐎(𝑡) = Re[⟚𝑆| exp(−𝑖𝐻𝑡)|𝑆⟩] 𝐎(𝑡) = 𝐎(𝑡) = 𝑁−1 ∑ 𝑛=0 𝑁−1 ∑ 𝑛=0 cos(𝐞 (𝑛) + 𝑡) cos(𝜖+𝜆−𝑛 + 𝑡). (4.12) (4.13) (4.14) (4.15) We observe that 𝐎(𝑡) is invariant under the rescaling 𝑡 → 𝜆+𝑡, thus it is a time fractal. We plot 𝐎(𝑡) in Fig. 4.1 for 𝐜0 = −1 and 𝑉0 = −14.2388293 so that 𝜆+ = 2. We also observe that 𝐎(𝑡) are closely related to the famous continous everywhere but differentiable nowhere Weierstrass function, 𝑀(𝑥) = ∞ ∑ 𝑛=0 𝑎𝑛 cos(2𝜋𝑏𝑛𝑥). (4.16) Fractal dimension is a way of quantifying the scaling of a system. Similarly to how a two- dimensional square has area proportional to the 2nd power of its side length and a 3-dimensional cube has volume proportional to the 3rd power of its side length, an object with fractal dimension 𝐷 has some size proportional to the 𝐷th dimension of some parameter. The Weierstrass function has fractal dimension 𝐷 = 2 + log 𝑎 log 𝑏 , (4.17) for several different definitions of fractal dimension [39]. Our amplitude 𝐎(𝑡), in the limit that 𝑁 → ∞, becomes the Weierstrass function where 𝑎 = 1 and 𝑏 = 𝜆+. This gives our amplitude a fractal dimension of 2. I was able to numerically estimate the fractal dimension to be about 1.9. 41 Figure 4.1 An example of time fractals. The amplitude 𝐎(𝑡) is displayed over the range from 𝑡 = 0 to 80 in the upper left, 𝑡 = 0 to 160 in the upper right, 𝑡 = 0 to 320 in the lower left, and 𝑡 = 0 to 640 in the lower right. When the time displayed doubles, the overall shape of the amplitude remains the same. 4.2 Projected Cooling The projected cooling algorithm, based on the principles of evaporative cooling, can be used on a quantum computer to prepare the localized ground state of any Hamiltonian with a translationally- invariant kinetic energy and interactions that vanish at long distances [26]. Suppose we have a one-dimensional chain of 2𝐿 + 1 qubits and a particle number-preserving Hamiltonian defined on it. We use the state |1⟩ to represent that a qubit is occupied by a particle. We define a compact region 𝜌 to be 2𝑅 + 1 qubits where 𝑅 ≪ 𝐿. We define a projection operator 𝑃 as the product of |0⟩⟚0| over all qubits outside 𝜌 (and the identity inside 𝜌). 𝑃 projects onto the subspace where all particles are within 𝜌. 42 Suppose that our Hamiltonian’s ground state |𝜓0⟩ is a localized bound state, and the only bound state of 𝐻. Let |𝜓𝐌⟩ be some initial state contained in 𝜌 that has nonzero overlap with |𝜓0⟩. The time evolution oeprator is 𝑈 (𝑡) = 𝑒−𝑖𝐻𝑡 . (4.18) In the limit 𝐿 → ∞, 𝑃𝑈 (𝑡)𝑃 has a stable fixed point proportional to 𝑃 |𝜓0⟩. This means that over time, the excited states leave the region 𝜌 and are killed off by the projection operator. So we are left with the ground state. The principle behind projected cooling can be illustrated by considering a Gaussian wave packet. Suppose that at time 0, we have a Gaussian wave packet centered on the origin with width √ 𝜎 = 𝑎. The standard formula for wave packet spreading gives that 𝜎(𝑡) = √ 𝑎2+(𝑡/𝑚)2 𝑎 . At time 𝑡, the fraction of the wave packet contained in our region 𝜌 is given by Performing a 𝑢-substitution with 𝑢 = 1 √ 2 𝑓 (𝑡) = ∫ 𝑅 (cid:17) 2 (cid:16) 𝑥 𝜎 (𝑡 ) 𝑒− 1 2 d𝑥. −𝑅 2𝜋 (cid:17), we obtain a result in terms of the error function, (4.19) 1 √ 𝜎(𝑡) (cid:16) 𝑥 𝜎(𝑡) 𝑓 (𝑡) = erf (cid:169) (cid:173) (cid:173) (cid:171) This allows us to normalize the Gaussian wave function truncated to within the region 𝜌: 2 𝑎2+(𝑡/𝑚)2 𝑎 (cid:170) (cid:174) (cid:174) (cid:172) √ . (4.20) 𝑅 𝜓(𝑥) = 0 √ 𝜋1/4 1 𝜎erf(𝑅/𝜎) 𝑒− 1 2 ( 𝑥 𝜎 )2 |𝑥| > 𝑅 . |𝑥| ≀ 𝑅 (4.21)    We compute the overlap of this wave function with the ground state wave function of a particle in an infinite square well, The overlap is 𝜓(𝑥) =    0 √ 1 𝑅 cos (cid:0) 𝜋 2𝑅 𝑥(cid:1) |𝑥| > 𝑅 . |𝑥| ≀ 𝑅 𝑒− 𝜋2 𝜎2 8𝑅2 𝜋1/4𝜎 𝑂 (𝑡) = (cid:16)erf (cid:104) 𝑅 (cid:105) √ 2𝜎 − 𝑖𝜋𝜎 + √ 2𝑅 2 √2𝑅𝜎erf (𝑅/𝜎) (cid:104)erf 𝑅 √ 2𝜎 43 (4.22) (4.23) (cid:105) (cid:17) + 𝑖𝜋𝜎 √ 2𝑅 2 . Figure 4.2 A comparison of my Gaussian wave packet model to projected cooling algorithm simulations. The ground state overlap computed in [26] is in blue, and the overlap between a truncated Gaussian wave packet and a particle-in-a-box wave function is in orange. In Fig. 4.2, I compare this model to the projected cooling algorithm simulations in Ref. [26]. I used a series expansion to eliminate the imaginary terms and fit my model to the simulated data. The shapes of the functions are in general agreement. 4.3 Rodeo Algorithm The rodeo algorithm is a method of eigenstate preparation and eigenvalue estimation on a quantum computer. It is exponentially faster than phase estimation and adiabatic evolution [35, 14]. Our initial work was with one-qubit Hamiltonians, but the Rodeo algorithm generalizes easily to larger Hamiltonians. One cycle of the (single-qubit) rodeo algorithm is illustrated in Fig. 4.3 . The object qubit 44 020406080100120140t0.20.40.60.81.0O(t) starts in some state |𝜓0⟩ with nonzero overlap with the ground state. Then 𝑁 cycles of the Rodeo algorithm are run, with 𝑡𝑘 chosen from a Gaussian distribution centered at 0 with a specified width 𝜎. Within each cycle, there is controlled time evolution of the Hamiltonian for time 𝑡𝑘 and a phase rotation applied to the ancilla qubit with a parameter 𝐞 that we choose. The probability that the ancilla qubit is still in the |0⟩ state after one cycle of the Rodeo Algorithm is 𝑃0(𝐞, 𝑡𝑘 ) = cos2 (cid:104) 𝑡𝑘 2 (𝐞obj − 𝐞) (cid:105) , (4.24) where 𝐞obj is the ground state energy of the object Hamiltonian. Performing 𝑁 cycles of the Rodeo Algorithm (each with a different 𝑡𝑘 but the same 𝐞) and marginalizing over all possible values of 𝑡𝑘 gives that the probability that every measurement of the ancilla qubit results in |0⟩ is 𝑃0𝑁 (𝐞) = (cid:34) 1 + 𝑒−(𝐞obj−𝐞)2𝜎2/2 2 (cid:35) 𝑁 . (4.25) We see that the probability goes exponentially to zero for 𝐞 ≠ 𝐞obj, and is one for 𝐞 = 𝐞obj. Intuitively, this is similar to interference in interferometry: when 𝐞 ≠ 𝐞obj there is destructive interference between the controlled time evolution and phase shift. Figure 4.3 One cycle of the rodeo algorithm. In general the object Hamiltonian can have multiple qubits. The object qubit is initialized to some state |𝜓0⟩, and the ancilla qubit is initialized to |0⟩. Since the probability of measuring the ancilla in the |0⟩ state decays exponentially, if we make a histogram of successful measurements for a variety of guesses 𝐞, we expect to find peaks at the eigenvalues of 𝐻obj. These peaks are robust to quantum computing errors: peaks will be lowered and the counts away from peaks will be increased, but the locations of the peaks will be unchanged. 45 Figure 4.4 Results from running the Rodeo Algorithm on IBM Q Casablanca. The eigenenergies can be extracted from the centers of the peaks. I determined that IBM’s at the time recent introduction of mid-circuit measurements allowed the Rodeo Algorithm to be run on actual quantum hardware without requiring an additional ancilla qubit for each cycle. We prepared the single-qubit Rodeo Algorithm and ran it on IBM Q Casablanca [35]. The results are shown in Fig. 4.4. We performed three scans over 𝐞, with increasing resolution as we located the peaks. The peaks are somewhat lower than predicted by theory and noiseless simulations, due to errors inherent in current quantum hardware. We were able to determine the energy eigenstates within a relative error of 0.08%. 4.4 Quantum Monte Carlo and Superfluidity I have also been involved with studies of superfluidity in neutron matter. I attempted to use quantum Monte Carlo methods to look for pairing effects indicative of simultaneous p-wave and s-wave superfluidity. I ran into some difficulties, as the particles appeared to enter bound states. Other people are continuing research into coexistence of superfluid states using auxiliary-field lattice simulations. 46 Figure 4.5 Results of Diffusion Monte Carlo simulation of fermions with an s-wave Pöschl-Teller interaction and a BCS wave function. There are pairs of up- and down-spin particles, with an additional up-spin for odd numbers of particles. We observe that completing a pair lowers the energy per particle. Superfluid systems exhibit pairing similar to Cooper pairs in superconductors, and quantum Monte Carlo simulations are able to observe this pairing [9]. Superfluidity in neutron matter is important because it occurs inside neutron stars. It is believed that there is s-wave pairing superfluidity in the low-density inner crust, and p-wave pairing in higher-density regions deeper inside neutron stars [31, 22, 50]. Neutron pairing is also relevant in neutron-rich nuclei [21]. We are interested in whether s-wave pairing and p-wave pairing can coexist, perhaps in an intermediate-density region inside neutron stars. We are able to observe s-wave pairing using a Pöschl-Teller interaction, 𝑉 (𝑟) = 1/cosh2(𝑟), in quantum Monte Carlo simulations, see Fig. 4.5. If we have only spin-up neutrons, we are able to observe 𝑝-wave pairing, shown in Fig. 4.6. To perform these calculations, we first run Variational Monte Carlo [2]. This uses the variational principle to tune a trial wave function and obtain an upper bound for the ground state energy. If we 47 Figure 4.6 Results of quantum Monte Carlo simulation of fermions with a p-wave Pöschl-Teller interaction and a Slater determinant wave function. All particles are spin-up. The plots show the differences between energies of adjacent numbers of particles. There appears to be pairing at large numbers of particles. The minima are at odd 𝑁 because the zero-momentum state is unpaired. We also observe shell closures at 𝑁 = 7 and 19. know the Hamiltonian and ground state wave function Κ0, we can calculate the energy: 𝐞0 = ∫ 𝑑Rι∗ 0 (R) ˆ𝐻Κ0(R) ∫ 𝑑R|Κ0(R)|2 . (4.26) Calculating this integral analytically or even numerically is impractical, as for 𝑁 particles there are 3𝑁 spatial coordinates, so we use Monte Carlo integration. We interpret the wave function as a probability distribution, and define a local energy So we can rewrite the energy as 𝑃(R) = ι∗(R)Κ(R) ∫ 𝑑R|Κ(R)|2 , 𝐞𝐿 (R) = ι−1(R) ˆ𝐻Κ(R). ∫ 𝐞 = 𝑑R𝐞𝐿 (R)𝑃(R). 48 (4.27) (4.28) (4.29) 2.55.07.510.012.515.017.520.0N0.50.00.51.01.52.0E [EF]ENEN1(Slater)EN.5EN1.5EN+1(Slater) If we have a set {R𝑖} distributed according to 𝑃(R), then by the central limit theorem 𝐞 ≈ 1 𝑀 𝑀 ∑ 𝑖=1 𝐞𝐿 (R𝑖). (4.30) We sample configurations R according to the probability distribution 𝑃(R) using the Metropolis algorithm. We construct a Markov chain with transition probabilities 𝑇 (R′ → R) that satisfy detailed balance, 𝑃(R)𝑇 (R′ → R) = 𝑃(R′)𝑇 (R → R′), (4.31) ensuring that the equilibrium of the Markov chain is 𝑃(R). If the probability of selecting a particular configuration update is Λ(R′ → R), then the acceptance rate (the probability that we accept the update) is 𝐎(R′ → R) = min (cid:20) 1, 𝑃(R)Λ(R → R′) 𝑃(R′)Λ(R′ → R) (cid:21) . (4.32) We can choose Λ(R → R′) to be symmetrical, such as by choosing a uniform or Gaussian distribution, causing the acceptance rate to be calculated solely from the wave function. We don’t know the exact ground state wave function, so we choose a wave function with parameters we can vary. Then by the variational principle, we obtain upper bounds of the ground state energy, 𝐞𝑉 = ∫ 𝑑Rι∗ 𝑇 (R) ˆ𝐻Κ𝑇 (R) ∫ 𝑑R|Κ𝑇 (R)|2 ≥ 𝐞0. (4.33) As a result of performing Variational Monte Carlo, we also obtain a set of particle configurations that samples the trial wave function. We use this as the initial state for Diffusion Monte Carlo (with importance sampling), a numerical method to solve the imaginary time Schrödinger equation which is a combination of a diffusion process and a branching process [37]. The main idea behind Diffusion Monte Carlo is the imaginary time evolution, Κ(𝜏) = 𝑒−( ˆ𝐻−𝐞𝑇 )𝜏Κ𝑇 ∑ = 𝑐𝑖𝑒−(𝐞𝑖−𝐞𝑇 )𝜏𝜓𝑖, (4.34) (4.35) where we expand in eigenstates of the Hamiltonian. If the trial energy 𝐞𝑇 ≈ 𝐞0, the ground state is projected out at large 𝜏, as the other eigenstates vanish. 𝑖 49 For each iteration of Diffusion Monte Carlo, the particle positions in each configuration are updated, and then the configurations are branched: configurations with low energies are duplicated, while configurations with high energies are removed. The position update includes a random term and a drift velocity term that moves the particles to regions with a large trial wave function. The moves are accepted or rejected using the Metropolis algorithm to satisfy detailed balance. The branching process duplicates and removes configurations according to their energies. Then the ground state energy is updated and the process is repeated. Diffusion Monte Carlo can in principle find the exact ground state energy. However, since fermion wave functions have both positive and negative regions, we avoid the fermion sign problem with the fixed node approximation which prohibits configurations from moving to a region with opposite sign. This means that the quality of the result depends on the choice of nodes of the trial wave function, and we can only obtain an upper bound for the ground state energy. 50 CHAPTER 5 CONCLUSIONS AND OUTLOOK There is a trade-off between number of qubits required and number of quantum gates required for fermion to qubit mappings. As we saw, for two-dimensional simulations, the 𝜂 𝜒 method requires 1.25 qubits per fermionic mode (lattice site), and up to 10 Pauli operators for each hopping term in the Hamiltonian. My composite lattice modification increases the number of qubits to 2 per mode, but decreases the maximum number of Pauli operators per hopping term to 4. The canonical circuit for exponentiation of a Pauli string of length 𝑛 requires 2(𝑛 − 1) two-qubit gates. I showed, using results from Ref. [42], that terms of the form (𝑋 𝑋 + 𝑌𝑌 )𝑍 𝑍, which are pairs of Pauli strings of length 𝑛, can be exponentiated using 2(𝑛 − 1) 2-qubit gates. This is half the number required for the canonical circuit. My motivation for obtaining these improvements is the limitations of current quantum com- puters. Gate error rates and decoherence time are two of the current limitations to quantum computation, and two-qubit gates have a larger computational cost for both of these than single- qubit gates. Thus I have sought to reduce the number of two-qubit gates required to simulate nuclear physics. Current nuclear lattice effective field theory simulations are limited by the amount of memory available on classical computers. I hope that my work will help enable larger systems to be studied using quantum computers once the hardware is sufficiently developed. 51 BIBLIOGRAPHY [1] M Bajdich et al. “Pfaffian pairing wave functions in electronic-structure quantum Monte Carlo simulations”. In: Physical review letters 96.13 (2006), p. 130201. [2] Michal Bajdich. Generalized Pairing Wave Functions and Nodal Properties for Electronic Structure Quantum Monte Carlo. 2007. arXiv: 0712.3066 [cond-mat.other]. url: https://arxiv.org/abs/0712.3066. [3] Max Bee-Lindgren et al. Controlled Gate Networks Applied to Eigenvalue Estimation. 2024. arXiv: 2208.13557 [quant-ph]. url: https://arxiv.org/abs/2208.13557. [4] P.Oscar Boykin et al. “A new universal and fault-tolerant quantum basis”. In: Information Processing Letters 75.3 (2000), pp. 101–107. issn: 0020-0190. doi: https://doi.org/ 10.1016/S0020-0190(00)00084-3. url: https://www.sciencedirect.com/science/article/pii/ S0020019000000843. [5] Sergey B Bravyi and Alexei Yu Kitaev. “Fermionic quantum computation”. In: Annals of Physics 298.1 (2002), pp. 210–226. [6] Jacob Bringewatt and Zohreh Davoudi. “Parallelization techniques for quantum simulation of fermionic systems”. In: Quantum 7 (2023), p. 975. [7] Nicholas Cariello. personal communication. Oct. 2024. [8] Davide Castelvecchi. “IBM releases first-ever 1,000-qubit quantum chip”. In: Nature 624.7991 (2023), pp. 238–238. [9] Soon-Yong Chang et al. “Quantum Monte Carlo studies of superfluid Fermi gases”. In: Physical Review A—Atomic, Molecular, and Optical Physics 70.4 (2004), p. 043602. [10] Yu-An Chen and Yijia Xu. “Equivalence between fermion-to-qubit mappings in two spatial dimensions”. In: PRX Quantum 4.1 (2023), p. 010326. [11] Riley W Chien and James D Whitfield. “Custom fermionic codes for quantum simulation”. In: arXiv preprint arXiv:2009.11860 (2020). [12] Andrew M Childs and Nathan Wiebe. “Hamiltonian simulation using linear combinations of unitary operations”. In: arXiv preprint arXiv:1202.5822 (2012). [13] Andrew M Childs et al. “Theory of trotter error with commutator scaling”. In: Physical Review X 11.1 (2021), p. 011020. [14] Kenneth Choi et al. “Rodeo algorithm for quantum computing”. In: Physical Review Letters 127.4 (2021), p. 040505. 52 [15] Richard Cleve et al. “Quantum algorithms revisited”. In: Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454.1969 (1998), pp. 339–354. [16] Matthew DeCross et al. “The computational power of random quantum circuits in arbitrary geometries”. In: arXiv preprint arXiv:2406.02501 (2024). [17] Charles Derby et al. “Compact fermion to qubit mappings”. In: Physical Review B 104.3 (July 2021). issn: 2469-9969. doi: 10.1103/physrevb.104.035118. url: http://dx.doi.org/ 10.1103/PhysRevB.104.035118. [18] Edward Farhi et al. Quantum Computation by Adiabatic Evolution. 2000. arXiv: quant-ph/ 0001106 [quant-ph]. url: https://arxiv.org/abs/quant-ph/0001106. [19] Richard P Feynman. “Simulating Physics with Computers”. In: International Journal of Theoretical Physics 21.6/7 (1982). [20] Keisuke Fujii et al. “Deep variational quantum eigensolver: a divide-and-conquer method for solving a larger problem with smaller size quantum computers”. In: PRX Quantum 3.1 (2022), p. 010346. [21] Stefano Gandolfi, Alexandros Gezerlis, and Joe Carlson. “Neutron matter from low to high density”. In: Annual Review of Nuclear and Particle Science 65.1 (2015), pp. 303–328. [22] A. Gezerlis, C. J. Pethick, and A. Schwenk. Pairing and superfluidity of nucleons in neutron stars. 2015. arXiv: 1406.6109 [nucl-th]. url: https://arxiv.org/abs/1406.6109. [23] Stuart Hadfield et al. “From the quantum approximate optimization algorithm to a quantum alternating operator ansatz”. In: Algorithms 12.2 (2019), p. 34. [24] Dean Lee. “Lattice simulations for few-and many-body systems”. In: Progress in Particle and Nuclear Physics 63.1 (2009), pp. 117–154. [25] Dean Lee. “Recent progress in nuclear lattice simulations”. In: Frontiers in Physics 8 (2020), p. 174. [26] Dean Lee et al. “Projected cooling algorithm for quantum computation”. In: Physics Letters B 807 (2020), p. 135536. [27] Dean Lee et al. “Time fractals and discrete scale invariance with trapped ions”. In: Physical Review A 100.1 (2019), p. 011403. [28] Kangle Li and Hoi Chun Po. “Higher-dimensional Jordan-Wigner transformation and aux- iliary Majorana fermions”. In: Phys. Rev. B 106 (11 Sept. 2022), p. 115109. doi: 10.1103/ PhysRevB.106.115109. url: https://link.aps.org/doi/10.1103/PhysRevB.106.115109. 53 [29] Jannes Nys and Giuseppe Carleo. “Quantum circuits for solving local fermion-to-qubit mappings”. issn: 2521-327X. doi: 10.22331/ q-2023-02-21-930. url: https://doi.org/10.22331/q-2023-02-21-930. In: Quantum 7 (Feb. 2023), p. 930. [30] Jannes Nys and Giuseppe Carleo. “Variational solutions to fermion-to-qubit mappings in two spatial dimensions”. In: Quantum 6 (2022), p. 833. [31] Georgios Palkanoglou and Alexandros Gezerlis. “Superfluid neutron matter with a twist”. In: Universe 7.2 (2021), p. 24. [32] Alberto Peruzzo et al. “A variational eigenvalue solver on a photonic quantum processor”. In: Nature communications 5.1 (2014), p. 4213. [33] Hoi Chun Po. Symmetric Jordan-Wigner transformation in higher dimensions. 2021. arXiv: 2107.10842 [cond-mat.str-el]. [34] John Preskill. (Aug. 2018), p. 79. //doi.org/10.22331/q-2018-08-06-79. “Quantum Computing in the NISQ era and beyond”. In: Quantum 2 issn: 2521-327X. doi: 10.22331/q-2018-08-06-79. url: https: [35] Zhengrong Qian et al. “Demonstration of the rodeo algorithm on a quantum computer”. In: arXiv preprint arXiv:2110.07747 (2021). [36] Qiskit contributors. Qiskit: An Open-source Framework for Quantum Computing. 2023. doi: 10.5281/zenodo.2573505. [37] Peter J Reynolds et al. “Fixed-node quantum Monte Carlo for moleculesa) b)”. In: The Journal of Chemical Physics 77.11 (1982), pp. 5593–5603. [38] Jacob T Seeley, Martin J Richard, and Peter J Love. “The Bravyi-Kitaev transformation for quantum computation of electronic structure”. In: The Journal of Chemical Physics 137.22 (2012). [39] Weixiao Shen. “Hausdorff dimension of the graphs of the classical Weierstrass functions”. In: Mathematische Zeitschrift 289 (2018), pp. 223–266. [40] Peter W. Shor. “Polynomial-Time Algorithms for Prime Factorization and Discrete Loga- rithms on a Quantum Computer”. In: SIAM Journal on Computing 26.5 (1997), pp. 1484– 1509. doi: 10.1137/S0097539795293172. [41] John C Slater. “The theory of complex spectra”. In: Physical review 34.10 (1929), p. 1293. [42] Adam Smith et al. “Simulating quantum many-body dynamics on a current digital quantum computer”. In: npj Quantum Information 5.1 (2019), p. 106. 54 [43] Matthias Steffen et al. “Quantum computing: An IBM perspective”. In: IBM Journal of Research and Development 55.5 (2011), pp. 13–1. [44] Mark Steudtner and Stephanie Wehner. “Quantum codes for quantum simulation of fermions on a square lattice of qubits”. In: Physical Review A 99.2 (2019), p. 022308. [45] Andrew Tranter et al. “The Bravyi–Kitaev transformation: Properties and applications”. In: International Journal of Quantum Chemistry 115.19 (2015), pp. 1431–1441. [46] Farrokh Vatan and Colin Williams. “Optimal quantum circuits for general two-qubit gates”. In: Physical Review A—Atomic, Molecular, and Optical Physics 69.3 (2004), p. 032315. [47] F Verstraete and J I Cirac. “Mapping local Hamiltonians of fermions to local Hamiltonians of spins”. In: Journal of Statistical Mechanics: Theory and Experiment 2005.09 (Sept. 2005), P09012. doi: 10.1088/1742-5468/2005/09/P09012. url: https://dx.doi.org/10. 1088/1742-5468/2005/09/P09012. [48] [49] James D Watson et al. “Quantum Algorithms for Simulating Nuclear Effective Field Theo- ries”. In: arXiv preprint arXiv:2312.05344 (2023). James D. Whitfield, Jacob Biamonte, and Alán Aspuru-Guzik. “Simulation of electronic structure Hamiltonians using quantum computers”. In: Molecular Physics 109.5 (2011), pp. 735–750. doi: 10.1080/00268976.2011.552441. [50] Shigehiro Yasui, Chandrasekhar Chatterjee, and Muneto Nitta. “Phase structure of neutron P 2 3 superfluids in strong magnetic fields in neutron stars”. In: Physical Review C 99.3 (2019), p. 035213. 55