Toward simulating Superstring/M-theory on a quantum computer

We present a novel framework for simulating matrix models on a quantum computer. Supersymmetric matrix models have natural applications to superstring/M-theory and gravitational physics, in an appropriate limit of parameters. Furthermore, for certain states in the Berenstein-Maldacena-Nastase (BMN) matrix model, several supersymmetric quantum field theories dual to superstring/M-theory can be realized on a quantum device. Our prescription consists of four steps: regularization of the Hilbert space, adiabatic state preparation, simulation of real-time dynamics, and measurements. Regularization is performed for the BMN matrix model with the introduction of energy cut-off via the truncation in the Fock space. We use the Wan-Kim algorithm for fast digital adiabatic state preparation to prepare the low-energy eigenstates of this model as well as thermofield double state. Then, we provide an explicit construction for simulating real-time dynamics utilizing techniques of block-encoding, qubitization, and quantum signal processing. Lastly, we present a set of measurements and experiments that can be carried out on a quantum computer to further our understanding of superstring/M-theory beyond analytic results.


Introduction
Quantum Field Theory (QFT) is the language of nature. In order to understand nature, we have to define QFT and solve it. In high energy physics, the lattice approach to QFT is a powerful conceptual and computational tool [1,2]. Supersymmetry (SUSY) is another important piece in modern theoretical physics. It may exist at a low-energy scale within reach by the LHC or next-generation particle accelerators, and at very least, it plays an important role in the holographic approach to quantum gravity. More specifically, via gauge/gravity duality [3], certain supersymmetric QFTs can give a nonperturbative formulation of superstring/M-theory. This gives us a strong motivation to define supersymmetric QFTs and solve them.
Typically, the traditional lattice QFT approach considers the Euclidean spacetime and uses the Markov Chain Monte Carlo simulation method. 1 Such approach is effective for various important problems, such as the derivation of spectrum of hadrons from QCD [6,7], determination of nuclear potential [8], and black hole thermodynamics in the holographic setup [9,10]. Still, there are many problems that cannot be accessed in this manner, most notably the real-time dynamics. Quantum simulation is a promising approach to such problems.
Despite a lot of effort and impressive progress, for many QFTs, the realizations on quantum computers remain challenging. The primary reason is that the lattice Hamiltonian is technically very complicated. Therefore, in this paper, we give an alternative approach, avoiding the use of a lattice. In order to explain the basic idea, let us recall a famous quote from Feynman -if you want to make a simulation of nature, you'd better make it quantum mechanical. An important idea implicit in this quote is that nature itself is a gigantic quantum computer, and systems that have natural physical realization can be simulated more easily. The lattice regularization is rather artificial. It may look natural to us humans, but it is safe to assume that nature is much smarter than us. Therefore, we should look for physical realization of QFTs in simpler quantum mechanical systems. In string theory, the system of D0branes and open strings can be described by quantum mechanics. D0-branes and open strings can have rich dynamics, and certain bound states in this system are equivalent to supersymmetric QFTs. By using such property, we can give physical realizations of those supersymmetric QFTs in a world described by certain quantum mechanics. As we can easily imagine, such physical realizations can be put on the quantum computer much more easily; essentially, we only have to put the quantum mechanics of D0-branes and open strings on a quantum computer. After that point, we only have to mimic what string/M-theory does. In some sense, it is similar to the Hamiltonian engineering for the analog quantum simulation; we realize a particular bound state of D0-branes and open strings, which is equivalent to the QFT we want to simulate, in a world described by the matrix model. This paper is organized as follows. In the rest of this section, we give several important remarks regarding the lattice regularization, in order to motivate the use of an alternative method introduced in this paper. In Section 2 and Section 3, we show how the matrix model and QFT can be realized in the Hamiltonian formulation. We introduce an explicit regularization scheme, such that it can be realized on a digital quantum computer with a large but finite number of qubits.
In Section 4, we present an explicit formalism for simulating regularized matrix model Hamiltonian on a quantum computer. We use the Wan-Kim fast adiabatic algorithm [11] to prepare ground state as well as thermofield double states on qubits. We then use block-encoding, qubitization, and quantum signal processing technique [12,13] to approximate the unitary real-time evolution. Lastly, we discuss several experiments and measurements that can be performed on this simulation of the matrix model to develop a deeper understanding of superstring/M-theory and holography beyond current analytic results.

Euclidean lattice: how it (sometimes) works on classical computers
When a quantum field theory in the Euclidean space is regularized on a lattice, it is important to keep the symmetries of the theory exactly at the regularized level. For example, Wilson's plaquette action [1] for the Yang-Mills theory preserves gauge symmetry, discrete translation (shift of one lattice unit), discrete rotation (90-degree rotation), parity, and charge conjugation. These exact symmetries control the radiative corrections such that the correct continuum limit is realized. Unless sufficiently many symmetries are preserved at the regularized level, the right continuum limit will not be obtained, because radiative corrections can break necessary symmetries. 2 One of the well-known cases is the chiral symmetry [14]. The Nielsen-Ninomiya no-go theorem [15] claims that the chiral symmetry cannot be preserved on the lattice with a few natural assumptions. For the vector-like models (i.e. left-handed and righthanded fermions appear in a pair, e.g., QCD), the overlap fermion [16] and domain-wall fermion [17] provide the ways to circumvent the Nielsen-Ninomiya theorem. However, a generic solution applicable to chiral gauge theories -including the standard model of particle physics, whose electroweak sector is chiral -is not known. Another wellknown case is supersymmetry on the lattice; because supersymmetry algebra contains the infinitesimal translation, which is broken on the lattice by definition, it is impossible to keep the entire supersymmetry algebra on the lattice. For (1 + 1)-and (2 + 1)dimensional theories, by keeping a part of supersymmetry, the right continuum limit can be taken [18][19][20][21][22][23][24][25][26]. But for (3 + 1) dimensions, no fine-tuning-free formulation is known. 3 The situation changes drastically for quantum mechanics, due to the lack of the ultraviolet divergence. Very often, naive regularizations that do not respect the symmetries lead to the correct continuum limit. Refs. [29,30] pointed out that this property is useful for the study of the supersymmetric matrix models. Such matrix models are important for quantum gravity via holography, and provided nontrivial tests of holographic duality at finite temperature and stringy level, which are out of reach with other approaches; see, e.g. Refs. [9,10,[31][32][33] for original references and Ref. [34] for a review.

Minkowski time and Hamiltonian formulation: why it is hard even on a quantum computer
The situation is similar in the Hamiltonian formulation. Let us consider the lattice Hamiltonian of pure Yang-Mills proposed by Kogut and Susskind [35]. It preserves various symmetries and hence leads to the desired continuum limit. Here, by continuum limit, we mean the continuum limit along the spatial dimensions; by definition, the time direction is continuous in the Hamiltonian formulation. On the other hand, for quantum mechanics, there is no need for the continuum limit in this sense, because there is no space by definition. When the Hamiltonian formulation is used on the digital quantum computer, yet another limit is needed: because we are considering the theories with bosonic degrees of freedom, whose Hilbert space is infinite-dimensional, the Hilbert space has to be truncated and expressed by using a finite number of qubits. There are two regularization parameters for QFT: the lattice spacing a and the dimension of the Hilbert space D = 2 nq , where n q is the number of qubits. The two-step limiting procedure is required: • For fixed lattice spacing a, we send n q → ∞, such that the correct lattice Hamiltonian acting on the infinite-dimensional Hilbert space is obtained; • Then we take the continuum limit, a → 0.
The first step is already nontrivial; for example, how can we express the unitary link variables by using qubits? We have to regularize the group manifold, which is doable but rather complicated; see, e.g., Ref. [36]. The second step is also highly nontrivial; actually, the situation can be worse than in the case of the Euclidean lattice, because now space and time are treated separately, and hence it is harder to keep large enough symmetry. For quantum mechanics, the second step is absent. The first step is also simplified because there is no dynamical gauge field. Gauge-singlet constraint is imposed on the states; the Hamiltonian only has to have the 'global' symmetry. In the class of theories we will consider, the dynamical variables are Hermitian matrices, which can be expressed as a collection of multiple harmonic oscillators interacting with each other in a certain manner. It allows us to use a simple truncation scheme of the Hilbert space.

QFT from Matrix Model: why it can work on quantum computer
The key idea we use in this paper is to embed the space to matrices, following concrete physical processes in string/M-theory. When we use lattice, the technical difficulty was that it is difficult to preserve sufficiently large symmetry at the discretized level. Hence we will borrow the idea from 'nature,' which is in this case string/M-theory. Regularizations utilizing actual physical setups in the matrix model, which have the origins in string/M-theory, respect necessary symmetries. Roughly speaking, the counterpart of the lattice volume is the matrix size N , and the detail of the lattice such as dimension corresponds to the choice of supersymmetric background in the matrix model. Hence the counterpart of the lattice-regularized QFT is the finite-N matrix model about a supersymmetric background.
Historically the first example of this kind of phenomenon is the Eguchi-Kawai equivalence [37]: at large N , gauge theory living at a point (essentially matrix model) is equivalent to infinite volume theory, if certain conditions are satisfied. The 'twisted' version of the Eguchi-Kawai equivalence [38,39] provides us with a natural way to embed the noncommutative space to matrices, which has a counterpart in string theory via the Myers effect [40]. Essentially three classes of quantum field theories can naturally be realized in matrix models: 4 supersymmetric gauge theories that naturally arise from string theory, gauge theories on noncommutative space, and large-N gauge theory. The counterpart of the continuum limit on lattice a → 0 is the appropriate large-N limit; the parameter-fine-tuning is not needed. The details of the construction of QFT from the matrix model will be explained in Sec. 3. Note that the limit of n q → ∞ is needed at each fixed N , but as we will see, this limit can be taken in a straightforward manner. See Sec. 2 for details.

Quantum Gravity in the Lab: matrix model on a quantum computer
A substantially large fraction of our motivation for the quantum simulation of the matrix model and supersymmetric QFT lies in quantum gravity. Quantum gravity in the lab [43,44] is a line of thinking which uses the holographic duality and experiments in the boundary QFT side to investigate hard problems in the bulk quantum gravity. For example, the quantum teleportation experiment can be used to test the existence of a wormhole [44]. Several other experiments have been carried out to test scrambling of quantum information and saturation of quantum chaos in physical models [45,46]. In addition to physically realizable models, quantum simulation on a digital quantum computer is also a possible avenue for progress. Recently, a quantum simulation proposal was made for the Sachdev-Ye-Kitaev (SYK) model [47,48] and sparse SYK model [49] that are good toy models for the holographic description of gravity [50].
In this paper, we take a similar approach and propose a framework for simulating a matrix quantum mechanics on a universal quantum computer. This theory is dual to superstring/M-theory and in the appropriate limit of parameters describes gravitational physics. Considering this duality, one can study various quantum gravitational phenomena in superstring/M-theory, for example, the presence of traversable wormholes, a saturation of quantum chaos, as well as the sub-AdS locality on this simulation models.
If the QFT side is sufficiently simple, it might be interesting also to address analog quantum simulation in the future. One could consider simulate matrix models in coldatomic physics, for instance, the Rydberg systems. Recently, there are some proposal of experimental implementations of some typical chaotic models, for instance Refs. [43,44,48,[51][52][53].

Regularizing the BMN matrix model
As we have briefly mentioned, we can realize various supersymmetric quantum field theories in terms of the matrix model. For example, we can realize 3d super Yang-Mills theory (SYM), 6d superconformal field theory (SCFT), and 4d SYM by taking suitable vacua in the BMN matrix model [54]. Hence giving an appropriate regularization scheme of the matrix model is a good starting point toward the realization of such QFTs on the quantum computer. In this section, we explain how to regularize the BMN matrix model in a suitable way to realize it on a quantum computer. In Sec. 3, we explain how QFT is realized by using the BMN matrix model.

Lagrangian formulation
Let us start with the Lagrangian formulation of the BMN matrix model [54]. All the dynamical variables of the BMN matrix model are N × N Hermitian matrices.
Although there are various ways 5 to construct the action of the BMN matrix model [54,58], here we skip their details and simply write the resulting Lagrangian. The theory is the U(N ) gauged supersymmetric matrix quantum mechanics whose field content can be interpreted as the dimensional reduction of (9 + 1)-dimensional super Yang-Mills theory: 6 • The gauge field A t .
• 16 adjoint fermions Ψ u with u = 1, · · · , 16. They are essentially the dimensional reduction of the Majorana-Weyl spinor in (9 + 1)-dimensions, which has 16 real degrees of freedom. We describe it as 16 complex fermions with a kind of reality condition called Majorana condition as explained later. Below we often write it as Ψ = (Ψ 1 , · · · , Ψ 16 ) T and do not explicitly write the indices.
Following the notation of [57], the Lagrangian is given by 7 where D t is a covariant derivative for the adjoint representation defined by The symbol ijk is the structure constant of SU(2), namely 123 = 231 = 312 = +1, 132 = 321 = 213 = −1. The matrices γ I are 16 × 16 real symmetric and traceless matrices satisfying 3) 5 1. Mass deformation of the Banks-Fishler-Shenker-Susskind (BFSS) matrix model [55] which was introduced as a matrix regularization of supermembrane in 11d [56]. Therefore the massless limit µ → 0 in (2.1) becomes the Lagrangian of the BFSS model. 2. Matrix regularization of supermembrane action in 11d pp-wave background [57]. 3. Dimensional reduction of the 4d N = 4 super Yang-Mills theory on R × S 3 along S 3 [58]. 6 The scalars can be interpreted as the dimensional reduction of the spatial component of the gauge field while the fermions are the dimensional reduction of the Majorana-Weyl spinor (gaugino) in (9 + 1)-dimensions, and hence, satisfy the Majorana condition. 7 A parameter R in Ref. [57] is related to our coupling constant by g 2 = R −3 .
which originally come from gamma matrices in higher dimensions 8 . There are various representations of the gamma matrices satisfying (2.3) but physical observables are independent of how to represent it. Therefore we can choose convenient representations depending on the problems under consideration.
In this paper, we take the following representation: where σ i are usual Pauli matrices, and g a are 4 × 4 matrices satisfying g a g b † + g b g a † = 2δ ab . As in the gamma matrices, there are various representations of g a but below, we do not use explicit forms 9 of g a . An advantage of this choice, particular for our purpose, is that the matrix γ 123 becomes diagonal and so does the fermion mass term Ψ † γ 123 Ψ. In our choice of the gamma matrices, the Majorana condition is written as The Majorana condition is solved as where I = 1, 2, 3, 4 and p = 1, 2. Thus, the Lagrangian in our choice of the gamma 8 In the (9 + 1)-dimensional perspective, the fermion originally had 32 complex components and we had ten 32 × 32 gamma matrices. The fermion is subject to Majorana and Weyl conditions, and finally has 16 real degrees of freedom. The Majorana condition for 32-component spinors isΨ ≡ Ψ † Γ 0 = Ψ T C, where C is charge conjugation matrix. The Weyl condition is the projection of the left-handed spinor. In the Lagrangian (2.1), the Weyl condition is already used while imposing the Majorana condition will be discussed later. 9 An explicit example is The BMN matrix model has various symmetries: • Translation symmetry along the time t.
• Supersymmetry: where is the 16 component Killing spinor satisfying • SO(3) global symmetry rotating the 3 scalars X i .
• SO(6) global symmetry rotating the 6 scalars X a .
Note also that the dimension of the parameters and fields are 10 Note that D t X I also transforms as Ω(D t X I )Ω −1 , and this fact makes the Lagrangian L invariant under the gauge transformation.

Hamiltonian formulation
Let us switch to the Hamiltonian formalism. We expand the matrices as where τ α is the generator of the U(N ) gauge group satisfying We take the temporal gauge A t = 0. The canonical conjugate momentum of X α I in this gauge is simply P α I = ∂ t X α I , while the conjugate of ψ is iψ † . Note that the conjugate momentum of the gauge field A t is zero since the Lagrangian does not contain ∂ t A t .
In the operator formalism, P I , X I , and ψ are promoted to the operators with the canonical (anti-)commutation relations: 14) The Hamiltonian is given bŷ Ip .
which satisfy the algebra whereM ij andM ab are the generators of SO(3) and SO(6) global symmetries defined byM respectively. It is important to note that the supercharge is gauge-invariant 20) and therefore supersymmetric condition is gauge-invariant.
In the anticommutation relation (2.18), we omitted a term proportional to the gauge generatorĜ. Such a term vanishes when acting on the gauge-singlet sector, but it can have a nontrivial consequence for the non-singlet sector. The Hamiltonian can be expressed without using the gauge generator asĤ = 1 32 {Q †Ip ,Q Ip }, and hence, the positive semi-definiteness holds in the non-singlet sector as well [59].

Regularization of the Hilbert space
The total Hilbert space of the BMN matrix model can be decomposed as where H X and H Ψ are subspaces associated with the scalars X I and fermions Ψ u , respectively. The dimension of H X is infinite because X I are bosonic, while H Ψ is finite-dimensional since it is associated with a finite number of fermions. Therefore we need to regularize H X in a certain way while we do not need a regularization for H Ψ . In this subsection we first regularize the subspace H X and then discuss the properties of the total Hilbert space after the regularization.

Fock basis of the bosonic part
First let us decompose the full Hamiltonian into the purely bosonic part and the other part: H =Ĥ X +Ĥ Ψ , (2.22) Here we focus on the purely bosonic partĤ X and further decompose it into free and interacting parts:Ĥ where ω I = for I = 4, 5, · · · , 9 (2.24) The free part is essentially 9N 2 harmonic oscillators described in the Fock basis, corresponding to α = 1, 2, · · · , N 2 and I = 1, 2, · · · , 9. The annihilation and creation operators are defined aŝ The free part of the Hamiltonian is essentially the number operator: Thus, the dimension of the full Hilbert space after the regularization is It is convenient to use the Fock basis regarding the fermion as well. We use the same notation |VAC free as before to denote the Fock vacuum both for the bosons and fermions, i.e.,Â α I |VAC free =ψ α Ip |VAC free = 0. (2.38) The minimum number of qubits needed for the regularization is 9N 2 log 2 Λ + 8N 2 . (Note that this is the number of logical qubits with a proper error correction.) In sec. 4, we will discuss how to realize the regularized theory in terms of qubits.

Sparseness of the Hamiltonian
In this section, we show that the matrix model Hamiltonian is very sparse for N, Λ 1 as long as potential is polynomial. In Sec. 4, we will see that such sparse Hamiltonian can be simulated efficiently.

Purely bosonic part
After the truncation (2.32), our Hamiltonian is Λ 9N 2 × Λ 9N 2 matrix and therefore has Λ 18N 2 elements. We are interested in how many elements are nonzero among them. We first discuss a rough estimate. For polynomial potential, 12 the Hamiltonian consists of a finite number of products of X Iα and P Iα whose number is polynomial in N . Noting that the operatorsX Iα =Â Iα +Â † , a finite product of X Iα and P Iα has only O(1) nonzero matrix elements in each row and column. Therefore, as long as we consider the polynomial potential, the number of nonzero matrix elements of the Hamiltonian is a polynomial of N among Λ 18N 2 elements. So the Hamiltonian in the Fock basis is very sparse. Now we make the estimate more precise. Let us start with the free partĤ free . From the expression (2.30), it is obvious that the free Hamiltonian has only 9N 2 nonzero elements. Therefore the free Hamiltonian is very sparse. Next, let us consider the interacting part. For our purpose, only the terms with the highest degree in the potential are relevant. Both for the BFSS and BMN matrix models, it is proportional to (2.39) We can easily see that each term has only 16 nonzero matrix elements at most in each row and column. So the problem is reduced to the counting of the number of the terms.
To do this, a little bit of group theory is needed. We normalize the generators of U(N ) as Tr (τ α τ β ) = δ αβ , Here α and β runs from 1 to N 2 . Instead of α, we can use p, q = 1, 2, · · · , N to label the generators, as τ ij They satisfy Tr(τ α τ β ) = δ αβ and α τ ij More precisely, we take the degree of the potential finite as N → ∞. 13 The results of this subsection are independent of the choice of the generators since it comes only from properties of the structure constants.
14 To derive the values of f αβγ , it is convenient to use N × N matrices M p,q whose (i, j)-component is (M p,q ) ij = δ pi δ qj . They satisfy M p,q M r,s = δ qr M p,s , and the U(N ) generators are written as The commutator of two generators can be expressed as From this, we can see that, among ∼ N 6 possible combinations of α, β and γ, only ∼ N 3 leads to nonzero f αβγ . The quartic interaction term is To see how many terms exist, we have to see how many combinations of α, β, γ and ρ give nonzero σ f αβσ f γρσ . It is ∼ N 4 , as one can check by using (2.41). Another way to understand that there are ∼ N 4 terms is to look at the definition of the trace, Tr(XŶẐŴ ) = N i,j,k,l=1X ijŶjkẐklŴli ; obviously, there are N 4 terms in the sum. Therefore, there are ∼ N 4 nonzero elements in each row and column.

Interaction between bosons and fermions
The interaction between fermion and boson leads to only ∼ N 3 terms at each row and column. Hence the leading contribution in the full Hamiltonian comes from Tr[X I , X J ] 2 , which gives ∼ N 4 terms at each row and column.

Remarks on cutoff dependence
We have truncated the Hilbert space by introducing the cutoff Λ in the harmonic oscillator basis. Needless to say, we have to take Λ sufficiently large in order to achieve a good approximation. Although the details of the cutoff effect depend on the details of problems under consideration, for reasonable physical setups of interest, our truncation prescription is valid, though it may or may not be optimal. For example, as a small perturbation about the ground states, we can imagine several incoming or outgoing objects (e.q. graviton or D-brane) described by the eigenvalues of matrices X I . If the eigenvalues take large values, then the free part of the Hamiltonian dominates the energies of those objects. As long as we consider the Hamiltonian time evolution, the energy is conserved, and hence those eigenvalues do not become arbitrarily large, which in turn means the arbitrary high-frequency modes in the Fock basis are not needed.
As we explain briefly in Appendix A, we could also use the coordinate basis for a regularization. The truncation errors in such a truncation scheme are discussed, e.g., in Refs. [60,61].

Remarks regarding gauge invariance
We are using the extended Hilbert space containing non-gauge-invariant states. The gauge transformation is generated byĜ α defined in eq. (2.16). Without the cutoff Λ, these generators commute with the Hamiltonian: Therefore, if an initial state |φ is gauge-invariant, the state remains gauge-invariant during the Hamiltonian time evolution One can construct gauge-invariant states as follows. From (2.16) and f αβγX β IP γ I = −if αβγÂ †β IÂ γ I , we can see that the vacuum of the free Hamiltonian is gauge-invariant: Then any state obtained by acting gauge-invariant operators such as Tr Â † Note that the gauge invariance can be broken at a finite cutoff due to the regularization effect. Technically this comes from the modification of the canonical commutation relation. For instance, neither the free Hamiltonian nor the interaction part commute with the gauge charge. Therefore the amount of breaking depends on how the states around the cutoff affect problems under consideration.
The gauge-invariant states span only a small fraction of the Hilbert space, and the unphysical, gauge-non-singlet states occupy the majority of the Hilbert space. Therefore, one might worry that accumulated simulation errors can spoil the gauge-singlet condition. One possible way to protect the gauge-singlet condition is to give a penalty to the non-singlet terms by adding a term proportional to αĜ 2 α to the Hamiltonian. In the case of the BMN matrix model, essentially the same thing might be happening automatically [59], as we will explain in Sec. 2.6. See also Ref. [62] which suggests that the non-singlet errors can be corrected rather easily, at least in the confining phase.

Maldacena-Milekhin proposal
Maldacena and Milekhin [59] argued that, at sufficiently large N and strong coupling where classical gravity is a good dual description, the non-singlet modes should be heavy and negligible. If it is true, once the initial condition is taken to be a singlet and the simulation is precise enough to keep the energy approximately constant, we do not have to do anything else to protect the singlet constraint.
This proposal can be tested on classical computers, by the lattice Monte Carlo simulation. Ref. [63] performed lattice simulation of the BFSS matrix model (µ → 0 limit of the BMN matrix model) and obtained the numerical results supporting the proposal.

QFT from the BMN matrix model
It is known that the BMN matrix model has a class of supersymmetric configurations called fuzzy spheres which preserves all the supersymmetries. By expanding the theory around fuzzy sphere backgrounds appropriately, we can obtain various supersymmetric QFTs which have origins from superstring/M-theory and are expected to be dual to some semiclassical gravities when the theories have large enough degrees of freedom. In this section, we give a brief review of the fuzzy sphere and discuss how it is realized on quantum computers. For trivial configurations for the fields except for X i , we can write the Lagrangian as This implies that we have the following class of classical solutions where J i is the generator of SU(2) algebra in N -dimensional representation (not necessarily irreducible), satisfying If we regard X i as a position, then the solution describes a space whose coordinates do not commute each other but obey the constraint X 2 i = (µ/3g) 2 Tr(J 2 i ) as in sphere. Therefore this solution describes a non-commutative version of the sphere and it is often called the fuzzy sphere solution. In addition, substituting the soliton (3.2) into the supersymmetric transformation (2.9), one can show that the fuzzy sphere solution preserves all the supersymmetries.
Note that the fuzzy sphere solution (3.2) is not invariant under generic gauge transformation, unless the representation of J i is N copies of one dimensional representation of SU(2) i.e., J i = 0 (we call this case "trivial"). For the generic representation of J i , the fuzzy sphere solution is invariant under a subset of U (N ) gauge transformation such that [ζ, J i ] = 0.

Fuzzy sphere in the Hilbert space
Let us see how the quantum states corresponding to the fuzzy sphere solutions (3.2) of the classical equations of motion, which preserves all supersymmetries, can be constructed in the Hamiltonian formalism. (A few ways to explicitly construct such states on quantum devices will be explained in Sec. 4.) Denoting such a state by |J i , the state has the following properties: • Invariance under supersymmetrŷ Although it is hard to construct such states analytically for strong coupling, we can explicitly construct it for weak coupling as follows. First, for the trivial case J i = 0, it is simply the free vacuum: Indeed one can show that this state satisfies the SUSY condition (3.4) up to O(g). In addition, this state is gauge-invariant and obviously satisfies (3.5).
In order to construct generic fuzzy-sphere states, it is convenient to redefine the operatorX i as as [64] Equivalently, we can use the translation operator e Then the Hamiltonian can be written as 15 H =Ĥ Note that the total HamiltonianĤ is unchanged and we have just decomposed the Hamiltonian in a different way from the trivial fuzzy sphere case. The "free part"Ĥ is quadratic in all the operators and therefore, we can rewrite it as a collection of simple harmonic oscillators after appropriate diagonalization of the mass terms. 16 Then we can construct the Fock space associated with the harmonic oscillators and and its Fock vacuum |J i ; VAC free is the fuzzy-sphere state in the weak-coupling limit: In Sec. 4, we will explain how the ground state of the Hamiltonian including the interaction term can be obtained on quantum device. Strictly speaking, the quadratic HamiltonianĤ free has zero modes, which corresponds to the gauge transformation of the fuzzy sphere, or in other words, the Nambu-Goldstone (NG) modes associated with the spontaneous breaking of the U(N ) symmetry due to the choice of a particular fuzzy-sphere configuration. Let us look into this point. 17 From the state |J i , we can obtain a gauge-invariant state simply by a 'symmetrization': where int depends explicitly on the representation of SU(2), and they are different from the 'free' and 'interaction' parts defined in (2.23). 16 To do this, it is most appropriate to expand the operators in terms of fuzzy sphere harmonics which is a non-commutative version of spherical harmonics (see, e.g., [65] and references therein).
Here we do not explicitly write their details. 17 When we discuss the realization of |J i on the quantum device in Sec. 4, we will give a prescription to control this flat direction.
Here the integral is over Haar measure of the U (N ) gauge group and N 1/2 is a normalization constant. The symmetrization commutes with the Hamiltonian since the Hamiltonian commutes withĜ α . Therefore, S(|J i ) is also the ground state. By construction, this state is gauge-invariant and hence satisfies the Gauss-law constraint. This can be seen explicitly by using the invariance of the Haar measure: Hence S(|J i ) is the gauge-invariant ground state. Now we have two states, |J i and S|J i . Apparently, |J i is much simpler when the realization on the quantum device is concerned. Hence, if possible, we would want to use |J i for the simulation; but does it make sense to use |J i ? A conjecture by Maldacena and Milekhin [59] has an interesting consequence regarding this point. According to their conjecture, at sufficiently large N and strong coupling, the energy of the non-singlet modes can be bounded from below. Below this lower bound, only the singlet modes can contribute to the dynamics, and hence, it does not matter whether we use |J i or S|J i . In some interesting cases (including the strongly-coupled QFT limits we consider below in Sec. 3.2), this lower bound is rather high, and we can study interesting quantum gravity problems by using |J i .
Even if the Maldacena-Milekhin proposal is correct, |J i and S|J i can give different results in certain parameter regions of interest. Therefore, preferably we should consider the singlet state S(|J i ). We will provide possible protocol for preparing the state S(|J i ) in Sec. 4.4. We will comment on the Maldacena-Milekhin proposal further in Sec. 3.4.

Mapping rule
We review the correspondence between the fuzzy-sphere states and QFTs.

3d maximally supersymmetric Yang-Mills theory (M2-branes and D2-branes)
The 3d maximally supersymmetric Yang-Mills theory (maximal SYM) describes the worldvolume theory of D2-branes in type IIA superstring [66,67]. The fuzzy sphere in the BMN matrix model can be regarded as a D2-brane made of D0-branes via the Myers effect [40]. Hence the 3d maximal SYM can be realized by taking (an appropriate) fuzzy-sphere configuration. At low energy, the same theory is expected to describe the M2-branes.
The precise construction is as follows [68]. Let us consider N 2 fuzzy spheres with spin s. The size of each fuzzy sphere is N 5 ≡ 2s + 1, and the size of the gauge group of the matrix model is N = N 2 N 5 . Hence the SU(2)-generators in the fuzzy sphere state take the following form: (3.14) Here J (s) i are the generators in the spin-s representation. This configuration can be interpreted as the N 2 -coincident spherical D2-branes. If we focus on the fluctuation about this configuration, then we can obtain the 3d U (N 2 ) maximal SYM on fuzzy sphere 18 . The ultraviolet momentum cutoff is proportional to µN 5 . The coupling constant of the 3d theory g 3d and the noncommutativity parameter 19 θ are related to g 1d and µ by 20 This is the standard relation between the matrix model and noncommutative space, which was first pointed out in Refs. [38,39] for the case of the fuzzy torus, in the context of the Eguchi-Kawai reduction [37]. In order to obtain the usual 3d SYM on noncompact and commutative space, we take three limits, namely • The continuum limit µN 5 → ∞, • Decompactification limit (flat limit) r S 2 = 3 µ → ∞, • Commutative limit θ → 0, while fixing the coupling constant and gauge group of the 3d theory, g 3d and N 2 . Note that one does not have to take the decompactification limit and the commutative limit, if one is interested in the theory at finite volume or with finite noncommutativity. 18 One can also realize the maximal 3d SYM in monopole background by choosing the representation of J i in a different way but here we have turned off the monopole background for simplicity of explanations. See e.g. [65] for details. 19 By using the standard embedding of the sphere to R 3 and zooming in the neighborhood of the north pole, we get two-dimensional noncommutative plane parametrized by x and y, with the noncommutative product x * y − y * x = iθ. The coordinate and matrices are related as follows. Firstly, gX i is identified with the gauge-covariant derivative. By writing gX i = µ 3 J i + ga i , the first term µ 3 J i is regarded as the momentum p i , or equivalently the derivative −∂ i , and the second term is regarded as the gauge field. At the north pole ( The coordinate x and y are defined as x = −θp 2 and x = θp 1 , such that [x, p 1 ] = [y, p 2 ] = i. 20 Here g 1d = g. We have added the subscript "1d" to emphasize that it is the gauge coupling in one dimension.
In this construction, the existence of maximal supersymmetry is crucial: in general, the commutative limit of QFT on noncommutative space is not the same as the corresponding theory on the commutative space, due to the UV/IR mixing [69], but this problem is absent for the theories with maximal supersymmetry [70,71].
The coupling constant g 2 3d has the dimension of mass. At a given energy scale ε, the effective, dimensionless coupling is g 2 3d /ε. M-theory and IIA string description is valid at g 2 3d /ε N −1/5 2 and N −1 , respectively [67]. If we are interested in the low-energy spectrum, the characteristic energy scale is the inverse of the radius of S 2 : ε ∼ µ. In the M-theory regime, the dual gravity description is the N 2 -coincident spherical M2-branes in the pp-wave geometry with the flux parameter µ, with tension The 6d N = (2, 0) SCFT is the worldvolume theory of M5-branes in M-theory. It is also expected to be dual to M-theory on AdS 7 × S 4 . Ref. [68] proposed that various nontrivial brane configurations can be described by the fuzzy-sphere vacua of the BMN matrix model. Among them, the M5-brane configuration is obtained in the following scaling limit in (3.14) [68,72,73]: The dual gravity description is the . The scalars X I in the matrix model with the canonical normalization (2.15), whose dimension is (mass) −1/2 , is obtained by multiplying µ 6N 2 to the coordinate in the gravity side: This parametrically large S 5 should be realized as the spherical distribution of the elements of six scalars X 4 , X 5 , · · · , X 9 .

4d N = 4 SYM in the large-N limit
The 4d N = 4 SYM describes the worldvolume theory of D3-branes in type IIB superstring. It also provides the canonical example of the AdS/CFT correspondence with type II superstring on AdS 5 ×S 5 . Around a certain concentric-fuzzy-sphere background, the BMN matrix model can describe the 4d N = 4 SYM on S 3 in the planar limit [65]. (Note that, unlike the M2/D2 and M5 theories, this construction is valid only in the planar sector.) The SU(2)-generators are chosen to be the following form 21 : 18) and the large-N limit is taken as k, n, T, n − T → ∞, λ 4d = 8π 2 g 2 1d k µn : fixed. (3.19) Note that the sum with respect to s is taken over the integer and half-integer. The radius of S 3 is given by 6/µ. An intuitive way to understand this construction is to regard S 3 as the S 1 -fibration of S 2 ; for each s, the 3d maximal SYM on S 2 is obtained as in Sec. 3.2.1, and the additional S 1 direction is generated via the summation over s interpreted as summing up KK momenta, in a way similar to the quenched Eguchi-Kawai reduction [82][83][84]. In this case, note that any decompactication limit is not required even if we are interested in flat noncompact space since R × S 3 is conformally equivalent to the flat space.

Comments on QFTs without string theory dual
It is straightforward to construct other matrix models, which admit the fuzzy sphere background; see, e.g., Ref. [85]. By using them, 3d QFTs on the noncommutative space can be regularized. Without supersymmetry, the fuzzy sphere background is usually unstable in the QFT limit, and hence, only supersymmetric theories admit simple constructions [86][87][88][89]. 22 Note that the commutative limit (θ → 0) of such noncommutative theory is different from the theory on the commutative space in general. In the opposite limit (θ → ∞), only the planar diagrams survive, and the large-N limit of the commutative theory is realized [38,39]. The construction of the large-N limit of the 4d N = 4 SYM [65,75]

Comments on the cutoff effect
In practical simulations, we shall approximate the fuzzy sphere state by imposing the truncation on the Hilbert space and the approximations depend on the cutoff Λ. Although it is difficult to estimate the cutoff effect precisely, one can estimate lower bound on Λ to get reasonable approximations by a simple argument below. Note that the bounds estimated below are very minimal requirement, and some amount of error can remain as long as Λ is large but finite. In this sense, our estimation is crude. Let us give a crude estimate of the cutoff effect in the fuzzy-sphere states. Here we only discuss the necessary conditions for the true ground state to be approximated by a certain linear combination of the Fock states in the truncated Hilbert space. In order for the ground state of the truncated Hamiltonian to be close to the true ground state, additional conditions should be required.
Let us first consider the non-gauge-invariant state |J i and take the gauge in which J 3 is diagonal. Then the diagonal elements of J 3 runs from −s to +s. Hence, depending on the value of the coupling g 1d , the flux parameter µ and the spin s, the fuzzy sphere can be rather large in the coordinate basis. Similar properties hold for other gauge choices as well. Therefore, in order to express the wave function properly, we have to be able to describe the functions which spreads from X α i ∼ − sµ g 1d to X α i ∼ + sµ g 1d , for all α = 1, 2, · · · , N 2 − 1. Similar consideration is needed for X α a and P α I as well. Then higher modes in the Fock basis, whose wave function in the coordinate basis is larger, are needed; the n-th excited mode of the harmonic basis is localized at X n µ , and linear combination of the modes at n < Λ can describe only the states localized at X Λ µ . Hence the following conditions are needed: • The width of the wave function of the n-th excited state is n|X 2 |n ∼ n µ . (Note that n|X|n = 0.) The size of the wave function of the Λ-th excited state ∼ Λ µ has to be larger than the typical value of X α I in the true ground state defined through the expectation value ofX 2 .
• Typical value of the momentum of the n-th excited state is n|P 2 |n ∼ √ nµ.
(Note that n|P |n = 0.) Hence √ Λµ has to be larger than the typical value of P α I .
In addition to these, let us require the following: • If we want some resolution to distinguish the possible values of X α I , then √ Λµ has to be smaller than that resolution scale. In order for the fuzzy sphere background to make sense, the resolution has to be finer than the difference between the values of the matrix entries of the fuzzy sphere, µ 3g 1d . Below, we consider the N 2 -coincident fuzzy sphere background (3.14). The spin is s = 1 2 (N 5 − 1), where N = N 2 N 5 . The largest value among the matrix entries We emphasize again that the three conditions listed above are very crude, and the bound obtained here is weaker than the requirement for a stronger condition that the ground state of the truncated Hamiltonian is close to the true ground state. Different realization of the fuzzy sphere background, which are related by the gauge transformation, can lead to different estimate. However it is not necessarily a bad news, in the following sense. Let us consider all possible realizations of the fuzzy sphere background and take the strongest bound; we use Λ crude,max to denote this. The requirement that the ground state of the truncated Hilbert space is close to the true ground state is a gauge-invariant notion, so there must be the value of cutoff appropriate for any realization of the fuzzy sphere; we call it Λ strict . Then Λ crude,max ≤ Λ strict .

Weak coupling
Eventually we want to consider the gauge-invariant state S(|J i ), but we start with a simpler case, the non-gauge-invariant state |J i . Furthermore take the specific realization of SU(2) in which J 3 is diagonal. The first condition for X 1,2,3 becomes Λ (3.20) Here we have used the relation (3.15), which relates the couplings of the 1d and 3d theories. The other two conditions are much weaker than (3.20). Note that the offdiagonal modes are heavy about this background. Therefore, it is likely that (3.20) can be a good estimate for the stronger requirement: the ground state of the truncated Hilbert space is close to the true ground state. As mentioned before, this stronger requirement is a gauge-invariant notion, and hence, we expect (3.20) gives a good estimate for the gauge-invariant state S(|J i ) as well (i.e., Λ crude and Λ strict can be roughly the same.) The wave function looks different depending on the embeddings of SU(2) to SU(N ). For example, we can perform a generic gauge transformation and make all the matrix entries to of order µ g 1d ×N 0 . If we use such embedding, much finer resolution (∼ µ g 1d N 5 ) is needed in order to distinguish the fuzzy sphere and other backgrounds. Hence the third condition matters and 1 √ Λµ µ g 1d N 5 is needed, which in turn becomes Λ g 2 1d N 2 5 µ 3 . This is weaker than (3.20); however note that, with such gauge choice, off-diagonal modes are not heavy and the expansion about this specific background cannot be truncated at low order; the condition (3.20) is required in order for the ground state of the truncated Hamiltonian to be sufficiently close to the true ground state.

M2-brane limit
At strong coupling, the M2-brane limit (N 2 , µ fixed, N 5 → ∞): By using (3.15) which relates the couplings of the 1d and 3d theories, we obtain µs 3g 1d ∼ √ N 5 6g 3d . It is reasonable to assume that the fluctuation is smaller than this classical value. (Otherwise, we cannot take the M2-brane limit.) Then the requirement is formally the same as the weak-coupling result (3.20), namely Λ µ µN 5 g 3d is required. Other conditions are much weaker.

M5-brane limit
According to Ref. [73], the eigenvalues of X 4,··· ,9 form S 5 , whose radius scales as Suppose we took the fuzzy sphere state |J i , which is not gauge-invariant. This choice appears to break the SU(N ) symmetry, and hence, we would expect the existence of the Nambu-Goldstone modes, which connect different realization of the fuzzy sphere. According to the Maldacena-Milekhin conjecture [59], it is not necessarily the case. In the matrix model, the spontaneous symmetry breaking can take place only in the strict large-N limit. In the M2/D2 or M5/NS5 limit, the fuzzy sphere configuration and the coupling constant are varied nontrivially with N , and hence, whether the Nambu-Goldstone mode can actually appear is a highly nontrivial issue. 23 According to the 23 We thank A. Milekhin for useful comments regarding this point.
Maldacena-Milekhin conjecture, the energy of the would-be Nambu-Goldstone modes are bounded from below, as If we are interested in the parameter region dual to weakly-coupled type IIA string or M-theory, we need to consider the energy scale much lower than λ 3d . There, such modes are negligible, and a sort of super-selection takes place. In the same manner, in the scaling region where the M5/NS5 brane theory appears (fixed N 5 , N 2 → ∞ with g 2 1d N 2 → ∞), the would-be Nambu-Goldstone modes are negligible. If the Maldacena-Milekhin conjecture is correct, then we do not have to prepare the gauge-invariant fuzzy sphere state for quantum simulation at strong coupling. This can simplify the state preparation significantly.

Realizations on quantum computer
So far, we have explained how certain supersymmetric quantum field theories can be formulated by using the matrix models, in such a way that they can be put on a universal quantum computer. We have not yet specified the detail of the implementation. In this section, we will give concrete protocols and show the potential benefit of using quantum simulation.

Encoding the bosonic part into qubits (compact mapping)
Now we wish to address the generic strategies we usually use to encode the Hamiltonian (see [91,92] for references). Firstly we focus on the bosonic part.
As we have seen, we regularize the Fock space by introducing the cutoff Λ to the excited modes of harmonic oscillators. For each harmonic oscillator, we assign K = log 2 Λ qubits. We use the compact mapping for the energy level j = 0, 1, · · · , Λ − 1, where we use the binary decomposition One could map the creation operator A † aŝ By using |j = |b K−1 |b K−2 . . . |b 0 and |j + 1 = b K−1 b K−2 . . . |b 0 , we can write |j + 1 j| as Note that each |b l b l | is just a Pauli spin operator: Therefore,Â † can be written as a linear combination of less than Λ 2 Pauli strings of at most length K (i.e., tensor products of K Pauli spin operators). 24 The same holds forÂ. InX α I =Â +Â † √ 2ω I , the non-Hermitian part inÂ andÂ † (odd number of iσ y ) cancel out and the number of Pauli strings becomes smaller compared toÂ andÂ † .
In our simulation problem, we prefer to use the compact mapping approach, since it has much fewer costs about Hilbert space, although it does not have a good locality. Still, however,X α I andP α I are K-local, where K = log 2 Λ and Λ is at most some powers of N in the situations under consideration; see Sec. 3.3. Therefore, the growth of K is slow (K ∼ log N ) and hence the lack of the locality may not be too problematic.

Encoding the fermionic part into qubits
The Jordan-Wigner transformation is a standard way to relate fermions and qubits. In order to expressΨ explicitly by using the Jordan-Wigner transformation, let us use a standard representation with the complex fermions shown in Sec. 2.2. The 16component Majorana fermion Ψ can be expressed by using four (I = 1, 2, 3, 4) 2component (p = 1, 2) complex fermions ψ Ip which satisfy the anticommutation relation {ψ †Ipα ,ψ β Jq } = δ I J δ p q δ αβ . To simplify the notation, we combine three indices to one index which runs from 1 to 8N 2 : {ψ † m ,ψ n } = δ mn (m, n = 1, 2, · · · , 8N 2 ). Then we can expressψ m acting on H Ψ aŝ This encoding is simple, but has a disadvantage: the operatorsψ is highly nonlocal. In the BMN matrix model, the maximum length of the Pauli strings appearing 24 Each |b l b l | contains at most two Pauli matrices, and l runs through K different values in each |j + 1 j|, so there are 2 K = Λ or less Pauli strings in each |j + 1 j|, and the maximum length is K. There are Λ − 1 different values of j, so the number of Pauli strings can be bounded by Λ(Λ − 1). More or less, the same bound can be obtained by noticing that the number of possible Pauli strings of length K or less, including the identity, is 4 K = Λ 2 .
in the realization ofψ is 8N 2 . Such long Pauli strings make the quantum simulation inefficient.
Alternatively, we can use the Bravyi-Kitaev transformation [93], which provides us with a good basis, which makes the length of the Pauli strings to scale log N . (For explicit form in terms of Pauli matrices, see, e.g., Ref. [94].) Whether the Bravyi-Kitaev basis is better than the Jordan-Wigner can depend on the value of N . At very large N , the Bravyi-Kitaev basis performs better.

Real-time dynamics: qubitization, quantum signal processing
For the real-time evolution, we consider the time-independent Hamiltonian. While we could apply various algorithms listed in Appendix F, here we consider another efficient algorithm, taking advantage of the encoding of matrix models. This algorithm, used in the context of quantum signal processing, has an explicit construction of oracles. Here, we will review an algorithm that is related to Refs. [13,47].
The Hamiltonians of the matrix model is written in terms of Pauli strings, The first step in the quantum simulation algorithm is the block-encoding of our Hamiltonian into a unitary matrix. We introduce ancilla states |i (i = 1, 2, · · · , L) and prepare a state |G defined by Then we can prepare a unitary operatorÛ acting on the Hilbert space times C L , which satisfiesĤ HereÎ is the identity operator acting on the Hilbert space. Such unitary operatorÛ can be constructed as a control-Π i operation for K-local Pauli stringΠ i , U |i |ψ = |i Π i |ψ . and by using it, we defineŴ =RÛ . where T n (·) stands for the n-th order Chebyshev polynomial of the first kind. (See the appendix E.1 for a proof of (4.14).) We are interested in approximating the time evolution operator e −iĤt . In fact, the approximation can be made with the help of the following expansion (Jacobi-Anger expansion) of the time evolution operator where J n is the Bessel function of the first kind. The right hand side can be handled by the quantum signal processing algorithm [13], namely quantum signal processing can be used as a black box to compute the Chebyshev polynomials efficiently when approximating e −iĤt up to an error for repetition n ∼ λt + log(1/ ) [13]. (See the appendix E for a review of quantum signal processing.) Thus, the total complexity of simulation e −iĤt is . 25 The value of ||α|| is dominated by the quartic interaction, which has O(N 4 ) combinations of color degrees of freedom. There is an overall factor g 2 , and the sum over the excitation levels gives a factor of the order j1,j2,j3,j4 √ j 1 j 2 j 3 j 4 ∼ Λ 3/2 4 = Λ 6 . Therefore, ||α|| ∼ g 2 N 4 Λ 6 . The cubic interaction gives a sub-leading correction of order µgN 3 Λ 9/2 . For more detailed discussions, we recommend the readers to read the appendix E and the original references [13,47]. We could regard this as an alternative algorithm for the real-time evolution of a fixed Hamiltonian, especially when we need to know the precise construction of the oracles. It will be interesting if one could also construct a time-dependent version of the algorithm.

Adiabatic state preparation
The simulation algorithms for the real-time dynamics have to be accompanied with the preparation methods for appropriate initial states. Furthermore the states themselves have rich information about the system. Hence, we will explain how the ground state can be constructed by using the adiabatic state preparation method. In this section, we present a novel application of the Wan-Kim algorithm [11] for fast digital state preparation that utilizes block encoding and the theoretical concept of quasi-adiabatic continuation. The interested reader can also look into alternative algorithms in Appendices D and D.2 that are based on time-dependent Hamiltonian simulation methods in Ref. [95].

Trivial vacuum
First, let us consider the trivial vacuum. In the weak-coupling limit, this vacuum is literally 'trivial': it is just the Fock vacuum, which is 'N fuzzy spheres with spin zero.' In the M5-brane limit (Sec. 3.2.2), it is expected that one M5-brane is described [68]. The nature of the wave function at strong coupling is expected to be not trivial at all, and we could see some crucial properties of the trivial vacuum by quantum simulation.
The full Hamiltonian of the BMN matrix model consists of the free part and the interaction part,Ĥ BMN =Ĥ free +Ĥ int (g), whereĤ int disappears when the coupling constant g is set to zero. The free partĤ free does not contain g. In the weak-coupling limit g → 0, the ground state is the Fock vacuum. Using the notation of Ref. [11] we introduce a notationĤ 0 =Ĥ free , (4.17) and use a parameter s ∈ [0, 1] to interpolateĤ 0 andĤ 1 aŝ To use the Wan-Kim algorithm, first we need to construct the block encoding ofĤ 0 , where β ≈ max{||Ĥ 0 ||, ||Ĥ 1 ||} and β ≈ ||Ĥ ||. For the BMN matrix model, β and β are also dominated by quartic interaction term (analogous to |α| in Eq. (4.16)) and thus, β, β ∼ g 2 N 4 Λ 6 . The cubic interaction gives a sub-leading correction of order µgN 3 Λ 9/2 .
In Appendix D, we discuss an alternative, more straightforward algorithm based on quantum simulation technique for time-dependent Hamiltonian evolution and the adiabatic theorem. The advantage of the Wan-Kim's fast digital algorithm is that it gives a poly-log scaling in 1/δ error, as opposed to polynomial scaling for the algorithm described in Appendix D.

Fuzzy-sphere vacuum
Next we consider the fuzzy-sphere states discussed in Sec. 3.1.2. The starting point is |J i ; VAC free , which is the ground state ofĤ free is quadratic, it is easy to determine the ground state in terms ofP I ,Ŷ i ,X a and the representation J i . But there is one issue, which is problematic in our specific regularization with the Fock basis:Ĥ (J i ) free has the flat direction, along which the ground state is zero momentum state which is not well described by using the truncated Fock space. Practically, we can lift the flat direction by adding a mass term proportional to i,α,β [Ĝ α ,Ŷ β i ] 2 . Let this modification be 'gauge fixing' term,Ĥ g.f. . ThenĤ free +Ĥ g.f. does not have flat directions and the ground state can be expressed as the Fock vacuum of this quadratic Hamiltonian. We will use this state, which we denote by |J i ; VAC g.f. , for the state preparation. In Appendix B we will show how |J i ; VAC g.f. can be expressed in the original basis used for the regularization. To obtain |J i , we perform the adiabatic state preparation by taking |J i ; VAC g.f. to be the initial state and Now we consider the gauge-invariant state S(|J i ). In order to eliminate potential NG modes (which can be lifted if the Maldacena-Milekhin conjecture is true), we can force the gauge-singlet constraint by adding a term proportional to αĜ 2 α to the BMN HamiltonianĤ. Then only the gauge-invariant fuzzy-sphere states can be the ground states. Hence it looks reasonable to useĤ 0 =Ĥ with a positive coefficient c.
Another possible option would be as follows. Once we prepare a gauge-invariant state close to S (|J i ) at weak coupling, the standard adiabatic state preparation can be used to go to strong coupling. It might be possible to obtain such state by starting with the Fock vacuum, which is gauge-invariant, and perform certain adiabatic state preparation compatible with the gauge invariance to obtain a gauge-invariant state close to S (|J i ), for example, by interpolatingĤ free and 26 for coincident fuzzy spheres with spin s. We can repeat the adiabatic state preparation once more to get the fuzzy-sphere state in the matrix model, by interpolating this Hamiltonian and the BMN Hamiltonian. In order to make sure that the gauge invariance is preserved, we can add a term like αĜ 2 α to the Hamiltonian, or we could keep doing measurement usingĜ α such that the charge will thus approximately be preserved during the time evolution.
The adiabatic state preparation can work when there ground states at s = 0 and s = 1 are smoothly connected, without a level-crossing. Further investigation is needed in order to check this property.

Thermofield double state (TFD)
Another important state is the thermofield double state, where we have two copies (left and right) of the same system with total Hamiltonian H =Ĥ L ⊗Î R +Î L ⊗Ĥ R , withĤ L =Ĥ R that have eigenvalues E i . There are numerous recent studies for CFT, SYK model, and generic chaotic systems [96][97][98][99] that show that, with good accuracy, TFD is realized as the ground state of the two independent copies of the system that are coupled aŝ The coupling term can be quite generic, and the only constraint is that operatorsÔ are local. The strength of the coupling g LR controls the temperature (1/β) of the TFD state. Analogous to the discussion in the previous section, we can construct oracles for Hamiltonian with left and right coupling and without it. 27 And again we can make use of Algorithm 1 and Theorem 1 of Ref. [11] with initial HamiltonianĤ L ⊗Î R +Î L ⊗Ĥ R , which has a simple ground state that is the tensor product of two ground states that we already know how to construct operationally. The TFD state, on the other hand, is the ground state of the target Hamiltonian Eq. (4.29) assuming the results of Refs. [96][97][98][99] will carry through to the matrix model Hamiltonian. The state preparation complexity of the thermofield double will be O( KLβ 2

Measuring quantum black holes on the quantum devices
As we have mentioned before, we could construct the 'trivial' vacuum and fuzzy sphere vacua in the BMN matrix model. The shape of the wave functions describing these vacua is already a highly nontrivial and interesting target of the quantum simulation. However, those are just the tip of the iceberg: quantum simulation could do much more than constructing those states.
In fact, one of our dreams is to probe the black hole dynamics using the quantum simulation of matrix models. For instance, let us consider physics near the trivial vacuum of the BMN matrix model, with sufficiently small µ, i.e., close to the BFSS matrix model. Suppose the ground state is prepared by using the adiabatic state preparation method. By performing a unitary transformation close to the identity, we can add a small amount of energy to the system. Then we can follow the unitary time evolution and see how the system thermalizes. We expect that the system thermalizes toward the Schwarzschild black hole in M-theory or black zero-brane in type IIA superstring theory, depending on the energy added to the system [67]. Because the quantum simulation allows us to access the quantum state, it might be possible to see how the black hole geometry is realized by the matrix degrees of freedom. One natural possibility is that the Schwarzschild black hole is realized as the partially-deconfined phase [100][101][102][103][104]. This possibility may be testable to some extent with the classical simulation as well, as demonstrated in Ref. [105] in a simpler matrix model, and hence, may serve as a benchmark for the power of the quantum simulation. Depending on the choice of µ and energy, the black hole can be unstable. It can be a resonance which eventually evaporates by emitting the Hawking radiation. It is extremely important to see such formation and evaporation of quantum black hole based on the first principle.
It is also interesting to consider QFT embedded in the matrix model and study the problems involving holographic scattering and bulk locality. Namely, when the semiclassical gravity dual ('bulk geometry') is expected, it is possible to perform the scattering experiments by shooting some excitations from the boundary towards the bulk. Revealing details of such experiments might tell us the existence of bulk locality, especially at the sub-AdS scale [106] (see also Ref. [107,108]). It might also be useful for the construction of the tensor network toy models of holography in the precision of sub-AdS [109,110].
Finally, we wish to mention more concrete proposals. The first example is a motion of D0-brane in the black zero-brane geometry [111]. Low-energy states above the trivial vacuum of the BFSS matrix model can describe the black zero-brane in type-IIA superstring theory [67]. (See Ref. [112] for a generalization to the BMN matrix model.) By exciting the (N, N )-components of matrices, we can obtain a state which is peaked around X I,ij = y I δ ij in the coordinate basis and P I,ij = q I δ ij in the momentum basis. Such a state admits an interpretation as a bound state of black zero-brane sitting at the origin of the bulk and a D0-brane located at y = (y 1 , · · · , y 9 ) moving with the momentum q = (q 1 , · · · , q 9 ). By monitoring how the peak moves in the coordinate basis and momentum basis via the Hamiltonian time evolution, and by comparing it with the dual gravity picture, we might be able to see how the bulk geometry is encoded in the matrix model. A natural possibility at strong coupling and large N is that the peak moves following the Dirac-Born-Infeld action for a probe D-brane [3].
Another interesting problem, which involves quantum chaos and the Lyapunov exponent. Given that we have an example of the vacua in the BMN matrix model or some other matrix models, it might be interesting to compute the Lyapunov exponents by constructing the following correlation functions on the quantum computer Ô 1 (t)Ô 2 (t )Ô 3 (t)Ô 4 (t ) (4.30) in the Heisenberg picture, where the operatorÔ's could be chosen as simple gaugeinvariant operators. This is the out-of-time-ordered correlator that could reveal the Lyapunov exponent. A quantum simulation protocol proposed in Ref. [113] can be useful for this purpose. Note that the Hamiltonian time evolution has to be controlled very precisely, that would be a challenging task without quantum error correction. In the strong coupling limit where Einstein gravity gives a precise dual description, the Lyapunov exponent should satisfy the Maldacena-Shenker-Stanford (MSS) bound [114]. It is interesting to consider the finite-oupling corrections away from the MSS bound which corresponds to the stringy corrections to the gravity dual 28 [114,[116][117][118].

Conclusion and discussion
In this paper, we discussed how matrix models (especially the BMN matrix model) and some classes of quantum field theories could be realized on a digital quantum computer. In Section 2 and Section 3, we illustrated how the matrix model could be realized in the Hamiltonian formulation. We introduced an explicit regularization scheme, such that it can be realized on a digital quantum computer with a large but finite number of qubits. Then in Section 4, we discussed the actual implementations. A very standard encoding prescription led to a rather simple form of the Hamiltonian, which allows us to use efficient simulation algorithms based on the block-encoding, qubitization, and quantum signal processing. The minimal number of qubits required to encode the BMN matrix model with U(N ) gauge group and cutoff Λ is 9N 2 log 2 Λ + 8N 2 and our protocol uses additional log 2 L qubits for L ∼ Λ 8 N 4 ancilla states for block-encoding and qubitization step. The circuit complexity for adiabatic state preparation via the fast digital Wan-Kim algorithm of Ref. [11] is and the circuit complexity for approximation time-independent Hamiltonian evolution unitary via qubitization and quantum signal processing Ref. [12,13] is Note that δ is the error in state preparation, ε is the error in unitary, and ∆ gap is set by the mass term ∆ gap ∼ µ at small coupling limit. 29 The cutoff Λ should be chosen appropriately depending on the physics under consideration, and in Sec. 3.3, we provided a rough discussion on several possible choices.
Here we comment on possible implications from this work, combined with quantum Church-Turing Thesis. The quantum Church-Turing Thesis states that any physical process that happens in the real world could be simulated in a quantum computer. We could write it in a more formal way: Any calculation that cannot be done efficiently by a quantum circuit cannot be done efficiently by any physical system consistent with the laws of physics. 30 Our work shows that some supersymmetric QFTs which have natural realizations in superstring/M-theory are simulatable, and hence, they are not excluded by the quantum Church-Turing Thesis. We do not yet know generic supersymmetric QFT, such as the minimal supersymmetric standard model, can be simulated efficiently. In principle, most theories can be simulated if the parameter-fine-tuning is allowed, and it would be important to understand the computational complexity of the parameter tuning.
In addition to the ones discussed in this paper, what kind of supersymmetric theories can be simulated without involving the parameter fine tuning on a digital quantum computer? Kaplan, Katz, and Unsal gave a regularization scheme of the spatial lattice with the continuous time for supersymmetric theories with four, eight, and sixteen supercharges in one, two, and three spatial dimensions [18], and showed that theories in one or two spatial dimensions do not require parameter fine tuning. It is likely that their lattice Hamiltonian can be simulated efficiently on a universal quantum computer. 31 By adding the flux deformation to the (1+1)-d theory and taking the fuzzy-sphere vacuum, the (3+1)-d theory might be obtained, as previously done for the theories on the two-dimensional Euclidean lattice [41,42,120]. In those theories, we can discuss the complexities of the calculations in a quantitative manner based on actual regularizations. It might be interesting to understand deeper about potential implications of supersymmetry in the context of the quantum Church-Turing Thesis, or discuss their relations between some recent discussions about quantum simulation capabilities and quantum black holes [121][122][123][124].
In this paper, we pointed out that the use of the equivalence between the matrix model and QFTs can simplify the implementation of the latter on a quantum computer. Let us take this opportunity to briefly discuss other 'non-lattice' approaches in the context of quantum simulation and quantum computation [125,126]. There also exists a similar approach to the method proposed in this paper about simulating quantum field theories using Hamiltonian simulation, which is called conformal truncation. The method of conformal truncation is proposed as an alternative of lattices for studying generic strongly-coupled quantum field theories: unlike simulating field theories by turning on couplings from free theories in the lattice, conformal truncation solves field theories from another side of the RG flow: turning on operators away from conformal field theories. The Hamiltonian we arrive at, in this case, might be non-local. One could discuss quantum simulation based on conformal truncation in Ref. [127]. One could also simulate quantum field theories using consistency relations existing already in quantum field theories. This is called the bootstrap approach. One example of the quantum setting for bootstrap problems is developed in Ref. [128], which depends on a theoretical speedup of semi-definite programming algorithms in a quantum computer.
The wave function satisfies which is equivalent tox Now let us make the truncation on the Fock space corresponding to take only the states |n with n = 0, 1, · · · , Λ − 1 as basis. With the truncation, the relation (A.5) for the position operatorx = 1 √ 2ω (a † + a) is no longer true for the highest energy state: Inspired by the relation (A.3) before the truncation, it would be natural to define a regularized version of the "position state" as However, this state is not an eigenstate of the operatorx in a precise sense: The second term is the deviation from the correct relation and irrelevant if we are interested in problems where the highest energy state |Λ − 1 is important.

Matrix model
The matrix model has 9N 2 harmonic oscillators labelled by I = 1, · · · , 9 and α = 1, · · · , N 2 . Before the truncation, we have the "coordinate basis" which satisfieŝ By taking the tensor product, we can define After the truncation, we can approximate the "coordinate basis" by replacing it with the state (A.7) for each I, α.

Regularization in the coordinate basis
A natural regularization scheme in the coordinate basis is to introduce a 'lattice', i.e., to restrict the values of X α I to be nδ X , where n is integer between ±n b . By sending δ X and n b to zero and infinity, respectively, in such a way that n b δ X becomes infinite, the original coordinate basis is reproduced. The momentum operatorP α I is approximated by the difference operator. Such a basis is used in Ref. [60] for scalar field theories. We do not expect a big difference from the Fock basis; the Hamiltonian can be expressed as a sum of the Pauli strings anyways, with more or less the same number of the terms.
B Construction of |J i ; VAC g.f.

We introduceâ
Hence the Fock vacuum ofâ Aâ A and that ofâ A are the same, in the sense Therefore,

the creation and annihilation operators in the original basis,Â andÂ
This completes the construction of |J i ; VAC g.f. :

C The Wan-Kim Algorithm
In Sec. 4, we used Algorithm 1 and Theorem 1 from Ref. [11] as a black box to construct the ground state of the matrix model. In this section, we give a brief summary of the Algorithm 1 and key ideas behind the proof of Theorem 1. Algorithm 1 prepares the ground state ofĤ 1 by simulating the adiabatic evolution viaĤ(s) = (1 − s)Ĥ 0 + sĤ 1 . The algorithm takes the low-level oracles that block-encode HamiltoniansĤ 0 , H 1 andĤ =Ĥ 1 −Ĥ 0 as inputs. Wan and Kim use the machinery of quasi-adiabatic continuation to provide a protocol that converges in polylogarithmic time in target state precision error, which is better than previously known adiabatic protocols. The goal of the construction is to approximate the entire adiabatic evolution unitaryÛ (s) generated by quasi-adiabatic continuation operator where W (t) is an odd function satisfying W (t) ≥ 0 at t ≥ 0. The main insight of the paper is finding a function W (t) that gives a good bound of the adiabatic error (|| |Ω(s) − U (s) |Ω(0) ||) and is easy to simulate digitally. Digital simulation complexity is closely related to the complexity of integrating function W (t) on a quantum computer. The function that satisfies both of these requirements is In Theorem 3, they prove that the adiabatic error is bounded as follows; where γ(s) denotes the gap of the Hamiltonian H(s), ∆ is additional control parameter and U ∆ (τ ) is the ordered exponential ofD(s). In Appendix C, the authors provide a low complexity circuit for performing the discrete integral of the function W (t). They define the log M -qubit state D Adiabatic state preperation, LCU decomposition, and oracles D.1 Naive adiabatic state preparation We start by reviewing how the simplest adiabatic state preparation works based on the adiabatic theorem. The Hamiltonian consists of the free part and the interaction part, whereĤ int disappears when the coupling constant g is set to zero. The free partĤ free does not contain g. For the adiabatic state preparation, the coupling g to depend on time t, and defineĤ(t) byĤ (t) =Ĥ free +Ĥ int (g(t)).

(D.2)
In the weak-coupling limit g → 0, we know the ground state precisely. 32 We take g(t = t init ) to be parametrically small, and take the initial state of the simulation to be the analytically-known ground state which we denote by |Ω(t init ) .
From t = t init to t = t fin , we gradually increase the coupling g(t), such that the final value g(t = t fin ) becomes the value we are interested. Then, if the change is sufficiently slow, the state at time t is the ground state ofĤ(t), due to the adiabatic theorem.
More precisely, the difference between the ground state ofĤ(t fin ) denoted by |Ω(t fin ) and the actual state |ψ(t fin ) is where ∆ gap is the scale of the mass gap of the Hamiltonian during the whole timedependent process, ∆t = t fin − t init , and s = (t − t i )/∆t. (We recommend a nice review [129] about this theorem.) Then the remaining problem is if we can perform the time-evolution with the timedependent HamiltonianĤ(t) on the quantum computer. For local Hamiltonians, a typical method is to use the product formula to decompose local terms and write a factorized product of the exponential (see one of the earliest papers [130] and a recent paper [131]). One of the most celebrated applications of this algorithm is the Jordan-Lee-Preskill algorithm, where they use it to study adiabatic state preparation in the λφ 4 theory in general dimensions; see also Refs. [60,[132][133][134][135][136][137][138][139]. 32 Note that we cannot take g to be exactly zero if we consider nontrivial a fuzzy sphere vacuum X α 1,2,3 |fuzzy sphere µ 3g J α 1,2,3 |fuzzy sphere .
The naive Trotterization is not suitable for the current setup, due to the lack of a manifestly local Hamiltonian in the computational basis. Thus, we consider alternative algorithms (a similar situation was discussed in the Hamiltonian truncation formalism of quantum field theories, see [127,140]). In quantum information science, there are alternative algorithms for non-local Hamiltonian evolution, mostly designed for quantum chemistry in the near-term or long-term simulation (for a review, see Ref. [141]).

D.1.1 Truncated Taylor-series method and linear combination of unitaries (LCU)
Here, we discuss a relatively simple implementation of the Hamiltonian simulation based on the truncated Taylor-series based on input models of the linear combination of unitaries (LCU), namely, the paper [142]. The discussion is very friendly to people who are not familiar with quantum simulation algorithms. We first explain the time-independent-Hamiltonian version and then modify it to the time-dependent-Hamiltonian version. That the Hamiltonian is easily expressed as a sum of Pauli strings (as we saw in Sec. 4.1) makes the simulation straightforward.

Time-independent-Hamiltonian version
In this method, the time evolution is written as and e −iĤt/r is approximated by truncating the Taylor expansion at some order, The value of positive integer r is chosen later, in such a way that a technical assumption (D.8) needed for an efficient method is satisfied.
In the BMN matrix model, the Hamiltonian is written as a sum of Pauli stringŝ Π i as eq.(4.7). By substituting (4.7) to (D.5) and then plugging it into (D.4), we can rewrite the write hand side of (D.4) as a sum of the products of the Pauli strings. The products of Pauli strings are again Pauli strings. Hence e −iĤt/r can be expressed as whereΠ's are again Pauli strings. Note that this decomposition depends on t, r and K. The values of β i 's can be evaluated without using a quantum computer, and they are used as a part of the inputs for the quantum simulation. We choose r so that the value of s defined by becomes 2: When K is sufficiently large, this is equivalent to t r L i=1 α i = log 2 with a good precision. This condition can always be satisfied, by allowing the identityÎ as one of Pauli stringP i i in (4.7) and tuning its coefficient α i . 33 In order to utilize this decomposition, we introduce ancilla states |i (i = 1, · · · , m), and define an operatorV acting on C m times the Hilbert space aŝ

Time-dependent-Hamiltonian version
The generalization to the time-dependent Hamiltonian is tedious but straightforward. e −iĤt/r should be replaced with the Dyson series T e −i nt/r (n−1)t/r dtĤ(t) n = 1, 2, · · · , r, (D. 13) and (D.5) is replaced by (D.14) Here T stands for the time ordering, i.e., operators at a later time come left. Because the time-dependence is only in the coupling constant g(t) in (D.2), this can be solved before using a quantum computer, and the same form as (D.6), with different values of the coefficients β i , can be obtained.
A full discussion about the time-dependent Hamiltonian simulation algorithms is included in [95], where the main theorem is given in Theorem 10'.

D.1.2 Utilizing the oracles and sparseness
The truncated Taylor series method explained in Sec. D.1.1 does not fully take advantage of the sparseness of the Hamiltonian; as we can see from Sec. 4.1, there are many cancellations among Pauli strings such that only a small number of nonzero entries remain, but if we treat each Pauli string separately as in (D.9), this cancellation is not utilized at all. By using another decomposition of the Hamiltonian, a more efficient simulation might be achieved. Then, in principle, how efficient can the simulation be?
In order to answer this question, one has to make some assumptions regarding the available oracles. An example of the oracles isV in (D.9); one assumes the existence of some oracles with which the sparseness can be fully utilized, and count the number of the gates and oracles necessary for the simulations.
Here, we discuss another input model that could manifest the sparseness. There has been substantial effort toward efficient quantum simulation algorithms utilizing the sparseness of the Hamiltonian. As one of such algorithms, we introduce the rescaled Dyson series method with sparse oracles designed in Refs. [95]. We will only provide the statement, and we provide some additional information in the appendix D.2 for selfcompleteness. The following algorithm takes advantage of sparsity, which is discussed before, pretty generic in matrix models.
We start with the definition of sparsity. For a given Hamiltonian H, the sparsity of the Hamiltonian d is the maximal number of non-zero entries it could have in any row or column. As we have seen in Sec. 2.4 the sparsity is d ∼ N 4 in the BMN matrix model.
Then the statement proven in Ref. [95], with some standard assumptions regarding the oracles (see the appendix D.2 for details), is  This theorem is the Theorem 10 of Ref. [95]. Here n q is the number of qubits in the Hilbert space, and H max,1 is the one-norm (defined in the L 1 space of the Hilbert space) for the maximal matrix element of the Hamiltonian, which is defined by 34 The notationÕ means that we are ignoring logarithmic factors. This algorithm is sufficient for us when doing time evolution for our adiabatic state preparation purpose. Finally, it is worth notice that although we only have logarithmic dependence about precision during time evolution, there is polynomial dependence of adiabatic errors when applying the algorithm towards adiabatic state preparation. The latter might be significantly improved according to the work by Wan and Kim [11], which we have discussed in the main text.

D.2 Hamiltonian simulation based on the rescaled Dyson series
Here we discuss some details about the algorithm appearing in Appendix. D.1.2. Let us start by repeating the definition of the sparsity: For the BMN matrix model in the Fock basis, the sparsity is d ∼ N 4 . Now, we start introducing our oracles. We assume that H(τ ) is at most d-sparse at any τ ∈ [0, T ]. (This is actually the case in Sec. 4.4.) We consider a set of basis states with four labels |τ, j, k, z = |τ |j |k |z , where τ is the (discretized) time, j and k are integers which run from 1 to dimH where H is the Hilbert space of the system under consideration (and hence the Hamiltonian is a dimH × dimH matrix), and z ∈ C is a complex number which is expressed by using binaries with some accuracy. We define two oracles O loc and O val to help us access the Hamiltonian. They are defined by 35 O loc (|j |s ) = |j |col(j, s) , (D. 18) and O val |τ, j, k, z = |τ, j, k, z ⊕ H jk (τ ) . (D.19) Here, the notation col(j, s) is used to denote the location at the j-th row and the s-th non-zero element in it. where ||H(τ )|| max ≡ max jk |H jk (τ )|. This is used to compute the max-norm. The latter requires another state |w , with w ∈ C. By using w = f (t) ≡ The cost and feasibility of the actual implementation of such oracles depend on the Hamiltonian and the architecture. Given that the matrix models have rather simple Hamiltonians, it would not be unrealistic to assume the existence of such oracles as a starting point for the discussion of efficient quantum simulations. For constructions of some oracles, see [143], G4 and [144], III A and III B.
The analysis of query complexity includes the determination of how many queries we address for given oracles, and moreover, how many additional gates we need in this process. By using the oracles defined above, the following theorem can be shown: Theorem 3. (Theorem 10 in [95]) Suppose that the time-dependent Hamiltonian H(τ ) acting on the 2 nq -dimensional Hilbert space consisting of n q qubits is d-sparse. Suppose also that there is an upper bound on the max-norm, H (f −1 (ς)) max , that is positive, and continuously differentiable. Then, there exists an algorithm such that the Hamiltonian evolution could be simulated using additional gates. ByÕ we mean an order estimate up to logarithmic corrections, and by we mean the given error.

E Quantum signal processing
Following Refs. [13,47], we review the quantum signal processing method, which was mentioned in Sec. 4.3.
E.1 Proof of (4.14) The n-th order Chebyshev polynomial of the first kind is defined by Let us prove (4.14) by using the mathematical induction. We can directly check (4.14) for n = 0 and n = 1. Suppose (4.14) holds for n − 1 and n. Then, for n + 1, we can check (4.14) as follows. Firstly, note that G|R = G|.

E.2 Quantum signal processing
By combining (4.14) and the Jacobi-Anger expansion (4.15), we obtain The ancilla state |G can easily be realized. Hence, if f (Ŵ ) can be realized efficiently, the Hamiltonian time evolution can be simulated. The quantum signal processing [13] provides us with a black box to construct f (Ŵ ) whenŴ is provided as an input.

F A short introduction on Hamiltonian simulation algorithms
In this appendix, we give a short introduction and overview of the existing algorithms about the Hamiltonian simulation. One could read Ref. [141] for some more detailed discussions.
As we mentioned before, perhaps the simplest way for simulating Hamiltonian evolution is through Trotterization, namely, to use the Lie-Trotter-Suzuki formula, or the product formula. The idea about Trotterization is that one could decompose the whole Hamiltonian towards local terms and simulate them separately for a relatively short time. This requires us to divide the total time towards short time steps. The more time steps we use, the more accurate result we obtain, but the more gates we need. There are many historical discussions along the line of product formulas, see, for instance, Refs. [130,131,[148][149][150]. An important type of improvement from the original product formula is to make use of randomization. See for instance, Refs. [151][152][153][154]. Moreover, one can also formulate a time-dependent version of the product formula algorithm. See, for instance, Refs. [155,156].
Beyond the naive application of product formulas, we might consider using some more advanced algorithms. Usually, those kinds of algorithms include some black box Hamiltonian input models, which we call oracles, and ask how many queries we need to access the oracles. This type of algorithms includes algorithms based on quantum walks [157,158], multiproduct formula [159,160], Taylor expansion [142,161], fractional-query models [145], Chebyshev polynomial approximations [162], qubitization [12,147], and quantum signal processing [13,163]. Many elements in the web of such algorithms are conceptually or technically related. In terms of input models, there are two common choices. One is through a linear combination of unitaries (LCU) by a decomposition of the Hamiltonian into a series of unitaries. While one could also construct oracles based on accessing matrix entries, which usually involves some assumptions about sparseness. For time-dependent situation, progresses are made in Refs. [142,146]. An alternative strategy is to construct oracles that encode the Hamiltonian as a block-diagonal element of an oracle unitary matrix, known as block encoding. See Refs. [12,147] and Ref. [11], where block encoding is used for time-independent simulation as well as digital adiabatic state preparation.
Those oracle-based methods are widely used in quantum computational chemistry [141]. For a general energy-level system, a helpful generic introduction is given in Ref. [164].
In this paper, we will mainly consider two different algorithms. For real-time dynamics with the static Hamiltonian, we construct a protocol based on block-encoding, qubitization, and quantum signal processing [13]. This example is particularly important since (similar to Ref. [47]) we have an explicit oracle construction for the matrix models. For the adiabatic state preparation, we will review and use the Wan-Kim algorithm [11]. The Hamiltonian simulation industry is under rapid construction, and some better algorithms may appear in the future. In Appendices D and D.2, we provide alternative algorithms based on LCU decomposition, rescaled Dyson series, and the naive adiabatic state preparation.