Mapping second-quantized bosonic modes onto qubits can be done using various encodings. Each of those have their pro et contra. Fewer qubits — more gates, and vice versa. Encoding an $N$-level bosonic mode with unary mapping requires $(N+1)$ qubits, with creation and annihilation operators being represented by $N$ Pauli operators. Using binary/gray code results in using $\log_2N$ qubits, with operators containing $N^2$ Pauli terms. A good overview is given here.
Now, let us move to practical applications and do some trotterization. Unlike fermions, not only we are going to lose precision due to the finiteness of the Trotter number — an additional problem with bosons is that their operator commutation relations cannot be satisfied exactly. Unlike the former, the latter circumstance poses a serious problem even for variational methods, for it causes WF leakage out of qubit subspace representing physical states.
Consider a Fock state with just two bosonic modes $|1;0\rangle$. For $N=1,2$ it is mapped onto qubits in the following ways: $$ \begin{array}{|l|l|l|l|l|} \hline & N=1 & N=1 & N=2 & N=2 \\ \hline & \text{Binary} & \text{Unary} & \text{Binary} & \text{Unary} \\ \hline |0;1\rangle & |01\rangle & |0110\rangle & |0010\rangle & |001010\rangle \\ \hline |1;0\rangle & |10\rangle & |1001\rangle & |0010\rangle & |010001\rangle \\ \hline \end{array} $$ For the unary mapping, I'm using the convention from (23) here.
Let us start from a qubit version of the state $|0;1\rangle$, and try to create its superpostion with $|1;0\rangle$. We could do it using $U_1=\exp(i(a_1^\dagger a_2 + a_2^\dagger a_1))$ (aka time evolution) or using $U_2=\exp(a_1^\dagger a_2 - a_2^\dagger a_1)$ (aka UCC). For our purposes, it does not really matter. For definiteness, let us choose $U_2$ (this would keep amplitudes real if not Qiskit...). For $N=1$ both encodings give identical results: $$ U_2 |01\rangle = 0.54 |01\rangle -0.84 |10\rangle $$ $$ U_2 |0110\rangle = 0.54 |0110\rangle -0.84 |1001\rangle $$ But take a look at what happens as we move to $N=2$: $$ U_2 |0010\rangle = (0.01942407020071564-0.01942407020071514j) |0000\rangle + (0.3018468699435287-0.30184686994353027j) |0010\rangle + (-0.4175689014626733+0.4175689014626758j) |0101\rangle + (-0.3664436305642885+0.36644363056428975j) |0111\rangle + (-0.24189290092458382+0.2418929009245852j) |1000\rangle + (-0.017433459057134393+0.017433459057134143j) |1010\rangle + (0.026209883033085912-0.026209883033085773j) |1101\rangle + (0.20090642541053194-0.2009064254105325j) |1111\rangle $$ $$ U_2 |001010\rangle = (-0.399515044839777-0.39951504483977407j) |001010\rangle + (-0.5803699021310753-0.5803699021310711j) |010001\rangle + (-0.02732498912217729-0.027324989122177044j) |010100\rangle + (0.0015213913160944605+0.0015213913160944115j) |100010\rangle + (0.05259767375695154+0.05259767375695122j) |111001\rangle + (-0.006559016720416016-0.00655901672041582j) |111100\rangle $$
See the difference? For the unary mapping, the dominant amplitudes are those from the table, corresponding to physical states. But for the binary mapping we are in a real trouble (same happens for the Gray code).
BTW, an interesting difference between binary/gray and unary mappings is the following: in the former case, all the qubit states represent some physical states, while in the latter case this is not true. This, however, by itself does not imply that one would benefit from using binary/gray encoding.
QUESTION
What are the known approaches to dealing with this problem? While for the unary mapping, we could probably tolerate the error, it seems that in the case of binary/Gray encodings the situation is pretty sad.
One obvious suggestion would be to somehow project the outcome of trotterization onto the "physical" subspace using ancillas but I don't how practical this can be.
UPDATE
Note that the situation with time evolution and VQE is quite different. For simulating time evolution, given enough gates, we could probably choose the time step to be small enough to make errors (arising due to the wave function leakage out of the physical subspace) arbitrary small. However, for VQE we want short circuits (of order of a single time evolution step), which implies having large parameters. Therefore, ONE CANNOT EFFICIENTLY SIMULATE BOSONS USING VQE AND UNITARY COUPLED CLUSTER, where by "efficiently" I mean "in a way that may potentially lead to quantum advantage". Therefore... we have to rely on sparse VQE!
