Efficient simulation of Gottesman-Kitaev-Preskill states with Gaussian circuits

We study the classical simulatability of Gottesman-Kitaev-Preskill (GKP) states in combination with arbitrary displacements, a large set of symplectic operations and homodyne measurements. For these types of circuits, neither continuous-variable theorems based on the non-negativity of quasi-probability distributions nor discrete-variable theorems such as the Gottesman-Knill theorem can be employed to assess the simulatability. We first develop a method to evaluate the probability density function corresponding to measuring a single GKP state in the position basis following arbitrary squeezing and a large set of rotations. This method involves evaluating a transformed Jacobi theta function using techniques from analytic number theory. We then use this result to identify two large classes of multimode circuits which are classically efficiently simulatable and are not contained by the GKP encoded Clifford group. Our results extend the set of circuits previously known to be classically efficiently simulatable.


Introduction
Identifying quantum computing architectures that are capable of yielding quantum advantage, in contrast to classically efficiently simulatable ones, is of paramount importance, both at the fundamental level, and to design useful quantum machines capable of surpassing classical capability for computing tasks [1].
Quantum computing architectures based on bosonic fields, to which quadrature operators Cameron Calcluth: calcluth@gmail.com with a continuous-variable spectrum can be associated, are attracting growing interest as they allow for implementing resource-efficient quantum error correction with the use of bosonic codes [2][3][4][5][6][7][8][9][10][11]. It has recently been shown that concatenating bosonic codes with qubit-types of codes increases the threshold of tolerable error probability in quantum error-correction, with respect to the use of only qubit error correction [12][13][14].
One of the most promising bosonic codes is the Gottesman-Kitaev-Preskill (GKP) code [2,5,8,10,11], which allows for correcting arbitrary errors. GKP codewords also have a highly negative Wigner function [2,15,16], which is a known necessary and quantifiable resource for universal quantum computation with continuous variables [17,18]. Furthermore, GKP codewords corresponding to the 0-logical state have been shown to be universal in the resource-theoretic sense [19, 20] -i.e., they promote a Gaussian set of operations, states and measurements, to faulttolerant universal quantum computation [21].
In contrast to these results, there are certain circuits involving GKP states which are efficiently simulatable. For instance, it is known that encoded stabilizer GKP states in combination with special Gaussian circuits, made of discrete displacements and encoded qubit and qudit Clifford operations, measured with homodyne detection, are classically efficiently simulatable [22]. Hence, the separation between circuits involving GKP states which can be simulated on a classical computer and those which cannot remains unclear.
In this work, we tackle the unexplored question: are GKP states in combination with arbitrary displacements, Gaussian gates and homodyne measurement universal, or classically efficiently simulatable? For these types of circuits, neither continuous-variable theorems based on the non-negativity of quasi-probability distribu-tions nor discrete-variable theorems such as the Gottesman-Knill theorem can be employed to assess the simulatability.
Here, we prove that a large class of such Gaussian circuits is classically efficiently simulatable. The operations that we consider are a large subset of the Gaussian operations, which are not contained within the continuous-variable Clifford group and are therefore not simulatable by the previous methods used in Ref. [22]. To highlight the distinction from these previous methods, we consider the case of single-mode and multimode measurements separately.
The circuits in Ref. [22] are simulatable when the operations consist of multimode Gaussian operations whereby the parameters of the symplectic operations are selected from a zero-measure set and the parameters of the displacements are selected from the set of rational numbers. In this work, when restricting to the measurement of a single mode, we demonstrate simulatability for multimode operations whereby the parameters of symplectic operations are chosen from a set that is dense in the reals. Furthermore, for singlemode measurements, simulation can be achieved in time which increases linearly with respect to the number of modes. Meanwhile, we extend the class of simulatable displacements from rational to all possible continuous displacements.
Furthermore, for multimode measurement after multimode Gaussian operations, we demonstrate that another large set of symplectic operations and all real displacements are simulatable. None of these sets is contained within the set of operations given in Ref. [22] and therefore this work further increases the known set of simulatable operations acting on 0-logical GKP states.
The simulation methods given in this paper are distinct from both stabilizer methods [23-31] and phase space methods [17,18,32,33]. Our method is based on deriving an expression for the probability density function (PDF) with respect to position for any Gaussian operation acting on any initial states.
When choosing 0-logical GKP states as input and restricting the parameters of the symplectic operations to a set which is dense in the reals, this PDF can be evaluated efficiently. The evaluation of the total PDF relies on the ability to evaluate the PDF of individual rotated and squeezed 0logical GKP states, which we demonstrate is pos-sible for certain rotation parameters. This calculation is performed by analyzing the rotated wave function and simplifying it using techniques from analytic number theory. In particular, we simplify a Jacobi theta function into an expression involving the quadratic Gauss sum which has a solution for specific parameters.
In Section 2 we present the general circuit class which we investigate throughout this work; we describe Gaussian operations and outline how to track the evolution of measurement operators in the Heisenberg picture. In Section 3 we present our first key result, i.e., an analytic expression for the PDF of a GKP state with respect to position, which has undergone an arbitrary Gaussian operation. Furthermore, we show that for angles of rotation selected from a set dense in the reals, the PDF can be reduced to a Dirac comb using theorems from analytic number theory. We use this result in Section 4 to show our second major contribution, namely that a large set of multimode Gaussian operations followed by a single-mode measurement are efficiently classically simulatable. Our final main finding presented in Section 5 demonstrates that we can extend these results to multimode measurements on the condition that the multimode operations are chosen from a further restricted set of Gaussian operations. In Section 6 we summarize the sets of operations which are simulatable in both the single-mode measurement case and the multimode measurement case and compare them with discrete-variable simulation techniques. In Section 7 we provide an analytic expression for the rotated and squeezed realistic GKP state. We also provide a brief overview of the challenges of extending our results to the simulation of circuits involving realistic GKP states. Finally, our concluding remarks are presented in Section 8.

The circuits considered
We consider bosonic systems consisting of an arbitrary number n of bosonic modes, each characterized by canonical positionq k and momentum p k operators, with [q k ,p l ] = iδ k,l . In particular, we focus on the class of circuits schematized in Fig. 1.
The input states are 0-logical square GKP states 1 which have a wave function in position representation given by [2] all in all, the input state can hence be compactly indicated by the tensor product over all the n input modes initialized in the 0-logical GKP state, |0 GKP = ⊗ n j=1 |0 GKP j . The density matrix of the initial state is given bŷ  The evolutions that we consider are a subset of all Gaussian operations on n modes. We will define the class considered in Sections 4 and 5. Measurement occurs through homodyne detection, say corresponding to the measurement of the position quadraturesq of a single mode (Section 4) or all modes (Section 5). In the following, we review the evolution of the quadrature operators 1 The GKP kets, as defined in Eq. (1), do not formally represent proper quantum states as they are not normalizable, and as such do not belong to the Hilbert space associated with a bosonic mode. However, in the remainder of this paper, we will nonetheless adopt the colloquial expression GKP states, for the sake of brevity.
under Gaussian evolutions which we will later use to determine the PDF of Gaussian circuits with input GKP states.

Tracking evolutions under Gaussian operations
Analogously to the Clifford group in discretevariable quantum computation, Gaussian operations preserve the continuous-variable Pauli (or Heisenberg-Weyl) group HW(n), which is the group generated by the set {e ic jpj , e −id jqj : j ∈ Z n c j , d j ∈ R} [24]. They can be constructed by selecting and applying any number of the following generators in any order [24,34,35]: Our approach to tackle the computation of the PDF of GKP states after Gaussian operations starts with computing the Heisenberg evolution of the position operatorsq j and their conjugate momentum operatorsp j under the action of a unitary operatorÛ . The generator gates in Eq. (3) have the effect of transforming the quadrature operatorsq j andp j linearly. Displacements e ic jqj have a non-trivial effect only on the momentum of the j-th mode, given bŷ Similarly rotationsR(θ j ) = e iθ j (q 2 The Fourier transform can be defined as the rotation with angle π/2, i.e. F =R(π/2) = e iπ(q 2 The Fourier transform acting on encoded GKP qubits corresponds to the Hadamard gate. Squeezing operationsŜ(s j ) = e −i ln s j (q jpj +p jqj )/2 have the effect Finally, the SUM gateĈ X = e −iq jpk , where j is the control mode and k is the target mode, has the effectq j →q j ,p j →p j −p k q k →q j +q k ,p k →p k . (8) As can be seen by these Heisenberg evolutions, we can track the evolution of the measurement operator associated with the modes of interest using a linear equation of the form i.e., we need to track 2n+1 real variables for each mode we would like to measure. To measure n modes we need to track (2n + 1)n real variables. Similarly to Ref.
[24], for any finite precision, this can be done efficiently as long as the number of gates is polynomial in the number of modes.

Decomposition of Gaussian operations
We will now analyze the effect of any quadratic Hamiltonian on the measurement modeq j . We first introduce the vector of operators [36, 37] r = q 1 , . . . ,q n ,p 1 , . . . ,p n T (10) and the group of symplectic matrices where Any unitary operator generated by a quadratic HamiltonianĤ =r T Hr/2 can be associated to a symplectic matrix M aŝ This operation acting on the operatorŝ r j →Û † Mr jÛM can be expressed in terms of a mapping of the vector of operatorŝ Any linear displacement can be achieved using an operatorDr parameterized by a vector of 2n real numbersr corresponding to the magnitude of displacement in position and momentum. Its effectD † rrDr on the vector of operators can be expressed aŝ r →r +r. (16) Utilizing this formalism, unitary operations corresponding to any quadratic Hamiltonian can always be decomposed as a Gaussian operation of the formÛ which acts on the quadraturesr aŝ We will denoteÛ (r, M ) =Û for simplicity, and will useÛ throughout the rest of this work. The matrix M can be written in block form as The effect ofÛ on the Heisenberg measurement modes can be evaluated as linear equations as in Eq. (9). Therefore, we see that any quantum computation composed of Gaussian operations in any order followed by measurements inq j can be decomposed in terms of a symplectic operation, followed by linear displacements in position with parameters and the block matrices Note that the values r i for i ≤ n can be any real value since they correspond to displacements in momentum, which commutes with the position measurement operator. In the following sections, we will show how to simulate a restricted set of Gaussian operations acting on 0-logical GKP states, followed by measurements inq j .

PDF of a Gaussian-transformed GKP state
The PDF with respect to position of a 0-logical GKP state following a general single-mode Gaussian operation consists of a Gaussian operator applied to an infinite sum of Dirac delta peaks. The Gaussian operator will introduce non-trivial coefficients at each of these peaks. It is known that for a restricted set of Gaussian operations -e.g., those which can be described in terms of GKP-encoded single-qubit Clifford operations -it is possible to analytically evaluate the wave function, and therefore also the PDF of a transformed GKP state [2]. However, to the best of our knowledge, no general result is known for arbitrary Gaussian operations.
In this section, we provide the evaluation of the PDF of a 0-logical GKP state that has undergone a generic Gaussian transformation -which generally cannot be described in terms of GKPencoded qubit Clifford operations -in the position basis. We note that we can decompose any single-mode symplectic operation using the Iwasawa decomposition [38] in terms of a shear operation followed by a squeezing operation and then a rotation. Furthermore, when measuring in the position basis, we can discard the initial shear operation, as it has no effect on the position quadrature. Therefore, the effect of all singlemode Gaussian operations acting on the position quadrature can be expressed in terms of a displacement, squeezing and a rotation. Since the displacement has the effect of a displacement of the position variable in the PDF, we can, without loss of generality, ignore the effect of displacement in the following calculation.
The wave function of the rotated and squeezed 0-logical GKP state can be expressed in terms of these single-mode operations as Furthermore, we can isolate the squeezing parameter s, which has the effect of rescaling the wave function of a rotated 0-logical GKP state, The full calculations detailed in Appendix A.1 show that the PDF of a rotated 0-logical GKP state measured in the position basis can be written compactly as being proportional to where the Jacobi theta function ϑ(ζ; τ ) is defined as We consider two cases for the rotation angle θ, for which the PDF can be evaluated analytically. For case 1 we consider rotation angles θ for which cot θ = u/v where u ∈ Z, v ∈ Z odd . For case 2 we consider rotations angles where θ mod π = 0. These cases correspond to selecting θ ∈ Θ where the set Θ is defined as and Q (2) is formally defined as the localization of the integers Z at the prime ideal 2Z [39]. First, we summarize the calculations of case 1, where the rotation angles can be expressed as cot θ = u/v. By expanding the theta function in terms of u, v, it is possible to express it in terms of a Dirac comb and a quadratic Gauss sum. The Gauss sum is given by (27) and the theta function can be written in terms of the Gauss sum as The delta function evaluates to 0 except at values of ζ equal to n /v so the theta function can be simplified to Using results from analytic number theory we identify that for the case of θ we consider, the Gauss sum evaluates to a constant. This implies that the norm of the theta function is proportional to a Dirac comb. Hence, the PDF for this case can be evaluated to where we have ignored the normalization constant as in Ref. [2]. Formally, this PDF is not normalizable, however, we can still interpret it as a probability density function in that it informs us that the measured position value x will be given by any integer m multiple of √ π sin θ/v with equal probability.
In the second case, we identify that a rotation of θ = 0 gives the identity, and a rotation by θ = π is equivalent to applying the Fourier transform twice, which is also the identity. Therefore, we can immediately write the PDF as We can unify both of these cases in a single PDF, with separation between peaks ∆ as calculated in Appendix A.2, yielding where gcd(u, v) = 1. The PDF of the rotated and squeezed mode can therefore be expressed as where we again ignore the normalization constant.
4 Simulatability of Gaussian circuits with single-mode measurement  For a restricted class of the circuits shown in Fig. 1, where the restriction arises from the two cases of allowed angles in Eq. (26), it is possible to simulate the outcome of the measurement. Specifically, we can consider a restricted class of multimode Gaussian operations A acting on 0logical GKP states followed by measurement of a single mode in the position quadrature, see Fig. 2.
We wish to establish whether the PDF of the corresponding measurement outcomes where x 1 are the eigenvalues of the position operator on the first modeq 1 |x 1 = x 1 |x 1 and ρ =Ûρ 0Û † , can be calculated efficiently classically. This refers to the notion of strong simulatability, whereby the PDF of the outcomes is calculated in a time which scales at most polynomially both with the number of modes and the number of digits of precision [40][41][42].
We begin by detailing the set of Gaussian operations A that are simulatable by our technique. We will then show that the PDF of the multimode operations in A can be decomposed into the evaluation of an integral involving functions which correspond to the single-mode PDF of a rotated and squeezed 0-logical GKP state, derived in Section 3. The restricted set of multimode Gaussian operations we consider ensures that each singlemode PDF satisfies the constraints defined in the previous section.
The restricted class of operations that we consider, defined as is the direct product 2 of the Heisenberg-Weyl group of all real displacements in phase space and the set of symplectic matrices M ∈ RSp(2n, R). The set RSp(2n, R) consists of all elements M ∈ Sp(2n, R) such that one of the following cases is satisfied where A, B are the block matrices of the symplectic matrix M as defined in Eq. (19). 2 Note that in the definition of the Gaussian operations, the semi-direct product of the Heisenberg-Weyl group and the symplectic operations provides a complete description of all possible Gaussian operations since both the set of Heisenberg-Weyl operations and the set of symplectic operations are groups [24]. However, in the restricted class of circuits we must use the direct product × since the restricted set of symplectic operations is not a group. This is equivalent to saying that the operations are selected from any combination of the two sets.
In other words, the class A contains operations of the form ofÛ in Eq. (17) where (r, M ) is selected from the set of all possible phase space displacementsr ∈ HW(n) and symplectic operations parameterized by the symplectic matrix M are such that the matrix satisfies the constraints in Eq. (37). We would like to demonstrate that when the operationsÛ are selected from the class A, acting on 0-logical GKP states, followed by single-mode homodyne measurement ofq, the circuit is simulatable.
To calculate the PDF of the single-mode measurement circuit which is shown in Fig. 2, we first track the evolution of the position quadrature in the Heisenberg picture, as in Eq. (9), which involves a summation over the n modes of terms of the form where we have introduced the re-scaled parameters s i and θ i such that a i = s i cos θ i and b i = −s i sin θ i , and where we have defined the rotated quadratureẑ i . This is equivalent to performing a squeezingŜ(s i ) operation, followed by a rotationR(θ i ), where s i ∈ R and θ i ∈ [0, 2π). The coefficients a i , b i , s i , θ i depend on the mode tracked but we will omit the mode index in this section because we are only tracking a single mode, which we chose to be the first mode without loss of generality.
By evolving the measurement operator, rather than the states, we can express the PDF in Eq. (35) as where the Heisenberg measurement operator can be calculated according to Eq. (9) andρ 0 is given in terms of Eq. (2). Using Eq. (38) we can then write each term of the sum as the Heisenberg operator for a transformed mode, which allows us to calculate the PDF of the measurement ofQ 1 . We can express the projection operator as a Dirac delta function, which allows us to write the PDF as This can be evaluated by inserting the identity over all the transformed modes, where z is the n-vector of the integration variables z j . This provides a PDF in terms of the wave functions of the individual transformed modes, where we have identified the PDF of individual transformed modes as Inserting the PDF of the individual modes given by Eq. (34) we find that the PDF of the measurement can be written compactly in terms of the parameters s i and ∆ i derived from the original symplectic matrix, Note that this is possible when the symplectic matrix satisfies Eq. (37), in virtue of Section 3. Further calculation details are provided in Appendix B.1. The simulation of sampling from this circuit can be performed by generating n random integers m i and returning The Gaussian operations selected from the class A can be parameterized in terms of the symplectic matrix elements A 1,i and B 1,i . From these variables, we can calculate each ∆ i using arithmetic operations and the greatest common divisor following Eq. (33) and Eq. (37). The squeezing parameters s i can be evaluated with trigonometry from Eq. (38). Furthermore, each of the variables ∆ i , s i can be calculated independently. Therefore, the total number of these independent calculations scales linearly with the number of modes n. Hence, we conclude that the PDF of the circuit can be calculated in linear time with respect to the number of modes. Therefore, the strong simulation of the class of circuits shown in Fig. (3) can be performed in linear time and hence exists within complexity class P.
For an example of a multimode circuit, involving the SUM and Fourier transform gates, which is simulatable by our technique, refer to Appendix B.2. A circuit diagram with an explicit construction of a general circuit is given in Appendix B.3.
These results are analogous to an input-GKP version of the results of Ref.
[24], which show simulatability for Gaussian circuits and Gaussian measurements for position eigenstates. However, we stress that the proof techniques used for the derivation of our results are different and that the input GKP states are non-Gaussian and highly Wigner negative [2,15,16], in contrast to the Gaussian input states of Ref. [24].
Practical tools for simulating Gaussian circuits with input GKP states were developed in Ref. [43], based on representing states as linear combinations of Gaussian states. Note however that their method is not computationally efficient for the circuits considered here, unlike ours. Tracking the evolution of the field quadratures is also used in the context of noise propagation under the twirling approximation to determine the fault-tolerance threshold for surface-GKP codes [13,14,44].
Finally, note that the 0-logical GKP state is possible to simulate using our technique due to the fact that its wave function can be written as a Dirac comb which does not have any positiondependent coefficients. We could also use other non-magic GKP states as input. For example, the state |1 GKP and the encoded X L basis state |+ GKP both satisfy the property that the peaks do not have position-dependent coefficients. Note that the encoded X L basis state |− GKP has an alternating sign for each peak, however, it can be generated using a momentum displacement applied to |+ GKP and therefore is simulatable. If we instead consider the magic T-logical GKP state,

Simulatability of Gaussian circuits with multimode measurement
|0 GKP q n Figure 3: Schematics of the multimode measurement circuits considered. As in Fig. 2, the input states are 0logical GKP stabilizer states. The operations considered belong to a further reduced set of Gaussian operations. Homodyne detection of multiple modes follows, corresponding e.g. to the measurement of the quadratureŝ q 1 , . . . ,q n .
If we make an additional restriction to the set of symplectic matrices that we wish to simulate, we can generalize the results of the previous section to multimode measurement. We will denote this set of matrices DSp(2n, R). The set includes matrices of the form Eq. (19) which can be decomposed to whereÃ,C are n×n matrices such thatÃ is nonsingular and bothÃ andC TÃ are symmetric 3 . The allowed angles θ = (θ 1 , . . . , θ n ) T are selected 3 A 2n-dimensional symplectic matrix can always be decomposed as [38] We make a restriction by choosing X = diag(cos θ) and Y = diag(sin θ). We also note the connection to a similar decomposition given in Ref. [46,47] which could also be used to define a more restricted class of simulatable operations.
from the set θ j ∈ Θ as defined in Eq. (26). Equivalently, we could restrict the symplectic matrices in Eq. (19) such that the block components satisfy A =Ãdiag(cos θ) and B =Ãdiag(sin θ) for some symmetricÃ.
When including the Heisenberg-Weyl operations we denote the full set of simulatable operations as We will now demonstrate that when the opera-tionsÛ are selected from the class B, acting on 0-logical GKP states, followed by multimode homodyne measurement ofq, the circuit is simulatable. The class B contains operations of the form U in Eq. (17) where (r, M ) is selected from the set of possible phase space displacementsr ∈ HW(n) and symplectic operations given in Eq. (48). In Section 6 we will show that the class B is contained in the class A defined in Eq. (36). The PDF of a general stateρ measured in the position basis in all modes is given by where and x = (x 1 , . . . , x n ). Similarly to the singlemode case, rather than evolving the state we can evolve the position measurement operators. By analyzing the decomposed matrix in Eq. (48) we can identify two independent operations which transform the measurement quadratures in a convenient way. The action of the second matrix in Eq. (48) can be expressed as a series of singlemode rotationŝ which rotate each mode independently according We then inspect the unitary operation corresponding to the first matrix withÃ in block (1, 1) which is given byÛÃ and has the effect where a (j) i is the element in the (j, i) position of the matrixÃ.
Eqs. (53) and (54) allow us to express the Heisenberg measurement operators aŝ We can also include arbitrary displacements in position by including a displacement operation e −ip·c before measurement, which results in the total Heisenberg measurement operator This can be expressed in terms of the individual rotated modes aŝ where, as said, a Writing the PDF in terms of the Heisenbergevolved measurement operators and using the cyclic property of the trace we find where we use the same notation as Eq. (51) to denote the simultaneous eigenkets of the operatorŝ Q j over x j , for all j ∈ {1, . . . , n}. Expressing the projection operators as Dirac delta functions, we have We can then substitute the expression for the Heisenberg evolved measurement operators, Eq. (57), to find Inserting the identity Eq. (42) and evaluating the integral over the resulting delta functions results in the final PDF The full calculation details of the PDF are given in Appendix C.1. If the symplectic operation is given in terms of the decomposition of Eq. (48) then it is possible to calculate all n of the parameters ∆ i in linear time by the same considerations as reported in Section 4. The PDF involves a n × n matrix of coefficients a (j) i and so evaluating the PDF involves a total of n 2 coefficients of this form. We have therefore demonstrated an efficient method for evaluating the PDF for any circuit of the type given in Fig. 3.
In summary, due to the form of the matrix decomposition in Eq. (48), we can simulate any circuit with modes initialized as 0-logical GKP states acted upon by parallel single-mode rotations parameterized by θ j ∈ Θ, followed by any symplectic operation parametrized by a symplectic matrix which when written in block form, as in Eq. (19), has block component B = 0. The restriction to generate a symplectic matrix such that the block component B = 0 is equivalent to ensuring that the operation does not transform q j →p k or vice versa, for any j, k. The SUM gate satisfies these criteria. Therefore, we can construct arbitrary examples of circuits that are restricted to the form in Eq. (48) by first choosing single-mode rotations parameterized by θ j ∈ Θ, followed by any combination of squeezing operations and SUM gates. An example simulation of a Gaussian operation acting on GKP states is given in Appendix C.2. A circuit diagram of a general circuit chosen from the set of allowed operations B is given in Appendix C.3.
Similar to the result given in Ref.
[24], it is clear that certain feed-forward operations can be efficiently simulated too, when they can be substituted with Gaussian two-mode operations followed by deferred measurements. In particular, Pauli operators implemented after classical feedforward can still be simulated efficiently.
It may also be possible to provide an alternative proof of the simulatability of these circuits by adapting the formalism developed in Ref.

Summary of simulatable operations
We have defined two new restricted classes of Gaussian operations, given in Eq. (36) and (49), in terms of sets of restricted classes of symplectic operations, RSp(2n, R) and DSp(2n, R). Opera-tions selected from the restricted class A, applied to 0-logical GKP states, followed by measurement of a single mode are simulatable. Operations selected from the restricted class B, applied to 0logical GKP states, followed by measurement of multiple modes are simulatable.
In Ref.
[22] it is proven that using the Gottesman-Knill theorem, one can simulate any continuous-variable quantum computation initialized to 0-logical GKP eigenstates, acted upon by encoded Clifford group operations C d for GKPencoded qudits in arbitrary dimension d, and measured with homodyne measurement.

HW(n)[Sp(2n, R)]
A B C d Figure 4: The classes of circuits we consider are displayed as a Venn diagram. We know that all A, B, C d are contained within the set of Gaussian operations HW(n)[Sp(2n, R)]. We also know that the logical encoded Clifford group does not completely contain, nor is completely contained by A or B. Finally, we know that B is contained by A. Note that the size of each of these regions in the diagram is arbitrary. It is however meant to represent the fact that the parameters of the symplectic and displacement operations are selected respectively from a zero-measure set and the set of rational numbers, for C d , and from a dense set and all reals for A.
The relationships between the classes A, B, C d demonstrate the power of the new simulation techniques outlined in this paper. First, we know that all A and B and C d are contained within the full set of Gaussian operations We also note that The relationship of Eq. (63) can be understood by considering that any symplectic matrix in DSp(2n, R), which satisfies the constraints of Eq. (48), must also satisfy the constraints of RSp(2n, R), given in Eq. (37). This is because the constraints of DSp(2n, R) given for the first row of the symplectic matrix are equivalent to all of the constraints of RSp(2n, R). Furthermore, the relationship of Eq. (64) can be seen immediately by considering that, for example, arbitrary displacements, which are contained in A and B, are not contained within C d . These relationships are presented as a Venn diagram in Fig. 4. Eq. (64) shows that there exists a large class of circuits which previously were not known to be simulatable which are simulatable by the methods detailed in this paper.
We will now compare the results of this work with those of Ref. [22]. It is always possible to define Gaussian operations in the form of Eq. (17) in terms of the real parameters of the symplectic matrix and displacements. The methods described in Ref.
[22] allow us to simulate only those parameters selected from a zero-measure set for the symplectic matrices and all rational displacements. In contrast, the methods described in this work demonstrate that all circuits selected from A followed by single-mode measurement are simulatable. The parameters of A constitute a set which is dense on the reals and thus significantly expands the class of simulatable Gaussian operations. The proof of this fact is given in Appendix D. Note, however, that the density of the parameters defining the operators should be considered a mathematical property characterizing the class of simulatable operations. The density does not imply that it is possible to simulate the PDF of operations outside of the set, i.e. it does not mean that operations in the closure of the set are necessarily simulatable.
Furthermore, in Appendix E we show that which informs us that neither A nor B completely contain the Clifford group operations which are simulatable by previous methods. The classes described in this section are those that we have shown to be simulatable in our work and are not necessarily fundamental. It is possible that there exists a new class of simulatable continuous-variable operations acting on 0-logical GKP states, e.g. D, which contains all of A, B, C d , which remains undiscovered.

Realistic GKP states
Here we briefly sketch how to extend the analysis that we have performed to the case of realistic GKP states. In principle, our simulation method could be extended to such case. We show that an analytical expression for the PDF of a rotated and squeezed single-mode GKP state, analogous to Eq. (34), can also be derived for the case of realistic GKP states. However, due to the more complex form of this PDF, which displays Gaussian peaks rather than a Dirac comb, we cannot directly apply the same derivation as before to evaluate the PDF after multiple modes have undergone symplectic operations.
The rotated and squeezed wave function of a realistic GKP state is given by which can -as in the case of infinitely squeezed GKP states, see Eq. (23) -be expressed in terms of the wave function ψ ∆ θ (x) of the rotated GKP state, i.e.
In Appendix F we demonstrate that it is possible to write this expression analytically in the form which is defined in terms of the constants and results in a single-mode rotated realistic GKP state having a PDF of the form The PDF of a circuit with input realistic GKP states followed by Gaussian operations and homodyne measurements is given by Solving this expression analytically would require solving integrals over a product of Jacobi theta functions, which hinders us from using the same techniques used for the ideal case. We are unaware of a method to evaluate this expression. We provide a more detailed overview of the challenges of evaluating the PDF for realistic GKP states, for both single-mode and multimode measurement, in Appendix F.

Conclusion
In this work, we have shown that large classes of Gaussian operations, in combination with input 0-logical GKP states and homodyne measurements, are classically efficiently simulatable. The proof structure used to prove our result is based on decompositions of the symplectic matrix and directly calculating the corresponding PDF, hence introducing very different methods compared to the ones used to prove the simulatability of a restricted class of Gaussian circuits in Ref. [22]. As a matter of fact, compared to C d , for single-mode measurements we define a larger class A of simulatable operations which includes symplectic operations with parameters which are dense in the reals and all continuous displacements. Meanwhile, for multimode measurements, we define a new set B, not contained by C d which is also simulatable. Identifying circuits which are simulatable is, beyond fundamentally relevant, also useful for the development of quantum computers, providing a possibility for benchmarking and verification. Indeed a quantum computer can be initially programmed to solve problems which are simulatable with classical computers, as a test to identify whether the computation result is accurate [48].
Our proof of efficient simulatability only holds for the case of infinitely squeezed 0-logical GKP states, and we have sketched the difficulties encountered when trying to extend this result to the case of realistic, i.e. finitely squeezed, GKP states [2]. Establishing conclusively whether realistic GKP states are also simulatable following Gaussian operations is hence left for future work.
We acknowledge useful discussions with Laura García-Álvarez and Juani Bermejo-Vega. G. F. and C. C. acknowledge support from the VR (Swedish Research Council) Grant QuACVA. G. F. and C. C. acknowledge support from the Wallenberg Center for Quantum Technology (WACQT). We begin by analyzing a simplified version of the circuit in Fig. 2 where we have one mode andÛ is the unitary operator corresponding to a single-mode symplectic operation. This would correspond to a circuit of the form given in Fig. 5. We will then generalize to include displacement in phase space (i.e. displacement inq orp). The general form of the transformed quadrature under symplectic operations for an individual mode is given byÛ †qÛ = aq + bp. (72) If we are given or can calculate the real coefficients a, b ∈ R in Eq. (72), then it is possible to reconstruct a single-mode operationÛ which could generate such a transformation. By using only a certain set of operations we can calculate a decomposition ofÛ which will help us to calculate the PDF after measurement of the position quadrature by applying that operation to the GKP state and calculating its wave function and hence also its PDF. We will begin by first deriving an appropriate decomposition, inspired by Ref.
[38] (but ignoring the shear gate which acts solely on thep quadrature), which can be applied to the symplectic operation for all values of a, b ∈ R.
Within the symplectic formalism, the rotation operatorR(θ) = e iθ(q 2 +p 2 )/2 has the effect of trans-|0 GKP Sp(2, R)q Figure 5: The measurement of a single mode following an arbitrary symplectic operation acting on the |0 GKP state can be tracked by considering the evolution of the position quadrature in the Heisenberg representation.
forming the mode quadratures as while the squeezing operatorŜ(s) = e −i ln s 2 (qp+pq) has the effect of Then, we can calculate the wave function of the transformed mode as We can interpret calculating each summand as a wave mechanics problem with Hamiltonian H = 1 2 (q 2 +p 2 ) parameterized with θ. This becomes which is the propagator for the simple harmonic oscillator [49] and can be evaluated to By completing the square, we can rewrite the propagator in the form which gives the expression for the wave function where we have used the definition of the theta function,i.e.
This gives a succinct analytic expression for the PDF of measuring theq quadrature of a GKP state under any arbitrary single-mode symplectic operation, In order to incorporate phase space displacements we note that any Gaussian operation involving a symplectic operationÛ Therefore, any Gaussian operation acting on a single mode, for which we will measure theq quadrature, can be decomposed asÛ

=DrŜ(s)R(θ). (87)
The PDF of measuring theq quadrature of a single-mode GKP state under any Gaussian operation can therefore be written as wherer 1 = c is the magnitude of displacement in position. The PDF involving the theta function in Eq. (88) can be evaluated for certain values of the angle θ, or equivalently certain values of τ . The angle θ cannot take values of the form kπ for integer k ∈ Z, but as we shall see later, these rotations are trivial to evaluate using a different method.
We now consider different cases where we can evaluate the PDF.
The theta function can be reduced to a sum over Gauss sums whenever cot θ = u/v is a fraction, i.e. a rational number. We also assume that the fraction u/v is written in its simplest form, i.e. gcd(u, v) = 1. If it is not then it should be reduced as such by first calculating gcd(u, v), which can be calculated in polynomial time with respect to the size of u, v [50]. Then, exp(πiτ ) = exp(2πiu/v) becomes a root of unity and so the summation repeats in cycles [51].
where we have identified the Gauss sum [52][53][54] By using the fact that the delta function evaluates to 0 for all values of ζ not equal to n /v for some n ∈ Z, we can rewrite the theta function as where the Gauss sum needs only to be identified for certain values of ζ, rather than the entire real axis. For these values of ζ, the Gauss sum is given by Since we have confined to the case that gcd(u, v) = 1 and v ∈ Z odd then we immediately have the result that is a fixed constant dependent only on u, v [53], where u v denotes the Jacobi symbol, a generalization of the Legendre symbol.
Then, we can substitute this into Eq. (91) We can reinsert the theta function into the probability in Eq. (88) to get (Case 2: θ = 0, π) When θ = 0, π the probability function appears to break down because 1/ sin θ is undefined. However, we can clearly identify that θ = 0 is simply the identity. For θ = π, this is equivalent to applying the Fourier transform twice, which on the |0 GKP state has the effect of the identity. In both cases we obtain

A.2 Evaluation of the separation between the peaks
We can now encompass both cases above, and write a general PDF for all rotation angles θ ∈ Θ that we define in Eq. (26), in terms of the separation between peaks. Note that the squeezing parameter corresponds to a rescaling of the wave function in position space, and we can trivially write the squeezed and rotated GKP PDF in terms of the rotated PDF, i.e.
Ignoring the normalization constant, as in [2], we can write a formula for the PDF of a general rotated mode Note that in the special case when θ = π/2 we have cot(π/2) = 0, so we identify u = 0 and v = 1, as the simplest form of the fraction u/v, to arrive at ∆ = 1/s. This coincides with what we would expect when performing a Fourier transform (equivalent to a π/2 rotation) on the input state in the logical basis F |0 GKP = |+ GKP = m |q = m √ π followed by a measurement.

B.1 Calculating the PDF
For the class of circuits that we consider in Section 4, namely those belonging to the class A defined in the main text, the PDF is given by whereρ is the state of all of the modes after the multimode operationÛ (i.e.ρ =Ûρ 0Û † , whereρ 0 is the density matrix of n modes prepared in the 0-logical GKP state |0 GKP ), and |q 1 = x 1 is the position eigenket for the first mode. Using the cyclic property of the trace, and then evaluating it as a trace over all the other modes in the position basis, we find that the PDF can be written as The PDF can then be evaluated by inserting the identity over all rotated quadraturesẑ θ i i which is given by and by writing the projection operators as a Dirac delta function [36], Using the expression for the Heisenberg evolved operator Eq. (40) allows us to evaluate the PDF as Then inserting the PDF of the single-mode rotated 0-logical GKP state given in Eq. (98), we can express the PDF of the single-mode measurement on a multimode circuit as where we have the ∆ j specified by Eq. (99) for each mode, i.e.

B.2 Simple Example
In this section, we provide a simple example to illustrate how our decomposition methods works. We would like to simulate the following operationŝ acting on GKP input states, i.e.Û For simplicity in this example case, these operations are chosen to be Clifford operations in terms of the encoded qubit. By considering the logical operations only, we know that the result of this circuit should be so that measurement of the quadratureq on a single mode (e.g. the first mode) will give a PDF Using our technique we can track the operations on the quadraturê Next, analyzing each quadrature we see that for which we can identify θ 1 = 0 s 1 = 1, From these values, we can calculate the ∆ j from Eq. (106) as This results in the following PDF: In this example, a further simplification can be made by reparameterizing the variables m 2 → m 2 − 2m 1 such that which informs us that the variable x 1 will be measured as some randomly selected integer multiple of √ π. This agrees with the qubit-based simulation at the beginning of this subsection since we measure, with equal probability, an outcome which would correspond to a |0 GKP or |1 GKP state.

B.3 Simulatable operations in single-mode measurement circuits
The class A of operations, for which single-mode measurements are possible to simulate, can be interpreted in two ways. In both of these interpretations we focus on the allowed symplectic matrices because we can simulate all displacements in phase space and insert them freely without changing the structure of the symplectic matrix [37], see Eq. (86).
The first interpretation can be seen directly from the calculation of the PDF in the Heisenberg picture in Eq. (104). Namely, we can understand the measurement mode as undergoing entangling SUM operations which transform it toq 1 →q 1 + · · · +q n . This Heisenberg-evolved measurement operator is then evaluated on the rotated and squeezed GKP state. We can therefore picture the class of simulatable circuits as those whereby the input GKP states undergo rotations by angles θ j ∈ Θ and squeezing operations followed by SUM gates from each mode onto the measured mode, and by an arbitrary symplectic operation on the remaining n − 1 modes. This is depicted in Fig. 6. Figure 6: The class of simulatable operations A for single-mode measurement as understood from Eq. (104). These are circuits initiated with 0-logical GKP states, acted upon by single-mode rotation operations with angles θ j ∈ Θ followed by arbitrary single-mode squeezing operations and SUM operations. Following this, the first mode is measured and any additional Gaussian operation may be applied to the remainder modes. Intermediate phase space displacements may be included at any point in the circuit.
The second interpretation relies on the decomposition of the restricted set of allowed matrices given in Eq. (37), and on the circuit interpretation of this class that will be provided in Appendix C.3. Indeed a symplectic matrix in class A must have elements of A 1,i = 0 or B 1,i = 0 or A 1,i /B 1,i = u i /v i where u i ∈ Z and v i ∈ Z odd . We can ensure that this is satisfied by using a symplectic matrix selected from the class B of simulatable operations for multimode measurement, followed by an arbitrary symplectic matrix applied to all but the measured mode. This can be seen by analyzing the symplectic matrix given in Eq. (48), i.e., which gives the block matrices A =Ãdiag(cos θ) and B =Ãdiag(sin θ). Inspecting the elements of the first row, we have For a given column i of the block matrices, if any ofÃ 1,i or cos θ i or sin θ i is zero then the i-th condition on the elements of the first row of the symplectic matrix is satisfied trivially because it ensures that A 1,i = 0 or B 1,i = 0. Otherwise, we have which will trivially satisfy the constraint if θ i is selected from Θ in Eq. (26). A further arbitrary symplectic matrix Sp(2(n − 1), R) can be applied to all but the measured mode. When applied to the symplectic matrix given in Eq. (116), it will have the effect of the identity on the measurement mode, while transforming the coefficients of the remaining modes. The first matrix on the left hand side of Eq. (116) corresponds to the lens subgroup [38] of symplectic matrices, T . The second matrix on the left-hand side of Eq. (116) corresponds to the intersection of the general linear group GL(n, R) and the subset of Hermitian positive definite symplectic matrices Π(n) in Sp(2n, R), yielding GL(n, R) ∩ Π(n) [38]. As we show in Appendix C.3, additional squeezing operations can be added, preserving the structure of the group GL(n, R) ∩ Π(n). We also provide a more detailed analysis of these sets of operations in Appendix C.3. We can therefore represent the set of allowed circuits as those depicted in Fig. 7.
|0 GKP R (θ n )Ŝ(s n ) Figure 7: The class of simulatable operations A for single-mode measurement as understood from Eq. (116). The circuit has 0-logical GKP states as input followed by single-mode rotations with rotation angles θ j ∈ Θ. This is followed by arbitrary single-mode squeezing and then operations selected from GL(n, R) ∩ Π(n) and T respectively. The single mode is measured and arbitrary symplectic operations may be performed on the remaining modes. Arbitrary displacements may be inserted at any intermediate step of the circuit.

C Simulation of multimode measurements for Gaussian operations in B
C.1 Calculating the PDF To calculate the PDF of measuring quadraturesq 1 , . . . ,q n for the class of operations B considered in Section 4, we will again utilize the Heisenberg picture formalism to track the measurement operators and then find an expression for the PDF in terms of rotated single-mode GKP states. The PDF of the measurement in the multimode case can be evaluated analogously to in the singlemode case by inserting the identity over all rotated quadraturesẑ θ i i , The PDF can then be evaluated in terms of the Heisenberg evolved measurement operators Eq. (57) as (120) We can then use to get We again have the ∆ j specified by Eq. (99) for each mode, i.e.

C.2 Simple Example
In this section, we expand the simple example of the previous appendix to illustrate how our simulation technique works for multimode measurement. We would again like to simulate the following operationŝ acting on GKP input states. Note that the operation in Eq. (124) belongs to classes A and B simultaneously. Using our technique we can again track the Heisenberg evolution of the measurement operators, however, this time we will track both quadratures: The symplectic matrix of this operation can be represented as a rotation of the first quadrature, followed by a mode mixing operation. The rotation is parametarized by θ = (π/2, 0) for which we can build the matrix diag(cos θ) diag(sin θ) −diag(sin θ) diag(cos θ) (127) with diag(cos θ) = diag(0, 1) and diag(sin θ) = diag (1, 0). The mode mixing operation is the C X operation which is in the block diagonal form Therefore, we can evaluate the PDF, which is given by Eq. (122) as where we evaluate the parameters ∆ i according to Eq. (123) as Therefore, the PDF can be written explicitly as

C.3 Simulatable operations in multimode measurement circuits
In this subsection, we express the allowed operations in circuit class B in terms of circuit diagrams. Similarly to Appendix B.3 we can analyze the decomposition of the allowed symplectic matrix. Note that we can freely insert displacements in phase space anywhere in the circuit, as these preserve the structure of the symplectic matrix, as in Eq. (86) [37]. The class B is composed of symplectic operations which can be decomposed as where the angles θ i ∈ Θ, as defined in Eq. (26). The first matrix corresponds to the lens subgroup [38] of symplectic matrices, T . These are transformations which takê q →q (135) The matrices of the form correspond to the intersection of the general linear group GL(n, R) and the subset of Hermitian positive definite symplectic matrices Π(n) in Sp(2n, R), yielding GL(n, R) ∩ Π(n) [38]. This describes transformations of the formq whereÃ is symmetric. We are unaware of a method to decompose this set of operations into single and two-mode operations. However, we can immediately notice that combining two operations of this form in sequence will result in a new matrix which also conforms to this structure. Formally this can be expressed as which can be confirmed by analyzing the product of two arbitrary matrices whereÃ,Ã are non-singular n × n matrices: The set of single-mode squeezing operations can be expressed as where s is a vector of non-zero values of squeezing and s + is the element-wise inverse of s. We can therefore combine single-mode squeezing operations with any operation selected from GL(n, R) ∩ Π(n) while maintaining that the corresponding symplectic matrix belongs to the same set. We can therefore consider that the simulatable circuits in class B are those for which there exists a circuit of the form shown in Fig. 8. Comparing this figure with Fig. 7 makes clear the interpretation of the class of circuits corresponding to the larger set A as the same operations in B, followed by arbitrary symplectic operations on the last n − 1 modes.

D Density of the set of simulatable operations in class A
In this Appendix, we will demonstrate that the circuit class A contains symplectic matrices which have parameters which are dense in the reals. Using the decomposition in Section 2 we have shown that we can construct a symplectic matrix M with a free choice parameters A 1,i = a i ∈ R and B 1,i = b i ∈ R, for i ∈ {1, . . . , n}. However, in order to simulate the measurement of the output mode we must make |0 GKP R (θ n )Ŝ(s n )qn Figure 8: The circuit diagram of arbitrary circuits selected from the simulatable class B. The first action of the circuit on the modes is a rotation of each mode with angles θ i ∈ Θ as defined in Eq. (26), which corresponds to the third matrix in Eq. (133). This is followed by arbitrary squeezing of each mode, and can be trivially included as a result of Eq. (140). The next operation is selected from GL(n, R) ∩ Π(n) which corresponds to the second matrix in Eq. (133). The final operation is selected from T and corresponds to the third matrix in Eq. (133). Finally, the modes are measured in the position basis with homodyne detection. We note that arbitrary phase space displacements can be inserted into the circuit due to the fact that they preserve the structure of the symplectic matrix, in accordance with Eq. (86).
a further restriction on these variables. We know that for a symplectic matrix to exist in RSp(2n, R) it must have either one of A 1,i = 0 or B 1,i = 0 or A 1,i /B 1,i ∈ Q (2) . By introducing a new set which is contained by Q (2) but also contains fewer elements that Q (2) , we can create arbitrary symplectic matrices which will always be contained within RSp(2n, R). If we select the matrix elements a i , b i ∈ Q odd we can guarantee that the matrix will be contained in RSp(2n, R), since the division of one element in Q odd by another element in Q odd , will always give an element in Q (2) . Next, we can show that the set from which we choose the matrix elements Q odd is dense on the reals. Formally we can prove this following the technique of Ref. [55] to prove the density of the rational numbers on the reals. The density of a set χ on the reals formally means that for any x, y ∈ R where x < y then there exists α ∈ χ such that x < α < y.
First we recall the demonstration that the rationals are dense on the reals, following the proof given by Ref. [55]. This means that for any x < y with x, y ∈ R we must have some q ∈ Q such that x < q < y.
This means that for any x < y such that x, y ∈ R we can define α ∈ Z and β ∈ N, where N = {1, 2, 3, . . . } are the natural numbers, not including 0. Now we will follow the same proof technique to show that the set Q odd is dense on the reals.
This means that for any x < y such that x, y ∈ R we can define α ∈ Z and β ∈ N such that x < 2α + 1 2β + 1 < y (148) i.e. there exists a number r ∈ Q odd such that x < r < y. Therefore, for any two symplectic matrices which have the effect on the output mode as M :q 1 → a 1q1 + b 1p1 + · · · + a nqn + b npn (149) and M :q 1 → a 1q 1 + b 1p 1 + · · · + a nq n + b np n (150) where a i , b i , a i , b i ∈ R, it is always possible to select parametersā i ,b i ∈ Q odd such that wherebyM :q 1 →ā 1q1 +b 1p1 + · · · +ā nqn +b npn (152) is simulatable by our technique.

E The Clifford group is not contained in the set of simulatable operations
We show using a simple example that the Clifford group is not contained within RSp(2n, R). We will use the Fourier transformF which transforms the quadratures according to Eq. (6). We will also use the phase gateP 1 = e iq 2 1 /2 [2] which transforms the quadratures aŝ q 1 →q 1 p 1 →q 1 +p 1 .
However, even with this approximation, it remains challenging to evaluate the PDF of the circuit even in the case of single-mode measurement. Note that in the limit ∆ = 0 we can identify that these parameters will reach which returns the infinite-squeezing GKP state PDF given in Eq. (84), |ψ ∆→0 θ,s (x)| 2 ∝ |ϑ(ζ = −x csc θ √ πs , τ = 2 cot θ)| 2 .

F.2 Single-mode measurement
We showcase the complexity of evaluating the PDF of a single-mode measurement result by providing an example of a calculation involving two modes. The PDF is given by PDF(Q 1 = x 1 ) = dz 1 dz 2 δ (s 1 z 1 + s 2 z 2 + c − x 1 ) 0 GKP ẑ θ 1 1 = z 1 2 0 GKP ẑ θ 2 2 = z 2 2 = dz 1 dz 2 δ ((s 1 z 1 + c − x 1 )/s 2 + z 2 ) 0 GKP ẑ θ 1 Note that in the case of infinite squeezing, the PDF of the rotated and squeezed GKP state is a Dirac comb, i.e. Eq. (34), which is a summation over delta functions. When inserted into the multimode PDF, the latter yields a product of summations of delta functions in terms of z 1 , . . . , z n , which results in a function that can be integrated over these variables, leading to Eq. (45). However, in the case of finite-squeezing, the rotated and squeezed single-mode PDF cannot be expressed in terms of delta functions, which limits the possibility of analytically solving this expression. Using the fact that the PDF of the rotated and squeezed GKP state is given by Eq. (169), we have PDF(Q 1 = x 1 ) = dz 1 e z 2 1 (γ+γ * ) |ϑ(ζ = z 1 η, τ )| 2 e ((x 1 −s 1 z 1 −c)/s 2 ) 2 (γ+γ * ) |ϑ(ζ = (x 1 − s 1 z 1 − c)/s 2 η, τ )| 2 . (179) We are not aware of any theorem which rules out the possibility to evaluate such a function in general. However, we are also unable to provide a general analytic solution of this integral.