Sigma models on quantum computers

We formulate a discretization of sigma models suitable for simulation by quantum computers. Space is substituted by a lattice, as usually done in lattice field theory, while the target space (a sphere) is replaced by the"fuzzy sphere", a construction well known from non-commutative geometry. Contrary to more naive discretizations of the sphere, in this construction the exact $O(3)$ symmetry is maintained, which suggests that the discretized model is in the same universality class as the continuum model. That would allow for continuum results to be obtained for very rough discretizations of the target space as long as the space discretization is made fine enough. The cost of performing time-evolution, measured as the number of CNOT operations necessary, is $12 L T/\Delta t $, where $L$ is the number of spatial sites, $T$ the maximum time extent and $\Delta t$ the time spacing.


I. INTRODUCTION
The advent of quantum computers opens up a new method to attack several physics problems which have, up to now, remained intractable. Perhaps the most interesting of those is the numerical treatment of manybody/field theories theories with sign problems. In particular, the nonperturbative calculation of real time observables, where very little progress has been made up to now [1][2][3], is an obvious target for quantum computation. Of course, the hope of attacking these problems hinges on being able to formulate quantum field theories in a way suitable for quantum computers. This topic is still in its infancy. The naive expectation is that fermionic fields can be more easily implemented in quantum computers as a qubit can encode the presence or absence of a fermion in a given state. This is born out by the few existing calculations that have been performed on quantum comptuers [4][5][6]. Bosonic fields are not so simply implemented. The attempts made up to now involve either eliminating the bosonic fields using some special property of the model or truncating the occupation number at any given site [7][8][9][10][11][12][13][14][15][16][17][18]. The situation is analogous to the early days of (classical) computing in field theory. Classical bits also seem more amenable at describing fermionic than bosonic fields as the cost of storing and manipulating reasonable approximations to real numbers was too high to be practical in the early days. There were at the time several attempts at substituting bosonic continuous field values by a finite set of values [19][20][21][22][23]. In all these schemes the symmetry of the model is reduced by the discretization of the bosonic fields.
When discretization reduces the symmetry it is unclear whether, in the spacetime continuum limit, the original model is recovered. For instance, the non-linear sigma model in one spatial dimension with fields taking values on a sphere was studied in the approximation where the sphere is substituted by the vertices of a platonic solid. It seems to still be controversial whether the dodecahedron model is in the same universality class as the original spherical model [24][25][26][27][28][29]. In the case of gauge theories the question, at least for abelian theories, was settled long ago: abelian gauge theories with any finite discrete group Z N are not in the same universality class as the U (1) model and do not approach the U (1) gauge theory as the spacetime continuum limit is taken [30][31][32].
This suggests that, to obtain the right continuum limit, we should construct a scheme where all the symmetries of the original model are maintained while discretizing the field variables in order to make the Hilbert space to have finite (and hopefully small) dimension. The topic of this paper is to present such a formulation. It is based on a well known construction in non-commutative geometry (the "fuzzy sphere" [33,34]) that has been used before in the study of (super)-membranes. The resulting system can be simulated on a quantum computer with two qubits per spatial site. We implement our simulation scheme on a simulated quantum computer and verify it produces the right results. The number of gates required is of the order of ∼ 12L(T /∆t), where L is the number of spatial sites, T the maximum time extent and ∆t the time discretization step.

II. SIGMA MODEL ON THE FUZZY SPHERE
The O(3) sigma-model is defined on a discretized space by the Hamiltonian H = r g 2 2 π(r) 2 + 1 2g 2 ∆x 2 (n(r + 1) − n(r)) 2 = r g 2 2 π(r) 2 + 1 g 2 ∆x 2 (1 − n(r + 1) · n(r)) , (1)  where n is a unit three-dimensional vector, π 2 is the Laplace-Beltrami operator on S 2 , and the sum runs over the L spatial lattice sites. The global symmetry n(r) → O ·n(r), where O is an orthogonal matrix, is evident. This model is asymptotically free in one spatial dimension [35].
In the Hamiltonian formalism, the wave function is a function of L copies of the sphere S 2 , ψ(n 1 , · · · , n L ). The Hilbert space is infinite dimensional even for L = 1. We will approximate this model by substituting the target space (the sphere) by the "fuzzy sphere" [33,34]. The fuzzy sphere is not defined as a subset of points of the sphere; instead, it is the functions on the sphere that are substituted by elements of a finite dimensional Hilbert space. Let us demonstrate the construction first in the L = 1 case. The wave function of the system is a function of n and can be expanded as with the constraint n i n i = 1. In the fuzzy sphere regularization we substitute this Hilbert space by the space of matrices where . Thus the infinite dimensional Hilbert space of functions on the sphere is substituted by a space of dimension (2j + 1) 2 . Since the space is finite dimensional it can be informally thought of as being the space of functions defined on a space with a finite number of points, the fuzzy sphere. In the j = 1/2 case, the space of matrices Ψ is four-dimensional and the fuzzy sphere can be informally thought of as a 4-point discretization of the sphere. However, the "points" of the fuzzy sphere are "spread out" and choosing them does not break rotation symmetry. Notice that the fuzzy sphere is not defined as a subset of points of the sphere. It is the algebra of functions on the sphere that is deformed into a (non-commutative) algebra given by matrix multiplication. Still, the fuzzy sphere is an approximation to the sphere in the sense that operation defined on the sphere can be approximated by equivalent constructions on the fuzzy sphere. For instance, the norm in the sphere and on the fuzzy sphere satisfy: We refer to the references [34] for a discussion of the standard geometrical constructs of the sphere framed in terms of the algebra of functions (and their extensions to the fuzzy sphere). In particular, the Hamiltonian of the sigma model with one spatial site is simply the Laplacian on the sphere. (This describes the quantum mechanics of a free particle on the sphere.) with κ a normalization factor. The eigenvalues of the Laplacian operator on the sphere are l(l + 1) for l = 0, 1, . . . with mulitplicities 2l + 1. When κ = j(j + 1) the spectrum of its fuzzy version H 0 is exactly the same but truncated to its lowest (2j + 1) 2 values. This is in contrast to other discretizations of the sphere where the lowest eigenvalues are reproduced only approximately.
Notice also that, as stressed before and contrary to other discretizations of the sphere, the fuzzy Laplacian H 0 has From now on we will work with j = 1/2 so the dimension of the Hilbert space at each site is 4. A convenient basis for this space is: In a system with L > 1 spatial sites, the Hilbert space of the system is the tensor product of L single-site Hilbert spaces. The generic wavefunction can be written as The kinetic term H 0 = n H 0 (r) of the Hamiltonian is the sum of Eq. (5) acting on the Hilbert space of each site and it is diagonal in the basis |a L−1 , . . . , a 0 (for a single site operator A, A(r) ≡ 1 1 ⊗L−r−1 ⊗ A ⊗ 1 1 ⊗r−1 denotes the operator acting on site r.) In the T basis the kinetic term is represented by a sum of similar tensor products of 4 × 4 matrices with The eigenvalues of H 0 are then E 0 a L−1 ...a0 = g 2 r χ(a r ), with χ(a) = 1 − δ a0 , that is, the kinetic energy of site r is 0 if a r = 0 or g 2 otherwise.
The "interaction" term arises from expanding the nearest-neighbor interaction term −n(r + 1) · n(r) in Eq. (1): They involve only two neighbor sites at a time (one link). In the T basis the J k operators are represented by the following 4 × 4 matrices (j k ) ij ≡ T i |J k | T j : where 1 1 is the two-dimensional identity matrix and σ's are Pauli matrices. The interaction term H Ik (r + 1, r) in the T basis is the matrix h Ik (r + 1, r) = −(κ/g 2 ∆x 2 )(j k ) r+1 (j k ) r . By this we mean that the element a L−1 , . . . , a 0 H Ik a L−1 , . . . , a 0 is h Ik (r + 1, r) i,j where i and j are the numbers that have the representation in basis 4 a L−1 , . . . , a 0 and a L−1 , . . . , a 0 respectively (the matrix indices run from 0 to 4 L − 1.)

III. IMPLEMENTATION OF TIME EVOLUTION AND ESTIMATE OF RESOURCES
The implementation of the time evolution of the model in terms of quantum gates starts by splitting the time evolution over a number of smaller steps ∆t = t/N and each time step using the Suzuki-Trotter formula. Each time step is further split into the evolution due to the four parts of The state of the system is time-evolved by first applying the kinetic term e −i∆tH 0 site-by-site; the site order does not matter since the H 0 (r) for different r commute. We follow the kinetic term with the first interaction term e −i∆tH I1 . This evolution is done link-by-link, and again, the order does not matter, as all H Ik (r, r + 1) commute with each other. We use periodic boundary conditions, so the evolution of the link from r = L − 1 to r = L is followed by evolution of a link from r = L to r = 1. The link-by-link evolution is repeated for e −i∆tH I2 and e −i∆tH I3 . This completes one time step of the total time evolution. The whole process is then repeated N times. The local four-dimensional Hilbert space is encoded by two qubits so with the pair of qubits (for instance, q 1 q 0 ) being the binary digits of the value of the corresponding index a. For instance, q 1 = q 0 = 0 corresponds to a 0 = 0 and q 1 = q 0 = 1 corresponds to a 0 = 3. In this basis, the kinetic term evolution at each site corresponds to the matrix where from now on we set g 2 = 1 and κ/g 2 ∆x 2 = 1 for simplicity. The circuit implementing the kinetic term evolu-tion is depicted on the left of Fig. (1). Since HS † σ 2 SH = σ 3 the interacting term H I1 = 1 3 1 1 ⊗ σ 2 ⊗ 1 1 ⊗ σ 2 is related by a similarity transformation to 1 3 1 1 ⊗ σ 3 ⊗ 1 1 ⊗ σ 3 , with the change of basis given by single qubit operations (here h is the Hadamard and S the phase one-qubit gates.) Similarly, since Hσ 1 H = σ 3 , H I2 and H I3 are related to 1 3 σ 3 ⊗ σ 3 ⊗ σ 3 ⊗ σ 3 , via similarity transformations involving only single qubit operations. We now use the fact that a controlled-NOT (CNOT) gate implements a similarity transformation that takes σ 3 ⊗ σ 3 into σ 3 ⊗ 1 1 and that exp[iθσ 3 ⊗ 1 1] is simply a rotation on the left qubit. For H I1 we apply this for q 2 and q 0 qubit pair, whereas for the other two terms we have to apply the CNOT transformation on the (q 1 , q 0 ) and (q 3 , q 1 ) pairs to reduce σ 3 ⊗ σ 3 ⊗ σ 3 ⊗ σ 3 to σ 3 ⊗ 1 1 ⊗ σ 3 ⊗ 1 1, and then a CNOT tranformation on the (q 3 , q 1 ) pair to reduce the exponentiation to a single qubit rotation. The quantum circuits implementing these operations are depicted in Fig. (1).
The counting of quantum gates is, of course, dependent on the instruction set available on the hardware. However, the difficulty in the hardware implementation of two-qubit gates makes it unlikely that any two-qubit gate besides the CNOT gate will be available in hardware. CNOT gates are also, by far, the ones likely to generate decoherence. Thus, we will only count the number of CNOT gates required in our implementation. The kinetic term circuit has 2 CNOT gates per site (notice that a controlled-U 1 operation requires two CNOT gates to be implemented in usual architectures). The H I2 and H I3 link terms seems to require 6 CNOT gates per site. However, since we apply e iδtH I2 on every link the gates shown in the dashed box in Fig. (1) cancel between adjacent links and do not have to be applied. The result is that only 4 CNOT gates per link are required. A similar thing happens to the H I3 links. Finally, H I1 requires 2 CNOT gates, for a total of 12 CNOT gates per site (for periodic boundary conditions, where there are as many links as sites). Since T /∆t steps are needed for a total time evolution T , 36(T /∆t) CNOT gates are required to implement U (T ) in our three-site model. This gate depth renders the model inaccessible to the current generation of processors. As a proof-of-principle, however, we run our algorithm in a quantum computer simulator (QISKIT [37,38]) for L = 3 and the results are shown in Fig. (2). The results in the figure were obtained with ∆t = 0.2, ∆x = 1, g 2 = 1 and show, for illustrative purposes, the probabilities of finding, as a function of time step, the states |000000 , |000001 and |111111 starting with the initial state with equal amplitudes of all elements of the basis |q 5 q 4 q 3 q 2 q 1 q 0 obtained by applying a Hadamard transformation on the state |000000 . In the same figure we show the exact time evolution and "Trotterized" evolution by multiplying the appropriate 4 3 × 4 3 matrices. The error bars reflect the expected variance from quantum mechanical measurement.
It is important to stress that every step in our construction can be easily can be carried out in much bigger lattices and even in more spatial dimensions. The size of the blocks of time evolution -involving at most 4 qubits at a time -are independent of the system size. Also, the time evolution due to the kinetic term can be applied simultaneously to all sites and the evolution due to the hopping term can be applied simultaneously to half of the links at once. The method is essentially unchanged as the number of spatial dimensions is increased. Unfortunately, in the implementation in terms of quantum circuits we found, the number of CNOT gates is a little too large for current quantum computers available to us. Our attempts at running it on the IBM's ibmqx4 machine resulted in mostly noise.

IV. CONCLUSIONS AND PROSPECTS
One of the issues to be faced on the road to using quantum computers in quantum field theory is the pres-ence of bosonic fields. The Hilbert space of a bosonic theory has an infinite number of dimensions per spatial site, while that for a fermionic theory is finite dimensional. On the other hand, the Hilbert space describing a quantum computer with a finite number of registers is finite dimensional. Thus, even after discretizing space, some further truncation of the field space will be required [5,7]. We propose a method to accomplish this while, at the same time, preserving the O(3) symmetry of the theory.  1), are likely to approach the same continuum limit as ∆t, ∆x → 0. In fact, the continuum limit of the sigma model is obtained by tuning ∆t, ∆x → 0 and g 2 is such a way as to keep physical quantities (mass gap, scattering amplitudes) fixed in physical units (perturbation theory indicates that the model is asymptotically free so the correct scaling is g 2 ∼ −1/ log(∆x) [35]). In this limit, details of the Hamiltonian become irrelevant and any other Hamiltonian with the same field content and symmetries, on account of universality, give rise to the same continuum limit [39]. More precisely, any other operators, consistent with the O(3) symmetry, is of higher dimension and, presumably, irrelevant in the continuum limit. The reasonable assumption of universality can be checked in a classical calculation. The "Trotterized" time evolution operator (Eq. (9)) corresponds to an action discretized in both time and space. A Monte Carlo calculation using this action (analytically continued to imaginary time) can demonstrate whether the fuzzy model is indeed in the same universality class as the sigma model and has, therefore, the same ∆t, ∆x → 0 limit.
Among theories of physical significance, bosonic fields also appear in principal chiral models (as, for instance, in low energy QCD) and gauge theories. These bosonic fields take values on group manifolds. A slight modification of the scheme proposed in the present paper can also be used in these cases, but it is somewhat more involved. A full account of these extensions will appear separately.