Fast simulation of bosonic qubits via Gaussian functions in phase space

Bosonic qubits are a promising route to building fault-tolerant quantum computers on a variety of physical platforms. Studying the performance of bosonic qubits under realistic gates and measurements is challenging with existing analytical and numerical tools. We present a novel formalism for simulating classes of states that can be represented as linear combinations of Gaussian functions in phase space. This formalism allows us to analyze and simulate a wide class of non-Gaussian states, transformations and measurements. We demonstrate how useful classes of bosonic qubits -- Gottesman-Kitaev-Preskill (GKP), cat, and Fock states -- can be simulated using this formalism, opening the door to investigating the behaviour of bosonic qubits under Gaussian channels and measurements, non-Gaussian transformations such as those achieved via gate teleportation, and important non-Gaussian measurements such as threshold and photon-number detection. Our formalism enables simulating these situations with levels of accuracy that are not feasible with existing methods. Finally, we use a method informed by our formalism to simulate circuits critical to the study of fault-tolerant quantum computing with bosonic qubits but beyond the reach of existing techniques. Specifically, we examine how finite-energy GKP states transform under realistic qubit phase gates; interface with a CV cluster state; and transform under non-Clifford T gate teleportation using magic states. We implement our simulation method as a part of the open-source Strawberry Fields Python library.

While concrete quantum computing architectures based on bosonic qubits have been proposed, the analysis and simulation of these qubits is challenging because of the infinite-dimensional Hilbert space that they occupy. This impedes the development and implementation of these architectures since determining fault-tolerance thresholds and overheads is limited by our ability to simulate these physical systems in realistic situations. The current most flexible method for simulating bosonic qubits relies on the Fock basis. Simulations in the Fock basis can be cumbersome, especially for CV states with large energy; in particular, high-quality and therefore high-energy cat and GKP states require a high photonnumber cutoff, incurring large memory loads and pro- * eli@xanadu.ai, nicolas@xanadu.ai, ilan@xanadu.ai, antal@xanadu.ai, theodor@xanadu.ai, josh@xanadu.ai, krishna@xanadu.ai, guillaume@xanadu.ai, ish@xanadu.ai cessing times. Moreover, determining how states change under CV channels and measurements is computationally expensive in the Fock basis representation [17] as the energy of a given state can increase significantly under paradigmatic CV transformations such as squeezing and displacements.
Here, we overcome the challenge of studying bosonic qubits by introducing a novel formalism that enables their analysis and simulation. Specifically, we present a mathematical framework for simulating a class of CV states, transformations, and measurements using linear combinations of Gaussian functions in phase space. This framework allows us to simulate the transformation of useful bosonic qubits such as GKP, cat, and Fock states under Gaussian channels and measurements, as well as under a class of valuable non-Gaussian channels effected through gate teleportation. Motivated as a tool to facilitate the design of quantum computing architectures, our framework can model important sources of decoherence (such as optical loss in photonics or dissipation in superconducting cavities) as well as transformations and measurements that are readily-implementable in photonics, such as linear optics, squeezing operations, homodyne and photon-counting detection. We accomplish this by leveraging the most convenient aspects of the Gaussian CV formalism-namely, the ability to regard the transformation of a state as a transformation of means vectors and covariance matrices-while providing the capability to simulate non-Gaussian systems, which is necessary to the construction of a quantum computer [18].
Informed by our formalism, we provide a method for the fast and accurate simulation of bosonic qubits. We find the scaling of the memory and processing time arXiv:2103.05530v1 [quant-ph] 9 Mar 2021 required for our simulator are vastly more favourable than current state-of-the-art Fock basis simulators [19]. We use this to conduct an in-depth numerical study of bosonic qubits in useful physical circuits under inevitable physical imperfections, including loss in optical components; finite-energy effects in resource states; and finite squeezing in the ancillae of measurement-based squeezing operations, the workhorse of inline squeezing [20]. Specifically, we analyze GKP states in three situations that are relevant for quantum computation. First, we examine GKP states passing through a qubit phase gate (introduced in [11] and studied in Sec. II D of [21]), a Clifford gate typically used in the universal gate set for GKP qubits and whose CV implementation we simulate being performed with measurement-based squeezing. Second, we consider the teleportation of a GKP state into a CV cluster state, a scenario present in proposals for measurement-based quantum computation with GKP qubits [1,4,22,23]. Third, we study applying a qubit T gate to GKP states via gate teleportation with finiteenergy GKP magic states and realistic entangling operations, which is an important scenario because the T gate, in conjunction with experimentally-accessible Gaussian operations, is a standard prerequisite for unlocking universal computation with GKP qubits [11]. While gate teleportation is expected to be favourable over other methods [24], simulation of the technique in the presence of realistic gate noise, to the best of our understanding, has not been performed, likely due to numerical challenges. Thus, by enabling the study of situations that were hitherto intractable, our work provides a valuable toolkit for the analysis and simulation of quantum computation based on bosonic qubits.
The structure of this paper is as follows. In Section II we provide background on the Gaussian CV formalism and on bosonic qubits. In Section III we introduce our new formalism for simulating a wide class of states which can be expressed as linear combinations of Gaussian functions in phase space. We provide rules for how these states evolve under Gaussian and a class of non-Gaussian transformations and measurements. With the framework in hand, in Section IV we show how useful classes of bosonic qubits-GKP, cat and Fock statescan be written in our formalism. In Section V, we detail how the formalism can be used for modelling loss, measurement-based squeezing, and useful non-Gaussian GKP qubit operations. Given the formalism, in Section VI, we provide our simulation methods, along with a comparison to simulation techniques in the Fock basis. Our formalism and methods allow us to present results from novel simulations of bosonic qubits in Section VII. We discuss additional areas of application for our simulator, and open research problems in Section VIII.
For a more hands-on introduction to the simulations that are enabled by our formalism, we invite the reader to consult the open-source code available in Strawberry Fields [19,25] and an accompanying set of tutorials available online [26][27][28].

II. BACKGROUND
In this section we provide overviews and pointers to the relevant literature of three different threads that are unified in the later parts of the manuscript. In Section II A we provide a brief survey of the Gaussian formalism for CV quantum systems, including quantum phasespace, the symplectic formalism, and Gaussian states, channels, and measurements. In Section II B we provide an overview of Strawberry Fields, the programming library in which we implement the formalism and methods developed in the rest of the paper. Finally, in Section II C we review the GKP and cat qubit encodings, which leverage the large Hilbert space of a CV mode.

A. Continuous-Variables and the Gaussian Formalism
Multimode continuous-variable systems are best described using the canonical positionq j and momentum p j operators acting on the infinite-dimensional Hilbert space associated with the system. It is often convenient to group these operators, representing N modes, into a vector:ξ T = (q 1 ,p 1 ,q 2 ,p 2 , . . . ,q N ,p N ). (1) The commutation relations these operators satisfy can be expressed succinctly in terms of the symplectic form as [ξ j ,ξ k ] ≡ξ jξk −ξ kξj = i Ω j,k .
As for any quantum mechanical system, a complete description of its state can be obtained by specifying its density matrixρ. For CV systems it is often useful to introduce the characteristic function [29] χ(r;ρ) = tr D (r)ρ , whereD(r) = exp iξ T Ωr is the Weyl or displacement operator and r ∈ R 2N is a real vector in phase-space.

Gaussian States
We recall Gaussian states [29] as the ones whose characteristic function takes the form where we introduced the vector of means µ and covariance matrix Σ of the stateρ with elements where the expectation value of an operatorÂ is defined as A = tr(Âρ). (8) The covariance matrix of a valid quantum state, whether or not it is Gaussian, satisfies the uncertainty relation A different characterization of Gaussian pure states can be obtained by noting that they can be prepared by applying a Gaussian unitaryÛ to the multimode vacuum state |0 ⊗N ; that is, |ψ G =Û |0 ⊗N , where the unitary is generated by a Hamiltonian that is at most quadratic in the quadrature operators [30,31]. The vacuum state is the unique state that is mapped to zero by its respective destruction operator a j |0 j = 0,â j = 1 √ 2 (q j + ip j ), (10) and has vector of means and covariance matrix µ vac = 0 and Σ vac = 2 1.
Finally, we introduce the Wigner function, which is the Fourier transform of the characteristic function: The Wigner function of a Gaussian state is a Gaussian function of the phase-space variables ξ. Later, we introduce a class of states with Wigner functions that can be expressed as a linear combination of Gaussian functions in phase space.

Gaussian Transformations and Measurements
A Gaussian unitary transformationÛ is equivalent to a homogeneous linear phase-space transformation ξ → S T ξ, followed by a phase-space displacement ξ → ξ + d. Here the symplectic map S (satisfying SΩS T = Ω) takes χ(r;ρ) → χ(S T r;ρ), which, for Gaussian states, is equivalent to transforming the covariance matrix as Σ → SΣS T and mean as µ → Sµ [31]. Finally, the displacement transforms the mean as Sµ → Sµ + d. As examples, the single-mode displacement and squeezing operators are defined respectively bŷ but in phase-space are represented as S sq. = e −r 0 0 e r , d T sq. = 0, where in the last line we assumed for simplicity that ζ = ζ * = r is real. A Gaussian channel is a linear completely-positive trace-preserving map from Gaussian states to Gaussian states. Gaussian channels can be described by a pair of real matrices (X, Y ) with Y + i 2 Ω ≥ i 2 XΩX T [32]. The action of the Gaussian channel described on the characteristic function is χ(r;ρ) → χ (r;ρ) = χ(X r;ρ) exp − 1 2 r T Y r .
It follows that the transformation on the covariance matrix and mean is given by In addition to the mapping provided by (X, Y ), a displacement can always be added to a Gaussian channel and the transformation will remain deterministic. One can also consider the characteristic and Wigner functions of an arbitrary operatorÂ. This is the socalled Weyl transform of a general Hilbert space opera-torÂ, obtainable by replacingρ withÂ in Eq. (4) and then taking its Fourier transform as in Eq. (12). The expectation value of the operator can now be written Â = tr(ρÂ) = d 2N ξ W (ξ;ρ)W (ξ;Â). (16) We allow thatÂ is a measurement operator, in which case Â is the probability of the outcome associated with the operator. In the case the measurement operator describes the detection of a Gaussian state parametrized by r M and Σ M , this is referred to as a general-dyne measurement, and the Weyl transform matches the Wigner function for that state [32]. A homodyne measurement is a limiting case of general-dyne measurement corresponding to projection onto an eigenstate of a quadrature operator, which can be treated as an infinitely squeezed state along that quadrature. The probability distribution for the outcome r M of a general-dyne measurement on a Gaussian state, as can be deduced from Eq. (16), is itself a Gaussian distribution since the integration is between two Gaussian functions. Gaussian measurement motivates a transformation beyond Gaussian channels, namely conditional Gaussian dynamics, i.e., an update to a subset of modes of a multimode Gaussian state conditioned on the outcome of a Gaussian measurement that is performed on the remaining modes. Following [32], the covariance and mean of the multimode state can be written as: where A denotes the modes that will remain active and B those that will be measured. If the Weyl transform of the measurement operatorM corresponds to a Gaussian state with mean r M and covariance Σ M , then the partial FIG. 1. A circuit representation of the Choi-Jamiolkowki method for applying the most general Gaussian CP map to a state. Here, |0 denotes the vacuum state; the vertical line linking two modes with ellipses on the ends indicates the twomode squeezing operation (in this case, we take the limit as squeezing goes to infinity); andΦ(·) is a general Gaussian CP map applied to one half of the Choi state so that it can be teleported onto the target stateρ. Since the two-mode squeezing operation is a Gaussian transformation, this circuit can be treated with conditional Gaussian dynamics whenρ is a Gaussian state.
trace tr B (ρ ABM ) corresponds to a Gaussian integral in phase space over the quadrature variables of modes B, yielding a Gaussian state in modes A with the following covariances and means: Conditional Gaussian dynamics also opens the door to effecting the most general Gaussian completely positive maps via their Choi state representation [32]. A CP map on a state is related by the Choi-Jamiolkowski isomorphism to a Choi state, which is found by acting the CP map on one half of a maximally entangled state. For CV systems, the maximally entangled state is the Gaussian two-mode squeezed state in the limit of infinite squeezing; a two-mode squeezed state can be constructed by interfering q-and p-squeezed vacuum states at a 50:50 beam-splitter. Consider the two-mode Choi state (covariance matrix Σ C and mean µ C ) associated to the CP map. The CP map can then be teleported onto an arbitrary target Gaussian stateρ by projectingρ and the untouched half of the Choi state onto the maximally entangled state, which is again just the two-mode squeezed state. We summarize the circuit in Fig. 1 Since the Choi state and the target are both Gaussian, before the projection they can be expressed in the form of Eq. (17), with the off-diagonal matrices being 0 since the states have not yet interacted. Because the measurement is also onto a Gaussian state, the final output state has covariances and means of the form (18). For the exact mapping, one must take care to evaluate the infinite squeezing limit for the two-mode squeezed states in the final expression for the means and covariances. Ifρ is a multimode state on n modes, the same procedure can be applied using n copies of the maximally entangled state on 2n modes [32].
As we show in Section III, we can take inspiration from the Gaussian phase space mathematical framework we have reviewed to introduce a class of states and measurements that can be expressed as a linear combination of Gaussian functions in phase space, along with how such states transform. Importantly, we find in Section IV that common bosonic qubit encodings fall within this formalism. Next, we review the definitions and properties of those bosonic qubits.

B. Continuous-Variable Simulation with Strawberry Fields
Strawberry Fields (SF) is a full-stack Python library for programming, designing, simulating, and optimizing continuous-variable quantum optical circuits [19,33]. The library has a unified frontend that allows to write CV quantum circuits and programs at a high-level. Moreover, it allows users with basic knowledge about CV and quantum photonics to access an application layer that can be used to solve practical problems in graph theory, point processes and chemistry. The frontend also provides functionality for gate decomposition and program compilation and verification. Once programs are verified and compiled they are passed to a software backend or directly to cloud-available hardware [34]. The software or quantum photonic hardware can return a number of useful results to the frontend, including batches of samples, cost functions for further numerical optimization or representations of the quantum state. The frontend provides further functionality for exploration such as sample processing and plotting. The three software backends handle the actual simulation using different internal numerical representations of quantum states, each with their own unique advantages and weaknesses. The gaussian backend simulates Gaussian states undergoing Gaussian and non-Gaussian (threshold and photonnumber-resolving) measurements [35,36]. The fock and tf backends use a Fock basis truncation to represent CV quantum states and operations as high-dimensional tensors [17]. They differ in the tools they rely on for the numerical implementation of the tensor operations: the former uses the NumPy package [37], while the latter employs TensorFlow [38].
We implement the results of this manuscript as a new, fourth backend of SF [25]. This bosonic backend can be regarded as a generalization of the gaussian backend, and integrates directly into the rest of the SF stack. For a series of beginner to advanced tutorials implemented by the authors on using the new backed, see [26][27][28]. The bosonic backend implements much of the formalism and methods we discuss in Sections III through VI, enabling new simulation capabilities while benefiting from the unified high-level frontend functionality of the rest of the library.

C. Bosonic Qubits
Residing in a two-dimensional subspace of the infinitedimensional Hilbert space of a CV mode, the bosonic qubit is robust unit for quantum computation in platforms such as photonics. Several classes of bosonic qubits can moreover correct errors within the CV space, adding an additional level of protection against physical noise. Here, we review the definitions and main properties of two promising encodings: GKP and cat states. Understanding how these states behave in realistic settings is especially valuable as their use becomes more widespread in quantum technologies. As we show in Section IV, we are able to express these qubits as linear combinations of Gaussian functions in phase space, allowing us to model them under realistic conditions such as loss, finite-energy, and noisy gate teleportations.

GKP States
The ideal square-lattice GKP logical states are defined as infinite combs of Dirac delta functions spaced by 2 √ π in the position quadrature: where |· q denotes an eigenstate of the position quadrature and k denotes the logical value. Rectangular and hexagonal lattice encodings are related to the square lattice via symplectic transformations [11], a mapping that we show falls neatly within our formalism. One advantage of the GKP encoding is that Clifford gates and measurements correspond to Gaussian transformations [11], which are experimentally accessible in the photonics context, as we review in Appendix H. Pauli X and Z gates correspond to displacements by √ π along the q and p quadratures, respectively. The Hadamard gate is a rotation by π/2 in phase space. The qubit phase, CX, and CZ gates correspond to a CV quadratic phase, CX, and CZ gates, which are active Gaussian transformations, in the sense of requiring a squeezing component. In practice, (measurement-based) inline squeezers are challenging to implement but are nonetheless feasible and deterministic, as we discuss in Section V B. Pauli X and Z measurements correspond to homodyne measurements along q and p quadratures. A universal gate set can be completed with the qubit T gate; in the GKP encoding, this gate can be implemented through gate teleportation with a magic state, a process we review and align with our formalism in Section V C.
GKP states can correct small displacement errors in phase space, and and can reduce larger displacement errors to qubit-level Pauli errors. A qubit error-correction code concatenated with GKP states can then be used to correct these discrete errors [11]. In Fig. 2 we review the GKP error-correction circuit from [11], noting that various other decompositions of the circuit exist [39,40]. Briefly, two ancillary GKP states are entangled with the data mode to be corrected and measured with homodyne detectors; the outcomes of these measurements determine the displacement that is then applied to the data mode to correct for the error (up to a logical Pauli error). We Error syndrome measurement with normalizable GKP states following the approach from [11]. First, shifts in q are corrected: an encoded data qubit |ψ gkp and an ancilla |+ gkp are sent through a SUM gate, and |ψ gkp is displaced according to the result of a homodyne q measurement on the ancilla, mod √ π . A similar procedure follows for shifts in p.
discuss in Section V C how our formalism can treat this circuit.
To the chagrin of experimentalists, ideal GKP qubits have infinite energy; to the chagrin of theorists, we must consider their finite-energy, normalizable forms. One such form is obtained by replacing each Dirac delta with a Gaussian peak corresponding to a squeezed state of variance ∆ 2 /2, and then applying an overall Gaussian envelope of width 1/∆ 2 so that peaks further from the origin are suppressed [11]: This normalization process is not symmetric in phase space because the peaks are constrained to remain centred at the initial positions of the ideal state [41]. A related normalization process, which has the Fock damping operator E( ) = e − n applied to the ideal state, is symmetric since the number operatorn acts symmetrically in phase space; we denote such states as |k gkp . In [41] the authors provide a thorough review of the connections and mappings between these finite energy forms of GKP states, noting that they can be related to each other by a simple squeezing operation; in [42] the authors explore alternative normalization envelopes to Gaussians, demonstrating sufficient conditions for the normalization process to yield physical states. Yet another option for finite energy GKP states are comb states [43]; these correspond to taking only a finite superposition of evenlyweighted q-squeezed states centred at the location of the peaks in the ideal state.
There is some freedom in the choice of logical basis states. For example, for resource-efficient preparation of Pauli X eigenstates, one can identify them with coherent states |±α , which are approximately orthogonal as α increases in magnitude. Then, the Pauli Z gate is given by a rotation in phase space by π. Pauli Z measurements become photon number parity measurements, since the wavefunction for k = 0 (1) is symmetric (antisymmetric) and therefore only contains even (odd) photon numbers. A cat-qubit Bell state can be prepared by splitting a higherenergy cat state |k √ 2α cat at a 50:50 beam-splitter. Bell state measurements on cat qubits can be performed by interacting two states at a beam-splitter, then measuring photon number patterns at the output. The Pauli X gate, small single qubit rotations, and two-qubit entangling gates can all be applied deterministically via gate teleportation, conditional on the availability of cat Bell states [12]. In Section IV B, we show how cat states can be written as a linear combination of Gaussian functions in phase space, and in Section IV C we show the same for Fock states. This means that our formalism enables simulation of the states, teleportation-based gates and measurements required for quantum computation with cat states.

III. LINEAR COMBINATIONS OF GAUSSIANS IN PHASE SPACE
Having reviewed the relevant background material on CV Gaussian formalism in Section II A, we present a new formalism for simulating a wide class of CV states, transformations and measurements. Specifically, in Section III A, we introduce states with Wigner functions that can be expressed as a linear combination of Gaussian functions in phase space. Next, in Section III B, we provide the framework for describing Gaussian and a class of non-Gaussian measurements on the aforementioned states. Finally, in Section III C, we detail how the states in our formalism transform under deterministic and conditional Gaussian maps, as well as under a class of non-Gaussian transformations.

A. States in the Wigner Representation
In this work, we consider n-mode states whose Wigner function can be written as a linear combination of Gaussian functions in phase space: where M is the set of indices which can in general be multiparametered, ξ ∈ R 2n is the phase space variable for an n-mode CV quantum system, andρ is the corresponding density matrix operator in Hilbert space. Each Gaussian in the linear combination is associated with a weight c m , a 2n-dimensional mean µ m and a 2n × 2n covariance matrix Σ m , all complex-valued in general. The normalized multivariate Gaussian distribution G is defined as Ifρ is a physical density matrix-in particular, if it has unit trace-then its corresponding Wigner function is normalized through At this point, we do not restrict our means, covariances and weights to be real, as imaginary components from these quantities can be required to produce interference fringes and negativity in the Wigner function; however, Hermiticity of the density matrix implies, at least, that the total Wigner function is real. Furthermore, we require the real part of the covariance matrices to be positive-definite so that the distribution remains bounded; however, the covariance matrices need not always respect the uncertainty relation from Eq. (9).
As an example of a wide class of states that fall into this framework, consider a pure state that consists of a superposition of Gaussian pure states (i.e. displaced squeezed vacuum states): Here, γ n and ζ n are the complex displacement and squeezing parameters of the n th term in the superposition; in this contextD(γ) is the displacement operator andŜ(ζ) is the squeezing operator (cf. Eq. (13)). The density matrix for this state is a sum of terms of the form κ m κ * n |γ m , ζ m γ n , ζ n |. Since the Weyl transform is linear, the Wigner function for the state is a linear combination of functions in phase space, each associated with an operator |γ m , ζ m γ n , ζ n |. For m = n, the phase space functions are simply products of |κ n | 2 and the Wigner function for |γ n , ζ n , which is a Gaussian state and hence has a Gaussian Wigner function. In Appendix A, we prove that, for m = n, the phase space function is still Gaussian, albeit with complex weight, means and covariances. Thus, states of the form (26) can be expressed in the form (23).
As we discuss in Section IV, the representation in Eq. (  as well as superpositions of Fock states created by performing photon-number-resolving measurements on some modes of a multimode Gaussian state can be expressed in the form of Eq. (23). In addition to giving us the ability to write down states useful for quantum computing, our formalism is well-suited to subjecting such states to general Gaussian transformations and measurements, as well as certain non-Gaussian transformations via gate teleportation, as we explore in the next few sections.

B. Gaussian and a Class of Non-Gaussian Measurements
Given a state whose Wigner can be written as a linear combination of Gaussians in phase space, we now describe the formalism for Gaussian measurements on such states. A general-dyne Gaussian measurement on n modes is characterized by the 2n × 2n covariance matrix Σ M of the Gaussian state onto which one projects. The outcome of a general-dyne measurement is a point in phase space, r M [32]. Given a state that can be described by Eq. (23), the probability of outcome r M is A special case of general-dyne measurement is homodyne measurement. Without loss of generality, we consider a measurement of the q quadrature, for which and the measurement outcome becomes r M → q M . Since the q-homodyne distribution can be retrieved by integrating out the p quadrature, and since the Wigner function is a linear combination of Gaussians, we have where µ (q) m and Σ (qq) m denote the q-quadrature components of the means and covariances. Importantly, the distribution is yet another linear combination of Gaussians, this time of a variable in n rather than 2n dimensions. Later, in Section VI, we provide a tailored simulation method for sampling outcomes of Gaussian measurements from states in our formalism.
The mathematical method for calculating the probability distribution of a Gaussian measurement can be straightforwardly generalized to a class of non-Gaussian measurements for which the measurement operator can itself be represented in phase space as a linear combination of Gaussian functions. Assuming the Weyl transform of a measurement operatorM associated with outcome M is of the form the probability of obtaining M is simply given by: Note that, for measurement operators, we do not necessarily have the unit trace condition, so j∈J d j = 1 in general. In Section IV C, we show, as a pertinent example, how Fock states can be expressed as a linear combination of Gaussians in phase space, which means photon-number-resolving measurements can be described with this formalism.

C. Gaussian and a Class of Non-Gaussian Transformations in a Gaussian-Inspired Framework
As we saw in Section II A, Gaussian transformations in phase space map Gaussian states to Gaussian states. As we show next, these maps motivate a wider class of Gaussian and non-Gaussian transformations for states with Wigner functions that can be represented as linear combinations of Gaussian functions, as in Eq. (23).
The first class of transformations we consider are deterministic Gaussian completely positive and trace preserving (CPTP) maps, that is, transformations that are not conditioned on any probabilistic measurement outcomes. As we reviewed in Section II A, deterministic Gaussian transformations acting on an n-mode Gaussian state can be parametrized by two 2n × 2n matrices, X and Y , and a length-2n vector d that together transform the covariance matrix and mean of the Wigner function [32]. Since the mapping is linear over phase space variables, it generalizes straightforwardly for a linear combination of Gaussian functions in phase space: When X is a symplectic matrix and Y = 0, this corresponds to a Gaussian unitary transformation.
Deterministic Gaussian transformations modify neither the number nor the weighting of the peaks in the linear combination, which makes them easy to apply to states of the form (23). A wide class of CV operationsdisplacement, squeezing, rotation, and beam-splittersas well as common noise models-loss and Gaussian random displacements-fall under this umbrella. By contrast, simple transformations such as displacements and squeezing can quickly push Fock distributions beyond the energy cutoff, an important limitation of the Fock representation.
The second class of transformations we consider are conditional dynamics: how a measurement of some modes updates the remaining active modes. In this case, our formalism opens the door to a class of non-Gaussian transformations of the Wigner function; if we take a set of target modes of the form (23) and interact them with a set of non-Gaussian ancillae, also of the form (23), and then perform a non-Gaussian measurement on the ancillary modes of the form (30), then the effective transformation on the target modes is non-Gaussian in general. This conclusion applies even if the ancillary modes or the measurement-but not both-are Gaussian. Effecting a non-Gaussian transformation on a state through an interaction with non-Gaussian ancillary modes or through a non-Gaussian measurement has been studied extensively in the context of CV gate teleportation and state preparation [11,21,[44][45][46][47][48][49][50]. However, we show next that we can represent such non-Gaussian transformations in the spirit of conditional Gaussian dynamics reviewed in Section II A.
To understand these Gaussian-inspired but nonetheless non-Gaussian transformations, consider two sets of modes: the set A of active modes, described initially by a Wigner function of the form (23) with weights a , means µ and covariances Σ , for ∈ L; and the set B of modes that will eventually be measured, with corresponding parameters b k , µ k and Σ k , for k ∈ K. If the two sets of modes are entangled via a deterministic Gaussian transformation parametrized by (X, Y , d), the weights, means and covariances become: In turn, we can express the means and covariances as where A and B indicate the active and measured modes, respectively. Consider now a measurementM with outcome M on modes B with a phase space representation of the form from Eq. (30), parametrized by weights d j , means µ j and covariances Σ j , with j ∈ J . Since the partial trace is linear, the corresponding phase space integral is a sum of many partial Gaussian integrals. Thus the covariances and means of modes A update as [32]: Until this point, the update rules have been inspired by the conditional dynamics for Gaussian states. However, our consideration of multiple Gaussians instead of just one necessitates a novel rule: an expansion and reweighting of the peaks in phase space. While outcome M occurs with probability p(M ;ρ AB ), as in Eq. (31), each peak from modes B and from the measurement operator contribute differently to this probability, with a weight given by: Therefore, given the result M , the weights update as: where we include a normalization by all the new weights. As a result of this process, the total number of Gaussians in the Wigner function for modes A, as well as their associated weights, means and covariances, is been grown both by modes B and by the measurement M : (38) such that the final number of Gaussian functions required to describe mode A is the product of the initial number of functions in modes A and B, and in the Weyl transform ofM .
While we have only been tracking transformations of Gaussian functions in phase space, it is worth emphasizing that these conditional dynamics increase the number of Gaussian functions initially in modes A, thereby necessarily describing non-Gaussian transformations. The cost of modelling non-Gaussian transformations with our formalism is that the number of peaks one needs to consider can grow; still, we show that this is a worthwhile trade-off compared to alternative methods for simulating certain classes of bosonic qubits.
In the case that modes B and the measurement are truly Gaussian states, as opposed to linear combinations of Gaussian functions, then b k = d j = 1, and (L, K, J ) → L. This means the initial number of weights does not increase, and the initial means µ and covariances Σ follow the traditional Gaussian conditional update rule reviewed in Section II A: where Σ M and r M parametrize the Gaussian measurement as in Eq. (27). The peak reweighting from Eq. (37) is still required; however, it is simpler because the total number of weights does not increase and d j = 1. Finally, we point out that since we can implement conditional Gaussian dynamics in this framework, we can apply the most general Gaussian CP maps via the Choi-Jamiolkowski isomorphism using the circuit in Fig. 1.
In Section V, we explore how various transformations such as loss, Fock damping, measurement-based squeezing, and gate teleportation onto bosonic qubits can be described in this formalism. Before this, in the following section, we present the description common bosonic states as linear combinations of Gaussians in phase space.

IV. USEFUL FAMILIES OF STATES EXPRESSED AS LINEAR COMBINATIONS OF GAUSSIANS
Having provided a general formalism in Section III, we now demonstrate how bosonic qubits, namely GKP, cat and Fock states, can be written in the form of Eq. 23.

A. GKP States
Care must be taken when dealing with the finite-energy forms of GKP states. Before we tackle those, let us recall the Wigner representation of the ideal states.

Ideal GKP States
Any pure GKP state can be expressed in the Bloch sphere representation as The Wigner function of an ideal single-mode squarelattice GKP state can be expressed in terms of the pa-rameters of Eq. (23) in the following way [11]: The weights c m (θ, φ)-which encode the logical content of the state-are presented in Appendix B as given in Ref. [51]. Recall each Gaussian peak in the Wigner function is a function of Σ −1 m , so we cannot directly evaluate the limit as δ → 0+. The Gaussian distributions for the ideal GKP Wigner function are Dirac delta functions located at the lattice points enumerated in M. This means that the ideal GKP state has infinite energy and cannot be normalized, rendering it unphysical. However, the covariances matrices do not vary with the index and are proportional to the identity, a feature that is put to use in the following section. Note that GKP states corresponding to alternative lattice spacings, such as rectangular or hexagonal GKP states, can be related to square-lattice states by symplectic transformations in phase space [11]. As we showed in Section III C, symplectic transformations are straightforward to apply to states that can be expressed as linear combinations of Gaussians in phase space, so it is easy to transform between lattices in our framework. Additionally, GKP qudits simply correspond to a different linear combination of δ-functions in phase space, so our formalism can be employed to treat those states as well.

Finite-Energy GKP States
There are different ways to obtain finite-energy versions of the GKP states. A summary of these is presented in Fig. 1 of Ref. [21]. Here, we focus on a commonlyused model of finite-energy GKP states: a Fock damping operatorÊ ( ) = e − n , > 0, (42) applied to the ideal GKP state in Eq. (19). Becausê E( ) is a single non-unitary Kraus operator, it is nontrace-preserving, meaning one needs to carefully account for the normalization of the resulting states. A straightforward but somewhat lengthy calculation (expounded in Appendix C 1) shows that the Wigner function corresponding to theÊ( ) operator applied to an ideal state can be represented in the same form as Eq. (23). With the set M, ideal coefficients c m (θ, φ), and ideal means µ m given in Eq. (41), we have that: Here N is chosen such that m∈M c m ( ; θ, φ) = 1. In this representation, the weights and means are real. In practice, one can take a finite number of terms from Ma numerical cutoff-depending on the desired precision. The value of the cutoff would depend on the strengths of the weights and the normalization of the Gaussian distributions in the Wigner expansion. A simpler derivation motivated by the dilation of the operatorÊ( ) is provided in Appendix C 2. We compare the model for finite-energy GKP states we have derived in Eq. 43 to the approach from [22]. There, the Dirac delta spikes of the ideal GKP Wigner function are replaced with Gaussians of non-zero variance. While that approach can capture the broadening of peaks due to finite-energy and other noise effects, and while the covariance matrix of each peak can be updated as Gaussian channels are applied, the state does not have an envelope, so it still has infinite energy and remains unphysical. Moreover, in that model, the locations of peaks are not tracked, so their contraction towards the origin due to finite energy effects and any shifts under Gaussian channels are ignored. This can have an impact on the estimation of qubit-level errors. By contrast, the mathematical form for finite-energy GKP states provided here offers a rich noise model that captures these effects.
For an alternative phase-space representation of finiteenergy GKP states, one can applyÊ( ) to the wavefunction to obtain a superposition of squeezed states, and then use this form to calculate the Wigner function directly. This representation differs from the one previously presented in two ways. First, the weights and means of the resulting Wigner function are complex rather than real numbers. Second, the covariances of the individual Gaussian functions now respect the uncertainty relation. Since in this representation, each peak is directly mappable to squeezed states in the wave function, it can be more useful for converting effects observed at the wavefunction level to phase space transformations. The derivation of this alternative representation is provided in Appendix C 3.
Let us denote a logical GKP qubit by |ψ(a) = a 0 |0 gkp + a 1 |1 gkp . We find that the parameters from Eq. (23) become Here N (a, ) is an overall normalization constant such that m∈M c m = 1; it depends on the choice of the GKP state and the strength of the Fock-damping operator ap-plied to it. This time, the Gaussians all have the same covariance matrix, and this matrix respects the uncertainty relation [30]. While the coefficients can be complex, this only stems from the initial qubit coefficients a 0 , a 1 . Finally, we can see that the means µ m ( ) are complex; the resulting sinusoidal oscillations in phase space are what generate the interference fringes and negative regions of the Wigner function in this representation.

Squeezed Comb States
Squeezed comb states form another class of finiteenergy GKP states [43]. These states can be obtained from the ideal GKP states by applying a symmetric step function in the position basis to pick out a certain number of peaks and then replacing the infinitely squeezed position kets with their finitely squeezed analogues.
Mathematically, squeezed comb states can be written in the following way [43]: where andq n = −(N + 1) d 2 + nd. The normalizaton constant is given by We can then write Using the linearity of the mapping between density operators and Wigner functions in Eq. (12) and the results from Appendix A, we find the Wigner function has the canonical form (23) with For a representation of comb states in terms of realvalued Gaussian functions in phase space, see Appendix F.

B. Cat States
As is the case with squeezed comb states, cat states are linear superpositions of pure Gaussian states. We can write the density matrix of these states as Again by linearity of the mapping between density operators and Wigner functions, we see that the Wigner functions of the first two terms | ± α ±α| are Gaussian functions with the covariance matrix of the vacuum and displacement µ ± = ± √ 2 ( (α), (α)). Using the results from Appendix A we find the Wigner function of |α −α| and |−α α| to be complex-valued Gaussians with prefactor e −2|α| 2 , vacuum covariance matrix, and complexvalued vectors of means and where all the Gaussian functions have the same vacuum covariance matrix given in (11). Cat states can also be approximated, to arbitrary precision, by a linear combination of Gaussian functions where the weights, means, and covariance matrices are all real. The trade-off for this real-valued representation is an increase in the number of terms in the linear combination. The derivation, shown in Appendix E, leads to where α ∈ R =0 , c ± , µ ± and Σ ± are the same as in Eq. (55), and D is a positive real parameter controlling the precision of this approximation through O(e −2πD 2 ). The general case of a cat state with α = |α|e iθ can be understood as a θ-rotation in phase space of a cat state with parameter |α|; as rotations correspond to symplectic transformation, we know from Section III C that this is easily implementable. Alternatively, a generalization of the description presented in Eq. (56) can be found in Appendix E.

C. Fock States
In this section we consider the representation of Fock states as linear combinations of Gaussians. Compared with GKP and cat states, it is less clear whether Fock states can be can be written this way, since they cannot be represented as discrete coherent superpositions of Gaussian states. To obtain a representation for Fock states we use the idea of photon addition; that is, we study a quantum-optical circuit in which postselected heralded outcomes in certain ancillary modes allow us to apply the creation operation of a given mode to an arbitrary input state.
We first consider photon addition applied to the vacuum state for the generation of a single photon, as considered in Ref. [52][53][54][55]. The circuit for the addition of a single photon is shown in Fig. 4 (a), where the vertical line with ellipses on the ends corresponds to a two-mode squeezing operationŜ 0,1 (r) = exp r â † 0â † 1 −â 0â1 with squeezing parameter r 1. To see how it works, we simply need to calculate whereÎ j is the identity operation in the Hilbert space of mode j. If |ψ = |0 0 is the vacuum state of mode 0, then the state at the output is a single photon in this mode. Since we are working in the regime where r 1, to make progress we replace the photon-number-resolving detection |1 1| by its poor-man's version,Î − |0 0|, the threshold detection.
We can now write the probability of successful and failed heralding more formally as where |Ψ =Ŝ 0,1 (r) |0 0 0 1 andn = sinh 2 r is the mean photon number in either of the two modes, 0 or 1. The state conditioned on successful heralding with threshold detectors in mode 0 is which is a linear combination of density matrices for two Gaussian states, namely, a thermal state with mean photon numbern and the single-mode vacuum. A thermal state with mean photon number â †â = n =n is a mixed Gaussian state with zero mean and covariance matrix Σ = (n+ 1 2 )1 and can be expressed in the Fock basis as [29] ρ th = 1 1 +n wheren = [e − 1] −1 or, equivalently, e = 1 + 1/n. With this expression one can easily confirm thatρ 1 in Eq. (61) approaches a single photon in the limit r → 0. We generalize this scheme to an n-photon addition that, as schematically shown in Fig. 4 (b), gives a Fock state with n photons when applied to the vacuum. The mathematical details of the derivation are provided in Appendix D; the final result is that we can approximate any n-particle Fock state using the general notation of Eq. (23) with Although the equations above were derived with the assumption that the squeezing parameter r 1, and they depend explicitly on this parameter, as long as r < 1/ √ n, the expressions correspond to a physical state. Formally, it is only in the limit r → 0 that one recovers with perfect fidelity a Fock state. In Fig. 5 we study the behaviour of the (in)fidelity between an ideal Fock state and our approximation. For r ∼ 10 −2 one gets a fidelity of at least 99.9% for n = 1, 2, 3. It is worth adding that a path to approximating superpositions of Fock states as linear combinations of Gaussian states to high fidelity is accessible through our formalism. We recall Gaussian Boson Sampling-type state preparation devices, wherein Fock measurements on all but one mode of a multimode Gaussian state herald a superposition of Fock states in the remaining mode. The exact form of the superposition can be tuned using the means and covariance matrix of the multimode Gaussian, as well as the choice of which Fock measurements to postselect upon [21,49,50,56]. As the input state is Gaussian, and we can approximate Fock measurements to arbitrarily high accuracy as linear combinations of Gaussians, the output state can be calculated using the partial trace form of Eq. (16); this materializes as a linear combination of Gaussian integrals over all but the two phase space quadratures of the output mode, yielding a linear combination of Gaussian functions in phase space for that mode. Notably, this state preparation device can be used as a means to prepare bosonic qubits, including GKP and cat states.

V. USEFUL TRANSFORMATIONS AND MEASUREMENTS IN A GAUSSIAN-INSPIRED FRAMEWORK
Given our formalism from Section III, and examples from Section IV for writing bosonic qubits as linear combinations of Gaussians, we now examine how our formal-ism can be used to treat useful transformations for studying quantum computing with bosonic qubits, namely loss, measurement-based squeezing and non-Gaussian gate teleportation.

A. Loss Channel and Fock Damping
Loss is a physical process modeled as an interaction with a thermal environment through a beam-splitter transformation, resulting in a bosonic Gaussian channel. The loss channel is closely related to the additive random noise channel by composition with an amplifier channel [57][58][59][60]. The strength of the loss parameter of the channel is set by the beam-splitter angle. The pure loss channel is easily described in the form (32) with the matrix pair where η is sometimes referred to as the transmission or transmissivity parameter, and where the environment is assumed to have zero temperature. Thermal lossinteraction with a thermal environment with covariance matrix (n + 1 2 )1-can also be incorporated into this description by multiplying Y η by (2n + 1). The Kraus operators and other details of the pure loss channel are provided in great detail in [60]. Loss is the dominant imperfection in the photonics context.
The Fock damping or finite-energy operator is given by (42). It is valuable for converting infinite-energy states, such as ideal position and momentum eigenstates and GKP states, into their normalizable, finite-energy forms. The operator is directly linked to the loss channel as it forms a leading-order Kraus operator of the channel (see Eq. 4.6 of [60]), while maintaining the purity of the state. In [61], in the context of GKP states, it was shown that E( ) can be derived by passing a state through a beamsplitter of transmissivity cos θ = e − with an ancillary vacuum state, and postselecting the ancillary mode on vacuum (see also Appendix B1 of [21] for a simple derivation of this result).
As we show in Appendix G, the use of only Gaussian states and operations in derivingÊ( ) makes the operator fit neatly into the formalism for Gaussian transformations developed in Section III C. If the initial means and covariances of the mode are (µ m,0 , Σ m,0 ), then we can simply use Eq. (34) with: where S θ = cos θ1 sin θ1 − sin θ1 cos θ1 is the symplectic matrix for a beam-splitter assuming a mode-wise ordering (q 1 , p 1 , q 2 , p 2 ). From there one can proceed with the update rule in Eq. (39), as well as the per-peak reweighting in Eq. (37) with d j = 1, r M = 0, and Σ M = 1/2, since the projection is onto vacuum.

B. Squeezed Ancilla-Assisted Gates
While passive Gaussian transformations can be effected in optics using phase shifters and beam-splitters, and while squeezing applied to vacuum can be achieved with nonlinear cavity resonators [62] or waveguides [63], inline squeezing-squeezing applied directly to an arbitrary input state-poses a greater challenge. A feasible approach to inline squeezing is possible by consuming an ancillary squeezed vacuum state [20]. To effect squeezing in q (p), one first interferes the target state with a highly q-squeezed (p-squeezed) vacuum on a beamsplitter parametrized by angle θ. Next, one performs a phomodyne (q-homodyne) measurement on the ancillary mode with a detector of efficiency η, producing result p M . Finally, one applies a feedforward p-displacement (q-displacement) by η −1 tan θp M to the target mode. As we discuss next, the resulting average map effects a transformation equivalent to squeezing the target mode by cos θ in q, along with some noise dependent on the level of squeezing of the ancillary state and the efficiency of the homodyne detection. Squeezing along any quadrature can be achieved by preceding and following inline squeezing with rotations.
Since the ancillary state, the transformation between the target and ancilla, the measurement on the ancilla, and the feedforward are all Gaussian, inline squeezing using a squeezed ancillary state-commonly referred to as measurement-based squeezing-falls within the formalism developed in Section III C. In Appendix H, we derive the resulting transformation on the covariance matrices, means and modes of the linear combination of Gaussians in phase space that form the target state. On average (integrating over p M ), the map due to this type of squeezing is a deterministic Gaussian CPTP map parametrized by: where r is the squeezing parameter for the ancilla. The superscript q denotes squeezing in the q quadrature. This result almost matches the result from [20]; however, we avoid any extra displacement on the target state by using knowledge of the homodyne detection efficiency in the feedforward displacement.
Importantly, while measurement-based squeezing produces noisy squeezed states on average, in the single-shot case (not averaging over p M ), the covariances, means and weights of the Gaussian functions of the target state transform in a non-linear fashion that precludes description with a deterministic CPTP map. We show this fact in Appendix H. For initial means and covariances (µ m,0 , Σ m,0 ), we can set Eq. (34) to be: where S r,2 models squeezing of the ancillary vacuum, S θ is the beam-splitter symplectic, and (X η,2 , Y η,2 ) model inefficiency in the homodyne detector for the second mode. Next, we use the update rule in Eq. (39), as well as the per-peak reweighting in Eq. (37) with d j = 1, modelling p-homodyne measurements. Finally, to that result we apply a displacement in p-quadrature to the means by η −1 tan θp M . Inline squeezing is an especially valuable operation to be able to model, as it is required to perform the most general Gaussian unitary operations [64]. Moreover, for GKP states specifically, inline squeezing is required to perform certain Clifford qubit gates such as the phase gate and the controlled-NOT and controlled-Z gates, as we review in Fig. 18 of Appendix H. Using our formalism, we now can model realistic noise effects from inline squeezing applied to bosonic qubit states in both the average and single-shot case. This is useful, for example, when quantifying the build-up of noise in stitching together cluster states of bosonic qubits using active entangling operations [1,22], although there exist methods for tailored clusters that do-away with inline squeezing [65].

C. GKP T gate and Error Correction
Just as inline squeezing operations can be teleported onto a mode via the consumption of a squeezed vacuum ancilla, non-Gaussian gates can be effected via gate teleportation using non-Gaussian ancillary states or measurements [11,21,[44][45][46][47][48][49][50]. As a pertinent example, we consider the implementation of the T gate for GKP states. A single-qubit T gate is a non-Clifford gate that, in conjunction with Clifford gates, completes the set of universal operations. Physically, in the GKP encoding, it is implemented through a non-Gaussian transformation, the only one strictly required in the universal set. To effect a T gate via gate teleportation, an ancillary GKP magic state |M gkp = 1 √ 2 (e −iπ/8 |0 gkp + e iπ/8 |1 gkp ) is first rotated in phase space by π/2 and entangled with the target mode at a CZ gate. Next, a p-homodyne measurement is performed on the ancilla; the measurement outcome is rounded to the nearest n √ π , and the parity 6. Optical implementation of the GKP qubit T gate up to global phase, following the method from [11]. Here, in the ideal limit, |M gkp = 1 √ 2 (e −iπ/8 |0 gkp + e iπ/8 |1 gkp ). The GKP Hadamard gate (with operatorĤ) is a π/2 phase shift, and the GKP CZ gate is discussed in Appendix H. The output p0 of the p-homodyne measurement is processed by the function s(p0), which rounds the value to the nearest n √ π and returns the parity of n. If n is even (odd), no phase gate (a phase gateP (s) = e is q 2 /2 with s = 1) is applied.
of n determines whether a feedforward phase gate is applied on the target mode [11]. The non-Gaussianity of the gate enters via the ancillary GKP resource state, which has significant Wigner negativity [51]. See Fig. 6 for a summary, and Fig. 18 of Appendix H for decompositions of the CZ and phase gate.
We emphasize that the ancillary GKP magic state and the target GKP state can both be described as Wigner functions of the form (23). The CZ and phase gatesimplemented directly or via ancilla-assisted squeezingcan be described through an average map or as a singleshot transformation, as we explored generally in Section III C and specifically to these transformations in Appendix H. Finally, a p-homodyne measurement of the ancillary mode containing the magic state falls neatly into the formalism from Sections III B and III C. Unlike measurement-based squeezing, we cannot write the average transformation effected by this circuit as a Gaussian CPTP map on the original state, but we can simulate the single-shot transformation of the gate teleportation by selecting or sampling an outcome for the last homodyne measurement and using Eqs. (35) and (37).
While we can now implement single-shot non-Gaussian operations via gate teleportation, a simple average map for T gate teleportation is untenable for several reasons. First, the ancillary state is itself described by a linear combination of Gaussian functions in phase space, in contrast to the single Gaussian peak of the ancillary squeezed state used for inline squeezing. In the latter case, when the ancillary mode is added to the description, the weights {c m } in the Wigner function of Eq. (23) do not change, since they are each multiplied by 1, while new rows and columns for the additional mode are added to the covariances and means of the Gaussian functions. When the ancillary magic state is added, however, in addition to new rows and columns being added to the covariances and means, the weights and the set of indices M of the Wigner function also change, since we must consider all the cross-terms of different Gaussian peaks in the target and ancillary state. This means that the final resulting state after the teleportation operation has a different set of weights than what it started with, which cannot be captured by a simple average map on each Gaussian function. Second, when applying the feedforward phase gate, the homodyne data is processed by a highly non-linear function because of the binning and parity check, so there would again be no opportunity for the integration over homodyne outcomes to yield a simple Gaussian CPTP map as in Eq. (32).
Notice that the GKP error-correction circuit in Fig. 2 is just another example of a gate teleportation circuit that can be treated by our formalism. The initial ancillary modes are GKP states, and hence can be expressed as a linear combination of Gaussian functions; the entangling operations are Gaussian CX gates; and the homodyne measurements and feedforward displacements are also all Gaussian. Thus, single-shot error correction, with appropriate noise sources such as loss and Fock damping incorporated into the circuit, can be studied using our formalism.

VI. SIMULATION METHODS
Having provided our general formalism in Section III, and useful states and transformations in Sections IV and V, we now describe our simulation methods. In Section VI A, we summarize how to track the weights, means and covariance matrices associated with the linear combination of Gaussian functions in phase space that describes the state of the modes in a photonic circuit. In Section VI B, we detail a technique for sampling outcomes of Gaussian measurements on states in our formalism. Finally in Section VI C, we compare our simulation method to a state-of-the-art method the employs the Fock basis, demonstrating the advantage of our formalism and method as compared to this alternative leading technique.

A. Tracking Weights, Means and Covariances
Our method for simulating states and transformations from our formalism is relatively straightforward. It can be summarized as follows: 1. Initialize the weights, vectors of means, and covariance matrices for each mode at the start of the circuit. For example, Eq. (55) provides these parameters for the cat state.
2. Combine the initial weights, means and covariances for each mode into weights, means and covariances for the multimode state of the whole circuit as follows: This procedure can be applied recursively to build the weights, means and covariances for the full multimode initial state.
3. For each gate or measurement in the circuit, apply the associated transformation to the weights, means and covariances of the state. For example, to apply the loss channel, the matrices that parametrize the channel from Eq. (65) are used in Eq. (32). Alternatively, for a given measurement outcome on a subset of modes, the other modes are updated using Eqs. (35) and (37). 4. The output of the simulation consists of the weights, means and covariances of all the modes after the application of transformations and measurements in the circuit, along with any of the measurement outcomes collected along the way. Since these weights, means and covariances can be used to construct the Wigner function of the state as in Eq. (23), we have access to full state information.
For some states, such as GKP and cat states, all the Gaussian peaks in the linear combination representing that mode have the same covariance matrix. In those cases, to save memory, we need not initialize the same copy many times, and can simply track a single covariance matrix.
In the next part, we detail our method for sampling outcomes of Gaussian measurements for states in our formalism.

B. Sampling Outcomes of a Gaussian Measurement
To further the utility of our formalism and method for performing simulations, we now provide an algorithm for sampling outcomes of Gaussian measurements on states that can be expressed as linear combinations of Gaussian functions in phase space. In Section III B, we detailed in Eq. (27) how to calculate probability distributions for Gaussian measurements within our formalism. A general approach to sampling from non-trivial Wigner functions can be found in Appendix D of [53]; however, we find that a more tailored approach one can use is a rejectionsampling technique [66]. In [67], the authors consider a rejection-sampling method for a distribution composed of real-valued Gaussian functions with (positive or negative) real weights. Here, we extend the method to include Gaussian functions with complex weights and means. We summarize the technique in Algorithm 1. First, we separate the peaks into those which have negative weights and real-valued means M − = {m|c m < 0} ∩ {m| (µ m ) = 0}, and everything else {m ∈ M − }: By writing µ m = (µ m ) + i (µ m ) we can obtain This allows us to define an upper-bounding function to the distribution of interest: 1 Next, we note that g(r M ) is a scalar multiple of a probability distribution q (r M ): Importantly, one can sample from q (r M ) more straightforwardly: first, one samples a value m 0 ∈ M − according to the distribution p m =c m /N ; since the weights c m from g(r M ) are now all positive, we can associate them with probabilities after proper normalization by N . Next, one samples a phase space value r 0 according to G (µm 0 ),Σ M +Σm 0 (r). Given r 0 , one then samples y 0 uniformly form [0, g(r 0 )]. If y 0 ≤ p (r 0 ;ρ, Σ M ), then r 0 is kept as the sampled value of p (r M ;ρ, Σ M ); otherwise, it is rejected and the process is restarted until a sample is drawn. This sampling technique leverages the form of the Wigner function of our states: we use the fact that our distribution can be bounded above by a mixture of Gaussian functions with all-positive weights, and that it is computationally easy to sample a value from a single Gaussian distribution.
Algorithm 1 Simulating Gaussian measurement outcomes using rejection sampling Input: Weights cm, means µm, and covariances Σm of the initial state. Covariance ΣM of the measurement. Output: Sampled phase-space position outcome r0 With the simulation methods in hand, we now compare our simulation method to simulation methods in the Fock basis.

C. Comparison to Fock Basis Simulations
It is well-known that representing a Gaussian state in terms of its vector of means and its covariance matrix is a more efficient representation than writing the state in the Fock basis, especially considering the energy of the state increases under displacements and squeezing. Here, we draw a similar conclusion for the states we have studied under our formalism. In this and upcoming sections, we implement our simulation technique in the new bosonic backend of Strawberry Fields.
Consider that the Hamiltonian for a single mode is given byĤ = 1 2 (q 2 +p 2 ) = 1 2r 2 = (n + 1 2 ). This means that the Wigner function of a Fock state |n has an associated radius in phase space of roughly |r n | = (2n + 1), beyond which the function decays monotonically. This fact can also be demonstrated rigorously using an analytic representation of Fock Wigner functions and the properties of Laguerre polynomials [68,69]. Thus, to determine the Fock representation for a state which has a phase-space Gaussian peak at a point r 0 , a conservative estimate has that we would need a photon number of at least to reach the required radius in phase space. Furthermore, Fock states beyond this value may be necessary, for example, to shape the phase space peak for a desired level of squeezing. In the presence of realistic noise sources like loss, the state additionally becomes mixed, requiring a density matrix rather than a state vector representation in the Fock basis, squaring the number of elements which one needs to track. From these considerations, we see that for a state with a phase space peak at r 0 , the number of elements one needs to track in the Fock basis to have a faithful representation of the state scales like Let us compare this Fock-basis scaling to what we would obtain by expressing our states of interest as linear combinations of Gaussian functions in phase space. The trivial case is a Gaussian state; regardless of its position, orientation, and level of squeezing in phase space, one need only track a 2-component vector of means and a 2 × 2 covariance matrix. For an N mode Gaussian state, the number of elements to track scales like 4N 2 +2N -no exponential scaling since all the modes are represented by a single Gaussian function. Adding a mode with a Gaussian state to a series of other modes with states represented by a linear combination of Gaussian functions only increases the dimension of the covariance matrices and means by 2, but does not change the number of weights, means or covariances one needs to track. This is especially useful for the GKP encoding which, when combined with Gaussian states, enables universal quantum computation [70,71].
We have shown in Section IV B that single-mode twolobe cat states can be represented using one 2 × 2 covariance matrix, along with four weights and 2-component vectors of means-two real-valued and two complexvalued. Crucially, the size of the mathematical objects required for this representation is independent of the energy of the cat state associated with α, and is invariant under Gaussian transformations. While the covariance for N modes encoded as cat qubits scales like 4N 2 , the number of weights still scales exponentially, like 4 N , and the number of elements in the means grows like 2N (4 N ). This is still preferable to the Fock representation, for which the number of density matrix elements scales like |α| 4N ; notably, the scaling for the linear combination of Gaussians representation does not depend on the energy We entangle the states via a CV CZ gate and and subject them 1% loss in each mode; next, we measure the mode originally containing the GKP state in the p-basis, postselecting on p = 0. In the top panel, we perform the simulation using the fock backend in Strawberry Fields with a cutoff of 50 photons. In the bottom row, we employ the bosonic backend. We see that the bosonic backend can maintain a correct description of the state even as the number of photons becomes very high due to the large per-peak squeezing.
of the state or modifications under Gaussian transformations. This can be valuable since, for example, the error rates for some qubit gates on cat states can scale like 1/|α| [2], so to examine regimes of low error, one must increase the energy of the cat state. In Fig. 7, we examine the case of two cat states sent through a lossy beam-splitter. We trace out one of the modes and plot the Wigner function of the remaining mode for increasing values of energy, as parametrized by α. We compare the results using the fock backend of Strawberry Fields using a photon number cutoff of 50 photons per mode-going beyond this value saturates the memory of the standard desktop terminal used for simulation-and using the bosonic backend, representing states as linear combinations of Gaussian functions in phase space. We find that, for α = 2, the Fock representation still works well, albeit requiring more memory and running more slowly when simulating the lossy beam-splitter. For α = 4 and 6, the Fock representation with 50 photons quickly becomes insufficient. This is unsurprising: for these values of α, the action of the beam-splitter leads to Gaussian peaks in phase space at distances, respectively, of 8 √ and 12 √ from the origin, requiring, by the conservative estimate in Eq. (74), greater than 32 and 72 photons for each case. In fact, although n = 32 is comfortably within the cutoff of 50 photons, the Fock representation cannot accurately construct the Wigner function for α = 4, confirming that Eq. (74) is an underestimation of the required photon number cutoff. As seen in Section IV C, even Fock states of n photons can be approximated by n real-valued weights, n 2×2 covariance matrices, and one 2-component vector of means. This is perhaps more memory-intensive than representing a pure number state in the Fock basis; however, under Gaussian transformations such as displacements or squeezing, or the common loss channel, representing the state by a linear combination of Gaussians becomes advantageous.
Finally, while GKP states have an infinite number of peaks, we can consider only a finite number since the weights in Eq. (43) decay exponentially under the Fock damping operator. Since the Fock damping operator decays exponentially in Fock space, an analogous procedure can be executed to establish a reasonable photon number cutoff. If a single-mode GKP state is approximated by a finite number of peaks, with the furthest peak located near √ π 2 (k, ), then the number of peaks one needs to track scales like k 2 + 2 , and consequently so does the number of elements for the vectors of means; by contrast, only a single 2×2 covariance matrix is required. Compare this with the Fock representation, where, by Eq. (74), the number of density matrix elements scales like (k 2 + 2 ) 2 . For N modes encoded as GKP states, we track (k 2 + 2 ) N weights and means, and a single 2N × 2N covariance matrix, an improvement over the (k 2 + 2 ) 2N scaling for the number of density matrix elements in the Fock basis. As a demonstration of the advantage of the bosonic backend, in Fig. 8 we plot the Wigner function for the output mode of a simulated CV teleportation of a high-energy GKP state ( = 0.01, corresponding to 20 dB of per-peak squeezing) onto a p-squeezed state with 15 dB of squeezing using a lossy CZ gate, homodyne measurement, and feedforward displacement. We see that the fock backend is limited in the radius of phase space that it can accurately capture, while the bosonic backend produces all peaks correctly. Moreover, the all-Gaussian nature of the teleportation circuit allows for more efficient computation when representing the states as linear combinations of Gaussians over the Fock representation, which requires many tensor contractions for large density matrices.

VII. NUMERICAL SIMULATIONS
Having established our simulation method and its advantage over competing techniques in Section VI, we now put it to use in simulations of bosonic qubits, leveraging a numerical implementation of our method available through the bosonic backend of Strawberry Fields [25]. In Section VII A, we provide a basic test-run of our simulator, plotting Wigner functions of and sampling from bosonic qubits. In Section VII B, we provide novel simulations of bosonic qubits undergoing realistic gate implementations, with a focus on GKP states. For an advanced tutorial implemented by the authors on using the bosonic backend of Strawberry Fields to simulate realistic GKP qubit gates, see [28].

A. Basic Examples: Wigner Plots and Homodyne Sampling
As a first simple demonstration of simulations with our technique, in Fig. 9 we plot the Wigner functions for |0 gkp with = 0.1 (10 dB of squeezing per peak of the wavefunction), for |0 α cat with α = 2, and for the single-photon Fock state. For the GKP state, we can identify the following features: the Gaussian functions centred near integer and half-integer multiples of √ π; the positive and negative peaks determined by the weighting function; and the per-peak variance, which is smaller than that of the vacuum. For the cat state, we see a similar distribution to Fig. 3, where the interference fringes are produced by the Gaussians with complex weights and means. Finally, for the Fock state, we see how two Gaussians, both with zero means but with slightly different covariance matrices, can combine to form a rotationally symmetric Wigner function with a region of negativity in the centre.
Algorithm 1 offers a straightforward method for sampling the results of general-dyne measurements of states with Wigner functions that are expressed as linear combinations of Gaussian functions. We use the algorithm to simulate 2000 samples of q and p quadrature homodyne measurements of a |0 α cat with α = 2, and of a |0 gkp with = 0.1 (10 dB per peak), and plot the output in Fig. 10. We see that the results align with the asymptotic marginal distributions for these quadratures. The marginals can also be easily attained in our formalism: integrating out one quadrature for a linear combination of Gaussian functions amounts to dropping the entries associated with that quadrature from each mean and covariance matrix, as we saw in Eq. (29).
For the cat state, we see the q quadrature distribution is the same as for mixture of coherent states centred at ± √ 2 α, since the interference fringes cancel out, while they are prominent in the p quadrature. For the GKP state, we can clearly identify the peaks at even (all) integer multiples of √ π in q (p) quadrature.

Measurement-Based Squeezing
In Section V B and in Appendix H, we examine how to use our formalism to describe measurement-based squeezing, a feasible method for applying in-line squeezing operations. Later subsections examine gates that employ measurement-based squeezing in more complicated circuits with bosonic qubits, but we restrict ourselves to a simpler case to develop a clear a picture of the action of the transformation. In Fig. 11 we plot the the Wigner function for measurement-based squeezing applied to vacuum. We assume that the ancillary state has a fixed level of r = 1.2 (∼ 10.5 dB) of squeezing, and that the efficiency of the homodyne detection is 0.99. We consider three different target levels of squeezing, r target = 0.3, 1, and 2, to apply to the vacuum state, and calculate the Wigner functions in the case of ideal, direct inline squeezing; the average map of measurement-based squeezing; and a single-shot occurrence of measurementbased squeezing.
There are a few key observations we can make from the simulation. First, if r target is comparable to or greater than the level of resource squeezing r, then the variance in the q quadrature exceeds the desired value. This is to be expected: Eq. (67) tell us that, while the level of noise in the q quadrature is modulated by e −2r (constant across the simulations), it grows with r target . Next, we see that, in the average case, the level of anti-squeezing essentially matches the ideal case. This is due to the high efficiency of the homodyne measurement, which forces the level of p quadrature noise-proportional to 1−η η , per Eq. (67)-to be small. Additionally, the mean in the average case matches the ideal case, as guaranteed by Eq. (67). Finally, in the single-shot case, the location at which the state is centred along the p quadrature can vary on the order of the ideal anti-squeezing; moreover, the level of anti-squeezing in the single-shot case becomes significantly lower than ideal as r target becomes comparable to or greater than the level of resource squeezing r. This indicates that the accurate level of anti-squeezing and the correct mean observed on average is an artifact of the averaging process itself: the states that factor into the average have less anti-squeezing but a broad distribution of where they are centred along the p quadrature, resulting in a broadened average output state. The differences between ideal, average, and single-shot instances of measurement-based squeezing are crucial to understand for realistic implementations of bosonic codes; our formalism and simulations are valuable for such analysis. In the single photon case, negativity is produced by the difference of two zero-mean Gaussians with slightly different covariance matrices.  . For different levels of target squeezing on vacuum rtarget, a comparison of the states produced from ideal squeezing (red), the average map from measurement-based squeezing (green), and a single-shot occurrence of measurement-based squeezing (blue). The centre and axis lengths are set by the mean and variances of the output state. We simulate using an ancillary squeezed state of r = 1.2 (∼10.5 dB), and a homodyne detection efficiency of 0.99. A key takeaway is that the single-shot instances have less antisqueezing and a distribution of locations for where they can be centred in phase space, indicating the level of antisqueezing and the correct mean value observed in the average map are a result of the averaging procedure itself.

GKP Phase Gate
The phase gate is a pertinent example of a singlemode GKP gate that employs inline squeezing, an operation that may be performed, in practice, through measurement-based methods. In Appendix H 2, we provide details for such an optical implementation. Here, we present results of a simulation of a realistic phase gate applied to |+ gkp with = 0.1 (10 dB). In Fig. 12, we plot the Wigner functions of the average output states, varying the level of squeezing of the ancillary squeezed state and the efficiency of the homodyne detection. We compare the graphs to those for |+i gkp with = 0.1 (10 dB), which is the finite energy version of the ideal applicationP |+ gkp . We find several differences. First, while the finite-energy |+i gkp has a symmetric envelope in phase space, and each Gaussian peak has an isotropic spread about its mean, the application of a phase gate to a finite-energy |+ gkp results in an asymmetric envelope and anisotropic Gaussian peaks about each mean in phase space. This would be the case even if the phase gate was applied perfectly to a finite energy state: the finite-energy |+ gkp is non-periodic due to the envelope, and each peak has finite width; even a perfect phase gate would cause shearing effects in both the envelope and the individual peaks, an effect not seen in ideal states with perfect periodicity and infinitely narrow peaks. Second, as the quality of the squeezed ancilla resource and homodyne efficiency worsen, we find that the peaks in phase space broaden, albeit asymmetrically, decreasing the height of each peak.
The differences in the Wigner function do not tell the full story, however; we are also interested in how they affect readout of the state. In Fig. 13, we plot the marginal distribution along q − p and q + p of three of the states in Fig. 12. Recall that binning the outcome along q − p to the nearest n √ π and taking the parity of n corresponds, ideally, to the a measurement in the qubit Pauli Y basis; alternatively, one can bin the outcome along q + p to the nearest n √ π and take the parity of n + 1. This choice comes from the fact thatĤσ yĤ † = −σ y for qubits, wherê σ y is the Pauli Y operator andĤ is a qubit Hadamard gate operator, effected for GKP states by a CV phase space rotation by π/2. In practice, these two CV operators can be measured by performing homodyne detection along q±p √ 2 and then rescaling the outcome by √ 2. We find that the choice of readout quadrature affects the probability of obtaining the correct qubit result. Because the Wigner function is sheared in phase space, measuring along q − p results in a narrower envelope and narrower peaks. For the purpose of qubit readout, the peak width is what matters since the likelihood of falling outside a 0-bin does not depend on the number of peaks. The results are presented in Table I. There is a sizable drop in readout quality going from q − p to q + p, even though both represent the same ideal qubit measurement. This is further confirmation that the one-to-many mapping of p(reading out a qubit 0) Simulation parameters q − p q + p |+i gkp , = 0.  Fig. 13, we calculate the probability of correctly identifying a qubit 0 readout from a Pauli Y measurement, depending on whether the CV homodyne measurement is performed along q − p or q + p. q − p provides better readout, due to narrowed peaks in the marginal distribution.
qubit-to-CV operators affords GKP states an advantageous flexibility [21]. Using our simulator, we identify how to adapt the readout strategy based on the noise to obtain a more faithful measurement of the binary qubit outcomes.

GKP to CV Cluster Teleportation
CV cluster states-multimode Gaussian states constructed by stitching together momentum-squeezed states with CV CZ gates-are a valuable resource for measurement-based quantum computing with GKP states. In particular, a GKP state can be added to and teleported along the cluster using the same operations as would be used for a momentum-squeezed state [22]. The teleportation circuit begins with a GKP state and a momentum-squeezed state interacting via a CZ gate; next, a p-homodyne measurement is applied to the GKP mode, yielding outcome p 0 ; finally, a shift in q by −p 0 is applied to the remaining mode. Here, we seek to understand the dynamics of a realistic optical implementation of GKP teleportation into a CV cluster. As discussed in Appendix H 2, we break apart the CZ gate into beamsplitters and inline squeezing operations effected using measurement-based squeezing, and then investigate the effects of finite squeezing and loss in the teleportation circuit. The circuit requires three squeezed states: the ancilla mode onto which the GKP is teleported and the resource states for measurement-based squeezing within the decomposition of the CZ gate. We vary the level of squeezing of these Gaussian modes, while fixing the GKP state to be |0 gkp with = 0.1 (10 dB). Additionally, we vary loss within the circuit; we apply a loss channel of transmissivity η to each mode after initialization and at the output of each of the four beam-splitters in the circuit (two for the decomposition of the CZ gate and one for each inline squeezing operation). For the three homodyne detectors in the circuit, we fix an efficiency of η det = 99%.
For each value of squeezing and η, we run 500 simulations, the teleported state output by each run depending on the probabilistic feedforward value. Then, given the final state in each run, we measure the qubit Pauli X and Z operators, achieved by performing homodyne We compare this to a method for ideally preparing the same state by generating |+ gkp with = 0.1, then applying a GKP phase gateP = e i q 2 /2 . Since the phase gate requires inline squeezing (see Appendix H), we simulate measurement-based squeezing with various levels of ancillary squeezing and ancilla detector efficiency, which we label in (b)-(d). We see that phase gates apply a shearing effect to the envelope and peaks of the finite energy GKP state, with lower quality measurement-based squeezing adding additional broadening to the peaks. This simulation was performed with the bosonic backend of Strawberry Fields.
measurements along p and q, respectively, binning the outcome to the nearest n √ π, then taking the parity of n. Ideally, teleportation of |0 gkp onto a CV cluster yields |+ gkp , which would return a qubit 0 outcome 100% of the time when measuring Pauli X, and 50% of the time when measuring Pauli Z. In Fig. 14 (a), we plot-as a function of squeezing and loss-the mean probability of obtaining the outcome 0 from a Pauli X measurement. As expected, we see that high squeezing and low loss provide the best readout. Interestingly, we also find that above 10 dB of squeezing (an amount comparable to the per-peak squeezing of the teleported state), loss is the major determinant of readout quality for this Pauli measurement. We are also interested in how far individual runs can deviate from the mean probability of correct readout, so in Fig. 14 (b), we plot the spread of the probability from (a) as a function of squeezing and loss. We find that the spread is quite small (only a tenth of a percentage in the worst case), meaning that individual runs yield values close to the correct readout.
The story differs for the Pauli Z measurement. Fig. 14 (c), where we plot the spread for the probability of reading out a qubit 0 outcome in the Pauli Z measurement, shows that there can be significant deviation from the 50% mean value in the presence of finite squeezing and loss. While in the ideal case, every instance of the teleportation circuit yields a teleported state that has a 50-50 chance of generating outcomes 0 or 1 in a Pauli Z measurement, finite squeezing and loss can distort the readout odds: certain instances of the teleportation circuit (heralded by homodyne measurements of the ancillae) can cause one outcome to be favoured by several extra percentage points. Interestingly, in contrast to the Pauli X measurement, the spread of results in the Pauli Z measurement are more sensitive to the level of squeezing than to the level of loss.
While it is intuitive that low loss and high squeezing should yield qubit outcomes closer to ideal, it is interesting-for this choice of initial GKP states-that that loss (squeezing) is the greater noise source for Pauli X (Z) measurements. It is worth investigating, however, whether a different feedforward function or postselection of outcomes can improve the faithfulness of qubit readout, and whether the homodyne outcomes can be used to assign a confidence rating to the qubit readout of the teleported mode. Such strategies have been explored in the context of fault-tolerant quantum computing architectures with GKP qubits [1,10,72].

GKP T Gate Teleportation
While Clifford gates can be performed using Gaussian resources in the GKP encoding, one requires a non-Clifford-for GKP states, non-Gaussian-gate to achieve universality. As we discussed in Section V C, the qubit T gate can be applied via gate teleportation using a GKP magic state |M gkp . A realistic T gate teleportation is no small feat: an ancillary GKP state is required in addition to the three ancillary squeezed states for the CZ gates and the feedforward phase gate, and the multiple beam-splitters each incur loss. As far as we are aware, no investigation of realistic T gate application has been carried out, likely due to the difficulty of simulation. Let us close this gap.
First, in Fig. 15, we present the result of a T gate applied to |+ gkp with = 0.1 (10 dB); in the ideal limit ( → 0), this should return the state |M gkp . For sake of clarity, in this figure we assume that the CZ gate and the phase gate are applied perfectly, but employ a finite-energy resource state also with = 0.1 (10 dB). In (a) we present the Wigner function for a finite energy |M gkp to use for compari- Binning the result to the nearest n √ π and taking the parity of (a) n or (b) n + 1 effects a GKP qubit Y measurement. We present the marginals for Wigner functions (a), (b), and (d) from Fig. 12. While marginal distributions along either q − p or q + p could be sampled for a Pauli Y measurement, the narrower peaks along q − p yield higher likelihood of falling within the correct bin, emphasizing the value of tracking shearing effects to optimize readout fidelity.
son, while in (b) we provide the Wigner function for the output of the T gate circuit given a measurement outcome on the ancillary GKP mode that does not require the feedforward phase. In (c), we present the same as (b) but for an instance of the T gate circuit that does require feedforward phase. We see that the gate works to some extent: the positive and negative peaks in (b) and (c) are in the same places as positive and negative peaks in (a), and we recall from Section IV A that the logical information of the state is encoded in the weights. However, even with the CZ and feedforward phase gate implemented perfectly, we still see some differences between Ideally, the teleported state should be |+ gkp , with Pauli X (Z) measurement of the state yielding a qubit 0 outcome with 100% (50%) probability. However, with noisy teleportation circuits (loss and measurement-based squeezing in the CZ gate), each instance of the circuit, conditioned on the probabilistic value of the teleportation feedforward, yields a slightly different state. We simulate the noisy teleportation circuit 500 times. In (a) and (b), we plot the mean and the (log of the) standard deviation for the probability of a Pauli X measurement yielding 0, as a function of circuit loss and squeezing.
Having found the mean probability of a Pauli Z measurement yielding 0 to be roughly 50% across squeezing and loss levels, in (c) we provide the (log of the) standard deviation for the probability a Pauli Z measurement yields 0. Wigner function for T gate teleportation applied to |+ gkp with = 0.1 (10 dB) using the circuit from Fig. 6 and postselecting on an outcome that does not require the feedforward phase gate. (c) Output from the same circuit, but postselecting on an outcome that does require the feedforward phase gate. The gate teleportation circuit for the T gate with realistic components is involved, but now simulatable using the bosonic backend of Strawberry Fields.
the Wigner functions. First, notice that the p quadrature variance of each peak is now larger than the q quadrature variance, even though they began as symmetric. This is because the CZ gate adds the q variance of the resource magic state to the data state. Second, in (c), there is a shearing effect of both the overall envelope and the individual peaks due the phase gate, as also seen in Section VII B 2.
Next, we examine the effect of finite-energy magic states on the quality of the T gate teleportation, keeping the rest of the circuit ideal. We again simulate the T gate circuit applied to |+ gkp with = 0.1 (10 dB), this time varying the value of for |M gkp . For each , we perform 500 simulations, each time producing a different output state based on the probabilistic measurement of the ancilla. For the output state of each simulation, we calculate marginal distributions in phase space to determine the probability of reading out a qubit 0 outcome for the four Pauli measurements (three Pauli operators, and two ways of measuring Pauli Y). In Fig. 16, we plot how the mean probability of obtaining a qubit 0 readout value for the Pauli measurements varies as a function of the magic state's finite energy parameter . We additionally provide what the ideal probabilities for readout should be for the Pauli measurements. We find that, with smaller epsilon (increased per-peak squeezing), the statistics become more closely aligned with the ideal values, along with less variation in the output state of each instance of the circuit, with on the order of 11 dB seeming to be sufficient to achieve correct readout to within 1%. It is good news that the quality of magic state required is not significantly higher than that of the data state. We also note that measuring the Pauli Y operator by performing homodyne along q − p again yields more faithful readout of the qubit outcomes.
Finally, we perform the same simulation of a GKP  16. Mean probability of obtaining a qubit 0 outcome from different Pauli operator measurements after applying a T gate to |+ gkp with = 0.1 (10 dB). We vary the perpeak squeezing of the GKP magic resource state and see how well the qubit readout matches the ideal readout values for the Pauli operators (dashed lines). We see that a per-peak squeezing on the order of 11 dB (close to that of the data state) is sufficient to retrieve correct readout to within 1%. The simulation provides another example of the usefulness of having two methods for performing Pauli Y measurements, as seen in Fig. 13; homodyne detection and binning along q − p yields a more faithful outcome.
T gate applied to |+ gkp , this time implementing the CZ and phase gates using realistic components: lossy beam-splitters, inefficient homodyne measurements, and measurement-based squeezing. We fix the magic state to have = 0.1 (10 dB); we set the efficiency for the four Here, the gate teleportation is effected using imperfect optical elements (see Figs. 6 and 18 for circuit breakdowns). We fix the ancillary magic state to have the same as the data state; we employ 12 dB squeezed states for measurement-based squeezing, and set the homodyne efficiency to 99%. We see how readout changes as we vary η, which parametrizes the loss channel we apply after each of the five states are initialized and after each of the five beam-splitters.
homodyne measurements (three for measurement-based squeezing in the CZ and feedforward phase gate, and one on the ancillary magic state) to be 99%; and we choose a squeezing level of 12 dB for the three squeezed states used for measurement-based squeezing. We investigate how the Pauli operator readout varies with respect to a loss parameter η, which we apply in multiple places throughout the circuit: right after each state (GKP or squeezed) is initialized and after each beam-splitter. A few comments are in order. First, in the case of no loss, we already see that the readout is not as good as when the CZ and phase gates are applied perfectly; this is undoubtedly due to the change from idealized squeezing to measurement-based squeezing; even with squeezed resource states of 12 dB, the readout has room for improvement, motivating deployment of even more highly squeezed states. Second, as loss per optical element is increased, the readout quality quickly deteriorates. This is likely due to the number of loss channels the data state undergoes, since it must pass through four beam-splitters to perform the teleportation. Notably, only the Pauli Z readout is still close to the desired value, but this is coincidence: a broad, noisy distribution is essentially uniform over the bins applied to the quadrature outcomes, resulting in a 50% chance of reading out a logical 0 from a Pauli measurement. The sensitivity of the opticallyteleported T gate to loss, due to its many components, is additional motivation for pursuing passive optical im-plementations of computing with GKP states [65]; by potentially eliminating in-line squeezing, these architectures reduce opportunities for loss and noise brought in by measurement-based squeezing.

VIII. SUMMARY AND OPEN PROBLEMS
In this work, we have introduced a formalism for simulating continuous-variable quantum states. At its heart is a representation of these states as linear combinations of Gaussian functions in phase space. This novel framework can be used to analyze and simulate valuable classes of CV states, namely bosonic qubits like GKP, cat, and Fock states, under realistic transformations. Our mathematical framework inherits the simplicity and convenience of the Gaussian CV formalism, requiring only to track how transformations and measurements update the weights, means, and covariance matrices of the Gaussian functions in the linear combination. In addition to the straightforward inclusion of Gaussian channels in our formalism, a salient class of non-Gaussian transformations and measurements-photon-number-resolving measurements and gate teleportation for GKP and cat qubitscan also be simulated using our framework.
Enabled by our formalism, we provided simulation methods that outperformed state-of-the-art CV simulators that employ the Fock basis. Using our new method, we were able to perform novel simulations of bosonic qubits in realistic settings, informing future design decisions for quantum computers based on bosonic qubits. We focused on GKP qubits, examining how they transform under gates that leverage measurement-based squeezing, how they interface with CV clusters, and how they transform under non-Clifford gates implemented via gate teleportation. While we only investigated how the readout of Pauli operators was affected in these settings, we emphasize that our simulator outputs complete state information in the form of the elements required to construct the Wigner function, which can be used as input to other simulators or to more advanced analytical tools, such as the modular subsystem decomposition [23,73].
Numerical simulations leveraging our formalism and methods were performed in the new bosonic backend of Strawberry Fields, an existing feature-rich Python library for simulating CV optical circuits. The backend, implemented by the authors, is accessible through a userfriendly frontend of the public repository [25], along with some tutorials on its use [26-28]; we hope this encourages exploration by those interested in CV quantum information.
We anticipate several useful applications of our work. For one, our new formalism will prove valuable in the design of circuit modules for quantum computers that employ bosonic qubits, as it can identify, for example, where losses need to be reduced in a circuit, or what level of ancillary squeezing is required to attain a target gate fidelity. While our focus in the simulations has been on GKP qubits, similar analysis of realistic gate application for cat states can be performed. Second, it is known that in optical settings where deterministic high-order nonlinearities are currently lacking, bosonic qubits will need to be produced probabilistically in the near term. There are promising heralded state generation techniques that involve photon-number-resolving measurements on some modes of a multimode Gaussian state [21,74], a scenario which we have already discussed falls neatly within our formalism. Previous investigations and optimizations of this technique were limited by the speed and memory limitations of the Fock basis; a straightforward extension would be to use our new formalism to develop these bosonic state preparation methods even further. Yet another state preparation proposal for GKP states relies on interacting or "breeding" cat states with beam-splitters and homodyne measurements on some of the modes [75,76]. Since the circuits in those protocols consist of Gaussian transformations and measurements applied to cat states, our formalism applies here too.
Finally, there are several avenues for future research. First, it is worth investigating the connection between other analytical representations of CV states (such as Riemann-Theta functions for GKP qubits [41,70], or the shifted Fock basis representation of cat states [2]) and the ones we have presented here to identify any additional opportunities for faster or more accurate simulation. Second, for the class of states and maps that we have considered, it would be of interest to determine the exact mathematical conditions for physicality in terms of the parameters of the states (weights, vectors of means and covariances matrices). Moreover, while we considered a wide class of transformations which fall under this formalism, it would be valuable to determine the most general class of transformations which can be reduced to manipulating the weights, means and covariances of the Gaussian functions in the linear combination. Lastly, large-scale simulation of bosonic qubits under realistic noise models is challenging due to exponential scaling. Our formalism is a step towards numerically tractable models capturing some of those effects, either by developing suitable approximations of the states based on the mathematical framework provided, or by using the simulator to develop more detailed models for the evolution of qubit-level information or noise, then leveraging existing work in the field of quantum fault-tolerance for qubits.
are pure Gaussian states. To obtain the Wigner function recall that it can be calculated as (cf. Eq. 4.5.21 of Ref. [77]) where α, β ∈ C and we parametrize α = 1 √ 2 (q + ip). We now recall (cf. Eq. 3.6.30 of Ref. [77]) the composition law for displacement operatorsD and the expression for the vacuum to vacuum amplitude of a Gaussian transformation (cf. Appendix 1 of Ref. [17]) With these two expressions it is direct to write the integrand in Eq. (A2) as a Gaussian in the real and imaginary parts of β and perform the integral to find the Wigner function to be a normalized Gaussian in ξ = (q, p) with covariance matrix and (complex) means multiplied by the prefactor In the limit where δ → γ we find c = 1 and µ = √ 2 ( (γ), (γ)) as expected for a pure Gaussian state. Finally, note that so far we only considered outer productsÔ = |ψ φ| where the squeezing magnitude and direction are the same. One can consider the more general situation in which they are not. The final expression for the Weyl transform ofD(γ)Ŝ(re −iθ ) |0 0| S † (se −iφ )D † (δ) is a Gaussian in the real and imaginary parts of α, δ and γ given by

Appendix B: Coefficients of Ideal GKP
For ideal GKP qubits in the state |ψ = cos θ 2 |0 gkp + e −iφ sin θ 2 |1 gkp , the coefficients for the Wigner function are provided by [51]: for k mod 2 = 0, mod 2 = 0, This Appendix explicitly computes the Wigner function shown in eq. (43) using two distinct mathematical methods: direct calculation using the Moyal product in Section C 1, and the physical application of the Fock damping operator using an optical circuit in Section C 2.

Real Weights, Means and Covariances Via the Moyal Product
Recall that the Wigner function of the physical state W phys ψ (ξ) is given by [78]: where denotes the Moyal product, ψ represents the encoded qubit state, and W denotes the Weyl transform of the Fock damping operatorÊ( ) = e − n . The -product for a single mode is defined by where the arrow above the derivative symbols indicate to take the derivative to the function to the left or right-hand side. The -product is linear and associative, but not commutative. Since W ideal ψ takes the form we find where we have defined W ideal µ0 = δ (ξ − µ 0 ). Before calculating W phys µ0 we recall that the Fock damping operator is simply e − n = (1 + n )ρ thermal with n = (e − 1) −1 and ρ th is a thermal state (cf. Eq. (63)). We can then write the Weyl transform of the Fock damping operator as an unnormalized Gaussian With this result in hand we again use that the -product is associative, and we have the identity [79]: where we implicitly use the fact that A, B admit Fourier transforms. We first evaluate where N = 1/2 . Reorganizing the terms and performing the change of variable k → −k, we get The second integral is readily recognised as the characteristic function of a multivariate normal probability distribution, so we find Where we used the change of variable y → 2y / , k → 2k / . Performing the remaining integral, we find Defining the variables q = q + q 0 and p = p + p 0 and rearranging the Gaussians, we find Performing the integral and replacing x and p , we find Finally regrouping the Gaussians together, we find Replacing the values for N , noting that 2Σ / 2 + (2Σ ) with µ = 2e − 1+e −2 µ 0 .

Real Weights, Means and Covariances Via Optical Circuits
The index set, weights, means and covariances matrices for the ideal GKP state are provided in Eq. (41). Finite energy GKP states can be recovered by applying the Fock damping channelÊ( ) = e − n to the state, a channel which neatly fits into the Gaussian-inspired transformation framework introduced in Section III C. For the Fock damping channel with e − = cos θ, we derive the exact update rules for states in Appendix G. Thus, using Eq. (G1), the covariances of the ideal GKP under Fock damping become: while the means become: The re-weighting of the peak is provided by Eq. (G2): where N is an overall normalization. This completes the derivation.

Complex Weights and Means, and Real Covariances Via the Wavefunction to Wigner Transformation
The ideal GKP state for a single qubit is: Next consider application ofÊ( ) = e − n : Next we use Mehler's kernel [80]: where the second last line was obtained with Mehler's formula. In the last line: Therefore, we have that the (subnormalized) wavefunction in position quadrature is given by: Note that we can write the amplitude inside the sum as whereD andŜ are the single mode displacement and squeezing operators (cf. Eq. (13)) and |0 is the single mode vacuum. We can use these results to write the Hilbert space vector of a finite energy GKP state as showing that finite-energy GKP states can be written as linear combinations of single-mode pure Gaussian states like the ones discussed in Eq. (26).
Here, the result was obtained by performing the Gaussian integral.
Putting it all together, we find that the normalized Wigner function is a linear combination of Gaussian functions in phase space: with: where N (a, ) is an overall normalization chosen so that m∈M c m = 1, since the Gaussian functions are already normalized over phase space.
implies that the covariance matrix is invariant under any permutation that does not involve the first mode. Thus, for example, the following two projectors are equivalent when acted on the state |Ψ Ψ| One can state this assumption mathematically by writing the adjacency matrix (cf. Eq. 24 of Ref. [81]) of this state, from which the covariance matrix is [53] and S T is the permutation matrix that maps the vector (q 0 , . . . , q n , p 0 , . . . , p n ) into the vector (q 0 , p 0 , . . . , q n , p n ). Note that the vector of means of this state is µ = 0, and since we will be heralding using threshold detectors any heralded state will also have a zero mean vector. The exact permutation symmetry of the state allows us to group the 2 n terms in Eq. (D1) into n + 1 terms with binomial coefficients asQ Now consider the term with j vacuum detections. For this term we need to trace out n − j modes. One can write the covariance matrix of the j + 1 modes after tracing out in the ordering (q 0 , p 0 , q 1 , p 1 , . . . , q j , p j ) as where A is used to indicate the 0 th mode and B the remaining 1 through j modes and where the submatrices are . . .
and for which physicality alone entails r < 1/ √ n. In the last equation we used 1 j for the j × j identity matrix and Z = diag(1, −1).
We can now write the covariance matrix of the state after n − j modes are traced and j modes are projected onto vacuum. We use Eq. 5.142 from Ref. [32] For our case we can write the answer analytically Also note that the weight of this Gaussian function should be up to normalization, which is given by Thus we conclude that the coefficients are given by c j = p j /N n .

Appendix E: Cat States as Linear Combinations of Real-Valued Gaussian Functions
Equation (55) obviously contains complex terms. However, it is possible to approximate it by a real linear superposition with arbitrary precision. We give the explicit expression here.
The complex terms in Eq. (55) can be written as Recall that cos q can be expressed as [82] cos q = e with R(q) = O(e −2π 2 D ), with D a positive parameter controlling the accuracy of the approximation.
where µ α,vac m,m = Σ −1 α (Σ −1 α + V −1 vac ) −1 µ α m,m . While this general expression is quite cumbersome, it is straightforward to find the simpler form for the case where either (α) = 0 or (α) = 0 (the case of α = 0 being trivial). Following the same procedure in case where α ∈ R, we find e −2|α| 2 e −iφ G(ξ; V vac , µ z ) + c.c. ≈ √ π e π 2 D 4 4α (Σ α + v) where we have The Wigner function of the |0 state of the squeezed comb state is given by [43] W comb (q, p) = exp(−e −2r p 2 ) N π Regrouping exponential terms together, we find exp −e 2r (q − q m ) 2 + p 2 + e Here we show how a state with a Wigner function of the form in Eq. (23) is transformed by the Fock damping operator in Eq. (42). Recall that theÊ( ) operator can be viewed as arising from passing a state through a beamsplitter of transmissivity cos θ = e − with an ancillary mode in a vacuum state, then measuring vacuum on the ancillary mode [21,61].
Appendix H: Measurement-Based Gates that Employ Squeezed Ancillae In Table II, we provide a summary of the Gaussian CPTP maps corresponding to measurement-based squeezing, implemented using an ancillary squeezed state, a beam-splitter and homodyne measurement followed by a feedforward displacement.

Gate
Average Map q-squeezing X In Table III, we provide a summary of the map for a single-shot run of ancilla-assisted squeezing.  III. Update rules for each Gaussian peak in a linear combination from a single-shot run of ancilla-assisted squeezing in q quadrature. r is the squeezing level of the ancilla state, θ is the angle of the beam-splitter with cos θ = e −s defining the squeezing parameter s, η is the loss parameter of the homodyne detector, and pM is the p-homodyne outcome on the ancilla mode. To match the average case, one sets the feedforward operations to fx(pM ) = 0 and fp(pM ) = pM tan θ/ √ η, and integrate over pM .

Inline Squeezing Gate
Here we derive the map effected by measurement-based squeezing, as described in [20], and shown in Fig. 18 (a).
We see that we recover the partial trace condition, i.e. that Σ m,0 → c 2 θ Σ m,0 + s 2 θ Σ B,0 and µ m,0 → c θ µ m,0 , which is exactly what happens if we pass through a beam-splitter and trace out mode B.
Feedforward quadratic in p M Now we go back to a feedforward function quadratic in m. We then perform the necessary integrals in Eq. (H18) using Eq. (H19). We have that Let us further choose the linear feedforward to be g 1 = 0 and h 1 = η −1 tan θ, which is based on the suggestion in [20]. Then the output covariance matrix can be written as This is a Gaussian CPTP map, with the X part of the map corresponding to squeezing mode A in q by a factor cos θ. An analogous derivation is possible for squeezing in p by using a p-squeezed state, performing homodyne in q and feedforward displacement in p. Moreover, we can achieve complex squeezing value re iφ by placing the squeezing circuit between phase shifters.

Gates that Employ Inline Squeezing
Given our examination of inline squeezing using squeezed vacuum ancillae, we now summarize valuable CV gates which employ inline squeezing.  [20]. An ancillary position eigenstate is combined at a beam-splitter of angle θ with the target mode. The ancilla is then measured in p quadrature with efficiency η, and informing a feedforward q displacement of p0 tan θ/ √ η. The level of q-quadrature squeezing applied to the target mode is r = ln cos θ. In practice, the quality of measurement-based squeezing is dependent on the level of squeezing of the ancillary state, and loss in the beam-splitter. Squeezing in p can be achieved with a momentum eigenstate, a q quadrature measurement, and a p displacement. Complex squeezing parameters can be implemented by sandwiching the circuit between two phase shifters. (b) A CV phase gate applied using rotations and inline squeezing [19]. Here, r = arccosh 1 + s 2 4 , θ = arctan s 2 , and φ = − π 2 sign(s) − θ. (c) A CV CZ gate [19], with r = arcsinh− s 2 , φ = π 2 , sin 2θ = (cosh r) −1 , and cos 2θ = tanh r. For the GKP encoding, the qubit phase and CZ gates corresponds to s = 1. In practice, the inline squeezing in the phase and CZ gates can be implemented using measurement-based squeezing, and additional noise in the CZ gate can be modelled by lossy beam-splitters.
Phase gate The CV quadratic phase gateP (s) is just a squeezing gateŜ(re iφ ) composed with a rotation gatê R(θ). The ideal decomposition is given by [19,83] P (s) =R(θ)Ŝ(re iφ ), θ = tan −1 s 2 , φ = −sign(s) as shown in Fig. 18 (b), with values provided for the GKP qubit phase gate. Further noise could be assumed by using a lossy rotation gateL(η)R(θ). If we assume perfect rotation operations, then the performance of the phase gate is limited by the squeezing gate.
A simple noise model is to add single-mode losses to the beam-splitter outputsL(η) ⊗L(η)BS(θ). A CZ gate can be achieved by applied Fourier rotations on mode B both before and after the CX gate. The CZ gate is depicted in Fig. 18 (c) with specific values provided for the GKP qubit CZ gate. Amplifier channel Amplification could be a useful in tackling photon loss effects [15]. In the ideal average case, an amplifier modifies the covariance matrices and means in the following manner [32]: The ideal amplifier amplifies the input signal and adds a symmetric noise. To construct a realistic noise channel with κ = cosh r, we use its Stinespring dilation to get [60]: S 0,1 (r) =BS(π/4) † [Ŝ(r) ⊗Ŝ(−r)]BS(π/4), whereŜ 0,1 (r) is a two-mode squeezing operation. Here, lossy beam-splitters and the noise we examined for inline squeezing form the dominant sources of error.