Resource-efficient digital quantum simulation of $d$-level systems for photonic, vibrational, and spin-$s$ Hamiltonians
Nicolas P. D. Sawaya Tim Menke Thi Ha Kyaw Sonika Johri Alán Aspuru-Guzik Gian Giacomo Guerreschi























nicolas.sawaya@intel.com


Department of Physics, Harvard University, Cambridge, MA 02138, USA
Research Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA










gian.giacomo.guerreschi@intel.com





































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- systems. Here, we instead consider encodings of d-level (i.e. qudit) 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 patterns, bosonic truncation, and other properties of local operators. Finally, we obtain resource counts for five common Hamiltonian classes used in physics and chemistry, while modeling 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, although 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.


















































Resource-efficient digital quantum simulation of d-level systems for photonic, vibrational, and spin-s Hamiltonians
    Gian Giacomo Guerreschi
    December 30, 2023
===================================================================================================================

























§ 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- systems <cit.> or fermionic systems <cit.>. The latter case is relevant for simulations of chemical electronic structure <cit.>, nuclear structure <cit.>, and condensed matter physics <cit.>. 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.





However, for a large subset of quantum physics problems, important roles are played by components that are d-level particles (qudits) with d>2, including bosonic fundamental particles <cit.>, vibrational modes <cit.>, spin-s particles <cit.>, or electronic energy levels in molecules <cit.> 
and quantum dots <cit.>. 
Accordingly, several qubit-based quantum algorithms were recently developed for efficiently studying some such processes, including nuclear degrees of freedom in molecules <cit.>, the Holstein model <cit.> and quantum optics <cit.>.

















In principle, there are combinatorially many ways to map a quantum system to a set of qubits <cit.>.
Mapping a d-level system to a set of qubits 
may be done by assigning an integer to each of the d levels and then performing an integer-to-bit mapping.
Some consideration of d-level-to-qubit mappings has been published in the very recent literature, primarily for truncated bosonic degrees of freedom <cit.>, but this is still an unexplored area of theory especially in regards to determining which encodings are optimal for which problem instances. 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. 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 <ref>.
Similar considerations have been studied for fermionic mappings <cit.>.
For near- and 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 <ref>.
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 <ref> shows the relative two-qubit operation requirements for a set of five prominent physics and chemistry problems (defined in Supplementary Section VI). 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 <ref>.
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 <ref>. After outlining how to map d-level operators to qubit operators for any arbitrary encoding in Section <ref>, we consider resource counts for most standard local operators used in bosonic and spin-s Hamiltonians in Section <ref>. In Section <ref> we present circuits for converting between encodings and enumerate their resource counts. Finally, in Section <ref> we obtain relative resource estimates for various commonly studied Hamiltonians in theoretical physics and chemistry. Conclusions are summarized in Section <ref>. 
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.






§ RESULTS







 §.§ 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 
ℛ^enc(l), which is a bit string on N_q bits, x_N_q-1⋯ x_1 x_0. To specify the value (0 or 1) of a particular bit i in the encoding, we use the notation 
ℛ^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

    l ↦ x_0 2^0 + x_1 2^1 + x_2 2^2 + ...


This straight-forward mapping has been used for qubit-based quantum simulation of bosons previously <cit.>. 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 <cit.>.
Its defining feature is that the Hamming distance d_H between nearest-integer bit strings is always 1, formally d_H(ℛ^Gray(l),ℛ^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 <ref>).
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( ℛ^Gray(l) ℛ^Gray(l+1) )=2 to be used for pentadiagonal matrices like q̂^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^N_q states of the Hilbert space are possible.
In the unary encoding (also called one-hot <cit.>, single-excitation subspace <cit.>, or direct mapping <cit.>), 
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 <cit.> under different names.
The unary encoding makes less efficient use of quantum memory, but it will become clear below that it 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/glog_2 (g+1). 


Examples of the block unary encoding are given in Table <ref>.
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.






Bitmask subset. Because some encodings do not make use of the full Hilbert space, it will be useful to define the subset of bits that is necessary to determine each integer l.
For a given encoding we call this subset the bitmask subset, denoted C(l) where C_l ⊆{0,1,2,...,N_q-1}.
The bitmask subset for the SB and Gray encodings is always C^SB(l) = C^Gray(l) = {0,1,2,...,N_q-1}, since all N_q bits must be known to determine the integer value.
In the unary encoding, the bitmask subset is simply C^Unary(l)={l}, because if one knows that bit l is set to 1, then one knows the other bits are 0.
In the block unary code, the bitmask subset simply contains the bits that represent the current block.
For example, for the Block Unary Gray code with g=3 (see Table <ref>), C^BU[g=3](2) = {0,1} and C^BU[g=3](3) = {2,3}.































 §.§ Mapping d-level Matrix Operators to Qubits












Any operator for a d-level system can be written as


    Â = ∑_l,l'=0^d-1 a_l,l' | l  l' |,


where l and l' are integers labelling pairwise orthonormal quantum states. In this work we conceptualize the mappings primarily in terms of Fock-type encodings (or alternatively second quantization encodings) where each |l represents one level in the d-level system. However, the mapping procedure is identical to the one used for first quantization operators <cit.> that we briefly discuss. When dealing with bosonic degrees of freedom, one must choose an arbitrary level d at which to truncate, since in principle a bosonic mode may have an unbounded particle number. Choosing this cutoff such that truncation error is below a given threshold is an essential step that has been previously studied <cit.>, though it is beyond the scope of the current work.







In performing a mapping of any d-by-d matrix operator to a sum of Pauli strings, the following approach may be used.
For each term in the sum, one first assigns an integer to each level and then uses an arbitrary binary encoding ℛ to encode each integer:



    |l↦ |ℛ(l) = |ℛ(l;N_q-1)⋯|ℛ(l;1)|ℛ(l;0)↦ |x_N_q-1⋯ |x_0,   x_i ∈{0,1}.


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 . 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



    | l  l' | ↦⊗_i ∈ C(l) ∪ C(l') | x_i  x'_i |_i,



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 1| =  (  + i ) ≡σ̂^-,



    |1 0| =  (  - i ) ≡σ̂^+,



    |0 0| =  (I + ),



    |1 1| =  (I - ).



For a single term in equation (<ref>), the result is a sum of Pauli strings, 




    Â↦∑_k^P c_k⊗_j^N_qσ̂_kj,


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∈ {, , }∪{I}. Note that in this work the set of Pauli matrices is defined to exclude the identity.
















 §.§ Significance of the Hamming distance




It is useful to analyze encoding efficiency based on the Hamming distance between ℛ(l) and ℛ(l').
The Hamming distance, which we denote d_H^ℛ(l,l') ≡ d_H(ℛ(l),ℛ(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 (<ref>), 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 | ℛ(l) ℛ(l') |, one performs the following mapping:



    |x_0  x_0'|⊗⋯⊗|x_N  x_N'| ↦1/2^N (â_0+b̂_0) (â_1+b̂_1) ⋯ (â_N+b̂_N)



where subscripts denote qubit number, â∈{, I}, and b̂∈{± i,±} (according to equations (<ref>) through (<ref>)). Expanding the RHS of equation (<ref>) leads to an equation of form (<ref>).
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 <ref>, 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 |34| + |43| to a set of qubits.
In the SB encoding, equations (<ref>)–(<ref>) yield the following Pauli string representation:

    |34| + |43|  Std. Binary |011  100| + |100  011| 
         = 1/4 ( ^(2)^(1)^(0) + ^(2)^(1)^(0) + ^(2)^(1)^(0) - ^(2)^(1)^(0) ).
 

Using the Gray code, the Pauli string instead takes the form

    |34| + |43|  Gray |010  110| + |110  010| 
       = 1/4 ( - ^(2)^(1)^(0) + ^(2)^(0) - ^(2)^(1) + ^(2) ).
 


The Hamming distance between ℛ(3) and ℛ(4) is d_H^SB=3 in the former case and d_H^Gray=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.











































 §.§ 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 operator n̂_d=3 = diag [0,1,2] ↦3/2 I - ^(1) - ^(2). If one begins with n̂_d=3^2 = diag [0,1,4], the encoded operator is 


    n̂_d=3^2 = diag [0,1,4]  Unary5/2 I - ^(1) - 2 ^(2).


If one instead squares the already-encoded Pauli operator for n̂_d=3, this yields

    (3/2 I - ^(1) - ^(2))^2 = 7/2 I - 3/2^(1) + ^(1)^(2) - 3^(2).


Superscripts denote qubit number. Pauli operators (<ref>) and (<ref>) behave identically on the subspace of the unary encoding, though operator (<ref>) 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.










 §.§ Trotterization and gate count upper bounds







Hamiltonian simulation often consists of implementing the unitary operator


    Û(t) = exp(-i Ĥ t)


for some user-defined time-independent Hamiltonian Ĥ, where t is the evolution time and we have set ħ=1. Any Hamiltonian can be expressed as a sum of local Pauli strings such that 


    Ĥ = ∑_k c_k ⊗_j σ̂_kj = ∑_k ĥ_k,


which takes the same form as equation (<ref>) with the Pauli strings and their coefficients compacted into terms {ĥ_k}. In practice, Hamiltonian simulation can be performed using a Suzuki-Trotter decomposition


    Û(t) = exp(-i t ∑_k ĥ_k)
         ≈(∏_k exp(-i ĥ_k t/η) )^η
         = Ũ(t),


where the expression is exact in the limit of large η or small t <cit.>. The numerical studies of this work consider the encoding-dependent resource counts for equation (<ref>), for a subset of prominent physics problems. We focus on determining the resources required for simulating a single Trotter step.












There are several variations and extensions to the Hamiltonian simulation approach of equation (<ref>), including higher-order Suzuki-Trotter methods <cit.>, the Taylor  series algorithm <cit.>, quantum signal processing <cit.>, and schemes based on randomization <cit.>. Notably for the current work, recent results suggest that simple first-order Trotterization will have lower error for near- and medium-term hardware <cit.>, 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.
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(-i t ĥ_k)=exp(-i t c_k ⊗_j σ̂_kj) may be implemented using the well-known CNOT staircase quantum circuit shown in Figure <ref>.
If a qubit j is acted on by  or , additional single-qubit gates are placed on qubit j as shown in the figure.
H ≡ i R_x+z(π) denotes the Hadamard gate that changes between the Z- and X-basis, and i R_x+y(π) converts between the Z- and Y-basis.
To exponentiate a single Pauli string, the number of CNOT (CX) gates required is 

    n_cx(p) = 2(p-1),


where p is the number of Pauli operators ({,,}, excluding I) in the term to be exponentiated.












It is instructive to calculate upper bounds for entangling gate counts. Consider a simple Hermitian term α + α^* and a single diagonal element . 
Here, d_H denotes d_H(ℛ(l),ℛ(l')), which will depend on the chosen encoding, and d_H=0 in the case of a diagonal element. 
We define K as |C^enc(l) ∪ C^enc(l')|, the number of qubits in the relevant bitmask subsets. As we are considering products of two-term sums (equations (<ref>)–(<ref>)), the distribution of Pauli strings can be analyzed in terms of binomial coefficients. In the Supplementary Section I we show that 

    n_cx,UB[d_H(l,l'),K] = ∑_p=2^K 2^d_HK-d_H  p-d_H (2p-2),

where UB denotes upper bound and the  factor is not present for diagonal terms. Some resulting upper bounds for particular Hamming distances are



    n_cx,UB[d_H=0,K] = (2^K K-2(2^K) - 2)  
    
    n_cx,UB[d_H=1,K] =  (2^K K-2^K)   
    
    n_cx,UB[d_H=2,K] =  (2^K K) .



Again, the above expressions are for a single Hermitian element pair or a single diagonal term. In practice, because substantial gate cancellation is possible once the quantum circuit has been compiled, these upper bounds are not always directly applicable when choosing an encoding. 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.

Below, we study a common sparsity pattern, a tridiagonal real matrix operator B̂ with zeros on the diagonal, i.e. with matrix structure i | B̂ | j = ∑_k B_k ( δ_i,kδ_j,k+1 + δ_i,k+1δ_j,k). 
This is the sparsity pattern of several commonly used d-level operators, such as the bosonic position operator q̂.











In Table <ref> we show analytical upper bounds for three different levels of sparsity, derived in Supplementary Section I. We consider a single Hermitian pair, the B̂ operator with O(d) nonzero entries, and a dense matrix operator with O(d^2) nonzero entries. In Supplementary Section I we show that the upper bound of entangling gate counts for an arbitrary operator on K qubits is O(K4^K). Because K = log_2 d in the compact codes, this means that the compact codes will never require more than O(d^2 log d) entangling gates.



An important consequence is that, as the matrix density increases, the comparative advantage of the unary encoding decreases. The compact codes' upper bound both for a B̂ and for a fully dense operator are both O(d^2 log d), since this is the maximum upper bound. The unary encoding's upper bound of O(d^2) for fully dense operator is only slightly lower. And because actual gate counts for various matrix instances will be less than these upper bounds, it appears possible that compact codes might often be superior for dense matrices in both qubit count and gate count. However, because the most commonly used quantum operators tend to have O(d) density, it is likely that unary will most often be superior in gate counts, at least for Hamiltonian simulation. Many exceptions to these trends are shown in Section <ref>.



As a more concrete demonstration of typical operator scaling, we calculated numerical upper bounds for B̂ with increasing d. Qubit counts and upper bounds for B̂ are shown in Figure <ref>, as a function of d for different encodings. We first encoded the entire operator B̂ into a sum of Pauli strings before collecting and cancelling terms, leading to some favorable cancellations. 
Then we applied equation (<ref>).
There is roughly an inverted relationship between the qubit counts and the operation counts, because sparser encodings like the unary and block unary have smaller bitmask subsets but require more total qubits. 



The differing gate count upper bounds between the SB and gray encodings (Figure <ref>) are explained by Hamming distances.
Because all non-zero terms in B̂ 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 (<ref>) 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 cancelling of Pauli terms and circuit optimization.



































































































































 §.§ 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 







    Ô_l,l = ∑_i=1^log_2d k_i R^SB(l;i)



where R^SB(l;i) ∈{0,1}. 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


    n̂ = diag [0, 1, 2, 3, ⋯]


and any linear combination a n̂ + b I 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, 


    n̂ = 2^0 ^(0) + 2^1 ^(1) + 2^2 ^(2) + ⋯


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 operator n̂ with at least d=11, it is most efficient simply to implement the SB encodings of n̂ with d=16 instead. A simple example illustrates this point. The standard operator for the bosonic number operator with truncation d=3 is


    n̂_d=3 = diag [0,1,2] Std. Binary3/4I + 1/4^(0) - 3/4^(0)^(1) - 1/4^(1)


while that for d=4 is


    n̂_d=4 = diag [0,1,2,3] Std. Binary3/2I + 1/2^(0) - ^(1)



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 <ref>, where it is clear that the SB mapping is superior when d is a power of 2. The right panel gives gate counts for operators such as n̂^2, where it is again advangatgeous for d to be a power of 2, though entangling gates are still required. As is also clear from the right panel, the square of a DES or DBD operator is in general not DBD.













 §.§ Local operators
























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 operator q̂=1/√(2) (â_j^† + â_j) is tridiagonal with zeros on the diagonal:


    q̂ = 1/√(2)[    0    1    0    0    …;    1    0 √(2)    0    …;    0 √(2)    0 √(3)    …;    0    0 √(3)    0    …;    ⋮    ⋮    ⋮    ⋮    ⋱ ]





This means that the square of q̂, often used in vibrational and bosonic Hamiltonians, is pentadiagonal but with zeros for terms where |l-l'|=1:


    q̂^2 = [      1      0 √(1·2)      0      0      …;      0      3      0 √(2·3)      0      …; √(1·2)      0      5      0 √(3·4)      …;      0 √(2·3)      0      7      0      …;      0      0 √(3·4)      0      9      …;      ⋮      ⋮      ⋮      ⋮      ⋮      ⋱ ].














Notably, |l-l'| 
for non-zero entries is either 0 or 2, making the Gray code less useful for this operator. The momentum operator p̂=i/√(2) (â_j^† - â_j) and its square p̂^2 have the same sparsity patterns as q̂ and q̂^2 respectively. The number operator n̂=â^†â of equation (<ref>) is diagonal and DBD, which leads to efficient SB mappings as discussed in Section <ref>. Finally, we study the two-site bosonic interaction operator 


    â_i ^†â_j + â_i â_j ^†.
















































In order to consider spin Hamiltonians such as 
Heisenberg models<cit.>, we encode 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 <cit.>


    l |  | l'   = ħ(s+1-l)δ_l,l'



    l |  | l'   = ħ/2(δ_l,l'+1 + δ_l+1,l') √((s+1)(l+l'-1)-ll')



    l |  | l'   = iħ/2(δ_l,l'+1 - δ_l+1,l') √((s+1)(l+l'-1)-ll')


where δ_α,β is the Kronecker delta. The  are DBD operators, while  and  are tridiagonal with zeros on the diagonal, the same sparsity pattern as bosonic p̂ and q̂ 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 <cit.> 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 as

    X̃_FQ = ∑_i=0^N_x-1 x_i | x_i  x_i |

and N_x is the number of discrete position points such that


    x_i = (i-N_x/2)Δ, i ∈ [0,N_x-1]

where Δ is chosen such that the desired highest-order Hermite-Gauss function is contained within (x_0,x_N_x-1). Advantages and disadvanteges are discussed in the Supplementary Section V. 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 that X̃_FQ is a DBD operator.















Quantum circuits for approximating the exponential of each operator and for each d were compiled and then optimized using the procedure given in the Supplementary Section II. 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 operators q̂, q̂^2, and  are plotted in Figure <ref>. 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 <cit.>. 



Comparing the tridiagonal operator q̂ with the upper bounds given in Figure <ref> 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 ) show the Gray encoding outperforming SB, though after optimization the advantage is less pronounced. In contrast, for the pentadiagonal q̂^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 of q̂^2, given in equation (<ref>). Also notable is the apparent dip in operation count at d=8, due to the fact that the diagonal of q̂^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 implementing q̂^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 for q̂ and d=4 and 8 for q̂^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  are plotted in Figure <ref>. 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 <ref>. Because  is diagonal, the unary always requires just d single-qubit rotations and no entangling gates. The two-particle operator  displays similar trends to .


As expected, the Gray code is usually superior to SB for the tri-diagonal , 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 low-d 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 <ref>) 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. 
In such cases (q̂ and  for d=9; n̂^2 ∼Ô^2 in Figure <ref> for several values;  and  for s=2), there may be a particular hardware budget would require this encoding for its particular memory/operation trade-off. Such highly specific hardware budgets seem 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 <ref>. As will be discussed in Section <ref>, it is often beneficial to Trotterize different parts of the Hamiltonian in different encodings, if the cost difference outweighs the overhead of conversion.



































 §.§ 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 operator n̂ is usually more efficient in SB, while the hopping term b̂_i ^†b̂_i+1 + h.c. is usually more efficient in the Gray encoding (see Figure <ref>). Here we show that the cost of converting from one encoding to another is often substantially less than the difference in resource efficiency between two encodings, which means that it can be advantageous to continually be 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 <ref>. 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 <cit.> as shown in Figure <ref>.











































































































































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 <cit.>, in Figure <ref> 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 <ref>, we provide the converter circuit resource count for a general d-level truncated quantum system, following Figure <ref>. Resource counts assume a decomposition of CSWAP into Clifford+T gates <cit.>. A general algorithm to build the converter circuit can be seen in the Algorithm <ref>.
· 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 <ref> and example circuits for d=5 and 7 are given in the Supplementary Section III.

















    
    





 
 
 
 
 



    



 

  
    
  
    
  
        
   
        
   
    
  

















































For completeness, we constructed a circuit for converting between SB and block unary, for BU^SB_g=3, shown in the Supplementary Section IV. We do not further analyze BU conversions, as block unary is expected to have limited utility, and even then it will usually be the case that decoherence times are too low to allow for conversions (see Section <ref>).




























 §.§ 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 <cit.>, multidimensional molecular Franck-Condon factors <cit.>, a spin-s transverse-field Heisenberg model <cit.>, and simulating Boson sampling <cit.> 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. Supplementary Section VI gives a more thorough overview of these problems.

Assuming that d remains constant as the particle number increases, 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 digital quantum simulation. In Figure <ref>, an example of encoding scheme (v) is shown. These encoding schemes do not directly correspond to the `scenarios' discussed below; the scenarios denote 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 again 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 include n̂, n̂^2, and q̂^2 (Figures <ref> and <ref>). 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 <ref>). 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 <ref>. 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 Supplementary Section VII. 


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 <ref>, 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  and  in Figure <ref>.

Scenario C applies when unary is the superior encoding and it is still considered the best encoding even if one repeatedly unravels and compacts to preserve memory (as in Figure <ref>). 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 <ref>), 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 Supplementary Section VII. 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.







































































§ DISCUSSION



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. Notably, 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. To the best of our knowledge, the Gray code had not been previously used in Hamiltonian simulation.

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,” each of which points to a different optimal encoding and simulation strategy. The simulation scenario depends on which encodings require the fewest operations, on whether interconverting between mappings is worth the additional cost, and on 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 simulation strategy 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 <cit.>. Additionally, we assumed all-to-all connectivity in this work, which tends to be a feature of ion trap quantum computers <cit.>.
But other quantum hardware types require one to consider the topology of the qubit connections and implementation of SWAP gates <cit.>, 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.












 








§ DATA AVAILABILITY:
 The data that support the findings of this study are available upon reasonable request.





§ ACKNOWLEDGEMENTS


A.A.-G., T. H. K. and  T. M. acknowledge funding from Intel  Research and Dr. Anders G. Frøseth. A.A.-G. also acknowledges support from the Vannevar Bush Faculty Fellowship and the Canada 150  Research Chairs Program. The authors declare that there are no competing interests.



§ AUTHOR CONTRIBUTIONS

N.P.D.S. conceived of and designed the study, developed the theory, and generated and analyzed the data. 
T.M. developed encoding conversion schemes, analyzed data, and created conceptual figures. 
T.H.K. developed encoding conversion schemes. 
S.J. optimized all quantum circuits. 
A.A.-G. supervised work on encoding conversions. 
G.G.G. conceived of the encoding conversions section and developed theory. 
N.P.D.S., G.G.G., T.M., and T.H.K. wrote the manuscript. 
All authors edited the manuscript.



§ COMPETING INTERESTS

The authors declare no competing interests.




























§ FIGURE LEGENDS




FIG. 1. 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 of q̂^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.


FIG. 2. 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. The five Hamiltonians are 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 <ref>.


FIG. 3. Canonical quantum circuits used to exponentiate Pauli strings on a universal quantum computer. One needs 2(p-1) two-qubit gates for such an operation, where p 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.



FIG. 4. Numerical upper bounds for resource counts of implementing one Suzuki-Trotter step of a d-by-d real Hermitian matrix operator B̂, where B̂ is tridiagonal with zeros on the diagonal. Top: Qubit counts for mappings considered in this work. BU^Gray_g stands 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 of B̂. This is the sparsity pattern of canonical bosonic position and momentum operators as well as the S_x 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 (<ref>). Notably,  encodings with higher qubit counts tend to have lower upper bounds for gate counts, and vice-versa.


FIG. 5. Entangling gates for a single Suzuki-Trotter step of an arbitrary diagonal evenly-spaces (DES) operator Ô (left), and its square (right). Gate counts are for optimized quantum circuits. DES operators are a subset of the diagonal binary-decomposable (DBD) operator class. Because it is diagonal, the unary code always requires only single-qubit operations. When log_2 d is an integer, the SB code requires no entangling gates and just log_2 d single-qubit operations, making it the most efficient encoding (in terms of both qubits and operations). DBD operators are a common operator class, encompassing e.g. the bosonic number operator n̂ and the spin-s operator Ŝ_z.


FIG. 6. CNOT gate counts for optimized circuits of the bosonic position operator (first row), position operator squared (second row), and two-site bosonic interaction operator (third row). Plots on the left correspond to enlargements of the dotted boxes in the plots on the right. There are several notable trends and anomalies: (a) Though unary usually requires the fewest operations as d increases, there are several operators for which the Gray or SB code is more efficient than the unary in both qubit and operation count. This occurs most pronouncedly at values such as d=4, 7, and 8. (b) The Gray code is usually more efficient than SB even after circuit optimization, especially for operators composed of tridiagonal operators (top and bottom rows). (c) The Gray code's advantage is either less pronounced or disappears for q̂^2, because q̂^2 is a pentadiagonal operator, for which the unity Hamming distance of the Gray code is less useful. (d) The reduction in operation count for q̂^2 at values such as d=8, 24, or 32 occur because the diagonal of q̂^2 is DBD.


FIG. 7. CNOT gate counts of optimized Suzuki-Trotter circuits, for approximating the exponentials of the spin-s operators shown. There is no clean overall trend for  and  (except that Gray tends to out-perform SB), highlighting the need to study encodings thoroughly for each new use case. Notably, because Ŝ_z is diagonal binary-decomposable, for values of s= 3/2 , 7/2 (which are 4- and 8-level system respectively) the SB code requires both fewest operations and fewest qubits. 



FIG. 8. Schematic showing the utility of swapping between encodings. 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 . 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 <ref>. 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.


FIG. 9. Quantum circuit for SB to Gray conversion.


FIG. 10. Circuit to convert the SB representation into the unary representation.
Every CSWAP is accompanied by a CNOT.
Modifications are required when d is not a power of two, as discussed in the main text and Supplementary Section III.

























