Resource-efficient digital quantum simulation of $d$-level systems for photonic, vibrational, and spin-$s$ Hamiltonians

Simulation of quantum systems is expected to be one of the most important applications of quantum computing, with much of the theoretical work so far having focused on fermionic and spin-half systems. Here, we instead consider encodings of $d$-level quantum operators into multi-qubit operators, studying resource requirements for approximating operator exponentials by Trotterization. We primarily focus on spin-$s$ and truncated bosonic operators in second quantization, observing desirable properties for approaches based on the Gray code, which to our knowledge has not been used in this context previously. After outlining a methodology for implementing an arbitrary encoding, we investigate the interplay between Hamming distances, sparsity structures, bosonic truncation, and other properties of local operators. Finally, we obtain resource counts for five common Hamiltonian classes used in physics and chemistry, counts that include the possibility of converting between encodings within a Trotter step. The most efficient encoding choice is heavily dependent on the application and highly sensitive to $d$, though clear trends are present. These operation count reductions are relevant for running algorithms on near-term quantum hardware, because the savings effectively decrease the required circuit depth. Results and procedures outlined in this work may be useful for simulating a broad class of Hamiltonians on qubit-based digital quantum computers.

Simulation of quantum systems is expected to be one of the most important applications of quantum computing, with much of the theoretical work so far having focused on fermionic and spin- 1 2 systems. Here, we instead consider encodings of d-level quantum operators into multi-qubit operators, studying resource requirements for approximating operator exponentials by Trotterization. We primarily focus on spin-s and truncated bosonic operators in second quantization, observing desirable properties for approaches based on the Gray code, which to our knowledge has not been used in this context previously. After outlining a methodology for implementing an arbitrary encoding, we investigate the interplay between Hamming distances, sparsity structures, bosonic truncation, and other properties of local operators. Finally, we obtain resource counts for five common Hamiltonian classes used in physics and chemistry, counts that include the possibility of converting between encodings within a Trotter step. The most efficient encoding choice is heavily dependent on the application and highly sensitive to d, though clear trends are present. These operation count reductions are relevant for running algorithms on near-term quantum hardware, because the savings effectively decrease the required circuit depth. Results and procedures outlined in this work may be useful for simulating a broad class of Hamiltonians on qubit-based digital quantum computers.

I. INTRODUCTION
Simulating quantum physics will likely be one of the first practical applications of quantum computers. In simulating the many body problem, most algorithmic progress so far has focused on systems with binary degrees of freedom, e.g. spin- 1 2 systems [1,2] or fermionic systems [3,4]. The latter case is relevant to simulations for chemical electronic structure [5,6], nuclear structure [7], and condensed matter physics [8]. This focus on binary degrees of freedom seems to be a natural development, partly because qubit-based quantum computation is the most widespread model used in theory, experiment, and the nascent quantum industry.
In principle, there are combinatorially many ways to map a quantum system to a set of qubits [23,24]. Mapping a dlevel system to a set of qubits may be done by assigning an integer to each of the d integer levels and then performing an integer-to-bit mapping. d-level-to-qubit mappings have been used in the very recent literature, primarily for truncated bosonic degrees of freedom [14,[17][18][19][20][21]25]. Still, unexplored areas include a consideration of more than the two henceforth proposed encodings, as well as the question of which encodings are optimal for which problems. The purpose of this work is both to provide a complete yet flexible framework for the mappings, and to analyze several encodings (both newly proposed herein and previously proposed) for a widely used set of operations and Hamiltonians. Especially for near-term noisy quantum hardware, gate counts and qubit counts will be limited. In principle, these constraints can be used to approximate a hardware budget for a set of hardware and a particular Hamiltonian simulation problem. For example, if one wants to simulate a collection of N bosons on a small quantum computer, the decoherence time and gate errors will constrain the allowed number of gates, while the total number of qubits will constrain the qubit count per boson. In this schematic, we show two arbitrary hardware budgets for Trotterizing the exponential ofq 2 for one boson with truncation d = 5. In device A, both the Gray and standard binary encodings are satisfactory, but the unary code requires too many qubits. However, because device B allows for more qubits but fewer operations, the unary code is sufficient while the former two encodings require too many operations. This highlights the need for considering multiple encodings, as an encoding that is best for one type of hardware is not necessarily universally superior.
Using an arbitrary selection of parameters for common physics and chemistry Hamiltonians, we have plotted the comparative computational costs required for first-order Trotterization. Costs are reported in terms of number of two-qubit entangling gates, relative to the cost of standard binary (SB). The three encodings shown here-standard binary, Gray code, and unary-are defined in the text. Explicit expressions are given in the Supporting Information for the five Hamiltonians: the Bose-Hubbard model, one-dimensional quantum harmonic oscillator (QHO), Franck-Condon calculation, boson sampling, and spin-s Heisenberg model. The optimal encodings are sensitive both to the Hamiltonian class and the number of levels d (determined by bosonic truncation or by the spin value s). In some cases, it is best to stay in a particular encoding for the duration of the simulation. Other times, it is worth bearing the resource cost of converting between encodings, because it saves on total operations. Still other times, the decision to save operations by converting between encodings will depend on whether available hardware is gate count limited or qubit count limited. Four Scenarios, A through D, are discussed in Section VI.
This aids in determining which mappings are more efficient for particular operators and specific hardware, including near-term intermediate-scale quantum (NISQ) devices. When choosing which encoding to use for a given problem, it is conceptually useful to think in terms of a "hardware budget," as shown in Figure 1. Similar considerations have been studied for fermionic mappings [26]. For nearand intermediate-term hardware, one will often have stringent resource constraints in terms of both qubit count and gate count. Imagine that one plans to perform Hamiltonian simulation for some N -particle system. Using some set of criteria for acceptable error and other parameters, one can in principle work backwards to determine how much of a quantum resource is available for each operation. This quantity would be different for each device. Perhaps one quantum computer would allow for more qubits but another allows for more operations, as in Figure 1. Because different encodings yield differing resource requirements, considering multiple encodings may be essential for determining whether the available resources are sufficient.
Here we briefly summarize our results for resource comparisons of real Hamiltonian problems, in order to highlight the utility of encoding analyses and to demonstrate the ultimate practical objective of this work. Figure 2 shows the relative two-qubit operation requirements for a set of five prominent physics and chemistry problems (defined in Section VI and the Supporting Information). All comparisons are made within a given Hamiltonian. Our investigations revealed a somewhat rich interplay between qubit counts, operation counts, encodings, and conversions. The difficulty in a priori predicting the optimal encoding scheme suggests that sophisticated compilation procedures, for automatically choosing and converting between multiple encodings, will play a large role in future quantum simulation efforts for d-level systems.
We note that the optimal encoding schemes have differing characteristics, all of which are present in Figure 2. The results for each Hamiltonian can be categorized as one of four scenarios. Scenario A: the optimal choice is either standard binary (SB) only or Gray only, with no benefit from converting between encodings (Bose-Hubbard d = 4; 1D QHO d = 4). Scenario B : the optimal choice is to convert between SB and Gray, in order to perform different local operators in different encodings (Heisenberg s = 7 2 ; Franck-Condon d = 4). Scenarios A and B are notable because they require both the fewest operations and the fewest qubits, as there is no benefit to expanding into the qubit-hungry unary encoding. Scenario C : unary-only is the optimal choice, and saving memory by compacting the data back to Gray or SB is still cheaper than remaining in the latter encodings (Bose-Hubbard d = 10; boson sampling d = 10; Franck-Condon d = 10). Scenario D: unary is the optimal encoding, but only if the data remains in high-qubit-count unary for the duration of the calculation (1D QHO d = 10; Heisenberg s = 2). The optimal encoding choice is highly sensitive to both the Hamiltonian class and the truncation value d. This suggests the need to perform an encoding-based analysis for any new digital quantum simulation of d-level particles.
Throughout this work, we study four encoding types: unary (also called one-hot), standard binary, the Gray code, and a new class of encodings we name block unary, all defined in Section II. After outlining how to map d-level operators to qubit operators for any arbitrary encoding in Section III, we consider resource counts for most standard local operators used in bosonic and spin-s Hamiltonians in Section IV. In Section V we present circuits for converting between encodings and enumerate their resource counts. Finally, in Section VI we obtain relative resource estimates for various commonly studied Hamiltonians in theoretical physics and chemistry. Conclusions are summarized in Section VII. This work highlights how the careful choice of the encoding scheme can greatly reduce the resource requirements when simulating a system of d-level particles or modes.

II. BINARY ENCODINGS
We begin by giving definitions for various integer to binary encodings. In this work, we use the term binary encoding to refer to any code that represents an arbitrary integer as a set of ordered binary numbers, not just the familiar base-two numbering system. These encodings can be used to represent states regardless of the type of system or basis we are considering.
For a given encoding enc, each integer l has some qubit (i.e. binary) representation, denoted R enc (l), which is a bit string on N q bits, x Nq−1 · · · x 1 x 0 . To specify the value (0 or 1) of a particular bit i in the encoding, we use the notation R enc (l; i).
Standard binary. We refer to the familiar base-two numbering system as the standard binary (SB) encoding, such that an integer l is represented by This straight-forward mapping has been used for qubit-based quantum simulation of bosons previously [14,17,18]. The SB mapping uses N q = log 2 d qubits when the range of integers under consideration is {0, 1, 2, . . . , d − 1}, where · is the ceiling function.
Gray code. In principle there are combinatorically many one-to-one mappings between a set of integers and a set of N q -bit strings. One mapping from classical information theory with particularly useful properties is called the Gray code or the reflected binary code [27]. Its defining feature is that the Hamming distance d H between nearestinteger bit strings is always 1, formally d H (R Gray (l), R Gray (l + 1)) = 1. The Hamming distance counts the number of mismatched bits between two bit strings. In other words, moving between adjacent integers requires only one bit flip (see Table I). As will be seen below, this encoding is especially favorable for tridiagonal operators with zero diagonals, since all non-zero elements |l l | then have hamming distance one. As far as we are aware, the Gray code has not been previously proposed for use in the simulation of bosons and other d-level systems on a qubit-based quantum computer. This encoding inspires the possibility of having a specialized encoding for many possible matrix sparsity patterns (for example a code for which d H = d H (R Gray (l)R Gray (l + 1)) = 2 to be used for pentadiagonal matrices likeq 2 ), but we do not consider this possibility here. Throughout this work we refer to the SB and Gray codes as compact encodings because they make use of the full Hilbert space of the qubits used.
Unary encoding. Mappings that do not use all 2 Nq states of the Hilbert space are possible. In the unary encoding (also called one-hot [28]), the number of qubits required is N q = d. Therefore, only an exponentially small subspace of the qubits' Hilbert space is used. The relationship between the unencoded and encoded ordered sets is ( Previous proposals for bosonic simulation on a universal quantum computer have used this encoding [17,18,[28][29][30], under different names. The unary encoding makes less efficient use of quantum memory, but it will become clear below that this mapping usually allows for fewer quantum operations. Block unary encoding. One can interpolate between the two extremes, using less than the full Hilbert space but more of the space than the unary code uses. In some limited instances this allows one to make tunable trade-offs between required qubit counts and required operation counts, which may be especially useful for mapping physical problems to the specific hardware budget of a particular near-term intermediate scale (NISQ) device. In this work, we introduce a class of such encodings that we call block unary.
The block unary code is parameterized by choosing an arbitrary compact encoding (e.g. SB or Gray code) and a local parameter g that determines the size of each block. It can be viewed as an extension of the unary code, where each digit (block) ranges from 0 to g. Within each block, the local encoding is used to represent the local digit. The number of qubits required is N q = d g log 2 (g + 1) . Examples of the block unary encoding are given in Table II. We use the notation BU enc g to define block unary with a particular pair of parameters. For a transition within a particular block, the number of operations is similar to a compact code with d = g + 1. For elements that move between different blocks, the transition will be conditional on all qubits in both blocks.
For codes using less than the full Hilbert space of the qubits (e.g. unary and block unary), some qubits can be safely ignored for a given element |l l |. This is because operations on these excluded qubits will not affect the manifold on which the problem is encoded. Therefore, for each element |l l | a mapping needs to consider only the bitmask subsets of the two integers, ignoring other bits. The operator on the qubit space is then where C(l) ∪ C(l ) is the union of the bitmask subsets of the two integers and the subscripts i denote qubit number. One then converts each qubit-local term |x j x j | to qubit operators using the following four expressions: |0 0| = 1 2 (I +σ z ), For a single term in equation (3), the result is a sum of Pauli strings, where P is the number of Pauli strings in the sum, c k is a coefficient for each Pauli term, and every operator is either a Pauli matrix or the identity:σ kj ∈ {σ x ,σ y ,σ z } ∪ {I}. Note that in this work the set of Pauli matrices is defined to exclude the identity.

B. Significance of the Hamming distance
It is useful to analyze encoding efficiency based on the Hamming distance between R(l) and R(l ). The Hamming distance, which we denote d R H (l, l ) ≡ d H (R(l), R(l )), is defined as the number of unequal bits between two bit strings of equal length. The important observation is that, for a given element a l,l |l l | in equation (3), the average length of the Pauli strings increase as the Hamming distance increases, where length is defined as the number of Pauli operators (excluding identity) in the term. In this subsection, for simplicity we at first assume that the bitmask subset is C(l) = {0, 1, ..., N q − 1}, implying that we are using a compact code such as Gray or SB. But we note that these Hamming distance considerations are relevant to all encodings. In the case of a non-compact encoding, one would consider only the union of the bitmask subsets. To clarify this result, consider the following. For an arbitrary element |l l | written in binary form |R(l) R(l )|, one performs the following mapping: where subscripts denote qubit number,â ∈ {σ x , I}, andb ∈ {±iσ y , ±σ z } (according to equations (6) through (9)). Expanding the RHS of equation (11) leads to an equation of form (10). For any mismatched qubit j, i.e. any qubit for which x j = x j , the clause (â j +b j ) contains two Pauli operators. For matched qubits (x j = x j ), the clause (â j +b j ) instead has one identity and one Pauli operator. Hence more matched bits lead to more identity operators in expression 11, leading to fewer Pauli operators in the final sum of Pauli strings. It follows that the number of non-identity operators can be reduced by having a smaller Hamming distance. This is relevant because Hamiltonians with more Pauli operators require more quantum operations to implement.
Consider the illustrative example of mapping the Hermitian term |3 4| + |4 3| to a set of qubits. In the SB encoding, equations (6)-(9) yield the following Pauli string representation: Using the Gray code, the Pauli string instead takes the form The Hamming distance between R(3) and R(4) is d SB H = 3 in the former case and d Gray H = 1 in the latter. The result is that the Gray code has fewer Pauli operators per Pauli string, meaning that it can be implemented with fewer operations.
FIG. 3. Canonical quantum circuits used to exponentiate Pauli strings on a universal quantum computer. It is clear from inspection that one needs 2(j − 1) two-qubit gates for such an operation, where j is the number of Pauli operators in the term. When a product of many exponentials is used, as in the Suzuki-Trotter procedure, there tends to be significant gate cancellation.

C. Avoiding superfluous terms in non-compact codes
When implementing local products of operators using non-compact codes (unary or BU), one should multiply operators in the matrix representation before performing the encoding to qubits. If one instead first maps local operators to qubit operators, and then multiplies the operators, superfluous terms may result. For example, when implementing an arbitrary squared operatorÂ 2 , one should begin with the matrix representation ofÂ 2 instead of squaring the qubit representation ofÂ. To see this, consider the unary encoding of the square of a 3-level matrix z . If one begins withn 2 d=3 = diag[0, 1, 4], the encoded operator iŝ If one instead squares the already-encoded Pauli operator forn d=3 , this yields Superscripts denote qubit number. Pauli operators (14) and (15) behave identically on the subspace of the unary encoding, though operator (14) is less costly to implement. One might attempt to eliminate superfluous terms after the mapping is complete, but this is likely a hard problem. In principle it may require combinatorial effort to determine which combinations of operators leave the encoding space unaffected. Hence the most prudent strategy is to always perform as much multiplication as possible in the matrix representation. These considerations are irrelevant when using one of the compact encodings.

D. Trotterization and gate count upper bounds
Hamiltonian simulation often consists of implementing the unitary operator for some user-defined HamiltonianĤ, where t is the evolution time and we have set = 1. For the purposes of this study, we consider Hamiltonians that are sums of local Pauli strings such thatĤ can be expressed aŝ which takes the same form as equation (10) with the Pauli strings and their coefficients compacted into terms {ĥ k }.
In practice, Hamiltonian simulation can be performed using a Suzuki-Trotter decomposition where the expression is exact in the limit of large η or small t [1,31]. The numerical studies of this work consider the encoding-dependent resource counts for equation (18), for a subset of prominent physics problems. We focus on determining the resources required for simulating a single Trotter step, which is equivalent to simulating the quantum propagator for small time t.
Note that there are many variations and extensions to the Hamiltonian simulation approach of equation (18), including higher-order Suzuki-Trotter methods [31], the Taylor series algorithm [32], quantum signal processing [33], and schemes based on randomization [34,35]. Notably for the current work, numerical evidence suggests that simple first-order Trotterization is predicted to have lower error for near-and medium-term hardware [36], even if the other methods are asymptotically more efficient. Since {ĥ k } takes a different form depending on the chosen encoding, the resource counts as well as the error |Û (t) −Ũ (t)| will be different as well. We leave the study of numerical error for future work, as the goal of the current work is to introduce these mappings and to understand some trends in their resource requirements.
Each term exp(−itĥ k ) = exp(−itc k jσ kj ) may be implemented using the well-known CNOT staircase quantum circuit shown in Figure 3. If a qubit j is acted on byσ x orσ y , additional single-qubit gates are placed on qubit j as shown in the figure. H ≡ iR x+z (π) denotes the Hadamard gate that changes between the Z-and X-basis, and iR x+y (π) converts between the Z-and Y-basis.
To exponentiate a single Pauli string, the number of CNOT (CX) gates required is where p is the number of Pauli operators ({σ x ,σ y ,σ z }, excluding I) in the term to be exponentiated. We briefly consider upper bound gate counts for a simple Hermitian term α|l l | + α * |l l| and for a single diagonal element |l l|. In this subsection, d H denotes d H (R(l), R(l )), which will depend on the chosen encoding, and d H = 0 refers to a diagonal element. We define K as |C enc (l) ∪ C enc (l )|, the number of qubits in the relevant bitmask subset. Due to the fact that we are considering products of two-term sums (equations (6)-(9)), the distribution of Pauli strings can be analyzed in terms of binomial coefficients. Some resulting upper bounds (derived in the Supporting Information) are In practice, because substantial gate cancellation is possible once the quantum circuit has been compiled, these upper bounds are usually not directly applicable when considering all the elements in a matrix operator. However, the above expressions may find direct utility in limiting cases, and they demonstrate the basic relationship between Hamming distance, size of the bitmask union, and gate counts.
As an illustrative example, consider a tridiagonal real matrix operatorB with zeros on the diagonal, i.e. with . This is the sparsity pattern of several commonly used d-level operators, such as the bosonic position operatorq. Qubit counts and upper bounds forB are shown in Figure  4, as a function of d for different encodings. To calculate upper bounds forB, we first encode the entire operatorB into a sum of Pauli strings before collecting and cancelling terms. Then we apply equation (19). There is roughly an inverse relationship between the qubit counts and the operation counts, because sparser encodings like the unary and block unary have fewer-qubit bitmask sets but require more total qubits.
The differing gate count upper bounds between the SB and gray encodings ( Figure 4) are explained by Hamming distances. Because all non-zero terms inB have unity Hamming distance, upper bounds for the Gray code are substantially lower. The other notable trend is that the unary code has lower upper bounds, asymptotically, than the other codes. This can be explained using equation (20) by noting that K = 2 for any matrix element, while K = log 2 d for Gray and SB. In other words, K stays constant in the unary encoding, whereas in the compact codes K increases with d. Upper bounds for BU Gray g=3 are between the compact codes and the unary code, as this encoding has an intermediate value of K. Below we will see that, though these trends generally persist, they are less pronounced and less predictable after circuit optimization. for block unary where g is the size of the block. Asymptotically, the number of qubits scales logarithmically for the SB and Gray encodings, and linearly for the unary and block unary encodings. Bottom: Upper bounds of CNOT operation counts for implementing one Suzuki-Trotter step ofB. This is the sparsity pattern of canonical bosonic position and momentum operators as well as the Sx spin operators in spin-s systems. Upper bounds were calculated by mapping the full operator to a sum of weighted Pauli strings, combining terms, and then using equation (19). Note that encodings with higher qubit counts tend to have lower upper bounds for gate counts, and vice-versa.

E. Diagonal Binary-Decomposable Operators
An important class of operators to consider is those which we call diagonal binary-decomposable (DBD). We define DBD operators as being diagonal matrix operators for which the diagonal entries of the operator (Ô) may be expressed asÔ A common subclass of DBD is the set of diagonal operators containing evenly-spaced entries. We call these diagonal evenly-spaced (DES) operators. An example is the bosonic number operator and any linear combination an + bI where a and b are constants. If log 2 d is an integer, then the Pauli operator is simply the base-two numbering system with k i = 2 i , The DBD class of operators is notable because, when log 2 d is an integer, exactly implementing exp(−iθn) requires only log 2 d single-qubit rotations and no entangling gates.
An operator for which log 2 d is a non-integer will not lead to this favorable only-single-qubit decomposition. For example, theŜ z operator for a spin-s system is DBD, but the advantage appears for the SB mapping only when d = 2s + 1 is a power of 2, namely s ∈ {3/2, 7/2, . . . }. However, for other operators one may simply increase d without changing the simulation result. For example, if one is required to exponentiate a truncated bosonic operatorn with at least d = 11, it is most efficient simply to implement the SB encodings ofn with d = 16 instead. A simple example illustrates this point. The standard operator for the bosonic number operator with truncation d=3 iŝ while that for d=4 isn The latter operator (d = 4) is composed only of single-qubit operators but the former (d = 3) is not. Operations counts for CNOT gates are shown in Figure 5, where it is clear that the SB mapping is superior when d is a power of 2. The right panel, for example, which gives gate counts for operators such asn 2 , yields some advantages when d is a power of 2, but entangling gates are still required. Note that, as is also clear from the right panel, the square of a DES or DBD operator is not necessarily DBD.

A. Operator definitions
Though local d-level operators can in principle contain arbitrary terms and even be entirely dense (i.e. a molecule's electronic energy levels with non-zero transitions between each), in practice there is a small set of sparse bosonic and spin-s operators that are used most often. Here we summarize the set of d-level operators used in this study, where it is conceptually useful to explicitly write down some d-by-d matrix representations. Bosonic operators can be constructed from the well-known ladder operatorsâ andâ † , where (importantly for encoding considerations) all non-zero terms |l l | obey |l − l | = 1. The position operatorq = 1 √ 2 (â † j +â j ) is tridiagonal with zeros on the diagonal:q This means that the square ofq, often used in vibrational and bosonic Hamiltonians, is pentadiagonal but with zeros for terms where |l − l | = 1:  (which are 4-and 8-level system respectively), the SB code requires both fewest operations and fewest qubits.
Notably, |l − l | for non-zero entries is either 0 or 2, making the Gray code less useful for this operator. The momentum operatorp = i √ 2 (â † j −â j ) and its squarep 2 have the same sparsity patterns asq andq 2 respectively. The number operatorn =â †â of equation (22) is diagonal and DBD, which leads to efficient SB mappings as discussed in Section III E. Finally, we study the two-site bosonic interaction operator In order to consider spin Hamiltonians such as Heisenberg models [11,37], we study spin-s operators of arbitrary s, where the number of levels is d = 2s + 1. Matrix elements for transitions |l l | are defined as follows [38] l|Ŝ z |l = (s + 1 − l)δ l,l where δ α,β is the Kronecker delta. TheŜ z are DBD operators, whileŜ x andŜ y are tridiagonal with zeros on the diagonal, the same sparsity pattern as bosonicp andq operators. The local operators considered thus far are effectively second quantization operators-each ket tends to correspond to an eigenstate in an isolated d-level system. Also of note is a recently proposed approach [19,20] which maps bosons to qubits using the first quantized representation of the quantum harmonic oscillator. The original proposal maps Hermite-Gauss functions, the eigenfunctions of the quantum harmonic oscillator, into a discretized position space. The approximate position operator is defined asX 32) and N x is the number of discrete position points such that where ∆ is chosen such that the desired highest-order Hermite-Gauss function is contained within (x 0 , x Nx−1 ). Advantages and disadvanteges are discussed in the SI. We raise the possibility of using this approach partly to point out that a N x -by-N x matrix operator may be mapped to qubits using the exact same procedure as the other operators, with N x replacing d. Note thatX F Q is a DBD operator.

B. Numerical results
Quantum circuits for approximating the exponential of each operator and for each d were compiled and then optimized using the procedure given in the SI. The optimization consists of searching for and performing gate cancellations where possible. For instance, two adjacent CNOT gates or two adjacent Hadamard gates will cancel. Entangling gate counts for the optimized circuits of bosonic operatorsq,q 2 , andâ † iâ i+1 + h.c. are plotted in Figure 6. We place significant focus on smaller d values because they tend to be more common in physics simulation, but we note that applications requiring larger d values do exist, for example in vibronic simulations where occupation numbers can approach d = 70 [18].
Comparing the tridiagonal operatorq with the upper bounds given in Figure 4 demonstrates that the circuit optimization greatly reduces gate count for the compact codes and block unary, often by a factor of 2 to 3. On the other hand, the unary encoding effectively sees no improvement from optimization, though it remains the code with fewest entangling gates for a large subset of d values.
As was the case in the upper bound calculations, operators built from tridiagonal matrices (q andâ † iâ i+1 + h.c.) show the Gray encoding outperforming SB, though after optimization the advantage is less pronounced. In contrast, for the pentadiagonalq 2 , the Gray code outperforms SB asymptotically, while SB is better for lower d values (and lower d values are likely to be more common in relevant bosonic Hamiltonian simulations). The changed trend can be explained by noting that the unity Hamming distance of the Gray code is not as advantageous for the sparsity structure ofq 2 , given in equation (27). Also notable is the apparent dip in operation count at d=8, due to the fact that the diagonal ofq 2 is DBD.
Importantly, when mapping bosonic problems using compact encodings, it is sometimes the case that increasing the truncation value d is beneficial. For instance, suppose one knows one can safely truncate at d = 5 for a bosonic problem. When implementingq 2 , one would instead simply implement the operator for d = 8, as the number of gates decreases while the number of qubits remains the same. Note that this is not possible in spin-s particles, as it would cause leakage to unphysical states.
One of the more intriguing results is that the unary code is often inferior to the Gray or SB encodings. Pronounced examples of this inversion include d=4,7,8 forq and d=4 and 8 forq 2 , among others. This is notable because, for these values of d, Trotterizing the operator requires both fewer qubit and fewer operations if Gray or SB is used. These results are in contrast to the naively expected trend that there would be a more consistent trade-off between qubit count and operation count. Results for single-particle spin-s operatorsŜ x andŜ z , as well as interaction operatorŜ z are plotted in Figure  7. Notably, unlike bosonic Hamiltonians, the d values are not simulation parameters but are determined by s in the system we wish to simulate. The trends in spin operators tend to be more unruly than those in the bosonic operators.
Analogous to the bosonic case,Ŝ z is DBD and therefore SB requires only single-qubit gates when d=4,8 (s= 3 2 , 7 2 ) and no entangling gates. For these two values, SB uses both the fewest operations and the fewest qubits (fewer than unary). However, the Gray code is superior to SB for other values of s, the same behavior seen in the general DES matrix of Figure 5. BecauseŜ z is diagonal, the unary always requires just d single-qubit rotations and no entangling gates. The two-particle operatorŜ (i) zŜ (j) z displays similar trends toŜ z . As expected, the Gray code is usually superior to SB for the tri-diagonalŜ x , because of the unity Hamming distance between nearest levels. Unary is inferior in both gate count and qubit count for most values, a result highlighted earlier in bosonic operators.
For both bosonic and spin-s operators, we have until now omitted discussion of the Gray-based block unary encoding with parameter g = 3. There is never a case where this BU mapping is the sole encoding with the lowest entangling gate count. However, at least in principle, there may be limited cases where a particular hardware budget (Figure 1) dictates the need for a block unary encoding. A necessary condition for even considering the use of BU Gray g=3 is that its operation count is less than both compact codes, but more than unary. BU would have the opposite relation ( Figure  4) in regard to qubit count for most cases. In such cases (q andâ † iâ i+1 + h.c. for d=9;n 2 ∼Ô 2 in Figure 5 for several values;Ŝ (i) zŜ (j) z andŜ z for s=2), it might rarely be that a particular hardware budget would require this encoding for its particular memory/operation trade-off. Needless to say, such highly specific hardware budgets are unlikely to often appear.
Note that the results herein should generally be considered constant-factor savings, because in most relevant systems d does not increase with system size, i.e. with the number of particles. For the simulation of scientifically relevant quantum systems, Hamiltonians are composed of more than one simple operator. For such situations, one may calculate the overall cost within a given encoding, as we do in Section VI. As will be discussed in Section V, it is often beneficial to Trotterize different parts of the Hamiltonian in different encodings, if the cost difference outweighs the overhead of conversion. When extra memory resources are available and the unary code is the most efficient for implementing an operator, one may expand into the unary representation, perform the operation, and then compact the data back to SB or Gray. The example shown here is the bosonic interaction operator a † iâ i+1 + h.c.. This operator is present in bosonic Hamiltonians and for digital simulations of beamsplitters. For many values of d, implementing this operator in unary is much cheaper than implementing it in Gray or SB. When this strategy is worth the cost of conversion, the Hamiltonian simulation is in Scenario C discussed in Section VI. Hence each particle starts with log 2 d qubits, expands out to d qubits, and then compacts back. Whether this procedure leads to cost savings is heavily dependent on the problem and the parameters (Section VI).

V. CONVERSIONS BETWEEN ENCODINGS
It is often the case that different terms in a Hamiltonian are more efficiently simulated in different encodings. For example, in the Bose-Hubbard model, the number operatorn is usually more efficient in SB, while the hopping term b † ib i+1 + h.c. is usually more efficient in the Gray encoding (see Figure 6). The cost of converting from one encoding to another is often substantially less than the difference in resource efficiency between two encodings, which often results in an advantage for compacting and uncompacting the data. For example, if unary is the most efficient for implementing an operator, one may wish to compact the data between operations to save memory resources, as shown in Figure 8. In this section we give general quantum circuits and resource counts for converting between all encodings considered in this work. One can convert between the Gray and SB encodings by applying ( log 2 d − 1) CNOTs in sequential order [39] as shown in Figure 9. From here onward, d is the number of truncation considered in a bosonic mode and K is the number of qubits used in the SB or Gray encoding.  Algorithm 1 An algorithm to build a quantum circuit for SB-unary conversion given an arbitrary truncation d. For multi-qubit gates, the first argument specifies control qubit and the successive ones are targets.
function stdbin-unary-converter(d) Move positions of SB bits Last CNOT gate end for end function The conversion between unary and SB is more complex. The conversion may be especially relevant in a future fault-tolerant quantum computing era, when extra quantum resources are available, because the unary encoding becomes more beneficial as d increases and because the conversion cost is significant. Inspired by previous work [40], in Figure 10 we show an example case for converting from SB to unary when d = 16. A state is initially encoded in SB using qubits on the left, and the memory space is enlarged to include the number of qubits needed for unary. No ancilla qubits are required. As quantum circuits are reversible, unary-to-SB conversion follows by inversion of the circuit. In Table III, we provide the converter circuit resource count for a general d-level truncated quantum system, following Figure 10. Resource counts assume a decomposition of CSWAP into Clifford+T gates [41]. A general algorithm to build the converter circuit can be seen in the Algorithm 1. · and · are respectively the ceiling and floor functions. The validity of the conversion procedure is most easily shown by tracing a single unary state through the reverse algorithm. When d is not a power of 2, modifications are needed. These modifications are already accounted for in Algorithm 1 and example circuits for d =5 and 7 are given in the SI.
For completeness, we constructed a circuit for converting between SB and block unary, for BU SB g=3 , shown in the SI. We do not further analyze BU conversions further, as block unary is in general expected to have limited utility, and even then it will usually be the case that gate count limits are too low to allow for conversions (see Section IV).

VI. COMPOSITE SYSTEMS
Here we consider resource counts for simulating five physically and chemically relevant Hamiltonian systems. The Hamiltonians correspond to the shifted one-dimensional QHO, the Bose-Hubbard model [9,42,43], multidimensional molecular Franck-Condon factors [17,18,44], a spin-s transverse-field Heisenberg model [11,37], and simulating Boson sampling [45] on a digital quantum computer. The former four systems consist of an arbitrary number of d-level particles. For the Franck-Condon factors, the Duschinsky matrix is assumed to have a constant k = 4 nonzero entries per row. With the exception of the simple QHO, all of these problem classes would benefit from digital quantum simulation, because there are limits to the theoretical and practical questions that can be answered by classical computers. The SI gives a more thorough overview of these problems. The differences in resource counts between mappings are constant-factor savings that are independent of system size.
Using resource counts from the optimized circuits for Trotterizing individual operators and from the circuits for inter-encoding conversion, we calculated and compared the required two-qubit entangling gate counts for the selected composite Hamiltonians. We considered five encoding schemes: (i) SB only, (ii) Gray code only, (iii) unary only, (iv) allowing for conversion between SB and Gray, and (v) using all three while compacting to save memory. For (iv) and (v), the reported results include the cost of conversion. To the best of our knowledge, schemes (ii), (iv), and (v) are novel to the context of digital quantum simulation. In Figure 8, an example of encoding scheme (v) is shown. These encoding schemes do not directly map to the 'scenarios' discussed below, as the scenarios consider the optimal encoding scheme under different hardware budgets.
The result for (iv), the encoding scheme that combines both the SB and Gray codes, is reported only when it represents an improvement over both SB-only and Gray-only. We give results for (v), which compacts and uncompacts the qubits for unary computations, only when the unary code was the most efficient of the first four encoding schemes. We give all results in terms of resource counts relative to the SB mapping, noting that the relative resource requirements between encoding schemes are independent of system size (i.e. number of particles or modes).
For some local bosonic operators, the number of entangling gates is not a monotonically increasing function of d. Such operators includen,n 2 , andq 2 (Figures 5 and 6). Our numerics for the composite systems take this into account, increasing the cutoff d if it is beneficial. For instance, if d = 5 is a sufficient truncation and we implement q 2 , we use resource counts for d = 8, because this uses the same number of qubits but fewer operations ( Figure 6). This trick is not possible for the spin-s systems, where d is determined not by a sufficient truncation value but by the nature of the particle itself (its spin s).
A selection of resource comparisons is shown in Figure 2. We show results from d = 4 and 10 because they highlight the variety of rankings that occur, and demonstrate that the best encoding scheme can be highly sensitive to d even within the same Hamiltonian class. Numerical results up to d = 16 are given in the SI.
In terms of which encoding class should be used, the results can be categorized into four scenarios. Scenario A applies when using just one of the compact encodings (SB or Gray) is the best choice. Of the results shown in Figure  2, the Bose-Hubbard and 1D QHO models for d = 4 fit this description. The optimal choice is to stay in one of the compact encodings for the entire calculation, while never using more than log 2 d qubits per particle for the calculation.
Scenario B refers to Hamiltonians for which the optimal strategy is to use a compact amount of memory but to allow for conversion between Gray and SB. This includes the Heisenberg model for s = 7 2 and the Franck-Condon Hamiltonian for d = 4. Because the cost of conversion is very small, this scenario usually implies that at least one of the local operators in the Hamiltonian are optimal in Gray, at least one is optimal in SB, and none are optimal in unary. To take the Heisenberg model with s = 7 2 as an example, one can see this is the case by comparingŜ x and S (i) zŜ (j) z in Figure 7.
Scenario C applies when unary is the superior encoding and it is still considered the best encoding if one unravels and compacts to preserve memory (as in Figure 8). This latter trait is important because it means that, even including the substantial cost of SB-to-unary conversion, with a cost of approximately 9d entangling gates (Table III), it is still better to convert back and forth between unary and SB/Gray. This is true even when memory constraints require that one stores the information compactly for most of the time. This occurs with d = 10 for the Bose-Hubbard, boson sampling, and Franck-Condon Hamiltonians, all of which are bosonic problems.
Scenario D refers to cases where unary is the superior encoding, assuming that the information is not compacted back to SB in order to save memory. This scenario implies that, if one has the qubit space to stay in unary form for the entire calculation, unary is optimal. If one does not have the memory resources for this, it is best to simply perform all operations in Gray and/or SB. The reason for this discrepancy is that the cost of converting binary to unary is substantial, as mentioned above. This scenario applies to the 1D QHO for d = 10 as well as the Heisenberg model for s = 2.
Results for d values up to 16 are given in the Supporting Information. Note that the novel memory-efficient schemes (ii), (iv), and (v) did not lead to improvements in every case; in a minority of the problem instances we considered, the optimal encoding scheme was to use only SB or only unary. By memory-efficient we mean that the scheme stores the encoded subsystems in log 2 d qubits, with the possible exception of when they are being operated on. Comparing to the memory-inefficient unary-only scheme, our novel approaches reduced two-qubit entangling gate counts by up to 33%. Compared to the memory-efficient SB-only scheme, we observed gate count reductions of up to 49%. The latter case is more relevant when qubit count is a substantial constraint. These savings are especially important for running algorithms on near-term hardware, since by simply modifying the encoding procedure one can substantially decrease the effective circuit depth.

VII. CONCLUSIONS
After introducing a general framework for encoding d-level systems to multi-qubit operators, we have analyzed the utility and trade-offs of several integer-to-bit encodings for qubit-based Hamiltonian simulation. The mappings may be used for Hamiltonians built from subsystems of bosons, spin-s particles, molecular electronic energy levels, molecular vibrational modes, or other d-level subsystems. We analyzed the mappings primarily in terms of qubit counts and the number of entangling operations required to estimate the exponential of an operator.
Of the Gray and SB codes, we demonstrated that the Gray code tends to be more efficient for tridiagonal matrix operators, while SB tends to be superior for a common class of diagonal matrix operator. Importantly, we show that converting between encodings within a Suzuki-Trotter step often leads to savings. Interestingly, though the unary code tends to require more qubits but fewer operations, it is often the case that the SB or Gray code is more efficient both in terms of qubit counts and operation counts. Notably, to the best of our knowledge, the Gray code had not been previously used in this context.
We compared resource requirements between encodings for the following composite Hamiltonians: the Bose-Hubbard model, one-dimensional quantum harmonic oscillator, vibronic molecular Hamiltonian (i.e. Franck-Condon factors), spin-s Heisenberg model, and boson sampling. The optimal encoding, and whether it was beneficial to interconvert between encodings, was heavily dependent both on the Hamiltonian class and on the truncation level d for the particle. We placed optimal encoding strategies into four different "scenarios." Which scenario to implement depends on which encoding is optimal in terms of operations, whether interconverting between mappings is worth the additional cost, and qubit memory constraints. The many anomalies in our results highlight the need to perform an analysis of each new class of Hamiltonian simulation problem, determining numerically which encoding scheme is optimal before performing a simulation on real hardware.
There are several directions open for future research. First, there are ways to analyze resource requirements other than enumerating the entangling operations. For long-term error-corrected hardware, estimating T gate count may be most relevant. Additionally, we assumed all-to-all connectivity in this work, a feature of many ion trap quantum computers. But other quantum hardware types require one to consider the topology of the qubit connections and implementation of SWAP gates [46], a consideration that would modify the resource counts and may modify some trends observed here.
We envision that the methodology and results of this work will be helpful for both theorists and experimentalists in designing resource efficient approaches to quantum simulation of a broader set of physically and chemically relevant Hamiltonians.

I. UPPER BOUNDS FOR ENTANGLEMENT GATES
The number of required quantum operations is closely related to the the lengths of the Pauli strings in the encoded Hamiltonian. The distribution of lengths of Pauli strings may be analyzed in terms of binomial distributions, because the quantity to be analyzed is a product of two-term sums. This is because each |x i x i | i in the expression is mapped to a sum of exactly two one-qubit unitary operators, as outlined in the main text. We consider a simple real Hermitian term |l l | + |l l| where l = l . For such a term, it is straight-forward to show that the number f of occurrences of length-p Pauli strings is where K is |C enc (l) ∪ C enc (l )|, i.e. the number of qubits in the relevant bitmask subset, and d H = d H (R(l), R(l )) where R depends on the chosen encoding. For a diagonal term |l l|, the number of occurrences of length-p Pauli strings is instead simply This is because every term in equation (1) is |0 0| or |1 1|, which respectively map to 1 2 (I + σ z ) and 1 2 (I − σ z ). Hence all terms in the product have a single identity-consider e.g. (constant)×(I + σ z )(I − σ z ) · · · (I + σ z )-which leads to the Pauli string lengths having a binomial distribution.
It is conceptually useful to consider the other extreme of all bits being mismatched, i.e. when d H = K for a term |l l | + |l l|. In this case, because there are no identity matrices present in the original product of equation (1), there are 1 2 2 d H 0 0 = 2 d H −1 Pauli strings, all with length K. The reason for the 1 2 factor is that non-Hermitian terms cancel out when we consider the sum |l l | + |l l|. The non-Hermitian term |l l | would not have the 1 2 factor when considering distribution of Pauli string lengths, but we do not analyze such terms because the goal is to exponentiate Hermitian terms using quantum circuits.
Combining equation (2) with n cx (p) = 2(p − 1), an upper bound (denoted U B) for the number of CNOT gates is * nicolas.sawaya@intel.com † gian.giacomo.guerreschi@intel.com where the sum starts at p = 2 because CNOT gates are required only for p > 1. This is an upper bound because gate cancellations are possible once the circuit has been constructed.
Using known binomial identities given below, and again considering only the term |l l | + |l l| or the diagonal |l l|, this becomes The special case of d H = 0 represents a diagonal term |l l|. In the case of compact codes (standard binary and Gray) for which log 2 d is an integer, these expressions can be expressed in terms of the number of levels d as We derived formulas (5) from the following binomial identities Again, thus far we have considered upper bounds only for one term |l l |+|l l| or one element |l l|. We now briefly consider upper bounds for typical full matrix operators. Equation (6) implies that for each term |l l | + |l l|, n cx,U B is of order O(d log d). Hence, for the sparsity pattern ofB discussed in the main text (for which the number of non-zero entries of the operator is O(d)), the compact encodings (Gray and SB) yield an upper bound of O(d)O(d log d)=O(d 2 log d) entangling gates. On the other hand, equation (5) shows that for the unary and BU encodings, one term pair |l l | + |l l| has n cx,U B ∼ O(1), since the number of CNOT gates is independent of d. For sparsity patternB this gives upper bounds of n cx,U B ∼ O(d), substantially less asymptotically than the compact encodings. We remind the reader that the unary and BU encodings use fewer qubits.
For an operatorÂ corresponding to a dense matrix and represented with a compact encoding, the worst-case scaling would be O( In most cases it is not particularly useful to use these expressions for calculating an analytical upper bound for the full operator, except to explain general trends as was done in the main text. When including all terms in the matrix operator, there will be many repeated Pauli matrices whose coefficients will be summed. Therefore, estimating an upper bound by multiplying the number of terms inÂ by the values in equation (5) will often yield unreasonably large estimates for gate counts. Circuit optimizers, like the one used in this work, should be used for more accurate resource estimates.

II. CIRCUIT OPTIMIZATION
The number of required two-qubit entangling gates (CNOT gates) per Suzuki-Trotter step was determined as follows. Each d-level operator was converted to a Pauli operator using the procedure of the main text. Collecting coefficients and cancelling Pauli terms was done using the QubitOperator class of the python software packages FermiLib and OpenFermion [MSK + 17]. Pauli terms were listed in a pseudo-alphabetical order, such that all terms containing σ (0) x are listed first, followed by terms beginning in σ y , etc. Quantum circuits were first synthesized by converting each Pauli term to a CNOT ladder, in the order described. One can see by inspection that circuit reductions will often occur. Pairs of CNOT, Hadamard, or iR x+y (π) gates, which appear after placing CNOT ladders next to each other, can be eliminated as they cancel to unity.
The circuit optimization was performed as follows. In the optimizer, every gate is represented in a data structure that contains its attributes, e.g. its inverse gate, its commutation properties with other gates, and whether it is a rotation gate. For each gate, the optimizer looks for an opportunity to cancel it with its inverse or merge it with another rotation gate by commuting it as far forwards and backwards as possible. The optimizer also looks for a limited set of commonly occurring patterns that allow for gate reductions. For the circuits used in this study, this pattern-searching allows us to reduce some sets of three CNOT gates to two: any gate sequence CNOT(i 0 ,i 1 ) × CNOT(i 1 ,i 2 ) × CNOT(i 0 ,i 1 ) is changed to the equivalent CNOT(i 0 ,i 2 ) × CNOT(i 1 ,i 2 ). Merging and cancelling gates, as well as this pattern-searching, are performed for several passes until the circuit converges.
The choice of ordering for the Pauli terms affects the gate counts, as different orderings of CNOT ladders affect the presence of particular pairs of eliminable gates. Finding the absolute optimal ordering is a combinatorially hard problem and is not the focus of this work [BKM18]. However, to test the quality of the default ordering, we took the encoded bosonicq operator for the unary, gray, and standard binary codes, testing >900 random orderings for each. We found that the default ordering was superior to every random ordering. This result is intuitive, as orderings that match similar Pauli strings side-by-side will lead to more pairs of adjacent CNOT gates and adjacent single-qubit gates.

III. SB TO UNARY CONVERSION: GENERAL CASE
Here we discuss how to modify the SB-to-unary circuit shown in the main text when d is not a power of two. First, we shift the most significant SB qubit up to qubit d − 1 rather than keeping it at qubit 2 K − 1, where K is the number of qubits required by the standard binary encoding. Second, we change the control qubit of SWAP and CNOT gates from 2 K − 1 to d − 1, and remove all unnecessary qubits. Third, the target for the final CNOT is shifted up by the number of qubits that was removed. For instance, using d = 5, the binary base line is shifted from 7 to 4. We then remove qubits 4 through 6 and any associated gates. All SWAP and CNOT controlled by 7 become controlled by 4. Finally, the last CNOT is placed with control 4 and target 0. The latter target has index d − 1 − 2 K−1 . Examples for d = 5 and 7 are given in Figure 1. Algorithm 1 in the main text correctly implements these cases; the examples discussed here are shown only as a pedagogical aid.

IV. SB TO BU SB g=3 CONVERSION
To show that one may in principle convert between all encodings, we provide a circuit for converting between SB and BU SB g=3 with d = 12. We do not write down the general case for this circuit, because the rare hardware budgets for which BU might be useful would likely be cases for which the hardware's gate count constraints would not allow for the cost of conversion.

V. FIRST QUANTIZATION
First quantized operators can be used to represent any bosonic degree of freedom because first quantized QHO operators obey bosonic algebra. The utility of this approach comes partly from the fact that one can quickly move between the position and momentum bases by using the quantum Fourier transform, using a first quantized momentum operatorP F Q that is the direct analog ofX F Q . For instance, consider implementing one Trotter step of an operator proportional to the approximate quantum harmonic oscillator Hamiltoniañ One may first exponentiateX 2 , then apply the quantum fourier transform (QFT) to move to the momentum basis, then applyP 2 , and finally use QFT −1 to move back to the position basis. This allows for all the operations to be done in the diagonal basis except for QFT, which takes just O(N 2 q ) ∼ O(log 2 d) operations. The Trotterization of X F Q andP F Q requires only log 2 d single-qubit rotations and zero two-qubit gates, as explained in the main text, although exponentiating their squares does require entangling gates.
However, a significant computational hurdle to the first quantization method is initial state preparation. Preparing the ground state of the QHO, for example, requires that one prepares the state where f (x) would be a normalized Gaussian, the first Hermite-Gauss function. Macridin et al. discuss existing options for preparing such an arbitrary state [MSAH18], showing that a variational approach is more efficient than one might expect for preparing ground states on near-or medium-term quantum computers. We expect that most early simulations of bosonic degrees of freedom on qubit-based quantum computers will use the Fock-type mappings that are the focus of the current work, since preparing a Fock state in second quantization is trivially done by flipping individual bits. However, there may be instances in which Fock states are not being used, such as wave packet dynamics or other approaches, and it may be that there are other operator-based advantages in using the first quantization mapping. Note that the procedure for mappingX andP to qubits is no different than mapping any other matrix representation, with d replaced by N x .

VI. HAMILTONIAN DEFINITIONS
In this section we provide the Hamiltonians studied in this work and discuss their relevance. Bose-Hubbard model. The Bose-Hubbard model is a ubiquitous model in condensed matter physics, consisting of a set of bosonic modes that occupy discrete lattice sites. The phase diagram of the basic model consists of the Mott insulator and superfluid phases [FWGF89]. Adding additional terms such as increasing the interaction radius or introducing an effective magnetic field can give rise to more complex phases like a density wave phase and a supersolid phase [SJG10,LGHL18,LDML08]. The model has been experimentally realized in ultracold atomic lattices [BDZ08]. Recent experiments [TLF + 19, BSW + 19, CPI + 19] report evidence for the creation of one of the more exotic predicted phases, a supersolid, a state of matter where an ordered material can flow without friction.
Simulating Bose-Hubbard models on a quantum computer will help elucidate the properties of these phases and possibly lead to the discovery of previously unknown phases. Quantum simulation may be especially important for studying behavior near the boundaries of quantum phase transitions, a regime that is notoriously difficult to study numerically on classical computers [HW17].
Here, we consider the 1D Bose-Hubbard model with periodic boundary conditions, expressed aŝ where t is the hopping term, U is the on-site site energy, and µ is the chemical potential. The parameters dictating the size of the problem are the number of levels d per site as well as the number of particles N .
Shifted one-dimensional QHO. We consider this Hamiltonian because it is one of the simplest prototypes for encoding a bosonic system into qubits. The Hamiltonian is directly relevant to chemistry, where it is equivalent to calculating the Franck-Condon factors [MMS + 18, SH19] between two electronic transitions for a diatomic molecule. Franck-Condon factors determine the intensity of the transitions to vibrational eigenstates in the final electronic state. The Hamiltonian isĤ where ω is the QHO frequency and δ is the displacement (related to the commonly used Huang-Rhys factor). The parameter δ indirectly determines how many levels d will have appreciable intensity when considering overlaps with the unshifted QHO. Multidimensional Franck-Condon factors. In the field of physical chemistry, in order to calculate an absorption spectrum, energy transfer characteristics, or other properties for molecules larger than two atoms, Franck-Condon profiles must include transitions to all possible vibrational Fock states of non-negligible intensity. Each vibrational normal mode in the final electronic state is a separate bosonic mode, and there are M = 3N − 6(5) such modes in nonlinear (linear) molecules. One defines q s and p s to be vectors of position and momentum operators, respectively, where s is one of two electronic potential energy surfaces (PESs) A or B. See references [Huh11, HGP + 15, SH19] for more details. The unitary M -by-M Duschinsky matrix S transforms between the bases of the two PESs such that where Ω s = diag([ω s1 , ..., ω sM ]) In the multidimensional case, all four parameters δ, S, Ω A , and Ω B indirectly determine the bosonic truncation in each mode, out of which the dimensionless displacement δ usually has the greatest effect. Depending on these parameters, the required d in a given mode can be anywhere from 1 in the case of negligible shifts, to as high as 70 in molecules such as nitrogen dioxide [SH19]. Though we consider the harmonic case here, it is the anharmonic Hamiltonians, where higher-order Taylor series terms are added to (15), where a quantum computer will likely first outperform a classical one for small molecules [SH19, LKC06, HNRB10, MR15,PR17].
For real molecules, the Duschinsky matrix S tends to be sparse. In our resource count simulations, we assume that each row of S has a constant number of non-zero entries per row, k. We set k = 4 for the simulations discussed here. We additionally assume that each vibrational mode is modelled with the same number of levels d, though for most real molecules one would use a different d for each mode. We make the approximations in order to compare typical constant-factor savings between different mappings, independent of system size.
Spin-s Heisenberg model. Quantum Heisenberg models, commonly used to study theoretical questions in condensed matter physics and to model magnetism in real solids, consist of quantum spins placed on a lattice. Though most work has focused on spin-1 2 particles, some problems require the use of higher-spin systems. Interpreting nuclear magnetic resonance experiments [Lev08] can require modeling higher-spin systems. Notably, a Heisenberg model simulation using a high spin of spin-7 2 has been recently reported [LSGB + 16], in order to study the magnetic properties of the material GdRhIn 5 . There are many stable isotopes up with nuclear spin I = 7 2 , 4, 9 2 [Sto05], and with unstable isotopes reaching still higher I.
For this work, we consider a one-dimensional spin-s Heisenberg model with ferromagneticŜ zŜz interactions and a magnetic field in the transverse direction, described by Hamiltonian Boson sampling. There is strong evidence that simulating multi-mode linear optical devices is hard. The problem has been formalized and can be stated as follows: consider a linear optical device with M modes and in which N single photons are injected in the first N ≤ M modes. The photons are detected in the output modes with occupation numbers that vary from repetition to repetition of the experiment. The task is to sample from this output distribution and, under the assumption that N M and that the linear optical device has been chosen at random (see [AA14] for details), together with some common theoretical conjectures, it has been shown that such a task is classically hard to perform. Such devices are composed of two kinds of optical elements: single-mode phase shifters and two-mode beam splitters. Here we consider the possibility of implementing the same dynamics not with optical systems, but using qubits. This approach may find utility, for example in testing computational complexity conjectures, if advanced qubit-based quantum computers come to fruition before advanced dedicated bosonic devices. A recent work discussed mappings to qubits for simulating quantum optics [Sab19].
In contrast to the other examples in this section, the dynamics of boson sampling is described in a stroboscopic way with layer-by-layer unitary operations, instead of being described by the time evolution of a large Hamiltonian. Nevertheless, the results of the previous sections can be readily applied when we describe the unitary operations in terms of their generators, in practice providing the effective Hamiltonians governing the evolution in each optical element. For single-mode phase shifter one has:ĥ p.s. =â † iâ i (17) while for two-mode beam splitter one has:ĥ b.s. =â † iâ j +â † jâ i (18) with indices i, j labelling the different optical modes. Note that this implementation is different from the previous Hamiltonian simulation examples, because the goal is to implementÛ BosonSampling = exp(−iθ RĥR ) exp(−iθ R−1ĥR−1 ) · · · exp(−iθ 1ĥ1 ) (19) whereĥ 1 , · · · ,ĥ R are the ordered list of effective Hamiltonians for each gate considered. Effectively, we produce quantum circuits that, after encoding the problem to qubits, performs a first-order Trotterization on each separate unitary.