Experimental demonstration of Gaussian boson sampling with displacement

Gaussian boson sampling (GBS) is quantum sampling task in which one has to draw samples from the photon-number distribution of a large-dimensional nonclassical squeezed state of light. In an effort to make this task intractable for a classical computer, experiments building GBS machines have mainly focused on increasing the dimensionality and squeezing strength of the nonclassical light. However, no experiment has yet demonstrated the ability to displace the squeezed state in phase-space, which is generally required for practical applications of GBS. In this work, we build a GBS machine which achieves the displacement by injecting a laser beam alongside a two-mode squeezed vacuum state into a 15-mode interferometer. We focus on two new capabilities. Firstly, we use the displacement to reconstruct the multimode Gaussian state at the output of the interferometer. Our reconstruction technique is in situ and requires only three measurements settings regardless of the state dimension. Secondly, we study how the addition of classical laser light in our GBS machine affects the complexity of sampling its output photon statistics. We introduce and validate approximate semi-classical models which reduce the computational cost when a significant fraction of the detected light is classical.


I. INTRODUCTION
Several recent experiments employing large-scale quantum systems entered a complexity regime where they cannot currently be simulated on a classical computer [1][2][3][4]. These experiments reached a milestone on the path to using quantum systems for solving computational tasks of practical importance that are intractable for classical computers [5]. One of the approaches used to reach this milestone is called Gaussian boson sampling (GBS) and consists in injecting a large number of nonclassical squeezed states of light into a multiport interferometer [6]. Light at the output of the interferometer is generally in a complex entangled state owing to quantum interference. This output state is then measured using an array of single-photon detectors. The complexity of calculating the output light photon statistics scales with the number of interferometer modes and the number of detected photons [7]. Over the years, these numbers have been increasing in part thanks to improvements in the quality and brightness of the squeezed light sources [8][9][10], the efficiency and energy resolution of the detectors [11,12], as well as the development of scalable chipbased experiments [13][14][15]. The largest GBS experiment performed thus far employed 25 squeezed light sources and measured over 100 photons at the output of a 144mode interferometer [3].
Although a GBS machine is not a universal quantum computer, drawing samples from the output photon distribution has several potential applications, including calculating the vibronic spectra of molecules [16] and characterizing features of graphs [17][18][19][20]. These applications require the ability to control the squeezing * These two authors contributed equally. and displacement of the input light as well as program the interferometer transformation. Experiments have already demonstrated the ability to implement arbitrary transformations by using reconfigurable multiport interferometers [15,[21][22][23]. In fact, updating such transformations in a feedback loop based on measurement outcomes has been used for machine learning [24] and variational quantum algorithms [25]. However, no experiment has yet demonstrated the ability to displace the squeezed light in a GBS machine. Displacements can improve the graph classification accuracy of GBS [19] and are needed for calculating the vibronic spectra of real molecules [16,26,27]. Moreover, a GBS machine equipped with displacements provides a powerful quantum state engineering tool that can conditionally prepare arbitrary single-mode states [28][29][30][31].
The displacement can be achieved by interfering the nonclassical squeezed light with laser light, i.e. a coherent state. Ref. [32] experimentally demonstrated that the interference between these two states can be used for quantum-enhanced interferometry. Other experiments observed nonclassical features of displaced quantum states such as oscillations in the photon-number distribution [33,34] and micro-macro entanglement [35]. In the context of GBS, the addition of laser light provides an easy way to increase the photon detection rate and thus reduce statistical errors when sampling the output photon-number distribution. However, it also raises questions regarding the complexity of the GBS problem. If the laser is much brighter than the squeezed light, then one might be able to find a classical approximation for the output photon distribution that can be efficiently sampled using classical algorithms [36][37][38]. Thus, the GBS machine would no longer provide a quantum computational advantage. It is not well understood where this transition in complexity occurs in practice.
In this work, we build a GBS machine that samples Gaussian boson sampling problem considered. Squeezed vacua |ζi and a single coherent state |α are injected into a d-mode lossy interferometer described by a transfer matrix T . The coherent state displaces the phase-space distribution of the output state. from a 15-dimensional displaced Gaussian state. Our experiment employs a single source of squeezed vacuum and thus it can be readily simulated on a classical computer. Rather than aiming to achieve a quantum computational advantage, we explore two new capabilities enabled by the displacement. Firstly, we show that it can be used to determine the multimode Gaussian state at the output of the circuit, i.e. perform a high-dimensional quantum state reconstruction. Secondly, we study the complexity of simulating GBS machines in the presence of displacements. We introduce approximate semiclassical models that reduce the computational cost of simulations. Similar to a quantum-to-classical transition, we find that the validity of these models generally improves as we increase the displacement strength.

II. THEORY
Consider injecting squeezed vacua states |ζ i and a single-mode coherent state |α (α = |α|e iφ ) into separate ports of a d-mode lossy interferometer, as shown in Fig. 1. Light at the output of the interferometer is in a Gaussian state that can be fully described by a 2d×2d covariance matrix Σ and a displacement vector δ of length 2d: Here ν = (â 1 , ...,â d ,â † 1 , ...,â † d ) is a vector of boson annihilation and creation operators. The former quantity describes the squeezing and thermal occupation of each mode (and their correlations) after propagating the input squeezed light through the interferometer. The latter quantity provides the displacement amplitude of each mode and is determined by the evolution of the coherent state. We consider all losses to be part of the interferometer circuit, i.e. Σ and δ define the state just before being measured by ideal detectors. We use the convention that uppercase (lowercase) bold symbols are matrices (vectors). GBS experiments performed thus far have not employed displacements, i.e. δ = 0. In this case, the probability to obtain the detection pattern n = (n 1 , ..., n d ) is given by [6] pr(n) = p vac where p vac is the probability to measure vacuum in all output modes. We introduced the A matrix: with I d being the identity matrix of dimension d and Σ Q = Σ+I 2d /2 is the covariance matrix of the state's Qfunction. The submatrix A n is determined by repeating n i times the ith and (i + d)th row and column of A and thus its size (2N × 2N ) grows with the total number of photons detected, N = i n i . The quantity haf(A n ) is called the hafnian. It sums over the product of pairs of elements in A n chosen from the set of perfect matching permutations of 2N , which involves summing (2N − 1)!! terms. This exponential scaling is at the heart of the complexity of GBS.
In contrast, when δ = 0, the probability to obtain a particular detection pattern is given by the loop hafnian [39]: which contains additional terms compared to Eq. (2) as it now includes the different possible ways that the photons could have originated from both the squeezers and the coherent state. The submatrixÃ n is determined the same way as A n but its diagonal entries are replaced bỹ γ whose elements are given by repeating n i times the ith and (i + d)th elements of Expanding the loop hafnian, we find [40]: The submatrix A n−{j1,j2} is obtained by removing rows and columns numbered j 1 and j 2 from A n . The first term in Eq. (6) is the same as Eq. (2) and accounts for the probability that all N photons originated from the squeezers. In contrast, the last term contains no hafnians and accounts for the probability that all photons originated from the coherent state. The remaining terms account for the different possible ways that the N photons could have originated from both sources and the interference between these possibilities.

A. Reconstructing the multimode Gaussian state
In continuous-variable quantum state tomography [41], one can reconstruct the quantum state of an unknown signal by interfering it with a coherent state and measuring the output photon statistics. To extend this idea to multimode signals, one should interfere the coherent state with every signal mode and measure the joint photon statistics across all modes [42][43][44][45]. This is precisely what the GBS circuit in Fig. 1 achieves, which raises the question of whether it is possible to reconstruct the state at the output of the circuit directly from the measured photon statistics in this configuration. This reconstruction would provide a way to verify that the GBS machine has been properly programmed for a desired calculation.
We find that the ability to control the coherent state's phase φ can be used to determine the matrix A and the vector γ. The former quantity determines the covariance matrix [Eq. (3)] and the latter determines the displacement vector [Eq. (5)], thereby fully characterizing the multimode Gaussian state just before detection. The A matrix can be written in terms of four d × d blocks, B and C: Here B (C) is a symmetric (Hermitian) matrix describing the squeezed (thermal) part of the state, e.g. C = 0 for a pure state [6]. It suffices to measure single-photon outcomes, p j ≡ pr(0 1 , ..., 1 j , ..., 0 d )/p vac , as well as twophoton outcomes, p j,k ≡ pr(0 1 , ..., 1 j , ..., 1 k , ..., 0 d )/p vac , to determine γ, B, and C, as follows.
The single-photon probabilities measured with and without the coherent state blocked are (respectively) obtained by expanding Eq. (4): Eq. (8a) directly determines C j,j . These can then be used in Eq. (8b) to determine γ j , i.e. γ j = p j − C j,j . We assume that γ j is a real number and thus neglect the phase difference between output modes in the displacement vector. These phases do not affect the output photon statistics since photon counters are phase insensitive. Next, the two-photon probabilities measured with and without the coherent state blocked are (respectively) given by By scanning the phase φ of the coherent state, one can determine |B j,k | from the amplitude and arg(B j,k ) from the offset of the observed fringe in p j,k (i.e. the last term in Eq. (9b)). Once B j,k is determined, one can solve for |C j,k | using Eq. (9a) and subsequently Re(C j,k ) using Eq. (9b), thus only leaving an ambiguity in the sign of the imaginary part of C j,k . The threefold statistics can be measured to determine this sign, e.g. by employing an algorithm like maximum likelihood to minimize the distance between the measured threefolds and those calculated via Eq. (4). Alternatively, we show in Appendix C that one can inject the coherent state into a different input mode of the interferometer in order to determine the imaginary part of C j,k . Thus, a total of three measurement settings and the ability to scan the phase φ is required to fully determine the multimode Gaussian state in situ. These measurement settings are: (i) |α blocked, (ii) |α injected into a first input mode, (iii) |α injected into a second input mode. Since p j,k = p k,j , the reconstructed B (C) matrix is constrained to be symmetric (Hermitian), as required for a physical state. Uncertainties in the measured probabilities (e.g. from counting statistics and fitting) can be propagated in Eqs. (8) and (9) to determine the uncertainty on each matrix element. The efficiency of our reconstruction technique is partly due to the strong constraint imposed by assuming that the output state is Gaussian. While this is a natural assumption to make in the context of GBS, in practice, imperfections such as phase drifts in the displacement field can render the experimentally generated state non-Gaussian [46,47]. Provided that these imperfections are minor, the reconstructed Gaussian state provides a good approximation of the experimentally generated state in that it accurately reproduces its photon statistics, as we demonstrate further below.

B. The k-order classical approximation
The complexity of calculating pr(n) is determined by the largest hafnian in Eq. (6), whose size is determined by the number of detected photons, N [39,40,48,49]. Injecting bright coherent light into a GBS circuit is an easy way to increase the likelihood of detecting a large number of photons and thus effectively increase the N achievable in an experiment. At first glance, this does increase the difficulty of simulating the experiment because one cannot rule out the possibility that the N photons originated from the squeezers. Of course, the likelihood of this occurring depends on the relative amount of photons originating from the squeezed and coherent light, which is reflected in the weights of the different terms in Eq. (6). This leads us to considering an approximate model which ignores terms that are small when the coherent light is bright relative to the squeezers. The "k-order approximation" only keeps terms in Eq. (6) for whichγ j appears at least 2N − 2k times. Roughly speaking, this assumes that at most k photons originated from the squeezers. For example, if k = 0, we assume that all the photons came from the coherent state and only calculate which contains no hafnians, whereas for k = N , we calculate all the terms in Eq. (6). Intermediate k values reduce calculation times by ignoring terms containing the larger hafnians (see Appendix G). We test the validity of this korder approximation on our experimental results further below.

III. EXPERIMENT
Our experimental setup is shown in Fig. 2(a). A mode-locked fiber laser produces pulses of 100 fs duration at a repetition rate of 10 MHz and with a center wavelength of 1550 nm. The pulses are split into two paths. In the top path, we frequency-double the pulses in a periodically-poled lithium niobate crystal and subsequently pump a periodically-poled potassium titanyl phosphate (ppKTP) waveguide which produces degenerate two-mode squeezed vacuum through type-II spontaneous parametric down-conversion. In the bottom path, we prepare the coherent state. The three beams are coupled into polarization-maintaining single-mode fibers with efficiency η c ∼ 50% which are then coupled into a chip using grating couplers with efficiency η g ∼ 70%. The chip is made using silicon-on-insulator and contains a 15 × 15 network of directional couplers which comprise the interferometer. We discuss its characterization below. The propagation efficiency inside the 2-mm-long chip is η p ∼ 70%. Finally, the 15 output modes are detected using superconducting nanowire single-photon detectors. Since these are not photon-number-resolving detectors, we can only determine "collision-free" outcomes where n j ≤ 1 (see Appendix D). We adjust the output light's polarization using fiber polarization controllers in every mode to maximize the detection efficiency (η d ∼ 80%). The total end-to-end efficiency of the experiment is η tot = η c η 2 g η p η d ∼ 10%. The interference quality between the three beams depends on their modal purity and indistinguishability. We engineer the ppKTP source to be spectrally uncorrelated and find that the modal purity of the down-converted photons is 0.85(2) via a second-order autocorrelation measurement. Bandpass spectral filters are used to block the sinc-sidelobes from the down-converted spectra and to filter the classical beam. The temporal overlap between the three beams is controlled by two delay stages. The single-mode nature of the on-chip directional couplers ensured spatial and polarization overlap. Since the three beams are indistinguishable, one cannot discern whether photons detected at the output of the interferometer originated from the squeezer or the coherent state. Thus, the probability of detecting two or more photons depends on the relative phase between these two sources, φ. We observe φ drifting on timescales of a few seconds [orange dashed line, Fig. 2(b)] due to various instabilities in the lab. By monitoring the twofold detection rates in real time, we construct an error signal that we then use to control the voltage applied to a phase modulator in the coherent state path and lock φ to π/4 [blue line, Fig. 2(b)]. More details on the phase locking can be found in Appendix B. We measure a fringe visibility of 82(2)% for the two-photon interference signal obtained by combining pairs of photons from the squeezer and the coherent state on a balanced beam splitter (see Appendix A). This visibility provides a benchmark of the overall indistinguishability and modal purity of the three beams.
The on-chip interferometer is described by a complex 15 × 15 transfer matrix T . Each element |T ij |e iθij gives the probability amplitude |T ij | that a photon entering port i exits through port j, while θ ij is the corresponding phase. Both quantities depend on the reflectivities and phases of the directional couplers. The reflectivities are chosen to follow a Haar-random distribution while the phases are randomised due to the fabrication tolerance [14]. Since we fix the three input modes in our experiment, we only characterize the relevant 3 × 15 submatrix. The probabilities |T ij | 2 [ Fig. 2(c) top] are determined by injecting light into each input mode i, one at a time, and measuring the single-photon detection rates R ij at every output j. We then normalize the rates for each input and multiply the normalized rates by the overall efficiency of the experiment, i.e. |T ij | 2 = η tot R ij / j R ij . The phases θ ij [Fig. 2(c) bottom] are determined from the visibility of two-photon interference signals obtained by injecting two photons from the squeezer into each possible pair of inputs and recording the twofold rates. By scanning the temporal delay between the photons, we observe Hong-Ou-Mandel-type dips whose visibilities can be related to θ ij [50].

IV. RESULTS
We begin by demonstrating the characterization of the A matrix using two methods. The first "direct" method is outlined in Sec. II A and determines A directly from the measured photon statistics. The second "indirect" method calculates A by propagating a twomode squeezed vacuum and coherent state through the measured transfer matrix T using the Python libraries Strawberryfields [51] and TheWalrus [52]. For this second method, the squeezing parameter r and coherent state intensity |α| 2 were determined by estimating the average photon numbers before losses. Throughout the experiment, we fix the pump power at 1 mW and measured n PDC = 0.01 photons per pulse from the squeezer, thus r = arcsinh( n PDC /η tot ) ∼ 0.3. In Fig. 3, we show the A matrix reconstructed using the direct method with |α| 2 = 1.9. The diagonal elements are undetermined because we do not have number-resolving detectors and thus cannot measure p j,j [Eq. (9a)] or p j,j [Eq. (9b)]. Details on the state reconstruction technique are provided in Appendix C. Since the quantum state contains mostly vacuum, metrics such as the fidelity do not provide a sensitive comparison between the A matrices obtained from both methods. Instead, we calculate the output photon statistics of the two matrices using Eq. (4) and compare these to the experimentally measured statistics. We calculate pr(n i ) for all 455 N = 3 collision-free detection patterns n i and normalize the resulting distributions, i pr(n i ) = 1. The distance between the experimental and theory distributions can be computed using the total variation distance: which varies between 0 and 1. We find D = 0.033(1) (D = 0.0477 (7)) for the direct (indirect) method, thus showing that both methods correctly characterized A.
The main advantages of the direct method are that it requires only three measurement settings and it is in situ, i.e. it directly determines each element of A from the output statistics of the GBS machine. In contrast, the indirect method requires injecting single photons or coherent states into every pair of input modes to determine T , resulting in at least 2d − 1 measurement settings [50,53]. Then, one still has to measure the squeezing parameter of each input squeezed state and calculate how these transform under T in order to determine A. Next, with n PDC = 0.01 fixed, we increase the coherent state intensity such that its measured value n α = η tot |α| 2 varies from 0 to 2.2. We record the photon statistics for each value for one hour. In Fig. 4(a) and (b), we show all 1365 measured collision-free fourfold probabilities for n α = 0 and n α = 0.15, respectively. As before, we quantify the discrepancy between experiment and theory by calculating D. We plot D [red triangles, Fig. 4(c)] as a function of n α and find a mean of 0.04(3) with a maximum value of 0.123(4) occurring at n α = 0. The trend of increasing D for small n α is likely due to slight distinguishability between the downconverted modes that has a more pronounced effect when it is more probable that the photons originated from the squeezer. We observe a similar trend for the distances obtained with the twofold and threefold distributions. We find an average D of 0.030(6) and 0.030(10) with a maximum value of 0.0421(1) and 0.0477(7) (also occurring at n α = 0) for the twofold and threefold distributions, respectively.
We also use the data collected above to study the validity of various approximate models. We first consider the "classical" model devised in Ref. [38]. Its strategy is to determine the displaced squeezed thermal state having a classical quasiprobability distribution that best approximates the experimentally prepared state (see Appendix E). The resulting photon-number distribution can be efficiently sampled using classical algorithms [37]. For low n α , the classical model has a large D [black squares] and performs far worse than the quantum model. However, for large n α , we find that D for the classical model is nearly equal to the quantum model, thus showing the classical model is a valid approximation in this regime. This is expected because it is more likely to detect photons originating from the coherent state than the squeezer at large n α . The kink at α = 0 is likely an artifact of not including distinguishability in the model and thus not sampling from the optimal classical state, i.e. D can be further reduced for α > 0 (see Appendix E). Between the classical and quantum models, we can make use of semiclassical k-order approximations discussed in Sec. II B. As we increase n α , these approximations become better at modeling the data, as expected.
Increasing n α also increases the rate at which we obtain larger N detection events [ Fig. 5(a)]. With n α = 0, we measure fivefold events at a rate of roughly 10 −2 Hz, whereas this increases to 10 4 Hz with n α = 2.2. When gauging the ability of models to predict large N samples, it is more practical to use a method which does not require calculating the distance between the entire distributions as in Eq. (10). To this end, we perform a likelihood test which compares two models A and B via the likelihood ratio Suppose S = {n 1 , ..., n P } is a set of P measured samples. It follows from Bayes theorem that L < 1 occurs if pr(B|S) > pr(A|S), meaning it is more likely the samples came from the probability distribution of model B.
We fix model B to be the full quantum model without approximations and calculate L for various approximate models A. In Fig. 5(b), we plot L for a sample set S of P = 500 randomly chosen samples from the experimentally collected data containing at least four-photon detection events. We find that L of the k = 3, 4 models increases with n α , and thus, as before, the validity of these models is improving as it becomes more likely that the detected photons originated from the coherent state. We also show the trend of L with each new N = 7 sample [ Fig. 5(c),(d)]. Despite the larger number of photons detected, the k = 4 approximation still appears to be valid. However, unlike in Fig. 4(c), the k ≤ 2 approximations appear to be inadequate to model the N ≥ 4 data regardless of n α . This suggests that these approximate models cannot accurately predict higher-order correlations even at high n α .

V. CONCLUSIONS
We experimentally implemented a GBS machine that samples from a displaced nonclassical Gaussian state. We introduced and tested the validity of approximate semiclassical models that exploit the classical nature of the displacement to speed up calculations when this quan-tity is large relative to the squeezer strength. Moreover, we showed that the displacement field enables the reconstruction of the Gaussian state at the output of a GBS machine using only three measurement settings. The techniques introduced here will be useful for characterizing and validating large-scale GBS experiments. In particular, the ability to efficiently reconstruct the output Gaussian state can be used to verify that the degree of squeezing and displacement as well as the interferometer transformation has been correctly set for a desired calculation. Moreover, as with approximate models that exploit experimental imperfections in sources and detectors to speed up GBS calculations [38,[54][55][56], the k-order models we introduced exploit the classical contribution of the displacement field. These various models can be used together to better gauge the computational difficulty of sampling the output light distribution of a GBS machine.
We briefly comment on the prospect of using GBS with displacement for simulating molecular vibronic spectra. The required displacement energy varies widely depending on the molecule, and can even be significantly larger than the squeezed vacuum energy. For example, simulating the vibronic spectra of formic acid [16] (sulfur dioxide [26]) uses roughly 0.07 (0.014) photons from squeezers and a displacement of about 1.5 (1.6) photons, whereas certain transitions in tropoline requires only squeezing and no displacement [57]. Although these numbers were achievable in our setup, the covariance matrix and displacement of the output Gaussian state was fixed by the static interferometer. We provide a recipe to implement arbitrary transformations using a reconfigurable multiport interferometer in Appendix F, which could simulate many molecules in a single GBS machine. Moreover, losses will reduce the fidelity of the simulated spectra, but this can be partially mitigated by optimizing the displacement and squeezing [20,57]. Finally, we also note that GBS can inspire more efficient classical algorithms for calculating vibronic spectra [58,59]. In particular, our k-order approximations could be useful to simulate systems having a large displacement energy relative to the squeezing.

ACKNOWLEDGMENTS
We thank Renyou Ge and Xinlun Cai for fabricating the silicon chip. We also thank Jacob Bulmer and Gabriele Bressanini for their comments on the manuscript. This work is supported by the Engineering and Physical Sciences Research Council (P510257 and T001062), H2020 Marie Sklodowska-Curie Actions (846073), Samsung GRC, and the KIST Open Research Program.

Appendix A: Source characterization
We use two-photon interference between the squeezer and coherent state to obtain a benchmark of the overall quality of the indistinguishability and modal purity of the three modes [ Fig. 6(a),(b)]. We first combine the two down-converted modes on a balanced beam splitter (BS1). Because of Hong-Ou-Mandel interference, the down-converted photons bunch and thus we observe a dip in coincidences at the BS1 output of visibility V = 94(4)% [ Fig. 6(c)]. Consequently, light in the bottom output port of BS1 is approximately in a single-mode squeezed vacuum state |ζ . On BS2, we combine |ζ with a coherent state |α whose amplitude is set such that the two-photon probability is roughly equal to that of the squeezed vacuum state, i.e. | 2|ζ | 2 ≈ | 2|α | 2 . By measuring coincidence events at the output of BS2, we observe a two-photon interference signal with V = 82(2)% [ Fig. 6(d)]. For our experimental parameters (r ∼ 0.3, |α| 2 ∼ 0.3, η ∼ 0.4), we numerically calculated that the upper limit on this visibility is 94%. The ratio of our measured visibility to the ideal one is consistent with the modal purity 0.85(2) of the down-converted light determined via a second-order autocorrelation measurement.

Appendix B: Phase locking
We use twofold photon statistics to lock the phase φ between the squeezer and the coherent state. The probability to measure a photon in output mode j and k is given by Eq. (9b). The last term in this equation is an interference term that depends on φ. The visibility of this interference depends on the relative likelihood that the coherent state and squeezer each produced two pho- tons and that these photons exit the circuit in output modes j and k. We construct an error signal by heuristically choosing (i.e. those with a high interference visibility) pairs of j, k and add their respective twofold rates p j,k (φ). A subset of these rates is shown in Fig. 7(a). The rates are either correlated or anticorrelated with respect to one another depending on whether Hong-Ou-Mandel bunching or antibunching occurs, which depends on the internal phases of the interferometer. The error signal shown in Fig. 7(b) is obtained by summing these rates with the anticorrelated ones multiplied by −1. We use this error signal in a PID loop in order to control the voltage of the phase modulator and lock φ to π/4. The voltage is updated every 0.1s.

Appendix C: Gaussian state characterization
In Sec. II A, we showed that by controlling the phase of a coherent state injected in one input mode of the interferometer and measuring single-photon and two-photon detection probabilities, we can nearly fully characterize the output Gaussian state. The only missing quantity is the sign of the imaginary part of C j,k . Here we show how to determine this sign by injecting the coherent state into a second input mode.
With the coherent state injected in a first input mode, Eqs. (8b) and (9b) provide the single-photon and twophoton probabilities. In these equations, we assumed the elements of γ to be real valued, thus fixing a phase reference in the output of the interferometer. Injecting the coherent state into a different input mode, we obtain the analogous equations where the elements of µ are now complex valued. We can determine the absolute value |µ j | via Eq. (C1a) by using the already known C j,j . As before, the last term in Eq. (C1b) is an interference term leading to a fringe that can be observed by scanning the phase φ of the coherent state [ Fig. 8]. The phase offset of this fringe can be used to determine arg(µ k ) since arg(B j,k ) is already known and we are free to choose one of the output phases of µ, e.g. arg(µ 1 ) = 0. This assumes that the coherent state injected in the second input mode has the same phase φ as when injected in the first input. If instead there is an unknown offset between the two phases, one can set arg(µ 1 ) =φ and solve for this single unknown parameter by minimizing the distance for the threefold photon statistics [Eq. (10)]. The imaginary part of C j,k is then determined by where j,k = (Im(µ j )Re(µ k ) − Re(µ j )Im(µ k )) −1 and Re[µ j µ * k C j,k ] is obtained from Eq. (C1b). Fig. 8 shows an example of this procedure for a particular pair of modes, (j, k) = (5, 6). We first collect data while sweeping the phase modulator voltage [ Fig. 8(a)]. We then fit Eq. (9b) in 2π regions of the phase scan. The final fit is obtained by averaging over the five regions with the smallest fitting errors to minimize the effect of phase fluctuations and reduce Poissonian counting statistic errors [ Fig. 8(b)]. Fitting errors are propagated through Eqs. (8b) and (9b) to determine the uncertainty on the recovered matrix elements B j,k and C j,k .
In theory, the optimal displacement value maximizes the amplitude of the oscillating term in Eq. (9b). In practice, since we are using non-number-resolving detectors, we employ a weaker displacement of n α = 0.19 to reduce the effect of collisions (see Appendix D). Four output mode pairs produced near-zero twofold detection rates (i.e. about 1 per second) due to very low transmission probabilities through the interferometer. Fitting these rates with Eq. (9b) leads to a near zero |B j,k | and |C j,k |, as expected, but also leaves the phase of these matrix elements undetermined. Other effects such as instabilities in the phase φ or counting statistics errors can also hinder the fitting when the twofold rates are low. A detailed study of the robustness and limitations of our reconstruction method will be presented in a future work.
To resolve these issues here, we determine the phases that minimize the distance [Eq. (10)] for the threefold distribution using a numerical optimization algorithm. The error on the resulting phases is determined via a Monte Carlo approach: we run the optimization ten times using a different set of initial values for the phases on each run. The initial values are obtained by sampling from a Gaussian distribution of mean arg(A j,k ) and standard deviation given by the corresponding uncertainty. For the phases of the four elements that could not be retrieved with the direct inversion, we sampled from a uniform distribution between [−π, π).

Appendix D: Collisions
Our experiment employs "click" detectors that cannot resolve photon numbers. As such, events in which an output mode contained more than one photon, n j > 1, are convolved in the measured probabilities of the collisionfree events. This leads to an error in the estimate of the collision-free probabilities.
To estimate the relative size of errors caused by collisions, we calculate the probability of a collision-free event n using the loop Torontonian [60]: pr(n) = p vac × ltor(Ã n ). (D1) Unlike the loop hafnian [Eq. (4)], the loop Torontonian determines the photon statistics measured by click detectors, i.e. it convolves the probabilities of collision events with n j > 1. An implementation of Eq. (D1) and Eq. (4) can be found in the Python package TheWalrus [52]. We compute the distance D between the distributions obtained using the two equations for all collision-free N = 4 detection outcomes. We find that D increases with n α with a maximum of D = 0.024 occurring at n α = 2.2. Thus, the error caused by collisions is relatively small.

Appendix E: Classical model
The classical model, devised in Ref. [38], determines the displaced squeezed thermal state having a classical quasiprobability distribution with the highest fidelity (i.e. state overlap) to the experimentally prepared GBS state. One can then calculate its photon statistics using classical algorithms such as those presented in Ref. [37].
To determine this classical state, we follow Algorithm 1 given in the Supplementary Material of Ref. [38], which we reproduce here.
We begin by finding the classical squeezed thermal state that approximates the down-converted light given FIG. 9. Distances D of fourfold distributions obtained using different theory models. In contrast to Fig. 4(c), here the parameters of the classical model are optimized to minimize D.
the total end-to-end efficiency of our experiment, η. After losses, a squeezed vacuum state with squeezing parameter r is transformed to a squeezed thermal state whose covariance matrix is given by V = diag(a + , a − ) where a ± = ηe ±2r +(1−η). Using our experimental parameters (η ∼ 0.1, r ∼ 0.28), this covariance matrix is nonclassical since V − I 2 is not positive semidefinite. The closest classical state is a squeezed thermal state with squeezing parameter s and thermal occupation number n [61]: with s c = ln( √ a + a − ). We propagate two such squeezed thermal states and a coherent state of intensity |α| 2 through the interferometer using Strawberryfields [51]. Since our down-conversion source produces two-mode squeezed vacuum, we interfere both squeezed thermal states on a fictitious balanced beam splitter before the interferometer. For the sake of comparing different models in Fig. 4(c), we do not optimize the measured parameters |α| 2 , η, r to minimize the distance of the classical model, i.e. we use the same parameters for all models. To gauge the best possible performance of the classical model, we perform this optimization in Fig. 9. The distance is further reduced compared to Fig. 4(c) likely because distinguishability is not included in the model.

Appendix F: Arbitrary transformation
In our experiment, the output state's displacement and squeezing is fixed by the static interferometer. If one instead uses a reconfigurable multiport interferometer capable of implementing any d-dimensional unitary transformation, then a more general (d − 1)-dimensional Gaussian state can be prepared using the recipe shown in Fig. 10. Such a multiport interferometer contains at least d(d − 1)/2 tunable beam splitters [62]. The displacement operationD(δ j ) is achieved by combining each output mode j with the coherent state on a beam splitter of low reflectivity R j 1 and phase shift φ j [63]. This leaves (d − 1)(d − 2)/2 beam splitters for the squeezers, which can be used for an arbitrary (d − 1)-dimensional unitary. We also note that setting T to the identity and measuring coincidences between detector D d and each D j can be used to lock the phases of every squeezed vacuum to the coherent state using a procedure similar to that presented in Appendix B.
Appendix G: k-order run time In Fig. 11, we plot the run time of a typical calculation of the loop hafnian [Eq. (2)] and the k-order approximation which truncates this quantity at a certain k (see Sec. II B). The calculations are performed on a desktop machine with a 16-core 2.9 GHz CPU and 16 GB of memory. The loop hafnian is calculated using the TheWalrus [52] whereas the k-order approximation uses our own code (available upon request). While the former code has been well optimized [56], our k-order approximation implementation is likely not optimal and we anticipate that its run time can be further improved.