Effective descriptions of bosonic systems can be considered complete

TL;DR

This paper proves that simplified models of bosonic systems accurately capture their true physics, enabling reliable simulations and quantum computing applications.

Contribution

The authors rigorously show that finite-dimensional approximations and polynomial Hamiltonians fully capture bosonic system behavior.

Findings

Any physical bosonic unitary evolution can be approximated by a finite-dimensional evolution.

Finite-dimensional evolutions can be generated exactly by polynomial bosonic Hamiltonians.

Abstract

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: typically, truncating their infinite-dimensional state space or restricting their dynamics to a simple class of Hamiltonians, such as polynomials of canonical operators. However, many natural bosonic Hamiltonians do not belong to these simple classes, and some quantum effects harnessed by bosonic computers inherently require infinite-dimensional spaces. Can we trust that results obtained with such simplifying assumptions capture real effects? We solve this outstanding problem, showing that these effective descriptions do correctly capture the physics of bosonic systems. Our technical contributions are twofold: firstly, we prove that any physical bosonic unitary evolution can be accurately 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, provide universal methods for engineering bosonic quantum states and Hamiltonians, show that polynomial Hamiltonians generate universal gate sets for quantum computing over bosonic modes, and lead to a bosonic Solovay–Kitaev theorem.

Bosonic systems live in an infinite-dimensional space, and in order to be able to describe them one usually reduces it to an effective finite dimension or constrain the dynamics in some way, but it is not fully understood whether this captures all the physics at play. Here, the authors fill this gap showing that such representations can successfully and rigorously approximate bosonic physics.

Figures

Completeness of effective descriptions of bosonic systems.Top: physical unitary dynamics of bosonic systems can be cut off to an effective dimension (Theorem 1), leading to rigorous precision guarantees for simulating these systems. Bottom: effective polynomial Hamiltonians can realize any finite-dimensional Hamiltonian (Theorem 2), enabling universal quantum control with these Hamiltonians and opening the way to universal bosonic simulators. Middle: combining the previous results, we design universal gate sets for bosonic computations based on effective polynomial Hamiltonians (Theorem 3), wh

Full text

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 satisfying the canonical commutation relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\hat{q},\hat{p}]=i\hat{I}$$\end{document} (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). 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^8,9^, and to draw equivalences with spin systems^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 ϕ^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^. Furthermore, 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–18^. 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos \hat{\delta }$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\delta }$$\end{document} is a canonical phase-difference operator^19^. Here, the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos$$\end{document} cannot be replaced by a polynomial—such as its truncated Taylor expansion—because the operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\delta }$$\end{document} is unbounded^20^. 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 quantum control and was first formalized in this context as an open problem in ref. ^21^. Beyond this controllability problem, it is also highly non-trivial 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^22^. 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^23^. Such a theorem is missing in the general infinite-dimensional setting^24^ 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^25^. 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 (see Fig. 1). Our technical contributions are twofold:Fig. 1. Completeness of effective descriptions of bosonic systems.Top: physical unitary dynamics of bosonic systems can be cut off to an effective dimension (Theorem 1), leading to rigorous precision guarantees for simulating these systems. Bottom: effective polynomial Hamiltonians can realize any finite-dimensional Hamiltonian (Theorem 2), enabling universal quantum control with these Hamiltonians and opening the way to universal bosonic simulators. Middle: combining the previous results, we design universal gate sets for bosonic computations based on effective polynomial Hamiltonians (Theorem 3), which satisfy a Solovay--Kitaev theorem (Theorem 4).

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, thus also yielding an algorithm for the simulation of bosonic dynamics on digital quantum computers with explicit space complexity (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 quantum controllability of polynomial Hamiltonians in infinite dimension from ref.^21^ (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 ref.^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).

Results

Bosonic quantum information theory

Quantum information theory with bosonic systems typically involves modeling unbounded operators over infinite-dimensional Hilbert spaces with continuous spectrum, and is thus often referred to as continuous-variable quantum information theory^26–28^. 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}^{\otimes m}$$\end{document} , with each single-mode Hilbert space spanned by a countably infinite (Fock) basis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\{\left\vert n\right\rangle \}}_{n\in {\mathbb{N}}}$$\end{document} .

Creation and annihilation operators \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{a}}^{{\dagger} }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{a}$$\end{document} for each mode are defined by their action on the Fock basis as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{a}}^{{\dagger} }\left\vert n\right\rangle=\sqrt{n+1}\left\vert n+1\right\rangle$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{a}\left\vert n+1\right\rangle=\sqrt{n+1}\left\vert n\right\rangle$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{a}\left\vert 0\right\rangle=0$$\end{document} . The (unbounded) canonical bosonic operators are then given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{q}=\frac{1}{\sqrt{2}}(\hat{a}+{\hat{a}}^{{\dagger} })$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}=\frac{1}{i\sqrt{2}}(\hat{a}-{\hat{a}}^{{\dagger} })$$\end{document} and satisfy the canonical commutation relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\hat{q},\hat{p}]=i\hat{I}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{I}$$\end{document} is the identity operator. Polynomial Hamiltonians over m modes are the Hermitian operators of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\hat{q}}_{1},{\hat{p}}_{1},\ldots,{\hat{q}}_{m},{\hat{p}}_{m})$$\end{document} , where P is a polynomial and where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\hat{q}}_{1},{\hat{p}}_{1},\ldots,{\hat{q}}_{m},{\hat{p}}_{m})$$\end{document} are the canonical operators of the modes 1, …, m. The number operator is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{n}}={\hat{a}}^{{\dagger} }{\hat{a}}$$\end{document} , and its expectation value is the average particle number of a bosonic mode, which we refer to as energy hereafter. For all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\in {\mathbb{N}}$$\end{document} , we denote by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}: \!\!={{\rm{span}}}{\{\left\vert n\right\rangle \}}_{0\le n\le d}$$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{6}{\sum }_{n\ge 1}\frac{1}{\pi n}\left\vert n\right\rangle$$\end{document} . 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^20^—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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{S}}}\subset {{\mathcal{H}}}$$\end{document} (see the Methods section). Similarly, denoting by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}({{\mathcal{H}}})$$\end{document} the group of unitary operators over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} , we define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}({{\mathcal{S}}}) : \!\!=\left\{{\hat{U}}\in {{\mathcal{U}}({{\mathcal{H}}})\,| \,\forall \left\vert \psi \right\rangle \in {{\mathcal{S}}},\,{\hat{U}}\left\vert \psi \right\rangle \in {{\mathcal{S}}}}\right\}.$$\end{document}

These form a group of physical unitary operators, i.e., which map physical states to physical states. This set 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in {\mathbb{N}}$$\end{document} and all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}\in {{\mathcal{U}}}({{\mathcal{S}}})$$\end{document} ,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{\hat{U}}(n) : \!\!={\sup }_{\left\vert \psi \right\rangle \in {{{\mathcal{H}}}}_{n}}\left\langle \psi \right\vert {\hat{U}}^{{\dagger} }\hat{n}\hat{U}\left\vert \psi \right\rangle < +\infty .$$\end{document}

In that case, the supremum is a maximum and the quantity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{\hat{U}}(n)$$\end{document} can be thought of as the maximal amount of energy involved when implementing the unitary evolution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} of an initial state with a number of particles at most n. This is a slightly weaker property than energy-limited quantum dynamics^29^ which map state of finite energy to states of finite energy.

We denote by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D(\rho,\sigma )=\frac{1}{2}\parallel \rho -\sigma {\parallel }_{1}$$\end{document} the trace norm distance between density operators^22^, 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^30^. 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^31^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\parallel {{\mathcal{E}}}{\parallel }_{\diamond }^{E} : \!\!={\sup }_{{{\rm{Tr}}}[\rho (\hat{H}\otimes {\hat{I}})]\le E}\parallel ({{\mathcal{E}}}\otimes {{\rm{id}}})\rho {\parallel }_{1}.$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}$$\end{document} is a Hermitian-preserving map, ρ is a density operator over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}\otimes {{{\mathcal{H}}}}^{{\prime} }$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}^{{\prime} }$$\end{document} isomorphic to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} is a Hamiltonian over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} with ground state energy equal to 0 used to specify the energy bound E > 0. We employ this topology hereafter, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}=\hat{n}$$\end{document} the number operator.

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} be a single-mode unitary operator, let E≥0 be an energy parameter, and let ϵ > 0 be an approximation parameter. Let also \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\in {\mathbb{N}}$$\end{document} . The unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} has (E, ϵ)-approximate effective dimension d + 1 if there exists a unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} which approximates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} in the following sense: for any unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}^{{\prime} }$$\end{document} over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{span}}}{(\left\vert n\right\rangle )}_{n\ > \ d}$$\end{document} , denoting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{V}={\hat{V}}_{d}\oplus {\hat{V}}^{{\prime} }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{V}}}=\hat{V}{{\boldsymbol{\cdot }}}{\hat{V}}^{{\dagger} }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}=\hat{U}{{\boldsymbol{\cdot }}}{\hat{U}}^{{\dagger} }$$\end{document} ,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\parallel {{\mathcal{U}}}-{{\mathcal{V}}}{\parallel }_{\diamond }^{E}\le \epsilon .$$\end{document}

When that is the case, we say that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} is a (d + 1)-dimensional (E, ϵ)-approximation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} .

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}=\hat{U}{{\boldsymbol{\cdot }}}{\hat{U}}^{{\dagger} }$$\end{document} be a single-mode unitary channel, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}\in {{\mathcal{U}}}({{\mathcal{S}}})$$\end{document} a physical unitary operator. Let E≥0 be an energy parameter, and let ϵ > 0 be an approximation parameter. Then, there exists \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d={{\mathcal{O}}}\left(\frac{E}{{\epsilon }^{4}}{E}_{\hat{U}}\left(\frac{64E}{{\epsilon }^{2}}\right)\right)\in {\mathbb{N}}$$\end{document} such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} has (E, ϵ)-approximate effective dimension d + 1. Moreover, a finite-dimensional (E, ϵ)-approximation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} can be computed in polynomial time in d, E, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\epsilon }$$\end{document} .

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 the Supplementary Information for the proof.

A direct consequence of Theorem 1 is that the corresponding sequence of (E, ϵ)-approximations strongly converges to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} as ϵ 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^9^. 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{\hat{U}}$$\end{document} 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^25^. We give an explicit construction in the Supplementary Information showing that this growth may be arbitrarily fast.

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\in {\mathbb{N}}$$\end{document} and let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} be a Hermitian operator over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} . There exists a polynomial Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\hat{H}}(\hat{q},\hat{p})$$\end{document} of degree at most 3d over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\hat{H}}(\hat{q},\hat{p})=\hat{H}\oplus {\hat{H}}^{{\prime} },$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{H}}^{{\prime} }$$\end{document} is a Hermitian operator over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{span}}}{(\left\vert n\right\rangle )}_{n\ > \ d}$$\end{document} . In particular, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\hat{H}}(\hat{q},\hat{p})$$\end{document} generates an evolution in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} given for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert \psi \right\rangle \in {{{\mathcal{H}}}}_{d}$$\end{document} by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{i{P}_{\hat{H}}(\hat{q},\hat{p})}\left\vert \psi \right\rangle={e}^{i\hat{H}}\left\vert \psi \right\rangle,$$\end{document}

Moreover, the polynomial \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\hat{H}}$$\end{document} 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 the Supplementary Information for the multimode generalisation of Theorem 2 and its proof. In particular, defining the unitary evolution operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{i{P}_{\hat{H}}(\hat{q},\hat{p})}$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} is highly non-trivial: we restrict its action to states in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} to avoid such mathematical technicalities.

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 quantum controllability of polynomial Hamiltonians in infinite dimension). Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert \psi \right\rangle \in {{\mathcal{H}}}$$\end{document} . For all ϵ > 0, there exists a polynomial Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\epsilon }(\hat{q},\hat{p})$$\end{document} such that evolving the vacuum state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert 0\right\rangle$$\end{document} under the Schrödinger equation with Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\epsilon }(\hat{q},\hat{p})$$\end{document} for constant time yields a state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\psi }_{\epsilon }\right\rangle$$\end{document} satisfying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D(\left\vert \psi \right\rangle,\left\vert {\psi }_{\epsilon }\right\rangle )\le \epsilon$$\end{document} .

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

Proof

Since the state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert \psi \right\rangle$$\end{document} is normalized, there exists dϵ such it is ϵ-close in trace distance to the state truncated at dϵ and renormalized, which we denote by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\psi }_{\epsilon }\right\rangle$$\end{document} . Then, there exists a finite-dimensional unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{U}}_{\epsilon }={e}^{i{\hat{H}}_{\epsilon }}$$\end{document} over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{{d}_{\epsilon }}$$\end{document} which maps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert 0\right\rangle$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\psi }_{\epsilon }\right\rangle$$\end{document} . By Theorem 2, the polynomial Hamiltonial \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\epsilon }({\hat{q}}, {\hat{p}}) : \!\!={P}_{{\hat{H}}_{\epsilon }}({\hat{q}},{\hat{p}})$$\end{document} satisfies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\epsilon }({\hat{q}},{\hat{p}})={\hat{H}}_{\epsilon }\oplus {\hat{H}}_{\epsilon }^{{\prime} }$$\end{document} . Evolving the vacuum under the Schrödinger equation with Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{\epsilon }(\hat{q},\hat{p})$$\end{document} thus leads to the same state as if evolving the vacuum under the (dϵ + 1)-dimensional Schrödinger equation with Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{H}}_{\epsilon }$$\end{document} , which concludes the proof. □

Beyond quantum state preparation, Theorem 2 provides a universal method for exactly reproducing any finite-dimensional Hamiltonian using polynomial Hamiltonians. This method can be used to engineer universal bosonic simulators capable of emulating the unitary evolution of any discrete-variable 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 [ref. ^25^, Theorem 4.1].

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 a 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}=\hat{U}{{\boldsymbol{\cdot }}}{\hat{U}}^{{\dagger} }$$\end{document} be a single-mode unitary channel, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}\in {{\mathcal{U}}}({{\mathcal{S}}})$$\end{document} a physical unitary operator. Let E≥0 be an energy parameter, and let ϵ > 0 be an approximation parameter. Then, there exists \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d={{\mathcal{O}}}\left(\frac{E}{{\epsilon }^{4}}{E}_{\hat{U}}\left(\frac{64E}{{\epsilon }^{2}}\right)\right)\in {\mathbb{N}}$$\end{document} and a polynomial Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\hat{q},\hat{p})$$\end{document} of degree at most 3d such that, writing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{V}}}$$\end{document} the corresponding unitary evolution induced on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} ,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\parallel {{\mathcal{U}}}-{{\mathcal{V}}}{\parallel }_{\diamond }^{E}\le \epsilon .$$\end{document}

Moreover, the polynomial P can be computed from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} in polynomial time in d, E, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\epsilon }$$\end{document} .

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

As such, this shows that polynomial Hamiltonians generate truly universal gate sets and paves the way for a Solovay–Kitaev theorem^23^ 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, ϵ > 0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\ge \frac{64E}{{\epsilon }^{2}}\in {\mathbb{N}}$$\end{document} . Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{G}}}$$\end{document} be a finite set of unitary operators over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} generating a dense subset of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}({{{\mathcal{H}}}}_{d})$$\end{document} and let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{P}}}$$\end{document} be its realization with polynomial Hamiltonians from Theorem 2. There is a constant c such that for any physical unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}\in {{\mathcal{U}}}({{\mathcal{S}}})$$\end{document} with (E, ϵ)-approximate effective dimension d + 1, there exists a finite sequence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{V}$$\end{document} of gates from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{P}}}$$\end{document} of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{O}}}({\log }^{c}(1/\epsilon ))$$\end{document} and such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\parallel {{\mathcal{U}}}-{{\mathcal{V}}}{\parallel }_{\diamond }^{E}\le 2\epsilon$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}=\hat{U}{{\boldsymbol{\cdot }}}{\hat{U}}^{{\dagger} }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{V}}}=\hat{V}{{\boldsymbol{\cdot }}}{\hat{V}}^{{\dagger} }$$\end{document} .

Recall that Theorem 1 ensures the existence of an (E, ϵ)-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^23^ of dimension d + 1, and provides an explicit algorithm for producing the sequence of polynomial Hamiltonians, given the target unitary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}\in {{\mathcal{U}}}({{\mathcal{S}}})$$\end{document} . We refer to the Methods section for the main techniques and to the Supplementary Information for the proof.

Theorem 4 generalises the Gaussian Solovay–Kitaev theorem from ref. ^24^, thus providing an answer to [^25^, 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.

Discussion

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, summarised in Fig. 1: (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 quantum controllability of polynomial Hamiltonians^21^; (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^25^.

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^22,32^.

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

On a different note, it has been argued in several works^11,15,21^ 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^37^. 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^23^. 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^38^.

Methods

Preliminary material

CV quantum information

In continuous-variable quantum information^26–28^, information is encoded in the state of m bosonic modes, which are mathematically modelled by a Hilbert space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} spanned by a countably-infinite Fock basis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\{\left\vert n\right\rangle \}}_{n\in {\mathbb{N}}}$$\end{document} , i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sum }_{n\in {\mathbb{N}}}{c}_{n}\left\vert n\right\rangle \in {{\mathcal{H}}}$$\end{document} if and only if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sum }_{n\in {\mathbb{N}}}| {c}_{n}{| }^{2} < \infty$$\end{document} .

To a single bosonic mode, we associate the annihilation operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{a}\left\vert n\right\rangle=\sqrt{n}\left\vert n-1\right\rangle$$\end{document} and creation operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{a}}^{{\dagger} }\left\vert n\right\rangle=\sqrt{n+1}\left\vert n+1\right\rangle$$\end{document} which satisfy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\hat{a},{\hat{a}}^{{\dagger} }]=\hat{I}$$\end{document} . The canonical operators (also called position and momentum or quadrature operators) are defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{q}=(\hat{a}+{\hat{a}}^{{\dagger} })/\sqrt{2}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}=(\hat{a}-{\hat{a}}^{{\dagger} })/i\sqrt{2}$$\end{document} , 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\hat{q}}_{1},{\hat{q}}_{2},\ldots,{\hat{p}}_{1},{\hat{p}}_{2},\ldots )$$\end{document} . In^11^, an informal 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^26^.

Although it is a mathematically convenient setting for analysing bosonic systems, the Hilbert space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} 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^20^. The *Schwartz space * \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{S}}}$$\end{document} is defined as the dense subspace of elements \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sum }_{n\in {\mathbb{N}}}{c}_{n}\left\vert n\right\rangle \in {{\mathcal{H}}}$$\end{document} such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lim }_{n\to 0}{c}_{n}{n}^{k}=0$$\end{document} for any natural number k^39^. The Schwartz space is not a Hilbert space, but it retains many of the useful analytical properties of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} whilst eliminating these physical inconsistencies: every state in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{S}}}$$\end{document} has bounded energy, and there is a well-defined action of polynomials of the canonical operators on those states^21^. 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}({{\mathcal{H}}})$$\end{document} be the group of unitaries on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} , the group of physical unitaries is defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}({{\mathcal{S}}}) :\!\!=\{\hat{U}\in {{\mathcal{U}}}({{\mathcal{H}}})| \forall \left\vert \psi \right\rangle \in {{\mathcal{S}}}:\hat{U}\left\vert \psi \right\rangle \in {{\mathcal{S}}}\}.$$\end{document}

The elements of this group are closely related to the notion of Schwartz operators defined in^40^. We note that, in fact, all our results involving physical unitary operators (Theorems 1, 3 and 4) are valid under the weaker condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle n\right\vert {\hat{U}}^{{\dagger} }\hat{n}\hat{U}\left\vert n\right\rangle < +\infty$$\end{document} for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\in {\mathbb{N}}$$\end{document} , rather than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}\in {{\mathcal{U}}}({{\mathcal{S}}})$$\end{document} .

Distance measures

In making claims about approximating operators on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{S}}}$$\end{document} , we must be explicit about which distance metric we are using to measure the quality of our approximations. The diamond norm for completely-positive channels is a natural candidate. It is defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${| | {{\mathcal{E}}}| | }_{\diamond } : \!\!=\sup \{{| | ({{\mathcal{E}}}\otimes {{\rm{id}}})\rho | | }_{1}| \rho \,{{\rm{density}}}\,{{\rm{operator}}}\,{{\rm{on}}}\,{{\mathcal{H}}}\otimes {{\mathcal{H}}}\},$$\end{document}

where ∣∣ ⋅ ∣∣1 denotes the trace norm. However, it has some undesirable properties: for instance, if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} is any unbounded Hermitian operator, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{U}}}}_{H,s}:\rho \mapsto {e}^{is\hat{H}}\rho {e}^{-is\hat{H}}$$\end{document} is the associated one-parameter family of unitary channels, then there are density operators ρ such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${| | {{{\mathcal{U}}}}_{H,s}(\rho )-{{{\mathcal{U}}}}_{H,s+\epsilon }(\rho )| | }_{1}$$\end{document} does not go to 0 even as ϵ does, so that the diamond norm \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${| | {{{\mathcal{U}}}}_{H,s}-{{{\mathcal{U}}}}_{H,s+\epsilon }| | }_{\diamond }$$\end{document} also does not go to 0 [ref. ^30^, 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^30,31^. They are defined by restricting the supremum to states whose energy is bounded by a parameter E > 0:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\parallel {{\mathcal{E}}}{\parallel }_{\diamond }^{E} :={\sup }_{{{\rm{Tr}}}[\rho (\hat{H}\otimes \hat{I})]\le E}\parallel ({{\mathcal{E}}}\otimes {{\rm{id}}})\rho {\parallel }_{1},$$\end{document}

where now \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} is a reference energy operator. Mathematically, the energy-constrained diamond norms are related to the strong operator topology, which gives the correct notion of pointwise continuity in time for unitary evolutions:

Proposition 1

(ref. ^31^, Proposition 3). If, for any E > 0 and quantum channel \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}$$\end{document} , the sequence of quantum channels \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\{{{{\mathcal{E}}}}_{p}\}}_{p\in {\mathbb{N}}}$$\end{document} is such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lim }_{p\to \infty }{| | {{{\mathcal{E}}}}_{p}-{{\mathcal{E}}}| | }_{\diamond }^{E}=0$$\end{document} then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\{{{{\mathcal{E}}}}_{p}\}}_{p\in {\mathbb{N}}}$$\end{document} converges to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}$$\end{document} in the strong operator topology.

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} 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≥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^41^, Lemma 9.4.2). Consider a density operator ρ and a projector Π where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\le \Pi \le {\mathbb{I}}$$\end{document} . Suppose that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{Tr}}}\left(\Pi \rho \right)\ge 1-\epsilon,\quad \epsilon \in \left[0,1\right]$$\end{document} . Then,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${| | \rho -\Pi \rho \Pi | | }_{1}\le 2\sqrt{\epsilon }.$$\end{document}

Note that this is a slightly weaker result than the more general version in ref. ^41^, but it is sufficient for our proof. Earlier proofs of essentially the same result can be found in^42,43^. Furthermore, while we will use this version for the sake of simplicity, a slightly tighter bound was recently proven^44^ (see Lemma 6 therein). We use this lemma twice, for the projectors Π_M_ and Π_N_ onto the subspace of states with a number of particles at most M and N, respectively. Now, the restriction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Pi }_{N}\hat{U}{\Pi }_{N}$$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\in {{\mathbb{C}}}^{N\times M}$$\end{document} be a matrix with linearly independent columns. Then, there exists a unitary matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q\in {{\mathbb{C}}}^{N\times N}$$\end{document} and an upper triangular matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\in {{\mathbb{C}}}^{M\times M}$$\end{document} such that:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=Q\left(\begin{array}{c}R\\ 0\end{array}\right).$$\end{document}

We apply this decomposition to a column full rank submatrix of the cut-off operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Pi }_{N}\hat{U}{\Pi }_{N}$$\end{document} and choose the resulting Q as our finite-dimensional unitary approximation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} , 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^45^, 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.

Tools for proving Theorem 2

Theorem 2 seeks to reproduce a finite-dimensional unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}={e}^{i\hat{H}}$$\end{document} 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} is an unbounded Hamiltonian, we cannot expect the usual power series expansion of the exponential, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{it\hat{H}}={\sum }_{n\in {\mathbb{N}}}\frac{{(it\hat{H})}^{n}}{n!}$$\end{document} , to converge to an operator on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} . For sufficiently well-behaved unbounded operators, the exponential may be defined instead via the spectral theorem or the resolvent formalism^20,46^. 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^20^.

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\hat{q},\hat{p})={\Pi }_{N}\hat{H}{\Pi }_{N}\oplus (I-{\Pi }_{N}){\hat{H}}^{{\prime} }(I-{\Pi }_{N})$$\end{document} for some N≥0, where Π_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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\hat{q},\hat{p})$$\end{document} remains in the subspace of states with a number of particles at most N, and reproduces the finite-dimensional unitary evolution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e}^{i{\Pi }_{N}\hat{H}{\Pi }_{N}}$$\end{document} on that subspace.

The proof of Theorem 2 then proceeds by constructing such a block-diagonal polynomial Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\hat{q},\hat{p})={\Pi }_{N}\hat{H}{\Pi }_{N}\oplus (I-{\Pi }_{N}){\hat{H}}^{{\prime} }(I-{\Pi }_{N})$$\end{document} , given the description of a finite-dimensional Hermitian operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} . This construction is based on decomposing the finite-dimensional Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} as a sum of terms pairing only two basis states, say \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert j\right\rangle$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert k\right\rangle$$\end{document} . 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 pairs (x0, y0), …, (xn, yn), the Lagrange interpolation polynomial is defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L(x)=\sum _{j=0}^{n}{y}_{j}{\prod}_{\begin{array}{c}0\le i\le n\\ i\ne j\end{array}}\frac{x-{x}_{i}}{{x}_{j}-{x}_{i}}.$$\end{document}

It satisfies L(xi) = yi, for all 0 ≤ i ≤ n.

In particular, polynomial operators of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\hat{a}}^{{\dagger} }\hat{a})$$\end{document} with P a well-chosen Lagrange interpolation polynomial allow us to reproduce diagonal entries of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} , while operators of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\hat{a}}^{{\dagger} }\hat{a}){\hat{a}}^{k}+h.c.$$\end{document} 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.

Proof of Theorem 3

Theorem 3 deals with the second kind of effective description we consider, namely approximating a physical infinite-dimensional unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} 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, ϵ)-approximation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} of the original unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d={{\mathcal{O}}}\left(\frac{E}{{\epsilon }^{4}}{E}_{\hat{U}}\left(\frac{64E}{{\epsilon }^{2}}\right)\right)$$\end{document} , in time polynomial in d, E, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\epsilon }$$\end{document} . Then, we extract a finite-dimensional Hamiltonian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}$$\end{document} generating the unitary evolution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} . 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}={e}^{i\hat{H}}$$\end{document} , itself approximating the original unitary operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} .

Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}\in {{\mathcal{U}}}({{{\mathcal{H}}}}_{d})$$\end{document} is a unitary operator on a (d + 1)-dimensional Hilbert space, we can write \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}=WD{W}^{{\dagger} }$$\end{document} for some unitary matrix W and a diagonal matrix D, whose diagonal elements are of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({e}^{i{\lambda }_{n}})}_{n=0,\ldots,d}$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda_n \in \left[0,2\pi \right)$$\end{document} . Let Σ be the diagonal matrix whose diagonal elements are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\lambda }_{n})}_{n=0,\ldots,d}$$\end{document} , then Σ is Hermitian and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}=W{e}^{i\Sigma }{W}^{{\dagger} }={e}^{iW\Sigma {W}^{{\dagger} }}$$\end{document} . Since diagonalising D and the matrix multiplications with Σ can both be performed in time polynomial in d, we can extract the generator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{H}=W\Sigma {W}^{{\dagger} }$$\end{document} of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} efficiently, which concludes the proof of Theorem 3. □

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} of the physical unitary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} , guaranteed to exist by Theorem 1, the proof proceeds by synthesizing a short sequence of finite-dimensional gates approximating \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{V}}_{d}$$\end{document} over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} , via the Solovay–Kitaev theorem for qudits:

Theorem 5

(Qudit Solovay–Kitaev theorem^23^). Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{G}}}$$\end{document} be a finite set of unitary gates with determinant one and their inverses, generating a dense subset of S**U(d). There is a constant c such that for any U ∈ S**U(d), there exists a finite sequence S of gates from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{G}}}$$\end{document} of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{O}}}({\log }^{c}(\frac{1}{\epsilon }))$$\end{document} and such that d(U, S) < ϵ, 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} generated by polynomial Hamiltonians approximation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} . Finally, we leverage the physicality of the target unitary channel to lift the approximation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{U}$$\end{document} over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\mathcal{H}}}}_{d}$$\end{document} to an approximation over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{H}}}$$\end{document} . We refer to the Supplementary Information for a detailed proof.

Supplementary information

Supplementary Information Transparent Peer Review file