Can effective descriptions of bosonic systems be considered complete?

Francesco Arzani,¹ Robert I. Booth,²,³ and Ulysse Chabaud¹

¹DIENS, École Normale Supérieure, PSL University, CNRS, INRIA, France

²University of Edinburgh, United Kingdom

³University of Bristol, United Kingdom

Bosonic statistics give rise to remarkable phenomena, from the Hong-Ou-Mandel effect to Bose-Einstein condensation, with applications spanning fundamental science to quantum technologies. Modeling bosonic systems relies heavily on *effective descriptions*: typical examples include truncating their infinite-dimensional state space and restricting their dynamics to a simple class of Hamiltonians, such as polynomials of canonical operators, which are used to define quantum computing over bosonic modes. However, many natural bosonic Hamiltonians do not belong to this simple class, and some quantum effects harnessed by bosonic computers inherently require infinite-dimensional spaces, questioning the validity of such effective descriptions of bosonic systems. How can we trust results obtained with such simplifying assumptions to capture real effects?

Driven by the increasing importance of bosonic systems for quantum technologies, we solve this outstanding problem by showing that these effective descriptions do in fact capture the relevant physics of bosonic systems. Our technical contribution is twofold: firstly, we prove that any physical, bosonic unitary evolution can be strongly *approximated* by a finite-dimensional unitary evolution; secondly, we show that any finite-dimensional unitary evolution can be generated *exactly* by a bosonic Hamiltonian that is a polynomial of canonical operators. Beyond their fundamental significance, our results have implications for classical and quantum simulations of bosonic systems, they provide universal methods for engineering bosonic quantum states and Hamiltonians, they show that polynomial Hamiltonians do generate universal gate sets for quantum computing over bosonic modes, and they lead to an infinite-dimensional Solovay-Kitaev theorem.

# I. INTRODUCTION

Bosonic systems, often referred to as continuous-variable quantum systems in the context of quantum computing, encompass a number of promising setups for quantum information processing. These include photonic setups, which have led to some of the first demonstrations of quantum computational speedup [1–4]. Moreover, bosonic error-correcting codes [5], specifically those theorized in 2001 by Gottesman, Kitaev and Preskill (GKP) [6], have recently led to the first demonstration of quantum error-correction beyond the break-even point [7].

The description of bosonic systems requires a Hilbert space of infinite dimension, since it must support position and momentum operators verifying the canonical commutation relation $[\hat{q}, \hat{p}] = i\hat{I}$ [8]. While this large state space allows us to robustly encode quantum information as in GKP codes [6], it also makes modeling bosonic systems and computations challenging. This has motivated the use of mathematically simpler *effective descriptions*, aiming to capture their important features. Here we consider two widely adopted simplifications: *effective dimension*—simplifying the state space—and *effective Hamiltonians*—simplifying the dynamics.

Choosing an effective dimension for a bosonic system amounts to imposing a cut-off of its infinite-dimensional Hilbert space leading to a finite-dimensional one, based for instance on an energy bound. It allows us to simulate these systems on classical or quantum computers, up to a truncation error [9, 10].

Alternatively, one may restrict to a specific class of

effective Hamiltonians. A standard choice for bosonic computations is the set of polynomials in the canonical bosonic operators [11], which we refer to as *polynomial Hamiltonians* for brevity. Polynomial Hamiltonians are ubiquitous in physics, such as in the Bose-Hubbard model of interacting spinless bosons on a lattice [12], in the $\phi^4$ theory [13] in quantum field theory, or in the multipole expansion of the dielectric polarization in nonlinear optics [14], to a name a few. In particular, the standard model of quantum computations with continuous variables was defined by Lloyd and Braunstein based on unitary gates generated by polynomial Hamiltonians [11]. In particular, *universality* for continuous-variable quantum computing is defined as the ability of a gate set to approximate any polynomial Hamiltonian evolution [11, 15].

While these effective descriptions greatly simplify the modeling of bosonic computations, they also have drawbacks: recall that descriptions based on a finite effective dimension cannot reproduce the bosonic canonical commutation relation exactly. In addition, there are fundamental differences between the sets of quantum correlations in finite and infinite-dimensional spaces [16, 17]. These features suggest that some bosonic quantum effects require infinite-dimensional spaces and may not be accounted for when restricted to a finite effective dimension.

Similar doubts can be raised about the ability of effective Hamiltonians to faithfully describe bosonic phenomena, as many bosonic systems are not naturally described by polynomial Hamiltonians. For instance, Hamiltonians of Josephson junctions in superconducting systems include terms of the form $\cos \hat{\delta}$, where $\hat{\delta}$ is a canonical
phase-difference operator [18]. Here, the function cos cannot be replaced by a polynomial—such as its truncated Taylor expansion—because the operator $\hat{\delta}$ is unbounded [19]. Moreover, in the context of bosonic computations, it is in fact not known whether the definition of universality, based on the ability to reproduce any evolution generated by polynomial Hamiltonians [11], is sound: by this we mean that it is far from obvious that universal continuous-variable quantum computers as defined by Lloyd and Braunstein satisfy plausible requirements that one would expect by analogy with their discrete-variable counterparts. One example is universal state preparation, loosely defined as the ability to approximately prepare any quantum state from a given reference state: as far as we know, there could be inaccessible regions of the Hilbert space when one restricts to using gates generated by polynomial Hamiltonians. This can be phrased as a question about universal control and was first formalized in this context as an open problem in [20]. Beyond this universal controllability problem, it is also highly nontrivial whether any physical unitary evolution may be approximated to arbitrary precision by a sequence of gates generated by polynomial Hamiltonians. By contrast, the definition of universality for discrete-variable quantum gate sets is based on the property that they allow us to approximate any unitary operation to arbitrary precision [21]. What is more, the Solovay–Kitaev theorem guarantees that these finite-dimensional universal gate sets provide efficient approximations of any unitary operator on a constant number of subsystems [22]. Such a theorem is missing in the general infinite-dimensional setting [23] and would be fundamental for understanding the computational complexity of bosonic systems, because it would imply that the complexity of bosonic computations is independent of the choice of gate set [24]. A pre-requisite for this result is to ensure that the definition of universality for continuous-variable quantum computing is sound: at the moment, it is unclear whether effective polynomial Hamiltonians provide a solid foundation for modeling bosonic computations.

Are these drawbacks simply mathematical artifacts, or should they be taken into account when describing bosonic systems? In other terms: *Can effective descriptions of bosonic systems be considered complete?*

Here we provide a positive answer to this question. As we shall demonstrate, this is not only a non-trivial conceptual question, but also a highly fruitful one that leads to several applications. Our technical contribution is twofold:

Firstly, we show that any physical unitary evolution can be *approximated* to arbitrary precision by a finite-dimensional evolution. The proof is based on an operationally meaningful and mathematically rigorous notion of approximation. Moreover, the corresponding cut-off unitary evolution and its effective dimension may be computed explicitly—in particular, the effective dimension relates to the amount of energy needed to implement the evolution (Theorem 1).

Secondly, we show that polynomial Hamiltonians can realize any finite-dimensional Hamiltonian *exactly*. We provide a constructive proof, based on Lagrange interpolation polynomials (Theorem 2). Together with the previous result, this shows that polynomial Hamiltonians can generate any physical unitary evolution *approximately* to arbitrary precision.

From a conceptual standpoint, our first technical result identifies a clear trade-off between space and energy when simulating bosonic systems, be it with classical or finite-dimensional quantum methods. Our second technical result provides universal methods for state and operator engineering using polynomial Hamiltonians, which can be used to simulate the dynamics of complex quantum systems using bosonic architectures. Furthermore, it resolves the open question about the universal controllability of polynomial Hamiltonians in infinite dimension from [20] (Corollary 1). Finally, these two technical contributions together show that any physical unitary evolution of a bosonic system may be approximated to arbitrary precision by a sequence of gates generated by polynomial Hamiltonians in a mathematically rigorous sense (Theorem 3). This solves the universality problem in continuous-variable quantum computing, showing that the definition of universality from [11] is sound, and further leads to a Solovay–Kitaev theorem in infinite dimension for specific universal gate sets generated by polynomial Hamiltonians (Theorem 4).

The rest of the paper is structured as follows. We set the stage and introduce the necessary notation in Sec. II. We show how to approximate physical unitary evolutions over infinite-dimensional Hilbert spaces by finite-dimensional ones in Sec. III. In Sec. IV, we show how to realize finite-dimensional Hamiltonians using polynomial Hamiltonians. Sec. V is devoted to a discussion of the implications of our results for bosonic computations. We conclude in Sec. VI. The Methods section provides an overview of the technical tools used in the proofs of our main results, the details of which are deferred to the Supplementary Information.

## II. BOSONIC QUANTUM INFORMATION THEORY

Quantum information theory with bosonic systems typically involves modelling unbounded operators over infinite-dimensional Hilbert spaces with continuous spectrum, and is thus often referred to as continuous-variable quantum information theory [25–27]. In this framework, systems corresponding to individual quantum harmonic oscillators are called qumodes, or simply modes, and quantum states of *m*-mode bosonic systems are elements of an infinite-dimensional Hilbert space $\mathcal{H}^{\otimes m}$, with each single-mode Hilbert space spanned by a countably infinite (Fock) basis $\{|n\rangle\}_{n \in \mathbb{N}}$.

Creation and annihilation operators $\hat{a}^\dagger$ and $\hat{a}$ for each mode are defined by their action on the Fock basis
as $\hat{a}^\dagger|n\rangle = \sqrt{n+1}|n+1\rangle$, $\hat{a}|n+1\rangle = \sqrt{n+1}|n\rangle$ and $\hat{a}|0\rangle = 0$. The (unbounded) canonical bosonic operators are then given by $\hat{q} = \frac{1}{\sqrt{2}}(\hat{a}+\hat{a}^\dagger)$ and $\hat{p} = \frac{1}{i\sqrt{2}}(\hat{a}-\hat{a}^\dagger)$ and satisfy the canonical commutation relation $[\hat{q}, \hat{p}] = i\hat{I}$, where $\hat{I}$ is the identity operator. Polynomial Hamiltonians over $m$ modes are the Hermitian operators of the form $P(\hat{q}_1, \hat{p}_1, \dots, \hat{q}_m, \hat{p}_m)$, where $P$ is a polynomial and where $(\hat{q}_1, \hat{p}_1, \dots, \hat{q}_m, \hat{p}_m)$ are the canonical operators of the modes $1, \dots, m$. The number operator is defined as $\hat{n} = \hat{a}^\dagger\hat{a}$, and its expectation value is the average particle number of a bosonic mode, which we refer to as energy hereafter. For all $N \in \mathbb{N}$, we denote by $\mathcal{H}_d := \text{span}\{|n\rangle\}_{0 \le n \le d}$ the subspace of states with a number of particles at most $d$.

In some sense, not all states in the Hilbert space are valid physical states, even if they are normalized: for instance, some may have infinite energy, such as $\sqrt{6}\sum_{n>1} \frac{1}{\pi n}|n\rangle$. Similarly, there are normalized states with infinite fluctuations in either position or momentum or both. To avoid inconsistencies and ensure that the action of polynomial Hamiltonians leads to physical states, it is common to define a set of physical states as a dense subspace of the Hilbert space known as the Schwartz space [19]—informally, the set of states with bounded canonical operator moments, with a position wave function and all of its derivatives decaying sufficiently fast at infinity, which we denote by $S \subset \mathcal{H}$ (see the Methods section). Similarly, denoting by $\mathcal{U}(\mathcal{H})$ the group of unitary operators over $\mathcal{H}$, we define

$$ \mathcal{U}(S) := \{ \hat{U} \in \mathcal{U}(\mathcal{H}) \mid \forall |\psi\rangle \in S, \hat{U}|\psi\rangle \in S \}. \quad (1) $$

These form a group of physical unitary operators, i.e., which map physical states to physical states. This set contains in particular all products of unitary operators generated by analytic functions of canonical operations [20], including polynomial Hamiltonians, and can be readily extended to the multimode setting. Physical unitary operators have the property that they map bosonic states with a finite number of particles to states of finite energy: for all $n \in \mathbb{N}$ and all $\hat{U} \in \mathcal{U}(S)$,

$$ E_{\hat{U}}(n) := \sup_{|\psi\rangle \in \mathcal{H}_n} \langle \psi | \hat{U}^\dagger \hat{n} \hat{U} | \psi \rangle < +\infty. \quad (2) $$

In that case, the supremum is a maximum and the quantity $E_{\hat{U}}(n)$ can be thought of as the maximal amount of energy involved when implementing the unitary evolution $\hat{U}$ of an initial state with a number of particles at most $n$.

We denote by $D(\rho, \sigma) = \frac{1}{2}\|\rho - \sigma\|_1$ the trace norm distance between density operators [21], which is used to define the diamond norm for quantum maps. The diamond norm induces a distinguishability measure between maps that is often too stringent for quantum mechanical applications involving unbounded operators, as physical unitary channels may be always maximally separated for this norm, even though they may have similar effects on all states below a certain energy [28]. Instead, the

strong topology of pointwise convergence provides a suitable notion of continuity in infinite dimensions, which is implied by the closeness in energy-constrained diamond norm (see the Methods section), defined as [29]:

$$ \|\mathcal{E}\|_{\diamond}^{E} := \sup_{\text{Tr}[\rho(\hat{H}\otimes\hat{I})]\le E} \|(\mathcal{E} \otimes \text{id})\rho\|_{1}. \quad (3) $$

Here, $\mathcal{E}$ is a Hermitian-preserving map, $\rho$ is a density operator over $\mathcal{H} \otimes \mathcal{H}'$, with $\mathcal{H}'$ isomorphic to $\mathcal{H}$, and $\hat{H}$ is a Hamiltonian over $\mathcal{H}$ with ground state energy equal to 0 used to specify the energy bound $E > 0$. We employ this topology hereafter, with $\hat{H} = \hat{n}$ the number operator.

### III. EFFECTIVE DIMENSION OF PHYSICAL UNITARY EVOLUTIONS

Having introduced the necessary notation, we show in this section that any physical unitary channel can be approximated to arbitrary precision by a finite-dimensional unitary channel, where the quality of the approximation is measured by the energy-constrained diamond norm. This is captured by the notion of approximate effective dimension, which we define as follows:

**Definition 1 (Approximate effective dimension).** Let $\hat{U}$ be a single-mode unitary operator, let $E \ge 0$ be an energy parameter, and let $\epsilon > 0$ be an approximation parameter. Let also $d \in \mathbb{N}$. The unitary operator $\hat{U}$ has $(E, \epsilon)$-approximate effective dimension $d+1$ if there exists a unitary operator $\hat{V}_d$ over $\mathcal{H}_d$ which approximates $\hat{U}$ in the following sense: for any unitary operator $\hat{V}'$ over $\text{span}(|n\rangle)_{n>d}$, denoting $\hat{V} = \hat{V}_d \oplus \hat{V}'$, $\mathcal{V} = \hat{V} \cdot \hat{V}^\dagger$ and $\mathcal{U} = \hat{U} \cdot \hat{U}^\dagger$,

$$ \|\mathcal{U} - \mathcal{V}\|_{\diamond}^{E} \le \epsilon. \quad (4) $$

When that is the case, we say that $\hat{V}_d$ is a $(d+1)$-dimensional $(E, \epsilon)$-approximation of $\hat{U}$.

When a unitary dynamics has a finite approximate effective dimension, it can effectively be modeled by a finite-dimensional unitary evolution. We now give our first main result:

**Theorem 1 (Effective dimension of physical unitary channels).** Let $\mathcal{U} = \hat{U} \cdot \hat{U}^\dagger$ be a single-mode unitary channel, with $\hat{U} \in \mathcal{U}(S)$ a physical unitary operator. Let $E \ge 0$ be an energy parameter, and let $\epsilon > 0$ be an approximation parameter. Then, there exists $d = O(\frac{E}{\epsilon^4}E_{\hat{U}}(\frac{64E}{\epsilon^2})) \in \mathbb{N}$ such that $\hat{U}$ has $(E, \epsilon)$-approximate effective dimension $d+1$. Moreover, a finite-dimensional $(E, \epsilon)$-approximation of $\hat{U}$ can be computed in polynomial time in $d, E$, and $\frac{1}{\epsilon}$.

The proof of this result is based on constructing a good unitary operator approximation of a truncated unitary operator, while carefully controlling truncation parameters for the state space, before and after applying the
unitary channel, so that the truncated unitary operator is itself a good approximation of the original unitary operator. We refer to the Methods section for the main techniques and to Appendix A for the proof.

A direct consequence of Theorem 1 is that the corresponding sequence of $(E, \epsilon)$-approximations strongly converges to $\hat{U}$ as $\epsilon$ goes to 0 in trace distance (see the Methods section). From an operational standpoint, this implies that we can effectively truncate the Hilbert space of physical bosonic computations to a certain effective dimension, as long as the cut-off is high enough. As such, Theorem 1 provides rigorous bounds on the effective dimension required to simulate bosonic dynamics with provable accuracy [10]. In other words, there will be no catastrophic loss of information in the measurements statistics when restricting to the cut-off evolution: by the operational property of the trace distance, results of a cut-off bosonic computation will be indistinguishable from results of the original computation. However, this cut-off comes at a price: Theorem 1 identifies a direct relation between the effective dimension sufficient to embed a unitary evolution to a good precision and the energy $E_{\hat{U}}$ of the unitary evolution, as defined in Eq. (2). Such an energy bound may grow very fast for seemingly simple dynamics, such as alternating Gaussian and non-Gaussian unitary evolutions [24].

## IV. UNIVERSALITY OF POLYNOMIAL HAMILTONIANS

We have shown that bosonic unitary evolutions may be well-approximated by finite-dimensional ones. In this section, we give our second main result, namely that any finite-dimensional Hamiltonian generating such a finite-dimensional unitary evolution can be realized exactly by a polynomial Hamiltonian:

**Theorem 2 (Finite-dimensional universality of polynomial Hamiltonians).** Let $d \in \mathbb{N}$ and let $\hat{H}$ be a Hermitian operator over $\mathcal{H}_d$. There exists a polynomial Hamiltonian $P_{\hat{H}}(\hat{q}, \hat{p})$ of degree at most 3d over $\mathcal{H}$ such that

$$P_{\hat{H}}(\hat{q}, \hat{p}) = \hat{H} \oplus \hat{H}', \qquad (5)$$

where $\hat{H}'$ is a Hermitian operator over $\text{span}(|n\rangle)_{n>d}$. In particular, $e^{iP_{\hat{H}}(\hat{q}, \hat{p})} = e^{i\hat{H}} \oplus e^{i\hat{H}'}$ and

$$e^{iP_{\hat{H}}(\hat{q}, \hat{p})} |\psi\rangle = e^{i\hat{H}} |\psi\rangle, \qquad (6)$$

for all $|\psi\rangle \in \mathcal{H}_d$. Moreover, the polynomial $P_{\hat{H}}$ can be computed in polynomial time in $d$.

The proof of this result is based on the use of interpolation polynomials for reproducing exactly a target finite-dimensional operator on a subspace of the Fock basis, taking advantage of the sparsity of the canonical operators in that basis. We refer to the Methods section for the

main techniques and to Appendix B for the multimode generalisation of Theorem 2 and its proof.

An important consequence of this theorem is that polynomial Hamiltonians allow us to explore the full Hilbert space, i.e., there are no inaccessible regions of the Hilbert space when starting from a fixed reference state, such as the vacuum state, and restricting to using gates generated by polynomial Hamiltonians:

**Corollary 1 (Universal controllability of polynomial Hamiltonians in infinite dimension).** Let $|\psi\rangle \in \mathcal{H}$. For all $\epsilon > 0$, there exists a polynomial Hamiltonian $P_{\epsilon}(\hat{q}, \hat{p})$ such that evolving the vacuum state $|0\rangle$ under the Schrödinger equation with Hamiltonian $P_{\epsilon}(\hat{q}, \hat{p})$ for constant time yields a state $|\psi_{\epsilon}\rangle$ satisfying $D(|\psi\rangle, |\psi_{\epsilon}\rangle) \le \epsilon$.

Corollary 1 proves the universal controllability of polynomial Hamiltonians in infinite dimensions, resolving the open question in [20] and [24, Open question 3]. We give a quick proof in the single-mode case hereafter, the multimode case being analogous.

*Proof.* Since the state $|\psi\rangle$ is normalized, there exists $d_{\epsilon}$ such it is $\epsilon$-close in trace distance to the state truncated at $d_{\epsilon}$ and renormalized, which we denote by $|\psi_{\epsilon}\rangle$. Then, there exists a finite-dimensional unitary operator $\hat{U}_{\epsilon} = e^{i\hat{H}_{\epsilon}}$ over $\mathcal{H}_{d_{\epsilon}}$ which maps $|0\rangle$ to $|\psi_{\epsilon}\rangle$. By Theorem 2, the polynomial Hamiltonian $P_{\epsilon}(\hat{q}, \hat{p}) := P_{\hat{H}_{\epsilon}}(\hat{q}, \hat{p})$ satisfies $P_{\epsilon}(\hat{q}, \hat{p}) = \hat{H}_{\epsilon} \oplus \hat{H}_{\epsilon}'$. Evolving the vacuum under the Schrödinger equation with Hamiltonian $P_{\epsilon}(\hat{q}, \hat{p})$ thus leads to the same state as if evolving the vacuum under the $(d_{\epsilon} + 1)$-dimensional Schrödinger equation with Hamiltonian $\hat{H}_{\epsilon}$, which concludes the proof. □

Beyond quantum state engineering, Theorem 2 provides a universal method for engineering exactly any finite-dimensional Hamiltonian using polynomial Hamiltonians. This method can be used to engineer universal bosonic simulators capable of simulating the unitary evolution of any finite-dimensional quantum system.

Furthermore, Theorem 2 also implies that bosonic circuits composed of input vacuum state, unitary gates generated by polynomial Hamiltonians, and number measurements, can simulate universal qudit computations efficiently (without requiring feed-forward of measurement outcomes), generalizing a recent result for qubit computations [24, Theorem 4.1].

## V. EFFECTIVE BOSONIC COMPUTATIONS

In this section, we explore the combined consequences of Theorem 1 (effective dimensions) and Theorem 2 (effective Hamiltonians) for the validity of effective descriptions of bosonic computations. At the fundamental level, both results show that effective descriptions are capturing bosonic computations well, i.e., one can always truncate the dimension of physical unitary evolutions up to
an arbitrarily small error and use polynomial Hamiltonians to model finite-dimensional evolutions exactly. Together, these two results imply that polynomial Hamiltonians can generate any physical unitary evolution approximately to arbitrary precision:

**Theorem 3** (Infinite-dimensional universality of polynomial Hamiltonians). Let $\mathcal{U} = \hat{\mathcal{U}} \cdot \hat{U}^{\dagger}$ be a single-mode unitary channel, with $\hat{\mathcal{U}} \in \mathcal{U}(S)$ a physical unitary operator. Let $E \ge 0$ be an energy parameter, and let $\epsilon > 0$ be an approximation parameter. Then, there exists $d = \mathcal{O}(\frac{E}{\epsilon^4}E_{\hat{U}}(\frac{64E}{\epsilon^2})) \in \mathbb{N}$ and a polynomial Hamiltonian $P(\hat{q}, \hat{p})$ of degree at most 3d such that, writing $\mathcal{V} := e^{iP(\hat{q}, \hat{p})} \cdot e^{-iP(\hat{q}, \hat{p})}$,

$$ \|\mathcal{U} - \mathcal{V}\|_{\diamond}^{E} \le \epsilon. \qquad (7) $$

Moreover, the polynomial $P$ can be computed from $\hat{U}$ in polynomial time in $d$, $E$, and $\frac{1}{\epsilon}$.

As we explain in the Methods section, this result is a direct consequence of Theorem 1 and Theorem 2. Formally, Theorem 3 shows that the group of unitary gates generated by polynomial Hamiltonians is dense in the group of physical unitary operators in the strong operator topology (see the Methods section). This places the definition of universality for continuous-variable quantum gate sets, based on the ability to approximate to arbitrary precision any unitary evolution generated by a polynomial Hamiltonian [11], on an equal footing with its discrete-variable counterpart, based on the property to approximate to arbitrary precision any unitary operation [21].

As such, this shows that polynomial Hamiltonians generate truly universal gate sets and paves the way for a Solovay-Kitaev theorem [22] in the infinite-dimensional setting. We obtain a version of this theorem hereafter, for universal gate sets generated by polynomial Hamiltonians as in Theorem 2:

**Theorem 4** (Solovay-Kitaev theorem for polynomial Hamiltonians). Let $E > 0$, $\epsilon > 0$ and $d \ge \frac{64E}{\epsilon^2} \in \mathbb{N}$. Let $\mathcal{G}$ be a finite set of unitary operators over $\mathcal{H}_d$ generating a dense subset of $\mathcal{U}(\mathcal{H}_d)$ and let $\mathcal{P}$ be its realization with polynomial Hamiltonians from Theorem 2. There is a constant $c$ such that for any physical unitary operator $\hat{U} \in \mathcal{U}(S)$ with $(E, \epsilon)$-approximate effective dimension $d+1$, there exists a finite sequence $\hat{V}$ of gates from $\mathcal{P}$ of length $\mathcal{O}(\log^c(1/\epsilon))$ and such that $\|\mathcal{U} - \mathcal{V}\|_{\diamond}^{E} \le 2\epsilon$, where $\mathcal{U} = \hat{U} \cdot \hat{U}^{\dagger}$ and $\mathcal{V} = \hat{V} \cdot \hat{V}^{\dagger}$.

Recall that Theorem 1 ensures the existence of an $(E, \epsilon)$-approximate effective dimension for all physical unitary operators. The proof of Theorem 4 combines Theorems 1 and 2 with the Solovay-Kitaev theorem for qudits [22] of dimension $d+1$, and provides an explicit algorithm for producing the sequence of polynomial Hamiltonians, given the target unitary $\hat{U} \in \mathcal{U}(S)$. We refer to the Methods section for the main techniques and to Appendix C for the proof.

Theorem 4 generalises the Gaussian Solovay-Kitaev theorem from [23], thus providing an answer to [24, Open question 4] for the gate sets based on the polynomial Hamiltonians used in Theorem 2. In particular, it implies the computational equivalence of bosonic quantum circuits based on these universal gate sets.

## VI. SIGNIFICANCE AND OUTLOOK

We have considered effective descriptions of bosonic systems based on effective dimension and effective (polynomial) Hamiltonians, and rigorously proved that they faithfully reproduce the features of the original dynamics. Our findings have several far-reaching consequences for the simulation of bosonic systems and the complexity of bosonic computations: (i) it is always possible to truncate the infinite-dimensional Hilbert space of physical bosonic unitary dynamics to an effective dimension depending on the energy involved, with negligible loss of information (Theorem 1), leading to explicit bounds on the space complexity of the simulation of bosonic systems; (ii) polynomial Hamiltonians provide universal bosonic simulators of finite-dimensional quantum systems (Theorem 2) and enable universal bosonic quantum state engineering (Corollary 1), which solves a long-standing open problem on the universal controllability of polynomial Hamiltonians [20]; (iii) polynomial Hamiltonians generate universal gate sets in infinite dimensions (Theorem 3), which is of foundational importance for continuous-variable quantum computing, since it proves the validity of the definition of universality based on polynomial Hamiltonians. This was left open since the seminal paper of Lloyd and Braunstein [11], and places it on an equal footing with its discrete-variable counterpart; (iv) these universal gate sets satisfy an infinite-dimensional Solovay-Kitaev theorem, leading to bosonic circuits that are equivalent from a computational complexity standpoint (Theorem 4), thus providing a solid foundation for the complexity theory of bosonic computations [24].

Our results also motivate several new research directions. For instance, we expect that our Theorem 1 on effective dimensions of infinite-dimensional unitary dynamics may be generalized to multimode, non-unitary dynamics, as a consequence of Stinespring's dilation theorem [21].

Moreover, Theorem 2 together with Corollary 1 on the universality of polynomial Hamiltonians may enable the experimental implementation of new quantum simulators [30] and the preparation of new exotic quantum states using bosonic platforms, including superconducting [31] and photonic [32] setups. It is also an intriguing question whether the universality of polynomial Hamiltonians can be extended to hybrid boson-fermion and boson-qubit systems [33].

On a different note, it has been argued in several works [11, 15, 20] that Gaussian Hamiltonians (of degree
≤ 2) together with any single non-Gaussian polynomial Hamiltonian (of degree > 2) are sufficient to generate any polynomial Hamiltonian evolution approximately. However, a formal proof of this statement is missing, although Trotter–Suzuki formulas with rigorous error bounds have been recently uncovered in the infinite-dimensional setting [34]. Combined with our results, formalizing these polynomial gate set equivalences could lead to versions of Theorem 2 and Theorem 3 valid for any choice of non-Gaussian polynomial Hamiltonian. In particular, this may unlock an infinite-dimensional Solovay–Kitaev for all polynomial gate sets rather than only those used in our Theorem 2.

In this regard, it is well-known that the complexity of the qudit Solovay–Kitaev algorithm for finding an approximating gate sequence scales poorly with the qudit dimension [22]. This implies in turn that our Solovay–Kitaev theorem (Theorem 4) has high complexity based on the effective dimension of the target physical unitary operator. It would be interesting to develop a version achieving better efficiency in terms of energy.

In summary, we have shown that effective descriptions are sufficient to capture and reproduce bosonic evolutions, with significant implications for quantum simulation, quantum state engineering, and universal quantum computing. From a conceptual standpoint, our work indicates that bosonic phenomena which cannot be reproduced, even approximately, in finite dimensions must involve unphysical states of infinite energy. Beyond the framework of bosonic systems, our work supports the view that there is much to be learned by challenging the assumptions underlying our idealized models of physical reality [35].

## ACKNOWLEDGMENTS

The authors warmly acknowledge inspiring discussions with F. Grosshans and D. Markham. U.C. also acknowledges stimulating discussions with S. Mehraban, S. Lloyd, A. Winter, L. Lami, M. Rördam and P.-E. Emeriau. The authors also thank S. Mehraban for helpful comments on a previous version of this manuscript. F.A. is grateful to J. Eisert for insightful discussions. We acknowledge funding from the Plan France 2030 project NISQ2LSQ (ANR-22-PETQ-0006). U.C. acknowledges funding from the European Union's Horizon Europe Framework Programme (EIC Pathfinder Challenge project Veriqub) under Grant Agreement No. 101114899.

## METHODS

### A. Preliminary material

a. CV quantum information. In continuous-variable quantum information [25–27], information is encoded in

the state of *m* bosonic modes, which are mathematically modelled by a Hilbert space *H* spanned by a countably-infinite Fock basis {|*n*⟩}*,n∈N*, i.e., Σ*,n∈N* c*,n|*n*⟩ ∈ *H* if and only if Σ*,n∈N* |c*,n|2 < ∞.

To a single bosonic mode, we associate the annihilation operator $\hat{a} |n\rangle = \sqrt{n} |n-1\rangle$ and creation operator $\hat{a}^\dagger |n\rangle = \sqrt{n+1} |n+1\rangle$ which satisfy $[\hat{a}, \hat{a}^\dagger] = \hat{I}$. The canonical operators (also called position and momentum or quadrature operators) are defined as $\hat{q} = (\hat{a} + \hat{a}^\dagger)/\sqrt{2}$ and $\hat{p} = (\hat{a} - \hat{a}^\dagger)/i\sqrt{2}$, respectively.

Polynomials in the canonical operators are central in the theory of quantum computation over continuous variables, as introduced by Lloyd and Braunstein [11]. In particular, a universal continuous-variable quantum computer is defined as a device that can approximately generate any evolution over multimode bosonic systems that can be written as a (hermitian) polynomial $P(\hat{q}_1, \hat{q}_2, \dots, \hat{p}_1, \hat{p}_2, \dots)$. In [11], an inductive reasoning is used to argue that arbitrary polynomials can be obtained though “judicious commutation” of all Hamiltonians of degree less than or equal to two and any Hamiltonian of higher order. The former generate so-called Gaussian transformations [25].

Although it is a mathematically convenient setting for analysing bosonic systems, the Hilbert space $\mathcal{H}$ has some drawbacks from a physical perspective. For instance, it contains states of infinite energy and states for which the position and momentum operators are not well-defined [19]. The Schwartz space $\mathcal{S}$ is defined as the dense subspace of elements $\sum_{n \in \mathbb{N}} c_n |n\rangle \in \mathcal{H}$ such that $\lim_{n \to 0} c_n n^k = 0$ for any natural number $k$ [36]. The Schwartz space is not a Hilbert space, but it retains many of the useful analytical properties of $\mathcal{H}$ whilst eliminating these physical inconsistencies: every state in $\mathcal{S}$ has bounded energy, and there is a well-defined action of polynomials of the canonical operators on those states [20]. Now, there are unitaries that map physical (Schwartz) states to non-physical states (e.g., with infinite energy), so we exclude these unitaries and allow only those that map physical states to physical states. Letting $\mathcal{U}(\mathcal{H})$ be the group of unitaries on $\mathcal{H}$, the group of physical unitaries is defined as

$$ \mathcal{U}(\mathcal{S}) := \{ \hat{U} \in \mathcal{U}(\mathcal{H}) \mid \forall |\psi\rangle \in \mathcal{S} : \hat{U} |\psi\rangle \in \mathcal{S} \}. \quad (8) $$

We note that, in fact, all our results involving physical unitary operators (Theorems 1, 3 and 4) are valid under the weaker condition $\langle n | \hat{U}^\dagger \hat{n} \hat{U} | n \rangle < +\infty$ for all $n \in \mathbb{N}$, rather than $\hat{U} \in \mathcal{U}(\mathcal{S})$.

b. Distance measures. In making claims about approximating operators on $\mathcal{H}$ or $\mathcal{S}$, we must be explicit about which distance metric we are using to measure the quality of our approximations. The diamond norm on completely-positive channels is a natural candidate. It is defined as

$$ \|\mathcal{E}\|_{\diamond} := \sup\{\|(\mathcal{E} \otimes \text{id})\rho\|_{1} | \rho \text{ density operator on } \mathcal{H} \otimes \mathcal{H}\}, $$
where $|| \cdot ||_1$ denotes the trace norm. However, it has some undesirable properties: for instance, if $\hat{H}$ is any unbounded Hermitian operator, and $U_{H,s} : \rho \mapsto e^{is\hat{H}}\rho e^{-is\hat{H}}$ is the associated one-parameter family of unitary channels, then there are density operators $\rho$ such that $||U_{H,s}(\rho) - U_{H,s+\epsilon}(\rho)||_1$ does not go to 0 even as $\epsilon$ does, so that the diamond norm $||U_{H,s} - U_{H,s+\epsilon}||_\diamond$ also does not go to 0 [28, Proposition 2]. This clearly precludes the use of the diamond norm for measuring the distance between two Hamiltonian evolutions in operational settings.

Instead, the energy-constrained diamond norms give a family of distance measures which avoid this pitfall [28, 29]. They are defined by restricting the supremum to states whose energy is bounded by a parameter $E > 0$:

$$ \|\mathcal{E}\|_{\diamond}^{E} := \sup_{\mathrm{Tr}[\rho(\hat{H}\otimes\hat{I})]\le E} \|(\mathcal{E} \otimes \mathrm{id})\rho\|_{1}, \quad (9) $$

where now $\hat{H}$ is a reference energy operator. Mathematically, the energy-constrained diamond norms are related to the strong operator topology, which gives the correct notion of continuity in time for unitary evolutions:

**Proposition 1** ([29], Proposition 3). If, for any $E > 0$ and quantum channel $\mathcal{E}$, the sequence of quantum channels $\{\mathcal{E}_p\}_{p \in \mathbb{N}}$ is such that $\lim_{p \to \infty} ||\mathcal{E}_p - \mathcal{E}||_\diamond^E = 0$ then $\{\mathcal{E}_p\}_{p \in \mathbb{N}}$ converges to $\mathcal{E}$ in the strong operator topology.

## B. Tools for proving Theorem 1

Theorem 1 is concerned with the first kind of effective description we consider, namely replacing an infinite-dimensional unitary operator $\hat{U}$ with a finite-dimensional approximation, while precisely quantifying the error introduced in the process.

The proof hinges on two types of approximations: firstly, we show that the error measured using an energy-constrained diamond norm can be bounded by considering only a truncated Fock space, introducing a suitable truncation parameter $M$; secondly, we use the physicality of the channel to focus on its action on a larger finite-dimensional subspace, introducing another suitable cut-off parameter $N \ge M$ in Fock space. The error introduced in both steps can be bounded making use of the so-called *gentle measurement lemma*:

**Lemma 1** (Gentle measurement lemma [37], Lemma 9.4.2). Consider a density operator $\rho$ and a projector $\Pi$ where $0 \le \Pi \le \mathbb{I}$. Suppose that $\mathrm{Tr}(\Pi\rho) \ge 1 - \epsilon$, $\epsilon \in [0, 1]$. Then,

$$ ||\rho - \Pi\rho\Pi||_1 \le 2\sqrt{\epsilon}. \quad (10) $$

Note that this is a slightly weaker result than the more general version in [37], but it is sufficient for our proof. We use this lemma twice, for the projectors $\Pi_M$ and $\Pi_N$ onto the subspace of states with a number of particles

at most $M$ and $N$, respectively. Now, the restriction $\Pi_N \hat{U} \Pi_N$ of a unitary operator to a finite-dimensional subspace is not necessarily unitary. Therefore we further need to rely on a notion of “closest unitary operator” in order to construct a suitable unitary version of the cut-off of the original channel. This notion is provided by the QR factorization (related to the Gram–Schmidt process):

**Lemma 2** (QR factorization for rectangular complex matrices). Let $A \in \mathbb{C}^{N \times M}$ be a matrix with linearly independent columns. Then, there exists a unitary matrix $Q \in \mathbb{C}^{N \times N}$ and an upper triangular matrix $R \in \mathbb{C}^{M \times M}$ such that:

$$ A = Q \begin{pmatrix} R \\ 0 \end{pmatrix}. \quad (11) $$

We apply this decomposition to a column full rank sub-matrix of the cut-off operator $\Pi_N \hat{U} \Pi_N$ and choose the resulting $Q$ as our finite-dimensional unitary approximation of $\hat{U}$, carefully bounding the error introduced. Combining the approximation errors at each step completes the proof, yielding an effective dimension $d = N$. The complexity of computing the QR factorization (see for example [38], Section 5.2) for the cut-off dimension introduced determines the complexity of constructing the unitary cut-off version of the original channel, i.e., computing its effective description. We refer to the Supplementary Information for a detailed proof.

## C. Tools for proving Theorem 2

Theorem 2 seeks to reproduce a finite-dimensional unitary operator $\hat{U} = e^{i\hat{H}}$ with a unitary operator generated by a polynomial Hamiltonian.

Since we are specifically interested in Hamiltonian evolutions where the Hamiltonian is an *unbounded* operator, we must unavoidably consider some of the mathematical technicalities in manipulating such operators. In particular, if $\hat{H}$ is an unbounded Hamiltonian, we cannot expect the usual power series expansion of the exponential, $e^{it\hat{H}} = \sum_{n \in \mathbb{N}} \frac{(it\hat{H})^n}{n!}$, to converge to an operator on $\mathcal{H}$. For sufficiently well-behaved unbounded operators, the exponential may be defined instead via the spectral theorem or the resolvent formalism [19, 39]. However, there are still many cases in which defining unambiguously the unitary operator generated by a polynomial Hamiltonian is challenging, because the corresponding Schrödinger equation does not have solutions for most choices of initial conditions [19].

In the present case, we avoid these mathematical technicalities by focussing on polynomial Hamiltonians that are block-diagonal in the Fock basis, i.e., which take the form $P(\hat{q}, \hat{p}) = \Pi_N \hat{H} \Pi_N \oplus (I - \Pi_N) \hat{H}' (I - \Pi_N)$ for some $N \ge 0$, where $\Pi_N$ is the projector onto the subspace of states with a number of particles at most $N$. This ensures that, given an initial state with at most $N$ particles, the Schrödinger evolution with Hamiltonian $P(\hat{q}, \hat{p})$ remains
in the subspace of states with a number of particles at most $N$, and reproduces the finite-dimensional unitary evolution $e^{i\hat{H}\hat{\Pi}_N}$ on that subspace.

The proof of Theorem 2 then proceeds by constructing such a block-diagonal polynomial Hamiltonian $P(\hat{q}, \hat{p}) = \Pi_N \hat{H} \Pi_N \oplus (I - \Pi_N) \hat{H}' (I - \Pi_N)$, given the description of a finite-dimensional Hermitian operator $\hat{H}$. This construction is based on decomposing the finite-dimensional Hamiltonian $\hat{H}$ as a sum of terms pairing only two basis states, say $|j\rangle$ and $|k\rangle$. These terms can be mimicked by polynomial Hamiltonians that are block-diagonal in the Fock basis, constructed using the canonical operators and polynomial functions that are zero on all integers in a specific range, except for the relevant ones, depending on $j$ and $k$. This is achieved using Lagrange interpolation polynomials:

**Definition 2** (Lagrange interpolation polynomials). Given $n+1$ points $(x_0, y_0), \dots, (x_n, y_n)$, the Lagrange interpolation polynomial is defined as

$$L(x) = \sum_{j=0}^{n} y_j \prod_{\substack{0 \le i \le n \\ i \ne j}} \frac{x - x_i}{x_j - x_i} \quad (12)$$

It satisfies $L(x_i) = y_i$, for all $0 \le i \le n$.

In particular, polynomial operators of the form $P(\hat{a}^\dagger \hat{a})$ with $P$ a well-chosen Lagrange interpolation polynomial allow us to reproduce diagonal entries of $\hat{H}$, while operators of the form $P(\hat{a}^\dagger \hat{a}) \hat{a}^k + \text{h.c.}$ lead to off-diagonal entries. We refer to the Supplementary Information for a detailed proof, including a generalisation of the result to the multimode case.

### D. Proof of Theorem 3

Theorem 3 deals with the second kind of effective description we consider, namely approximating a physical infinite-dimensional unitary operator $\hat{U}$ with a unitary approximation generated by a polynomial Hamiltonian, while precisely quantifying the error introduced in the process.

*Proof.* The proof combines Theorem 1 and Theorem 2 as follows. We first use Theorem 1 to derive a $(d+1)$-dimensional $(E, \epsilon)$-approximation $\hat{V}_d$ of the original unitary operator $\hat{U}$, with $d = O(E/\epsilon^4 E_{\hat{U}} (\frac{64E}{\epsilon^2}))$, in time

polynomial in $d, E$, and $\frac{1}{\epsilon}$. Then, we extract a finite-dimensional Hamiltonian $\hat{H}$ generating the unitary evolution $\hat{V}_d$. We explain hereafter how this step can be performed efficiently. Finally, we use Theorem 2 to construct a polynomial Hamiltonian of degree at most $3d$ which generates a unitary evolution that reproduces the unitary operator $\hat{V}_d = e^{i\hat{H}}$, itself approximating the original unitary operator $\hat{U}$.

Since $\hat{V}_d \in \mathcal{U}(\mathcal{H}_d)$ is a unitary operator on a $(d+1)$-dimensional Hilbert space, we can write $\hat{V}_d = WDW^\dagger$ for some unitary matrix $W$ and a diagonal matrix $D$, whose diagonal elements are of the form $(e^{i\lambda_n})_{n=0,\dots,d}$ for $\lambda \in [0, 2\pi)$. Let $\Sigma$ be the diagonal matrix whose diagonal elements are $(\lambda_n)_{n=0,\dots,d}$, then $\Sigma$ is Hermitian and $\hat{V}_d = We^{i\Sigma W^\dagger} = e^{iW\Sigma W^\dagger}$. Since diagonalising $D$ and the matrix multiplications with $\Sigma$ can both be performed in time polynomial in $d$, we can extract the generator $\hat{H} = W\Sigma W^\dagger$ of $\hat{V}_d$ efficiently, which concludes the proof of Theorem 3. □

### E. Tools for proving Theorem 4

Theorem 4 provides a Solovay-Kitaev theorem for universal unitary gate sets generated by block-diagonal polynomial Hamiltonians constructed in Theorem 2.

Given a good finite-dimensional unitary approximation $\hat{V}_d$ of the physical unitary $\hat{U}$, guaranteed to exist by Theorem 1, the proof proceeds by synthesizing a short sequence of finite-dimensional gates approximating $\hat{V}_d$ over $\mathcal{H}_d$, via the Solovay-Kitaev theorem for qudits:

**Theorem 5** (Qudit Solovay-Kitaev theorem [22]). Let $\mathcal{G}$ be a finite set of unitary gates with determinant one and their inverses, generating a dense subset of $SU(d)$. There is a constant $c$ such that for any $U \in SU(d)$, there exists a finite sequence $S$ of gates from $\mathcal{G}$ of length $O(\log^c(\frac{1}{c}))$ and such that $d(U, S) < \epsilon$, where $d$ is the distance induced by the operator norm.

We then employ Theorem 2 to generate each finite-dimensional unitary gate in the sequence using polynomial Hamiltonians, obtaining a short sequence of unitary operators over $\mathcal{H}$ generated by polynomial Hamiltonians approximation $\hat{U}$ over $\mathcal{H}_d$. Finally, we leverage the physicality of the target unitary channel to lift the approximation of $\hat{U}$ over $\mathcal{H}_d$ to an approximation over $\mathcal{H}$. We refer to the Supplementary Information for a detailed proof.

[1] S. Aaronson and A. Arkhipov, Theory of Computing 9, 143 (2013).

[2] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Physical Review Letters 119, 170501 (2017).

[3] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C.

Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, et al., Science 370, 1460 (2020).

[4] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, A. E. Lita, T. Gerrits, S. W. Nam, V. D. Vaidya, M. Menotti, I. Dhand, Z. Vernon,
N. Quesada, and J. Lavoie, *Nature* **606**, 75 (2022).

[5] V. V. Albert, *arXiv preprint arXiv:2211.05714* (2022), 10.48550/arXiv quant-ph/0505030.

[6] D. Gottesman, A. Kitaev, and J. Preskill, *Phys. Rev. A* **64**, 012310 (2001).

[7] V. Sivak, A. Eickbusch, B. Royer, S. Singh, I. Tsioutsios, S. Ganjam, A. Miano, B. Brock, A. Ding, L. Frunzio, et al., *Nature* **616**, 50 (2023).

[8] It is well-known that such a relation cannot be satisfied by finite-dimensional linear operators, as can be seen by contradiction by taking the trace, giving zero for the left hand side and a non-zero value for the right hand side.

[9] S. P. Jordan, K. S. Lee, and J. Preskill, *Science* **336**, 1130 (2012).

[10] Y. Tong, V. V. Albert, J. R. McClean, J. Preskill, and Y. Su, *Quantum* **6**, 816 (2022).

[11] S. Lloyd and S. L. Braunstein, *Physical Review Letters* **82**, 1784 (1999).

[12] H. A. Gersch and G. C. Knollman, *Physical Review* **129**, 959 (1963).

[13] P. Ramond, *Field theory: a modern primer* (Routledge, 2020).

[14] L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University Press, 1995).

[15] S. Sefi and P. Van Loock, *Physical Review Letters* **107**, 170501 (2011).

[16] B. Tsirelson, *Problem statement from website of open problems at https://www.scienceDirect.com/doi/10.1016/j.science.2011.06.001* (2006).

[17] Z. Ji, A. Natarajan, T. Vidick, J. Wright, and H. Yuen, *Communications of the ACM* **64**, 131 (2021).

[18] J. M. Martinis, in *Les Houches*, Vol. 79 (Elsevier, 2004) pp. 487–520.

[19] B. C. Hall, *Quantum theory for mathematicians*, Vol. 267 (Springer, 2013).

[20] R.-B. Wu, T.-J. Tarn, and C.-W. Li, *Physical Review A* **73**, 012719 (2006).

[21] M. A. Nielsen and I. L. Chuang, *Quantum Computation and Quantum Information: 10th Anniversary Press* 3rd edition (1996). 10th ed. (Cambridge University Press, New York, NY, USA, 2011).

[22] C. M. Dawson and M. A. Nielsen, *arXiv preprint quant-ph/0505030* (2005), 10.48550/arXiv quant-ph/0505030.

[23] S. Xie, B. Datta, L. Lami, and C. Rouzé, *Physical Review Letters* **126**, 190504 (2021).

[24] U. Chabaud, M. Joseph, S. Mehra-ban, and A. Motamedi, *arXiv preprint arXiv:2410.04274* (2024), 10.48550/arXiv:2410.04274.

[25] S. L. Braunstein and P. Van Loock, *Reviews of modern physics* **77**, 513 (2005).

[26] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, *Reviews of Modern Physics* **84**, 621 (2012).

[27] G. Adesso, S. Ragy, and A. R. Lee, *Open Systems & Information Dynamics* **21**, 1440001 (2014).

[28] A. Winter, *arXiv preprint arXiv:1712.10267* (2017), 10.48550/arXiv:1712.10267.

[29] M. E. Shirokov, *arXiv preprint arXiv:1706.00361* (2017), 10.48550/arXiv:1706.00361.

[30] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, *Nature* **607**, 667 (2022).

[31] J. H. Busnaina, Z. Shi, A. McDonald, D. Dubyna, I. Nsanzineza, J. S. Hung, C. S. Chang, A. A. Clerk, and C. M. Wilson, *Nature Communications* **15**, 3065 (2024).

[32] C. Wang, F.-M. Liu, M.-C. Chen, H. Chen, X.-H. Zhao, C. Ying, Z.-X. Shang, J.-W. Wang, Y.-H. Huo, C.-Z. Peng, et al., *Science* **384**, 579 (2024).

[33] C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Demler, F. Grusdt, M. Greiner, and D. Greif, *Science* **365**, 1501 (2017).

[34] D. Burgarth, P. Facchi, A. Hahn, M. Johnsson, and K. Yuasa, *Physical Review Research* **6**, 043155 (2024).

[35] A. Einstein, B. Podolsky, and N. Rosen, *Physical Review* **47**, 777 (1935).

[36] L. Schwartz, in *Annales de l'Université de Grenoble*, Vol. 23 (1947) pp. 7–24.

[37] M. M. Wilde, *Quantum information theory* (Cambridge University Press, 2013).

[38] G. H. Gloub and C. F. Van Loan, Johns Hopkins University Press 3rd edition (1996).

[39] N. Dunford and J. T. Schwartz, *Linear operators, part 1: general theory*, Vol. 10 (John Wiley & Sons, 1988).
# Supplementary Information

## Appendix A: Proof of Theorem 1: Effective dimension of physical unitary channels

Hereafter, we denote by $\Pi_N = \sum_{n=0}^N |n\rangle\langle n|$ the projector onto the subspace of states with a number of particle at most $N$, for $N \ge 0$. We give a formal version of Theorem 1 from the main text:

**Theorem 1 (Effective dimension of physical unitary channels).** Let $\mathcal{U} = \hat{U} \cdot \hat{U}^\dagger$ be a unitary channel, with $\hat{U} \in \mathcal{U}(S)$ a physical unitary operator. Let $E \ge 0$ be an energy bound, and let $\epsilon > 0$ be an approximation parameter. Then, there exists $N = O(\frac{E}{\epsilon^4} E_{\hat{U}}(\frac{64E}{\epsilon^2})) \in \mathbb{N}$ and a unitary operator $\hat{V}_N$ over the range of $\Pi_N$ that approximates $\hat{U}$ in the following sense: for any unitary operator $\hat{V}'$ on $\text{Range}(I - \Pi_N)$, denoting $\hat{V} = \hat{V}_N \oplus \hat{V}'$ and $\mathcal{V} = \hat{V} \cdot \hat{V}^\dagger$,

$$ \|\mathcal{U} - \mathcal{V}\|_{\diamond}^{E} \le \epsilon, \qquad (A1) $$

where the energy-constrained diamond norm is defined with respect to the number Hamiltonian $\hat{n} = \hat{a}^\dagger\hat{a}$. Moreover, the cut-off operator $\hat{V}_N$ can be computed from $\hat{U}$ in time

$$ \mathcal{O}\left(\frac{N E^2}{\epsilon^4}\right) = \mathcal{O}\left(\frac{E^3}{\epsilon^8} E_{\hat{U}}\left(\frac{64E}{\epsilon^2}\right)\right). \qquad (A2) $$

*Proof sketch.* We focus on the single-mode case for simplicity. With the notations of the theorem, by [23, Theorem 1], it is equivalent to show that there exists $\hat{V}_N$ such that for all $\hat{V}'$,

$$ \sup_{\substack{|\psi\rangle \in \mathcal{H} \\ \langle\psi|\hat{n}|\psi\rangle \le E}} \|\hat{U}|\psi\rangle\langle\psi| \hat{U}^\dagger - (\hat{V}_N \oplus \hat{V}') |\psi\rangle\langle\psi| (\hat{V}_N^\dagger \oplus \hat{V}'^\dagger)\|_1 \le \epsilon. \qquad (A3) $$

As explained in the Methods section of the main text, the idea of the proof is to perform two successive cut-offs of the state space, before and after the unitary evolution, respectively parametrised by integers $M$ and $N$, and to carefully control the approximation at each step, followed by a rounding of the cut-off operator to the closest unitary operator. More precisely, the steps of the proof are as follows: let $N \ge M \ge 0$,

(i) show that the supremum over energy-constrained states is well approximated by a supremum over $\mathcal{H}_M$ when $M$ is large, using the gentle measurement lemma;

(ii) show that the action of $\hat{U}$ on $\mathcal{H}_M$ is well-approximated by that of $\Pi_N \hat{U} \Pi_N$ when $N$ is large, also using the gentle measurement lemma;

(iii) show that the action of $\Pi_N \hat{U} \Pi_N$ itself on $\mathcal{H}_M$ is well-approximated by that of a unitary operator $\hat{V}_N$ defined on $\mathcal{H}_N \supset \mathcal{H}_M$ when $N$ is large, using the Gram–Schmidt orthonormalisation procedure;

(iv) combine the successive approximations to conclude.

*Proof of Theorem 1.* (i) Let $\hat{V}$ be a unitary operator and let $|\psi\rangle = \sum_{n \ge 0} \psi_n |n\rangle \in \mathcal{H}$ such that $\langle\psi|\hat{n}|\psi\rangle \le E$. Denoting $\mathcal{U} = \hat{U} \cdot \hat{U}^\dagger$ and $\mathcal{V} = \hat{V} \cdot \hat{V}^\dagger$, we have

$$ \begin{aligned} \|(\mathcal{U} - \mathcal{V})(\psi)\|_1 &\le \|\hat{U}\Pi_M |\psi\rangle\langle\psi| \Pi_M \hat{U}^\dagger - \hat{V}\Pi_M |\psi\rangle\langle\psi| \Pi_M \hat{V}^\dagger\|_1 \\ &\quad + \|\hat{U}(|\psi\rangle\langle\psi| - \Pi_M |\psi\rangle\langle\psi| \Pi_M) \hat{U}^\dagger - \hat{V}(|\psi\rangle\langle\psi| - \Pi_M |\psi\rangle\langle\psi| \Pi_M) \hat{V}^\dagger\|_1 \\ &\le \|\hat{U}\Pi_M |\psi\rangle\langle\psi| \Pi_M \hat{U}^\dagger - \hat{V}\Pi_M |\psi\rangle\langle\psi| \Pi_M \hat{V}^\dagger\|_1 \\ &\quad + \|\hat{U}(|\psi\rangle\langle\psi| - \Pi_M |\psi\rangle\langle\psi| \Pi_M) \hat{U}^\dagger\|_1 \\ &\quad + \|\hat{V}(|\psi\rangle\langle\psi| - \Pi_M |\psi\rangle\langle\psi| \Pi_M) \hat{V}^\dagger\|_1 \\ &= \|(\mathcal{U} - \mathcal{V})(\Pi_M \psi \Pi_M)\|_1 + 2\| |\psi\rangle\langle\psi| - \Pi_M |\psi\rangle\langle\psi| \Pi_M \|_1 \\ &\le \|(\mathcal{U} - \mathcal{V})(\Pi_M \psi \Pi_M)\|_1 + 4\sqrt{\frac{E}{M}}, \end{aligned} \qquad (A4) $$

where we used the triangle inequality for the first two steps, the fact that the trace distance is unitarily invariant in the third step, and the gentle measurement lemma in the last step.
(ii) Given $\hat{U} \in \mathcal{U}(S)$ and $M \ge 0$, we now look for $\hat{V}$ such that $\|(\mathcal{U}-\mathcal{V})(\phi)\|_1$ is uniformly bounded over $\text{Range}(\Pi_M)$. With $N \ge M$, we pick a unitary operator of the form $\hat{V} = \hat{V}_N \oplus \hat{V}'$, where $\hat{V}_N$ is a unitary operator over $\mathcal{H}_N$ to be defined later and $\hat{V}'$ is arbitrary. For all $|\phi\rangle \in \mathcal{H}_M$, we write $|\phi_U\rangle := \hat{U}|\phi\rangle$ and $|\phi_V\rangle := \hat{V}|\phi\rangle = \hat{V}_N|\phi\rangle$. We have $\Pi_N |\phi_V\rangle\langle\phi_V| \Pi_N = |\phi_V\rangle\langle\phi_V|$, since $\hat{V} = \hat{V}_N \oplus \hat{V}'$ and $N \ge M$, and $\langle\phi_U|\hat{n}|\phi_U\rangle \le E_{\hat{U}}(M)$, where

$$E_{\hat{U}}(M) := \sup_{|\chi\rangle \in \mathcal{H}_M} \langle\chi|\hat{U}^{\dagger}\hat{n}\hat{U}|\chi\rangle, \qquad (A5)$$

which is a bounded quantity for each $M \ge 0$ since $\hat{U} \in \mathcal{U}(S)$ and $\mathcal{H}_M$ is compact. Hence, following the same steps as in Eq. (A4), for all $N \ge M$,

$$ \begin{aligned} \|(\mathcal{U} - \mathcal{V})(\phi)\|_1 &\le \|\Pi_N |\phi_U\rangle\langle\phi_U| \Pi_N - \Pi_N |\phi_V\rangle\langle\phi_V| \Pi_N\|_1 \\ &\quad + \| |\phi_U\rangle\langle\phi_U| - \Pi_N |\phi_U\rangle\langle\phi_U| \Pi_N\|_1 \\ &\quad + \| |\phi_V\rangle\langle\phi_V| - \Pi_N |\phi_V\rangle\langle\phi_V| \Pi_N\|_1 \\ &\le \|\Pi_N |\phi_U\rangle\langle\phi_U| \Pi_N - \Pi_N |\phi_V\rangle\langle\phi_V| \Pi_N\|_1 + 2\sqrt{\frac{E_{\hat{U}}(M)}{N}} \end{aligned} \qquad (A6) \\ = \|\Pi_N \hat{U} |\phi\rangle\langle\phi| \hat{U}^{\dagger} \Pi_N - \hat{V}_N |\phi\rangle\langle\phi| \hat{V}_N^{\dagger}\|_1 + 2\sqrt{\frac{E_{\hat{U}}(M)}{N}}. $$

where we used the triangle inequality for the first step, and the fact that $\Pi_N |\phi_V\rangle\langle\phi_V| \Pi_N = |\phi_V\rangle\langle\phi_V|$ and the gentle measurement lemma for the second step.

(iii) Given $\hat{U} \in \mathcal{U}(S)$ and $N \ge M \ge 0$, we now look for a unitary operator $\hat{V}_N$ over $\mathcal{H}_N$ such that $\|\Pi_N \hat{U} |\phi\rangle\langle\phi| \hat{U}^{\dagger} \Pi_N - \hat{V}_N |\phi\rangle\langle\phi| \hat{V}_N^{\dagger}\|_1$ is uniformly bounded over $\mathcal{H}_M$. We first derive useful properties of the column vectors $(|u_0\rangle, \dots, |u_M\rangle)$ of $\Pi_N \hat{U} \Pi_M$ (i.e., the first $M$ columns of $\Pi_N \hat{U} \Pi_M$), which are obtained in Eqs. (A11,A12) and Lemma 3.

Let us denote $u_{ij} := \langle i|\hat{U}|j\rangle$ the Fock basis coefficients of the unitary operator $\hat{U}$, for all $i, j \in \mathbb{N}$. For all $k = 0, \dots, M$, we have

$$ \begin{aligned} \langle u_k | u_k \rangle &= \sum_{n=0}^{N} |u_{nk}|^2 \\ &= 1 - \sum_{n>N} |u_{nk}|^2 \\ &\le 1, \end{aligned} \qquad (A7) $$

since $\hat{U}$ is unitary, and

$$ \begin{aligned} \langle u_k | u_k \rangle &= 1 - \sum_{n>N} |u_{nk}|^2 \\ &> 1 - \frac{1}{N} \sum_{n>N} n |u_{nk}|^2 \qquad (A8) \\ &\ge 1 - \frac{E_{\hat{U}}(M)}{N}, \end{aligned} $$
where we used the definition of $E_{\hat{U}}(M)$ from Eq. (A5) in the last line. Moreover, for all $0 \le i \ne j \le M$,

$$ 
\begin{aligned}
|\langle u_i | u_j \rangle| &= \left| \sum_{n=0}^{N} u_{ni}^* u_{nj} \right| \\
&= \left| \sum_{n>N} u_{ni}^* u_{nj} \right| \\
&\le \sqrt{\sum_{n>N} |u_{ni}|^2 \sum_{n>N} |u_{nj}|^2} \\
&\le \frac{1}{N} \sqrt{\sum_{n>N} n |u_{ni}|^2 \sum_{n>N} n |u_{nj}|^2} \\
&\le \frac{E_{\hat{U}}(M)}{N},
\end{aligned} 
\quad (A9) $$

where we used the fact that $\hat{U}$ is unitary in the second line, Cauchy-Schwarz inequality in the third line, and the definition of $E_{\hat{U}}(M)$ from Eq. (A5) in the last line.

Eqs. (A7), (A8) and (A9) imply that taking $N$ large enough ensures that the $M+1$ first column vectors ($|u_0\rangle, \dots, |u_M\rangle$) of the Fock basis matrix of $\Pi_N \hat{U} \Pi_N$ are close to being an orthonormal family. In particular, we define

$$ \delta := \frac{E_{\hat{U}}(M)}{N}, \quad (A10) $$

and assume in what follows that $N$ is chosen large enough so that $\delta < \frac{1}{2M}$. Then,

$$ |\langle u_i | u_j \rangle| \le \delta \quad 0 \le i \ne j \le M \quad (A11) $$

$$ 1 - \delta \le \langle u_k | u_k \rangle \le 1 \quad k = 0, \dots, M. \quad (A12) $$

We also rely on the following technical result, which says that each $|u_j\rangle$ is almost orthogonal to the subspace spanned by ($|u_0\rangle, \dots, |u_{i-1}\rangle$), for any $i \le j$:

**Lemma 3.** With $\delta > 0$ defined in Eq. (A10) and $\delta < \frac{1}{2M}$, let $P_k$ denotes the orthogonal projector onto $\text{span}(|u_0\rangle, \dots, |u_{k-1}\rangle)$ for all $k = 1, \dots, M$. For all $0 \le i \le j \le M$, $\langle u_j | P_i | u_j \rangle \le \delta$. As a result, $\Pi_N \hat{U} \Pi_M$ is column full rank.

*Proof.* For all $1 \le i \le j \le M$ we have $\langle u_j | P_i | u_j \rangle \le \langle u_j | P_j | u_j \rangle$, so it is sufficient to prove the result when $i = j$. For $j = 1, \dots, M$, defining the normalised projection

$$ |\phi_j\rangle := \frac{P_j |u_j\rangle}{\|P_j |u_j\rangle\|}, \quad (A13) $$

we have $\langle u_j | P_j | u_j \rangle = |\langle u_j | \phi_j \rangle|^2$. Moreover, $|\phi_j\rangle \in \text{span}(|u_0\rangle, \dots, |u_{j-1}\rangle)$ so we may write $|\phi_j\rangle := \sum_{i=1}^{j-1} \varphi_i |u_i\rangle$, and
we have

$$
\begin{align*}
1 &= \langle \phi_j | \phi_j \rangle \\
&= \sum_{k,l=0}^{j-1} \varphi_k^* \varphi_l \langle u_k | u_l \rangle \\
&= \sum_{i=0}^{j-1} |\varphi_i|^2 \langle u_i | u_i \rangle + 2 \Re \left( \sum_{0 \le k < l \le j-1} \varphi_k^* \varphi_l \langle u_k | u_l \rangle \right) \\
&\ge (1-\delta) \sum_{i=0}^{j-1} |\varphi_i|^2 - 2 \sum_{0 \le k < l \le j-1} |\varphi_k||\varphi_l|| \langle u_k | u_l \rangle | \tag{A14} \\
&\ge \sum_{i=0}^{j-1} |\varphi_i|^2 - \delta \sum_{i=0}^{j-1} |\varphi_i|^2 - 2\delta \sum_{0 \le k < l \le j-1} |\varphi_k||\varphi_l| \\
&= \sum_{i=0}^{j-1} |\varphi_i|^2 - \delta \left( \sum_{i=0}^{j-1} |\varphi_i| \right)^2 \\
&\ge (1-j\delta) \sum_{i=0}^{j-1} |\varphi_i|^2,
\end{align*}
$$

where we used Eq. (A12) in the third line, Eq. (A11) in the fourth line, and Cauchy–Schwarz inequality in the last line. This shows that

$$
\sum_{i=0}^{j-1} |\varphi_i|^2 \le \frac{1}{1-j\delta}, \qquad (A15)
$$

and we thus obtain

$$
\begin{align*}
\langle u_j | P_j | u_j \rangle &= \left| \sum_{i=0}^{j-1} \varphi_i \langle u_j | u_i \rangle \right|^2 \\
&\le \sum_{i=0}^{j-1} |\varphi_i|^2 \sum_{i=0}^{j-1} |\langle u_j | u_i \rangle|^2 \tag{A16} \\
&\le \frac{j\delta^2}{1-j\delta} \\
&\le \delta,
\end{align*}
$$

where we used Eqs. (A11) and (A15) in the third line and $j \le M$ and $\delta < \frac{1}{2M}$ in the last line. This completes the first part of the lemma.

For all $j = 1, \dots, M$, we obtain $\langle u_j | I - P_j | u_j \rangle = \langle u_j | u_j \rangle - \langle u_j | P_j | u_j \rangle \ge 1 - 2\delta$, so $|u_j| \notin \operatorname{span}(|u_0\rangle, \dots, |u_{j-1}\rangle)$ and $\Pi_N \hat{U} \Pi_M$ is column full rank. $\square$

Now that these properties are established, recall that we are looking for a unitary operator $\hat{V}_N$ over $\mathcal{H}_N$ such that $\|\Pi_N \hat{U} |\phi\rangle\langle\phi| \hat{U}^\dagger \Pi_N - \hat{V}_N |\phi\rangle\langle\phi| \hat{V}_N^\dagger\|_1$ is uniformly bounded over $\mathcal{H}_M$. To this end, we relate $\Pi_N \hat{U} \Pi_N$ to a unitary operator through the Gram–Schmidt process. Note that we can assume $\Pi_N \hat{U} \Pi_N$ to be full rank without loss of generality, because by Lemma 3, $\Pi_N \hat{U} \Pi_M$ is column full rank and can be completed to a full rank matrix on $\mathcal{H}_N$ which matches with $\Pi_N \hat{U} \Pi_N$ on $\mathcal{H}_N$.

We thus assume $\Pi_N \hat{U} \Pi_N$ to be full rank and write its QR decomposition $\Pi_N \hat{U} \Pi_N = \hat{Q}_N \hat{R}_N$ where $\hat{Q}_N$ is unitary over $\mathcal{H}_N$ and where $\hat{R}_N$ has an upper-triangular matrix in Fock basis given by:

$$
\begin{pmatrix}
\langle e_0 | u_0 \rangle & \langle e_0 | u_1 \rangle & \langle e_0 | u_2 \rangle & \cdots & \langle e_0 | u_M \rangle \\
0 & \langle e_1 | u_1 \rangle & \langle e_1 | u_2 \rangle & \cdots & \langle e_1 | u_N \rangle \\
0 & 0 & \langle e_2 | u_2 \rangle & \cdots & \langle e_2 | u_N \rangle \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & \langle e_N | u_N \rangle
\end{pmatrix},
\quad (A17)
$$
where $(|u_0\rangle, \dots, |u_N\rangle)$ are the column vectors of the Fock basis matrix of $\Pi_N \hat{U} \Pi_N$, and where $(|e_0\rangle, \dots, |e_N\rangle)$ are their orthonormalised version through the Gram-Schmidt process, forming the column vectors of the Fock basis matrix of $\hat{Q}_N$. Recall from Lemma 3 that for all $k = 1, \dots, N$, $P_k$ denotes the orthogonal projector onto $\text{span}(|u_0\rangle, \dots, |u_{k-1}\rangle)$. Let

$$ \begin{aligned} |v_0\rangle &:= |u_0\rangle \\ |v_k\rangle &:= (I - P_k) |u_k\rangle. \end{aligned} \qquad (A18) $$

The orthonormal vectors $(|e_0\rangle, \dots, |e_N\rangle)$ are then given by

$$ |e_k\rangle := \frac{|v_k\rangle}{\|v_k\|}, \qquad (A19) $$

for all $k = 0, \dots, N$.

We now set $\hat{V}_N = \hat{Q}_N$. For all $|\phi\rangle \in \mathcal{H}_M$, we have

$$ \begin{aligned} \|\Pi_N \hat{U} |\phi\rangle\langle\phi| \hat{U}^\dagger \Pi_N - \hat{V}_N |\phi\rangle\langle\phi| \hat{V}_N^\dagger \|_1 &= \|\Pi_N \hat{U} \Pi_N |\phi\rangle\langle\phi| \Pi_N \hat{U}^\dagger \Pi_N - \hat{V}_N |\phi\rangle\langle\phi| \hat{V}_N^\dagger \|_1 \\ &= \|\hat{Q}_N \hat{R}_N |\phi\rangle\langle\phi| \hat{R}_N^\dagger \hat{Q}_N^\dagger - \hat{Q}_N |\phi\rangle\langle\phi| \hat{Q}_N^\dagger \|_1 \\ &= \|\hat{R}_N |\phi\rangle\langle\phi| \hat{R}_N^\dagger - |\phi\rangle\langle\phi| \|_1 \\ &= \|\Pi_M \hat{R}_N \Pi_M |\phi\rangle\langle\phi| \Pi_M \hat{R}_N^\dagger \Pi_M - |\phi\rangle\langle\phi| \|_1, \end{aligned} \qquad (A20) $$

where we used the fact that the trace distance is unitarily invariant in the third line and the fact that the Fock basis matrix of $\hat{R}_N$ is upper-triangular in the last line.

To conclude, we show that this quantity can be made arbitrarily small uniformly over $\mathcal{H}_M$ for all $M$, by picking $N$ large enough. To do so, we bound the coefficients of the Fock basis matrix of $\Pi_M \hat{R}_N \Pi_M$. From Eqs. (A18) and (A19) we have

$$ \langle e_i | u_j \rangle = \frac{\langle v_i | u_j \rangle}{\sqrt{\langle v_i | v_i \rangle}}. \qquad (A21) $$

for all $0 \le i, j \le M$. In particular, for all $0 \le i < j \le M$, with the convention $P_0 = 0$,

$$ \begin{aligned} |\langle v_i | u_j \rangle| &= |\langle u_i | I - P_i | u_j \rangle| \\ &\le |\langle u_i | u_j \rangle| + |\langle u_i | P_i | u_j \rangle| \\ &\le \delta + \sqrt{\langle u_i | P_i | u_i \rangle \langle u_j | P_i | u_j \rangle} \\ &\le 2\delta, \end{aligned} \qquad (A22) $$

where we used the triangle inequality in the second line, Eq. (A11) in the third line, Cauchy-Schwarz inequality in the fourth line and Lemma 3 twice in the last line. With $\sqrt{\langle v_i | v_i \rangle} = \sqrt{\langle u_i | u_i \rangle - \langle u_i | P_i | u_i \rangle} \ge \sqrt{1 - 2\delta}$ by Eq. (A12) and Lemma 3, this yields

$$ \begin{aligned} |\langle e_i | u_j \rangle| &= \frac{|\langle v_i | u_j \rangle|}{\sqrt{\langle v_i | v_i \rangle}} \\ &\le \frac{2\delta}{\sqrt{1 - 2\delta}} \\ &\le 4\delta, \end{aligned} \qquad (A23) $$

where we used $\delta \le \frac{3}{8}$ in the last line. Similarly, for all $k = 0, \dots, M$,

$$ \begin{aligned} \langle e_k | u_k \rangle &= \frac{\langle v_k | u_k \rangle}{\sqrt{\langle v_k | v_k \rangle}} \\ &= \frac{\langle u_k | I - P_k | u_k \rangle}{\sqrt{\langle u_k | I - P_k | u_k \rangle}} \\ &= \sqrt{\langle u_k | I - P_k | u_k \rangle} \\ &= \sqrt{\langle u_k | u_k \rangle - \langle u_k | P_k | u_k \rangle}. \end{aligned} \qquad (A24) $$
Hence, $0 \le \langle e_k | u_k \rangle \le 1$ by Eq. (A12), and $\langle e_k | u_k \rangle \ge \sqrt{1-2\delta}$ by Eq. (A12) and Lemma 3, leading to

$$ \begin{aligned} & |\langle e_k | u_k \rangle - 1| \le 1 - \sqrt{1 - 2\delta} \\ & \le 2\delta, \end{aligned} \quad (A25) $$

where we used $\delta \le \frac{1}{2}$ in the last line.

To convert these bounds to a bound on the trace distance in Eq. (A20), we make use of the following technical result:

**Lemma 4 (Trace distance between unnormalised pure states).** Let $\hat{A}$ be a linear operator over $\mathcal{H}_M$ and let $|\phi\rangle \in \mathcal{H}_M$ be a normalised pure state. Then

$$ \| \hat{A} |\phi\rangle \langle \phi | \hat{A}^\dagger - |\phi\rangle \langle \phi | \|_1 = \sqrt{(1 + \langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle)^2 - 4 |\langle \phi | \hat{A} | \phi \rangle|^2}. \quad (A26) $$

*Proof.* We write $\hat{A} |\phi\rangle$ in a basis ($|\phi\rangle, |\phi^\dagger\rangle$), where $|\phi^\dagger\rangle$ is orthogonal to $|\phi\rangle$. The density matrix of $\hat{A} |\phi\rangle \langle \phi | \hat{A}^\dagger - |\phi\rangle \langle \phi|$ in that basis is given by

$$ \left( \begin{array}{c|c} |\langle \phi | \hat{A} | \phi \rangle|^2 - 1 & |\langle \phi | \hat{A} | \phi \rangle| \sqrt{\langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - |\langle \phi | \hat{A} | \phi \rangle|^2} \\ |\langle \phi | \hat{A} | \phi \rangle| \sqrt{\langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - |\langle \phi | \hat{A} | \phi \rangle|^2} & |\langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - |\langle \phi | \hat{A} | \phi \rangle|^2 \end{array} \right), \quad (A27) $$

with eigenvalues

$$ \lambda_{\pm} = \frac{1}{2} (\langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - 1) \pm \sqrt{\frac{1}{4} (\langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - 1)^2 + \langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - |\langle \phi | \hat{A} | \phi \rangle|^2}. \quad (A28) $$

Hence,

$$ \begin{aligned} \| \hat{A} |\phi\rangle \langle \phi | \hat{A}^\dagger - |\phi\rangle \langle \phi | \|_1 &= |\lambda_+| + |\lambda_-| \\ &= 2 \sqrt{\frac{1}{4} (\langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - 1)^2 + \langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle - |\langle \phi | \hat{A} | \phi \rangle|^2} \quad (A29) \\ &= \sqrt{(1 + \langle \phi | \hat{A}^\dagger \hat{A} | \phi \rangle)^2 - 4 |\langle \phi | \hat{A} | \phi \rangle|^2}. \end{aligned} $$

<figure><img src="image_1.png" /></figure>

We apply this lemma to $\hat{A} = \Pi_M \hat{R}_N \Pi_M$. By Eqs. (A23,A25),

$$ |\langle i | \hat{A} | j \rangle| \le 4\delta \quad 0 \le i \ne j \le M \quad (A30) $$

$$ 1 - 2\delta \le \langle k | \hat{A} | k \rangle \le 1 \quad k = 0, \dots, M. \quad (A31) $$

As a result, for $0 \le i < j \le N$,

$$ \begin{aligned} |\langle i | \hat{A}^\dagger \hat{A} | j \rangle| &= \left| \sum_{k=0}^{M} \langle i | \hat{A}^\dagger | k \rangle \langle k | \hat{A} | j \rangle \right| \\ &\le \sum_{k=0}^{i-1} |\langle i | \hat{A}^\dagger | k \rangle| |\langle k | \hat{A} | j \rangle| + |\langle i | \hat{A}^\dagger | i \rangle| |\langle i | \hat{A} | j \rangle| \quad (A32) \\ &\le 16N\delta^2 + 4\delta \\ &\le 12\delta, \end{aligned} $$

where we used $\delta < \frac{1}{2M}$ in the last line. Similarly, for all $k = 0, \dots, M$,

$$ \begin{aligned} |\langle k | \hat{A}^\dagger \hat{A} | k \rangle| &\le \sum_{l=0}^{k-1} |\langle k | \hat{A}^\dagger | l \rangle| |\langle l | \hat{A} | k \rangle| + |\langle k | \hat{A}^\dagger | k \rangle|^2 \\ &\le 16N\delta^2 + 1 \\ &\le 1 + 8\delta, \end{aligned} \quad (A33) $$
where we used $\delta < \frac{1}{2M}$ in the last line. For $|\phi\rangle := \sum_{k=0}^{M} \phi_k |k\rangle \in \mathcal{H}_M$, we thus have

$$ \begin{aligned} |\langle\phi|\hat{A}^\dagger\hat{A}|\phi\rangle| &\le \sum_{k=0}^{M} |\phi_k|^2 |\langle k|\hat{A}^\dagger\hat{A}|k\rangle| + \sum_{i \ne j} |\phi_i||\phi_j||\langle i|\hat{A}^\dagger\hat{A}|j\rangle| \\ &\le 1 + 8\delta + 12\delta \sum_{i \ne j} |\phi_i||\phi_j| \\ &\le 1 + (8 + 12N)\delta, \end{aligned} \quad (A34) $$

where we used Cauchy-Schwarz inequality in the last line. We also have

$$ \begin{aligned} |\langle\phi|\hat{A}|\phi\rangle| &= \left| \sum_{k=0}^{M} |\phi_k|^2 \langle k|\hat{A}|k\rangle + \sum_{i \ne j} \phi_i^* \phi_j \langle i|\hat{A}|j\rangle \right| \\ &\ge \sum_{k=0}^{M} |\phi_k|^2 \langle k|\hat{A}|k\rangle - \sum_{i \ne j} |\phi_i||\phi_j||\langle i|\hat{A}|j\rangle| \\ &\ge 1 - 2\delta - 4\delta \sum_{i \ne j} |\phi_i||\phi_j| \\ &\ge 1 - (2 + 4N)\delta, \end{aligned} \quad (A35) $$

where we used Cauchy-Schwarz inequality in the last line. With Lemma 4, this implies

$$ \begin{aligned} \| \hat{A} |\phi\rangle \langle\phi| \hat{A}^\dagger - |\phi\rangle \langle\phi| \|_1 &= \sqrt{(1 + \langle\phi|\hat{A}^\dagger\hat{A}|\phi\rangle)^2 - 4|\langle\phi|\hat{A}|\phi\rangle|^2} \\ &\le \sqrt{(2 + (8 + 12M)\delta)^2 - 4(1 - (2 + 4M)\delta)^2} \\ &= \sqrt{\delta}\sqrt{12(4 + 5M) + 16(3 + 8M + 6M^2)\delta} \\ &\le \sqrt{\delta}\sqrt{12(12 + 9M)}, \end{aligned} \quad (A36) $$

where we used $\delta < \frac{1}{2M}$ in the last line. Hence, with Eq. (A20) we finally obtain

$$ \| \Pi_N \hat{U} |\phi\rangle \langle\phi| \hat{U}^\dagger \Pi_N - \hat{V}_N |\phi\rangle \langle\phi| \hat{V}_N^\dagger \|_1 \le \sqrt{\delta} \sqrt{12(12 + 9M)}. \quad (A37) $$

(iv) Combining Eqs. (A4), (A6) and (A37) yields

$$ \begin{aligned} \| (\mathcal{U} - \mathcal{V})(\psi) \|_1 &\le 4\sqrt{\frac{E}{M}} + 2\sqrt{\frac{E_{\hat{U}}(M)}{N}} + \sqrt{\delta}\sqrt{12(12 + 9M)} \\ &= 4\sqrt{\frac{E}{M}} + \sqrt{\frac{E_{\hat{U}}(M)}{N}}(2 + \sqrt{12(12 + 9M)}). \end{aligned} \quad (A38) $$

where we have used $\delta = \frac{E_{\hat{U}}(M)}{N}$ and where $\mathcal{V} = \hat{V} \cdot \hat{V}^\dagger$ with $\hat{V} = \hat{V}_N \oplus \hat{V}'$, where $\hat{V}_N$ is a unitary operator over $\mathcal{H}_N$ obtained via the orthonormalisation of a column full rank completion of $\Pi_N \hat{U} \Pi_M$ and $\hat{V}'$ is any unitary operator over $\text{Range}(I - \Pi_N)$. Finally, choosing

$$ \begin{aligned} M &= \frac{64E}{\epsilon^2} \\ N &= \frac{4E_{\hat{U}}(M)(2 + \sqrt{12(12 + 9M)})^2}{\epsilon^2} = \mathcal{O}\left(\frac{E}{\epsilon^4} E_{\hat{U}}\left(\frac{64E}{\epsilon^2}\right)\right), \end{aligned} \quad (A39) $$

ensures $\delta = \frac{E_{\hat{U}}(M)}{N} < \frac{1}{2M}$ and

$$ \| (\mathcal{U} - \mathcal{V})(\psi) \|_1 \le \epsilon. \quad (A40) $$

The time complexity of computing $\hat{V}_N$ is given by the time complexity of the QR decomposition [38] of $\Pi_N \hat{U} \Pi_M$, i.e.,

$$ \mathcal{O}(NM^2) = \mathcal{O}\left(\frac{NE^2}{\epsilon^4}\right) = \mathcal{O}\left(\frac{E^3}{\epsilon^8} E_{\hat{U}}\left(\frac{64E}{\epsilon^2}\right)\right). \quad (A41) $$

□
Appendix B: Proof of Theorem 2: Finite-dimensional universality of polynomial Hamiltonians

Hereafter, we write $(\hat{q}, \hat{p}) = (\hat{q}_1, \hat{p}_1, \dots, \hat{q}_m, \hat{p}_m)$ for $m \ge 1$. We recall Theorem 2 from the main text, extended to the multimode setting:

**Theorem 2** (Finite-dimensional universality of polynomial Hamiltonians). Let $m \ge 1$, let $N_1, \dots, N_m \in \mathbb{N}$ and let $\hat{H}$ be a Hermitian operator over $\bigotimes_{k=1}^m \mathcal{H}_{N_k}$. There exists a polynomial Hamiltonian $P_{\hat{H}}(\hat{q}, \hat{p})$ of degree at most $3^m N_1 \cdots N_m$ over $\mathcal{H}^{\otimes m}$ such that

$$ P_{\hat{H}}(\hat{q}, \hat{p}) = \hat{H} \oplus \hat{H}', \qquad (B1) $$

where $\hat{H}'$ is a Hermitian operator over $\text{Range}(I - \bigotimes_{k=1}^m \Pi_{N_k})$. In particular,

$$ e^{iP_{\hat{H}}(\hat{q}, \hat{p})} = e^{i\hat{H}} \oplus e^{i\hat{H}'}, \qquad (B2) $$

and

$$ e^{iP_{\hat{H}}(\hat{q}, \hat{p})} |\psi\rangle = e^{i\hat{H}} |\psi\rangle, \qquad (B3) $$

for all $|\psi\rangle \in \bigotimes_{k=1}^m \mathcal{H}_{N_k}$. Moreover, the polynomial $P_{\hat{H}}$ can be computed efficiently in the size of $\hat{H}$.

*Proof sketch.* The proof proceeds by showing that any finite-dimensional linear operator can be written as the restriction of an infinite-dimensional operator that is a polynomial in canonical bosonic operators, with the property that it does not mix the finite-dimensional subspace with the rest of the Hilbert space. In other terms, in the single-mode case, if $\hat{A}$ is a linear operator acting on the finite-dimensional Hilbert space $\mathcal{H}_N$, we construct a polynomial $P_{\hat{A}}$ such that

$$ P_{\hat{A}}(\hat{q}, \hat{p}) = \Pi_N \hat{A} \Pi_N \oplus (I - \Pi_N) P_{\hat{A}}(\hat{q}, \hat{p}) (I - \Pi_N). \qquad (B4) $$

This construction is based on adding together $(N+1) \times (N+1)$ interpolation polynomials, each reproducing a single entry of the linear operator $\hat{A}$, taking advantage of the sparsity of the canonical operators in Fock basis. Using this result for $\hat{A}$ Hermitian and taking the operator exponential completes the proof in the single-mode case.

The proof then extends to the multimode setting, by taking tensor products and linear combinations of the single-mode polynomials constructed above.

*Proof of Theorem 2.* We first prove the result in the single-mode setting. Let $N \ge 0$ and let $\hat{A}$ be a linear operator on $\mathcal{H}_N$. We construct a polynomial $P_{\hat{A}}$ such that $P_{\hat{A}}(\hat{q}, \hat{p}) = \Pi_N \hat{A} \Pi_N + (I - \Pi_N) P_{\hat{A}}(\hat{q}, \hat{p}) (I - \Pi_N)$. To do so, we only need to treat the cases where $\hat{A} = |n\rangle\langle n|$, $\hat{A} = |n\rangle\langle n+k|$, and $\hat{A} = |n+k\rangle\langle n|$, for $n = 0, \dots, N$, $k = 1, \dots, N$ and $n+k \le N$, since the general case is obtained by taking linear combinations of these elementary cases.

For the case $\hat{A} = |n\rangle\langle n|$, for $n = 0, \dots, N$, we pick $P_{|n\rangle\langle n|}(\hat{q}, \hat{p}) = P_n(\hat{n})$, where $\hat{n} = \hat{a}^\dagger \hat{a} = \frac{1}{2}(\hat{q}^2 + \hat{p}^2 - \hat{I})$ and where $P_n$ is a polynomial satisfying $P_n(m) = \delta_{nm}$ for all integers $0 \le m \le N$. Using the Lagrange interpolation polynomial for these values, we thus define:

$$ P_n(X) := \prod_{\substack{k=0 \\ k \neq n}}^{N} \frac{X-k}{n-k}. \qquad (B5) $$

Then, for all $m \in \mathbb{N}$,

$$ \begin{aligned} P_{|n\rangle\langle n|}(\hat{q}, \hat{p}) |m\rangle &= P_n(\hat{n}) |m\rangle \\ &= \left( \prod_{\substack{k=0 \\ k \neq n}}^{N} \frac{\hat{n}-k}{n-k} \right) |m\rangle \\ &= \left( \prod_{\substack{k=0 \\ k \neq n}}^{N} \frac{m-k}{n-k} \right) |m\rangle \qquad (B6) \\ &= \begin{cases} |n\rangle & m=n \\ 0 & m \neq n \text{ and } m \le N \\ \left( \prod_{\substack{k=0 \\ k \neq n}}^{N} \frac{m-k}{n-k} \right) |m\rangle & m > N. \end{cases} \end{aligned} $$
By construction, $P_{|n\rangle\langle n|}(\hat{q}, \hat{p}) = \Pi_N |n\rangle\langle n| \Pi_N + (I - \Pi_N)P_{|n\rangle\langle n|}(\hat{q}, \hat{p})(I - \Pi_N)$, and $P_{|n\rangle\langle n|}$ has degree $2N$.

For the case $\hat{A} = |n\rangle\langle n+k|$, for $n = 0, \dots, N$, $k = 1, \dots, N$ and $n+k \le N$, we pick $P_{|n\rangle\langle n+k|}(\hat{q}, \hat{p}) = P_{n,k}(\hat{n})\hat{a}^k$, where $P_{n,k}$ is a polynomial to be determined. For all $i, j \in \mathbb{N}$,

$$ \begin{aligned} \langle i | P_{|n\rangle\langle n+k|}(\hat{q}, \hat{p}) | j \rangle &= \langle i | P_{n,k}(\hat{n}) \hat{a}^k | j \rangle \\ &= \begin{cases} P_{n,k}(i) \sqrt{\frac{(i+k)!}{i!}} & j = i+k \\ 0 & \text{otherwise.} \end{cases} \end{aligned} \quad (B7) $$

Hence, to ensure that $P_{|n\rangle\langle n+k|}(\hat{q}, \hat{p}) = \Pi_N |n\rangle\langle n+k| \Pi_N \oplus (I - \Pi_N)P_{|n\rangle\langle n+k|}(\hat{q}, \hat{p})(I - \Pi_N)$, we choose

$$ P_{n,k}(i) = \begin{cases} \sqrt{\frac{n!}{(n+k)!}} & i = n \\ 0 & i \neq n \text{ and } i \le N. \end{cases} \quad (B8) $$

Note that enforcing the last condition for $i \le N-k$ would be sufficient to strictly reproduce the action of the target operator on the finite-dimensional subspace Range($\Pi_N$), but enforcing this condition for $i \le N$ instead ensures that the corresponding infinite-dimensional operator also has the desired block-diagonal structure in Fock basis. Using the Lagrange interpolation polynomial for these values, we thus define:

$$ \begin{aligned} P_{n,k}(X) &:= \sqrt{\frac{n!}{(n+k)!}} \prod_{\substack{m=0 \\ m \neq n}}^{N} \frac{X-m}{n-m} \\ &= \sqrt{\frac{n!}{(n+k)!}} P_n(X). \end{aligned} \quad (B9) $$

This polynomial has degree $N$, so $P_{|n\rangle\langle n+k|}(\hat{q}, \hat{p}) = P_{n,k}(\hat{n})\hat{a}^k$ has degree $2N+k \le 3N$.

For the case $\hat{A} = |n+k\rangle\langle n|$, for $n = 0, \dots, N$, $k = 1, \dots, N$ and $n+k \le N$, we pick $P_{|n+k\rangle\langle n|}(\hat{q}, \hat{p}) = P_{|n\rangle\langle n+k|}(\hat{q}, \hat{p})^\dagger$, which also has degree $2N+k \le 3N$.

For the general case, we write $\hat{A} = \sum_{n=0}^N a_{nn} |n\rangle\langle n| + \sum_{0 \le i < j \le N} a_{ij} |i\rangle\langle j| + a_{ji} |j\rangle\langle i|$. Combining the previous cases, we obtain

$$ P_{|i\rangle\langle j|}(\hat{q}, \hat{p}) = \sqrt{\frac{\min(i, j)!}{\max(i, j)!}} \hat{a}^{\dagger(\max(i, j)-j)} P_{\min(i, j)}(\hat{n}) \hat{a}^{\max(i, j)-i}, \quad (B10) $$

where $P_n$ is defined in Eq. (B5), and we define

$$ \begin{aligned} P_{\hat{A}}(\hat{q}, \hat{p}) &:= \sum_{0 \le i,j \le N} a_{ij} P_{|i\rangle\langle j|}(\hat{q}, \hat{p}) \\ &= \sum_{n=0}^{N} a_{nn} P_n(\hat{n}) + \sum_{0 \le i < j \le N} a_{ij} \sqrt{\frac{i!}{j!}} P_i(\hat{n}) \hat{a}^{j-i} + a_{ji} \sqrt{\frac{j!}{i!}} \hat{a}^{\dagger(j-i)} P_i(\hat{n}), \end{aligned} \quad (B11) $$

which has degree less or equal to $3N$. By construction, we have

$$ P_{\hat{A}}(\hat{q}, \hat{p}) = \Pi_N \hat{A} \Pi_N + (I - \Pi_N) P_{\hat{A}}(\hat{q}, \hat{p})(I - \Pi_N), \quad (B12) $$

and $P_{\hat{A}}(\hat{q}, \hat{p})$ is Hermitian when $\hat{A}$ is Hermitian.

Let $\hat{U}$ be a unitary operator over $\mathcal{H}_M$. There exists a Hermitian operator $\hat{H}$ over $\mathcal{H}_M$ such that $\hat{U} = e^{i\hat{H}}$. Then, $P_{\hat{H}}(\hat{q}, \hat{p}) = \Pi_N \hat{H} \Pi_N + (I - \Pi_N) P_{\hat{H}}(\hat{q}, \hat{p})(I - \Pi_N)$ and

$$ \begin{aligned} e^{i P_{\hat{H}}(\hat{q}, \hat{p})} &= e^{i [\Pi_N \hat{H} \Pi_N + (I - \Pi_N) P_{\hat{H}}(\hat{q}, \hat{p})(I - \Pi_N)]} \\ &= e^{i \Pi_N \hat{H} \Pi_N} + e^{i (I - \Pi_N) P_{\hat{H}}(\hat{q}, \hat{p})(I - \Pi_N)} \\ &= \Pi_N e^{i \hat{H}} \Pi_N + (I - \Pi_N) e^{i (I - \Pi_N) P_{\hat{H}}(\hat{q}, \hat{p})(I - \Pi_N)} (I - \Pi_N) \\ &= \Pi_N \hat{U} \Pi_N + (I - \Pi_N) \hat{U}' (I - \Pi_N), \end{aligned} \quad (B13) $$
where we have set $\hat{U}' := e^{i(I-\Pi_N)P_{\hat{H}}}(\hat{q},\hat{p})(I-\Pi_N)$.

We now turn to the multimode setting. Let $m \ge 1$ be the number of modes, and let $\hat{A}$ be a linear operator over $\bigotimes_{k=1}^m \mathcal{H}_{N_k}$, for $N_1, \dots, N_m \in \mathbb{N}$. Writing $\hat{A} = \sum_{k=0}^m \sum_{i_k, j_k=0}^{N_k} a_{ij} |i\rangle\langle j|$ and $(\hat{q}, \hat{p}) = (\hat{q}_1, \hat{p}_1, \dots, \hat{q}_m, \hat{p}_m)$, we define

$$P_{\hat{A}}(\hat{q}, \hat{p}) := \sum_{k=0}^{m} \sum_{i_k, j_k=0}^{N_k} a_{ij} \bigotimes_{k=1}^{m} P_{|i_k\rangle\langle j_k|}(\hat{q}_k, \hat{p}_k), \quad (B14)$$

which is a polynomial of degree at most $3^m \times N_1 \cdots N_m$, where $P_{|i_k\rangle\langle j_k|}(\hat{q}_k, \hat{p}_k)$ is the single-mode polynomial defined in Eq. (B10). By Eq. (B12), this polynomial satisfies

$$P_{\hat{A}}(\hat{q}, \hat{p}) = \sum_{k=0}^{m} \sum_{i_k, j_k=0}^{N_k} a_{ij} \bigotimes_{k=1}^{m} [\Pi_{N_k} |i_k\rangle\langle j_k| \Pi_{N_k} + (I - \Pi_{N_k}) P_{|i_k\rangle\langle j_k|}(\hat{q}_k, \hat{p}_k)(I - \Pi_{N_k})]. \quad (B15)$$

Expanding the tensor product leads to a single product term of the form

$$\bigotimes_{k=1}^{m} (\Pi_{N_k} |i_k\rangle\langle j_k| \Pi_{N_k}) = \left(\bigotimes_{k=1}^{m} \Pi_{N_k}\right) |i\rangle\langle j| \left(\bigotimes_{k=1}^{m} \Pi_{N_k}\right), \quad (B16)$$

while all the other product terms contain at least one term of the form $(I - \Pi_{N_k}) P_{|i_k\rangle\langle j_k|}(\hat{q}_k, \hat{p}_k)(I - \Pi_{N_k})$, for some $k \in \{1, \dots, m\}$. As a consequence, these terms can be all be expressed as operators over $\text{Range}(I - \bigotimes_{k=1}^m \Pi_{N_k})$, and thus

$$\begin{align}
P_{\hat{A}}(\hat{q}, \hat{p}) &= \sum_{k=0}^{m} \sum_{i_k, j_k=0}^{N_k} a_{ij} \left( \bigotimes_{k=1}^{m} \Pi_{N_k} \right) |i\rangle\langle j| \left( \bigotimes_{k=1}^{m} \Pi_{N_k} \right) + \left( I - \bigotimes_{k=1}^{m} \Pi_{N_k} \right) P_{\hat{A}}(\hat{q}, \hat{p}) \left( I - \bigotimes_{k=1}^{m} \Pi_{N_k} \right) \nonumber \\
&= \left( \bigotimes_{k=1}^{m} \Pi_{N_k} \right) \sum_{k=0}^{m} \sum_{i_k, j_k=0}^{N_k} a_{ij} |i\rangle\langle j| \left( \bigotimes_{k=1}^{m} \Pi_{N_k} \right) + \left( I - \bigotimes_{k=1}^{m} \Pi_{N_k} \right) P_{\hat{A}}(\hat{q}, \hat{p}) \left( I - \bigotimes_{k=1}^{m} \Pi_{N_k} \right) \tag{B17} \\
&= \left( \bigotimes_{k=1}^{m} \Pi_{N_k} \right) \hat{A} \left( \bigotimes_{k=1}^{m} \Pi_{N_k} \right) + \left( I - \bigotimes_{k=1}^{m} \Pi_{N_k} \right) P_{\hat{A}}(\hat{q}, \hat{p}) \left( I - \bigotimes_{k=1}^{m} \Pi_{N_k} \right). \nonumber
\end{align}$$

Moreover, this operator is Hermitian when $\hat{A}$ is Hermitian. With the same derivation as in Eq. (B13) this concludes the proof. $\square$

## Appendix C: Proof of Theorem 4: Solovay-Kitaev theorem for polynomial Hamiltonians

We give a formal version of Theorem 4 from the main text:

**Theorem 4 (Solovay-Kitaev theorem for polynomial Hamiltonians).** Let $E > 0$, $\epsilon > 0$ and $N \ge \frac{64E}{\epsilon^2} \in \mathbb{N}$. Let $\mathcal{G}$ be a finite set of unitary operators over $\mathcal{H}_N$ generating a dense subset of $\mathcal{U}(\mathcal{H}_N)$ and let $\mathcal{P}$ be its realization with polynomial Hamiltonians from Theorem 2. There is a constant $c$ such that for any physical unitary operator $\hat{U} \in \mathcal{U}(\mathcal{S})$ with $(E, \epsilon)$-approximate effective dimension $N+1$, there exists a finite sequence $\hat{V}$ of gates from $\mathcal{P}$ of length $O(\log^c(1/\epsilon))$ and such that $\|\mathcal{U} - \mathcal{V}\|_{\diamond}^E \le 2\epsilon$, where $\mathcal{U} = \hat{U} \cdot \hat{U}^{\dagger}$ and $\mathcal{V} = \hat{V} \cdot \hat{V}^{\dagger}$.

*Proof.* The physical unitary operator $\hat{U}$ has $(E, \epsilon)$-approximate effective dimension $N+1$, so with Theorem 1 there exists a unitary operator $\hat{V}_N$ over $\mathcal{H}_N$ such that for any unitary operator $\hat{V}'$ over $\text{Range}(I - \Pi_N)$,

$$\sup_{\langle \psi | \hat{n} | \psi \rangle \le E} D[(\hat{V}_N \oplus \hat{V}') |\psi\rangle, \hat{U} |\psi\rangle] \le \frac{\epsilon}{2}, \quad (C1)$$

where we have expressed the energy-constrained diamond norm using the trace distance $D(\cdot, \cdot) = \frac{1}{2}\|\cdot - \cdot\|_1$ as in Eq. (A3).
The Solovay-Kitaev theorem for qudits [22, Theorem 1] of dimension $N + 1$ ensures that there is a constant $c$ such that there exists a finite sequence $\hat{S}$ of gates from $\mathcal{G}$ of length $\mathcal{O}(\log^c(1/\epsilon))$ and such that

$$ \sup_{|\phi\rangle \in \mathcal{H}_N} D(\hat{V}_N |\phi\rangle, \hat{S} |\phi\rangle) \le \frac{\epsilon}{4}, \qquad (C2) $$

and thus

$$ \begin{aligned} |\langle\phi|\hat{V}_N^\dagger \hat{S}|\phi\rangle| &= \sqrt{1 - D(\hat{V}_N |\phi\rangle, \hat{S} |\phi\rangle)} \\ &\ge \sqrt{1 - \frac{\epsilon^2}{16}}, \end{aligned} \qquad (C3) $$

for all $|\phi\rangle \in \mathcal{H}_N$.

Now for all $|\psi\rangle = \sum_{n \ge 0} \psi_n |n\rangle \in \mathcal{H}$ such that $\langle\psi|\hat{n}|\psi\rangle \le E$,

$$ \begin{aligned} \langle\psi|\Pi_N|\psi\rangle &= \sum_{n=0}^{N} |\psi_n|^2 \\ &= 1 - \sum_{n>N} |\psi_n|^2 \\ &\ge 1 - \sum_{n>N} \frac{n}{N} |\psi_n|^2 \qquad (C4) \\ &\ge 1 - \frac{1}{N} \sum_{n \ge 0} n |\psi_n|^2 \\ &\ge 1 - \frac{E}{N}. \end{aligned} $$

Moreover, for all $|\psi\rangle = \sum_{n \ge 0} \psi_n |n\rangle \in \mathcal{H}$ such that $\langle\psi|\hat{n}|\psi\rangle \le E$, and for any unitary operator $\hat{V}'$ on Range$(I - \Pi_N)$,

$$ \begin{aligned} D[(\hat{V}_N \oplus \hat{V}') |\psi\rangle, (\hat{S} \oplus \hat{V}') |\psi\rangle] &= \sqrt{1 - |\langle\psi|(\hat{V}_N \oplus \hat{V}')^\dagger(\hat{S} \oplus \hat{V}')|\psi\rangle|^2} \\ &= \sqrt{1 - |\langle\psi|(\hat{V}_N^\dagger \hat{S} \oplus \hat{I})|\psi\rangle|^2} \\ &= \sqrt{1 - |\langle\psi|\Pi_N \hat{V}_N^\dagger \hat{S} \Pi_N |\psi\rangle + \langle\psi|I - \Pi_N |\psi\rangle|^2} \\ &= \sqrt{1 - |\langle\psi|\Pi_N |\psi\rangle \langle\psi_N|\hat{V}_N^\dagger \hat{S} |\psi_N\rangle + 1 - \langle\psi|\Pi_N |\psi\rangle|^2}, \end{aligned} \qquad (C5) $$

where we have defined $|\psi_N\rangle := \Pi_N |\psi\rangle / \|\Pi_N |\psi\rangle\| \in \mathcal{H}_N$. Using the reverse triangle inequality $|a + b| \ge ||a| - |b||$ we obtain

$$ \begin{aligned} D[(\hat{V}_N \oplus \hat{V}') |\psi\rangle, (\hat{S} \oplus \hat{V}') |\psi\rangle] &\le \sqrt{1 - \left[\langle\psi|\Pi_N |\psi\rangle \left|\langle\psi_N|\hat{V}_N^\dagger \hat{S} |\psi_N\rangle\right| - (1 - \langle\psi|\Pi_N |\psi\rangle)\right]^2} \\ &\le \sqrt{1 - \left[\left(1 - \frac{E}{N}\right) \sqrt{1 - \frac{\epsilon^2}{16}} - \frac{E}{N}\right]^2} \\ &= \sqrt{\frac{\epsilon^2}{16} + \frac{2E}{N} \left(\sqrt{1 - \frac{\epsilon^2}{16}} + 1 - \frac{\epsilon^2}{16}\right) - \frac{E^2}{N^2} \left(1 + \sqrt{1 - \frac{\epsilon^2}{16}}\right)^2} \qquad (C6) \\ &\le \sqrt{\frac{\epsilon^2}{16} + \frac{4E}{N}} \\ &\le \frac{\epsilon}{4} + 2\sqrt{\frac{E}{N}} \\ &\le \frac{\epsilon}{2}, \end{aligned} $$
where we have used Eqs. (C3) and (C4) in the second line and $N \ge \frac{64E}{\epsilon^2}$ in the last line. Combining this with Eq. (C1) and the triangle inequality we obtain

$$ \sup_{\langle \psi | \hat{n} | \psi \rangle \le E} D[(\hat{S} \oplus \hat{V}') |\psi\rangle, \hat{U} |\psi\rangle] \le \epsilon, \qquad (C7) $$

for any unitary operator $\hat{V}'$ over Range($I - \Pi_N$).

Let $s$ denote the length of the sequence $\hat{S}$ given by the qudit Solovay-Kitaev theorem. For each unitary gate $\hat{G}_k \in \mathcal{G}$ over $\mathcal{H}_N$ appearing in the sequence $\hat{S} = \hat{G}_s \cdots \hat{G}_1$, there exists a Hermitian operator $\hat{H}_k$ over $\mathcal{H}_N$ such that $\hat{G}_k = e^{i\hat{H}_k}$, which can be computed efficiently in $N$, for instance by diagonalizing $\hat{G}_k$. For each gate, Theorem 2 provides a realization generated by a polynomial Hamiltonian $P_{\hat{H}_k}(\hat{q}, \hat{p})$ of degree $3N$, such that

$$ P_{\hat{H}_k}(\hat{q}, \hat{p}) = \hat{H}_k \oplus \hat{H}'_k, \qquad (C8) $$

where $\hat{H}'_k$ is a Hermitian operator over Range($I - \Pi_N$), and

$$ e^{iP_{\hat{H}_k}(\hat{q}, \hat{p})} = \hat{G}_k \oplus e^{i\hat{H}'_k}, \qquad (C9) $$

so that

$$ e^{iP_{\hat{H}_s}(\hat{q}, \hat{p})} \cdots e^{iP_{\hat{H}_1}(\hat{q}, \hat{p})} = \hat{S} \oplus \hat{V}', \qquad (C10) $$

where $\hat{V}'$ is a Hermitian operator over Range($I - \Pi_N$). With Eq. (C7) we finally obtain

$$ \sup_{\langle \psi | \hat{n} | \psi \rangle \le E} D[e^{iP_{\hat{H}_s}(\hat{q}, \hat{p})} \cdots e^{iP_{\hat{H}_1}(\hat{q}, \hat{p})} |\psi\rangle, \hat{U} |\psi\rangle] \le \epsilon, \qquad (C11) $$

which concludes the proof. □