Squeezing and multimode entanglement of surface acoustic wave phonons

Exploiting multiple modes in a quantum acoustic device could enable applications in quantum information in a hardware-efficient setup, including quantum simulation in a synthetic dimension and continuous-variable quantum computing with cluster states.We develop a multimode surface acoustic wave (SAW) resonator with a superconducting quantum interference device (SQUID) integrated in one of the Bragg reflectors. The interaction with the SQUID-shunted mirror gives rise to coupling between the more than 20 accessible resonator modes. We exploit this coupling to demonstrate two-mode squeezing of SAW phonons, as well as four-mode multipartite entanglement. Our results open avenues for continuous-variable quantum computing in a compact hybrid quantum system.


I. INTRODUCTION
Quantum computation and simulation show potential for tackling difficult computational problems by leveraging superposition and entanglement in engineered quantum devices. Although individual quantum systems can be controlled with excellent precision, scaling the hardware to the complexity required while maintaining sufficient control remains a challenging problem [1]. Most architectures proposed for quantum simulation and computation [2][3][4][5][6] use one circuit component for each node in the processor, leading to demanding hardware requirements for practical applications. It is therefore attractive to explore alternative approaches to quantum computing that provide for compact encoding and processing of quantum information.
In principle, the use of continuous variables (CV) allows for realizations of measurement-based quantum computing with frequency combs, requiring only a small number of coupled quantum systems [7,8].
This paradigm of quantum computing relies on entangling a large number of modes rather than qubits, and does not face fundamental restrictions preventing universality and fault-tolerance [9]. Much experimental progress in CV encoding of quantum information has been achieved in the domain of quantum optics [10][11][12][13][14], where optical parametric oscillators can be used to generate large cluster states, a type of multipartite entangled states providing the resource for CV quantum computation. With superconducting circuits, CV encoding has been pursued mainly for error-correction schemes on logical qubits encoded in many-photon superconducting cavity states [15,16]. While parametric devices are important for lownoise amplification [17], multimode measurement-based schemes for quantum computing at microwave frequencies have received relatively little attention [18][19][20]. The dominant approach to quantum computation with superconducting quantum circuits has been the gate-based quantum processor, with most effort expended on scaling up the number of physical qubits [21].
A limiting factor for realizing CV encoding in circuit quantum electrodynamic (QED) systems is the typically large electromagnetic mode spacing, making devices with a large number of accessible modes very long or difficult to design. On the other hand, microwave frequencies are amenable to digital signal processing and thereby a greater degree of programmable control than is currently possible in optical systems. The prospect of integrating superconducting qubits as a means of providing non-Gaussian operations necessary for quantum advantage in computation [22] is an additional strength of microwave circuits.
We demonstrate an approach towards realizing CV quantum computation based on cluster state generation and control in a multi-mode hybrid superconducting quantum acoustic device. The interaction between mechanical oscillators and superconducting circuits has been used to show quantum effects [23], including entanglement [24]. Surface acoustic wave (SAW) resonators support dense mode spectra with high Q-factors (> 10 5 ) and have been used in multi-mode experiments in the quantum regime [25][26][27]. Substantial progress has also been made in recent years in the controlled generation of non-classical phononic states [28][29][30], and applications as quantum random access memories have been proposed [31].
Here, we develop a multimode quantum acoustic device by integrating a superconducting quantum interference device (SQUID) into one of the Bragg reflectors of a SAW resonator. The SQUID inductance modulates the reflectivity of a unit cell in the mirror and hence the ef- The reflector on the right hand side has 500 unit cells of fingers alternatingly connected to either the top or bottom electrode, similarly to an IDT. As shown schematically in b), the electrodes are shunted by a SQUID which is modulated by a fluxline. The IDT at the center provides a single input and output port to the resonator. The unit cell of the IDT has a period of λIDT = 736 nm and a double finger structure to suppress mechanical reflections. c) False-color scanning electron micrograph of the right hand side reflector with the SQUID and on-chip fluxline.
fective length of the resonator. Due to the narrow free spectral range, the SQUID reflector gives rise to coupling of more than 20 modes. We exploit this coupling to generate two-mode squeezed states between phonons in different SAW modes with a parametric drive. Extending the pump scheme to four tones, we demonstrate multipartite entanglement between four acoustic modes. Our results suggest this quantum acoustic platform can be used to create highly entangled multimode states for CV quantum computing.

II. DEVICE DESIGN AND SETUP
The SAW resonator, shown in Fig. 1, is defined by two reflectors with the leading edges separated by a distance of 600 µm. The reflector on the left hand side has 1200 fingers all shorted together. On the right hand side, the reflector has an interdigitated structure, where fingers are connected to either the top or bottom electrode in an alternating pattern. The top and bottom electrodes each have N p = 500 fingers with an overlap of W = 100 µm, and are connected via a SQUID. An interdigitated transducer (IDT) centered between the reflectors provides a single port to the resonator. The port IDT has 75 periods and a double-finger structure to suppress mechanical reflections [32]. An on-chip fluxline is used to apply an RF flux through the SQUID. The IDT and reflectors are fabricated from aluminium on a gallium arsenide substrate. Due to the piezoelectric coupling, the SAW field inside the resonator induces an electric potential difference between the top and bottom electrodes, generating currents through the SQUID.
Configurations where an interdigitated reflector is shunted by a variable load impedance have been used for SAW-based sensors [33]. Here, the SQUID impedance provides a means of flux tuning the SAW resonator, as well as a cross-Kerr interaction between the modes. Integrating the SQUID makes the device a kind of acoustic analogue to the superconducting cavity-based Josephson parametric amplifier [17], where the short wavelength of sound allows for a much denser mode spacing than in the purely electromagnetic case. The SQUID reflector is equivalent to a dispersively coupled nonlinear resonator, as the interdigitated fingers give rise to a large capacitance connected in parallel with the SQUID inductance. We use this model to explain the effect of parametric modulation in this system.
The mode structure of the resonator is shown in Fig. 2. While the IDT is centered with respect to the leading edge of each reflector, the broken symmetry due to the SQUID allows coupling to both odd and even modes with a free spectral range of F SR = 2.3 MHz. The alternating pattern of even and odd modes is apparent in the external and internal quality factors Q c , Q i extracted from fits to reflection measurements. The IDT couples more effectively to the even modes, resulting in a lower Q c . The frequency dependence of the IDT and mirrors provide a bandwidth of around 40 MHz around the IDT center frequency where SAW modes are overcoupled.
The narrow free spectral range of the resonator allows for simultaneous measurement of the response at multiple resonances all multiplexed in a single channel. For this measurement we use a digital microwave measurement platform [34] to directly digitally synthesize and measure signals at multiple frequencies simultaneously without analog mixers for frequency conversion.

III. COUPLED PARAMETRIC RESONATOR INTERACTION
The electromagnetic mode of the mirror has a frequency ω LC which is parametrically modulated in time. The coupling to the SAW modes gives rise to the effective Hamiltonian (Appendix D) Here b j and b † j are the ladder operators for the SAW modes, while s(t) is a time-dependent factor determined by the flux pump amplitude and frequency. Assuming a uniform vacuum coupling strength g, the effective coupling ratesg j depend on the mirror and SAW mode resonance frequencies as The bare SAW frequencies ω j are renormalized due to the interaction toω The second term in Eq. (1) contains both beamsplitter and two-mode squeezing interactions. Depending on the modulation frequency of s(t), either interaction can be selected. A beamsplitting interaction may be implemented by a parametric drive close to the difference frequency of the SAW modes. Modulating near the sum frequency will induce two-mode squeezing between pairs of modes. In our experiment we modulate the magnetic flux through the SQUID without DC flux bias at the frequency of a SAW mode. The parabolic dependence of the frequency on flux implies the electromagnetic mirror mode, and hence the effective coupling rate, are modulated at twice the pump frequency s(t)g jgk ∼ (cos 2ω p t + 1)g jgk . This gives rise to two-mode squeezing as two photons from the pump tone are converted to one phonon each in modes symmetric around the pump.

IV. TWO-MODE SQUEEZING
To observe two-mode squeezing of the SAW field, we apply a parametric pump via the on-chip fluxline at the frequency of a SAW mode f p . We measure the output field from the IDT at SAW mode frequencies symmetrically around the pump f i,± , such that 2f p − f i,− − f i,+ = 0. The frequency configuration of the measurement is illustrated in Fig. 3a. To characterize the correlations, we obtain reference histograms of the I and Q quadratures with the pump turned off. To minimize the effect of slow drift in the experimental setup, the pump output is switched on and off at a rate of 2 Hz. The output signal is amplified using a travelling-wave parametric amplifier [35]. Data are collected over approximately 7 hours (3.5 hours each with the pump on and off). A quadrature rotation is applied to the measured data to compensate for slow phase drift in the experiment.
In Fig. 3b we show subtracted quadrature histograms measured simultaneously in four pairs of SAW modes. Histograms are generated from N = 1.25 · 10 6 points measured in each mode with the pump on (off). The unsqueezed histograms obtained with the pump off are then subtracted from those produced with the pump on. We observe squeezing below the pump off level in all four mode pairs, extending the two-mode squeezing effect to SAW phonon fields. The ellipticity, defined as the ratio of the squeezed and anti-squeezed axes is well above unity for all four mode pairs and shows a diminishing trend with increased detuning. The ratio of the standard deviation in the squeezed quadrature to the pump off case, given by has a value R p < 1 across the four mode pairs. The correlation ratios are plotted as a function of pump-probe detuning in Fig. 3c. As expected, no correlations are observed outside the two-mode squeezed pairs. From our analysis of mode correlations (see appendix F 1) we estimate that the squeezing is below the vacuum level, if the modes are cooled below an effective temperature of 80 mK. Even if our device is not perfectly thermalized to the 10 mK cryostat temperature, it is unlikely the effective phonon temperature should exceed this bound, leading to strong indication that we have observed squeezing below the vacuum in the phonon field. The ability to generate two-mode squeezing with a single pump tone is an important step towards multimode entanglement, as this can be achieved using multimode squeezing [36]. In the next section we present such an experiment with a multitone pump and calibrated measurement chain.

V. MULTIMODE ENTANGLEMENT
Following the two-mode squeezing measurement, we develop the experiment further to observe multipartite entanglement involving more SAW modes. For this purpose we use another device with reduced cross-talk between the pump line and IDT described in more detail in Appendix B. With a multi-tone modulation, the Hamiltonian of Eq. 1 can provide coupling between any pair of Dashed lines indicate probe frequencies for collecting quadrature noise data. b) Histograms of two-mode correlation data where the measurement with the pump turned off has been subtracted. The label I+ (I−) denotes the I quadrature of modes at frequencies above (below) the pump frequency. The label ∆ indicates the detuning of the modes from the pump. The histograms for all two-mode quadrature combinations for the ∆ = F SR case are shown in appendix G. c) Squeezing ellipticity and ratio of standard deviations of the noise in the squeezed quadrature relative to the pump-off case. The quantities are plotted as a function of detuning of the measured mode pair from the pump. All modes show squeezing below the pump off level, and the error bars are smaller than the plot markers. modes in the resonator. Instead of single pump tone, we now apply a regularly-spaced comb of pump frequencies containing up to four tones. Due to the slight deviation from equidistance in the SAW modes the resulting entanglement is restricted to a set of modes in the vicinity of the pumps. The sharp mode structure also implies the amplitude and phase of correlations are sensitive to the pump comb settings. This is apparent in the scattering response shown in Appendix H. For the two-mode squeezing measurement data we subtract the noise measured with the pump off. In order to establish multipartite entanglement we instead perform a calibration of the gain (Appendix F 2) and added noise of the measurement amplification chain. The calibration is based on Planck spectroscopy [37] and provides an estimate of the power level corresponding to vacuum fluctuations in the SAW modes, allowing us to test the measurement data for continuous-variable entanglement.
The Gaussian state of the probe modes is characterized by the quadrature covariance matrix. Drift and noise in the measurement setup can diminish mode correlations and render averaged covariance matrices unphysical. To mitigate this problem, we divide the 2.5-minute measurement into two-second intervals and perform a reconstruction [38] to ensure a physical state and test for entanglement on each interval separately.
The pump and measurement configuration for four modes using four pump tones is shown in Fig. 4a. Figure 4b shows the quadrature covariance matrix obtained in this measurement. Using the calibration reference, the measured amplitudes are scaled with the single photon energy and measurement bandwidth ∆ BW as where A ∈ {I, Q}. With this scaling the vacuum state is given by the identity matrix. The measured mode set extends beyond the four modes analyzed here, but we are not able to recover physical covariance matrices for all modes from our calibration of the amplifier gain and added noise. The histograms of the output quadratures for the f k modes are all measured in parallel. In this case the correlations are not restricted to pairwise two-mode squeezing, but all modes are mutually correlated. As shown in Fig. 4c, the correlation features are qualitatively captured by our theoretical model presented in Appendix E.
The multimode correlated state can be analyzed for entanglement. We evaluate the entanglement using a variant of negativity of partial transpositions [39]  ) in [40]. This test relies on violating the inequality which holds for separable states.
In computing the quantity E, the covariance matrix is rotated to eliminate correlations between I and Q quadratures and V II (V QQ ) denotes the submatrix containing the I − I (Q − Q) correlations. The vectors h and g are real-valued with lengths equal to the number of modes. The subscripted terms indicate elements (and their corresponding modes) of h and g belonging to the bipartition subsets I and J . One should consider h and g as the coefficients of a general test operator acting as our entanglement witness [41]. We are free to optimize h and g to maximize any violation of the inequality Eq. (7) for a given bipartition I and J . The entanglement measure is then obtained as a weighted mean over all 75 intervals within the full integration time. We estimate the standard deviation in the measurement and express entanglement in terms of the significance Σ w , given by The uncertainty σ is calculated as where the matrix elements σ ij are obtained by error propagation accounting for the uncertainty in the calibration gain and noise parameters as well as noise in the measurement (appendix F 2). The entanglement significance is computed for all bipartitions of the four-mode set. As shown in Table I, all bipartitions yield entanglement by at least 2.4 standard deviations. Negativity of the entanglement test for all bipartitions is a signature of full multipartite entanglement [7]. We observe that the entanglement significance is substantially higher for bipartitions where the modes 3 and 4 appear in separate sets. This is due to the imperfect alignment of the equidistant pump comb with the mode structure. In the measured covariance matrix in Fig. 4b, this leads to stronger correlations involving modes 3 and 4. More uniform correlations and enhanced entanglement significance can be obtained by optimizing the pump settings.
Tailoring the digitally synthesized microwave frequency pump spectrum is also a way to obtain different entanglement structures in this setup. The square lattice is one example of a cluster state that can be used to to implement universal quantum computation and can be generated from the Hamiltonian of Eq. 1 [36]. The digital control of the amplitude and phase of each pump tone also enables extending the measurement scheme to observe multipartite entangled states with sizes approaching the number of modes in the SAW resonator. Beyond the straight-forward generation of particular target cluster states, technical (and theoretical) challenges remain to overcome errors due to the finite squeezing and achieve fault-tolerance [42].

VI. CONCLUSIONS
We have demonstrated two-mode squeezing in a surface acoustic wave resonator, likely below the phononic vacuum level. Extending this scheme to a multitone pump spectrum and calibrated measurement chain, we observed fully inseparable multipartite entanglement of four resonator modes. The dense mode structure of the resonator enables multiplexing all modes to one measurement channel without analog frequency conversion. The small on-chip footprint of our device (< 0.2 mm 2 ) further contributes to scalability.
For the measurements presented here, the pump strengths are similar to the loss rates. This limits the amount of squeezing and correlations that can be induced. To enable stronger pumping and enhance the entanglement generation, a 3-wave mixing scheme could be adopted where the pumping occurs at around twice the mode frequencies. For the same drive amplitude, this yields stronger pumping as well as less effect of saturation in the parametric amplifier. A prospect for further development is using superconducting qubits to implement non-Gaussian operations on the resonator state such as the addition or subtraction of single phonons. Non-Gaussianity is important to many applications [43] and qubit-controlled operations on the resonator state are more readily implemented in a circuit QED setting than optical experiments where nonlinearities are typically weaker. As a step in this direction, the device used to measure multimode correlations has an integrated transmon qubit, although it was not used for this experiment.
Another promising application for this device is in quantum simulation using the resonator modes as lattice sites in a synthetic dimension. Modulating the reflector SQUID at a frequency corresponding to the the free spectral range will induce nearest-neighbour hopping of phonons, giving rise to an effective lattice Hamiltonian in a hardware-efficient way.

VII. ACKNOWLEDGEMENTS
We acknowledge IARPA and Lincoln Labs for providing the TWPA used in this experiment. We are grateful to G. The vacuum coupling strength between the LC mode of the SQUID mirror and the SAW modes is given by the overlap of the zero-point voltage fluctuations of the SAW mode with the charge fluctuations on the mirror fingers [25]. As any acoustic mode that is confined in the resonator is necessarily efficiently reflected by the mirror, we make the simplifying assumption that the SAW wavelength matches the mirror period for all modes. The amplitude of the voltage zero-point fluctuations is given by The piezoelectric coefficient e 14 and the dielectric constant , as well as the substrate density ρ and SAW velocity v SAW are material parameters, while A denotes the effective mode area. The charge fluctuations across the mirror fingers are where E L = (Φ 0 /2π) 2 /L J is the characteristic inductive energy, and E C = e 2 /(2C) sets the charging energy scale. The capacitance ratio β indicates the ratio of the mirror capacitance seen by the SAW modes to the total mirror capacitance. Because the SAW field decays exponentially into the mirror with a penetration depth L p , this ratio is approximately given by β = L p /L m , where L m is the total length of the mirror. This yields an approximate coupling strength With literature values for the material parameters [44,45] and the penetration depth L P estimated from the measured free spectral range, we obtain g ≈ 2π ·1.6 MHz. This weak vacuum coupling strength implies the dispersive interaction between SAW modes and the mirror are small compared to the linewidth, and parametric excitation thus relies on strong pumping.

Appendix B: Device B
The device used for the multimode entanglement experiments of Sec. V has slight design variations. The mirror edge separation is L e = 560 µm. To enhance the coupling of the electromagnetic mirror mode to SAW, the number of finger pairs in the mirror has been reduced to N p = 275. The estimated vacuum coupling (Eq. A3) is g = 2π · 1.6 MHz.
The total capacitance is C = 3.3 pF and the SQUID shunting the two electrodes has a total critical current of I = 190 nA (L J = 1.7 nH). The associated LC mode has a frequency of ω LC = 2π · 2.1 GHz. To enable further operations on the resonator state, a transmon qubit has been integrated between the port IDT and left hand mirror. While not used in the present experiment, qubit operations could be relevant to creating non-gaussian SAW states. A microscope image of the device is shown in Fig. 5. The signal output and input have band-pass filters to suppress harmonics and spurious tones. To avoid crosstalk from the flux pump line saturating the TWPA, we use a destructively interfering cancellation tone at the pump frequency. This cancellation tone, only used in the two-mode squeezing experiment, is output from a separate channel and combined with the TWPA pump. The input line to the IDT is used for characterizing the resonator reflection.

Appendix C: Experimental details
A schematic illustrating the measurement setup is shown in Fig. 6. The digital microwave measurement platform has 8 channels of high speed digital-to-analog (DAC) and analog-to-digital converters (ADC) serviced by a large field programmable gate array (FPGA), all synchronized to one stable clock [34]. With this setup we are able to digitally synthesize drive signals and digitize response signals in the band 2-4 GHz, without analog IQ mixers. We set the sampling clock at 4 GSamples/s for the ADCs, and at 5 GSamples/s for the DACs, resulting in a Nyquist frequency of 2 GHz and 2.5 GHz respectively. The SAW cavity modes are designed to fall in the band 3.8 -3.9 GHz, within the second Nyquist zone of the converters. Tones outside the second Nyquist zone are removed with external bandpass filters.

Appendix D: Parametric coupling model
In this model we consider the coupling induced between linear SAW modes by the common interaction with the LC resonance of the mirror under parametric modulation. If we indicate the LC mode and SAW modes by the a and b j operators, respectively, the total Hamiltonian is a sum of three terms, H = H 0 + V + D(t) given by where the driving term D(t) is the time-dependent part of the Josephson energy due to the flux pump Φ(t).
If |Φ(t)| Φ 0 , we can expand the driving term D(t) to second order in Φ(t) In the case for a single flux pump, the time-dependence is described by a sinusoidal function Φ(t) = Φ AC cos(ω p t + θ). Inserting this into Eq. (D4), we arrive at the expression of the drive term for a single flux pump at ω p where d = ω LC (πΦ AC /2Φ 0 ) 2 /2 is the effective pump amplitude. For multiple flux pumps at different frequencies, Φ(t) is instead a superposition of sinusoidal functions.
In the limit of weak interaction between the electromagnetic mirror mode and SAW, g |ω i − ω LC |, we perturbatively expand the Hamiltonian by applying the Schrieffer-Wolff transformation [46] (D7) and we use the form of the drive term D(t) given by Eq. D5. The commutator [S, V ] is therefore where we introduce the effective coupling ratesg j andḡ j : We also treat the pump-dependent term e −S D(t)e S perturbatively, by expanding up to second order in S according to The commutators are calculated to be [S, [S, D(t)]] = −2d (cos (2ω p t + 2θ) + 1) where we have neglected an unimportant constant term.
To summarize, we write down the Hamiltonian H where the last term produces squeezing and beamsplitter interactions. Finally, we perform a resonance approximation by dropping all rapidly oscillating terms. This leaves us with the HamiltonianH where the final sum is only over SAW-modes k satisfying the 4-wave mixing criterion ω k = 2ω p − ω j . For multiple pumps, this would result in several 4-wave mixing criteria (one for each pump) and thus couple each SAW mode to more modes. We also define the renormalized frequenciesω LC andω j as Appendix E: Calculating the theoretical covariance matrix Here we outline how to arrive at a covariance matrix from a system of Langevin equations. Calculating the Heisenberg equations of motion for the SAW modes b using HamiltonianH Eq. (D15), we arrive at a system of equations describing multiple parametrically coupled modes. Assuming identical external couplings γ and no internal losses, we arrive at [47] b where we define the complex parametric coupling jk to be The sum is made over all modes k satisfying all 4-wave mixing criteria. Since the FSR is small compared toω j and the detuning between the SAW modes and the mirror electromagnetic mode, we make the approximationg j ≈g k which corresponds to the idealized case with identical parametric couplings. Eq. (E1) is then insteaḋ Working in the frequency domain is more convenient. The Fourier transform of Eq. (E3) is where ∆ j = Ω j −ω n + 4| | + iγ/2. If the measurement frequency Ω j =ω j − 4| |, the expression simplifies to ∆ = iγ/2. For many pumps and modes, the system of equations can be conveniently summarized by a complex weighted directed graph [48]. We consider the case for the multipartite entanglement result of Fig. 4 and draw the corresponding graph in Fig. 7. The graph illustrates the mode-couplings and can also be associated with a mode-coupling matrix M . In the basis ofb = (b 1 , ..., b 4 which provides a complete description of the frequency domain expression by −iMb = √ γb in .
The covariance matrix V can be obtained from M via the scattering matrix S, given by The scattering matrix relates the incoming modes to the outgoing modesb out = Sb in according to input-output theory [47][48][49]. The covariance matrix of the incoming noise modes V in then transforms as [7] where the subscript IQ indicates the scattering matrix has been (linearly) transformed into the quadrature basis defined by Any Gaussian state is fully characterized by V out . In the case of internal losses γ int , we need to make some minor adjustments to the preceding method. The definition of ∆ j is adjusted to be ∆ j = Ω j −ω j + 4| | + iγ tot /2 with γ tot = γ +γ int . In addition, we introduce the diagonal matrix K int = γ int I which is used to define a scattering matrix for the loss channel as In the presence of internal losses, Eq. (E8) is instead replaced by The noise coming from the internal loss port is characterized by V loss . Typically, it is assumed to be identical to the incoming noise port V in = V loss . So far in this discussion, the pump modes have been ignored. These modes form a set of correlated modes which is separate from the probe modes and can therefore be ignored in the analysis. We also note that due to the restricted probe mode set, the highest-frequency pump tone used in the measurement does not contribute to the measured correlations.
A covariance matrix calculated theoretically from Eqs. (E5-E10) is shown in Fig. 4c. We use a uniform pump strength of | | = 30 kHz and equal external and internal loss rates of γ = γ ext = 20 kHz. Although simplified, this configuration corresponds approximately to that of the multimode entanglement experiment and the theoretical covariance matrix qualitatively reproduces the features of the measured data shown in Fig. 4b.

Appendix F: Amplifier Gain and added noise
We model the effect of amplification on the covariance matrix according to [7] where T = √ G I and N = (G − 1)(2n + 1) I. Our amplifier chain is characterized by an effective amplitude gain √ G and effective added mean photon number n. The covariance matrix measured after amplification is given bỹ V , while V represents the quantum state. Thus a good estimate of √ G and n would allow us to reconstruct the quantum statistics.

Two-mode squeezed state
There are different methods to calibrate gain and added noise, which typically require some form of calibrated noise source [20] or a temperature sweep [37,50]. A method to roughly estimate the gain of the amplification chain using only our device is by measuring crosscorrelations in a two-mode squeezed state. Assuming the added noise is thermal, the cross-correlations corresponding to squeezing are independent of the added noise but not the amplifier gain. We quantify these correlations by adding the relevant elements of the covariance matrix as We measure C by applying a flux pump at 3.8732 GHz, while measuring the noise in a pair of neighbouring modes. As illustrated in Fig. 8, the pump frequency is placed directly on a SAW mode and probe frequencies are swept across neighbouring SAW modes, always keeping them strictly symmetric with respect to the pump to satisfy the 4-wave mixing criterion. In Fig. 8, C is plotted as a function of the detuning between the pump and the probe frequencies. The probe frequency sweep results in a Lorentzian-like shape of C as a function of detuning, where the peak occurs when both probes are located within their respective SAW modes.
The C lineshape is calculated by deriving the covariance matrixṼ for two coupled SAW modes, according to the method outlined in Appendix E. More specifically, the mode-coupling matrix M for two modes is graphically represented in Fig 9. The outgoing noise is then fully characterized by the 4 × 4-matrix V out , which we can find by following Eq. (E6) -(E10). Amplification is taken into account by substituting V → V out in Eq. (F1). The resulting covariance matrixṼ is used to calculate C. Given the resonance frequencies of the SAW modes along with their linewidths and assumingg 1 =g 2 , we are left with three unknown parameters: the gain G, parametric coupling and the effective phonon temperature T eff . If we fix the phonon temperature, a fit to C will give us the gain G and the parametric coupling . An example fit is shown as a solid line in Fig. 8, with an assumed phonon temperature of T eff = 30 mK. This yields an estimate of the gain to be G ≈ 80 dB, which lies within range of our expectations.
This method is not a substitute for proper gain and noise calibration procedures. However, we will use this method to make an estimate on the maximal phonon temperature possible for the two-mode squeezing in Fig.  3 to be a signature of entanglement. The procedure con-sists of essentially three steps: extract G by fitting to C, estimate the added noise N and finally reconstruct the original quantum statistics from data in Fig. 3 according to Eq. (F1).
After the gain G is extracted from fitting to C, the amplifier noise is estimated by solving for n in Eq. (F1), by replacingṼ by the pump off statistics V off while V is substituted by a thermal state with the corresponding temperature T eff . Finally, solving for n with these assumptions yields n ≈ 0.08. This added noise value is a very low estimate. Using a higher added noise value during reconstruction of V however, will result in more squeezing.
With the gain G and noise n, one can attempt reconstructing the pre-amplified covariance matrix from the two-mode squeezing data in Fig. 3. To determine whether the reconstructed state is entangled, we apply the partial positive transpose (PPT) criterion [7,39]. If the reconstructed two-mode squeezed state is labelled V TMS , the PPT criteria states that if we do a partial transposition: then a necessary and sufficient condition for separability for bipartite Gaussian states is that the matrix H =V TMS + iΩ is positive semidefinite. Ω is the symplectic matrix, defined as Thus we test for entanglement by calculating the eigenvalues λ of H and checking whether the smallest eigenvalue λ min is negative. Note that for Ω and Λ we are assuming the covariance matrix is in the basis (I 1 , Q 1 , I 2 , Q 2 ). However, the eigenvalue λ min depends on the phonon temperature, since it affects our estimate of G and n. We take this into account by calculating the value of λ min at various phonon temperatures T eff , presented in Fig. 10. Accordingly, we observe that entanglement persists for phonon temperatures up to 63 mK. This should be compared to the mixing chamber temperature of roughly 10 mK and previous experiments estimating the effective SAW phonon temperature to 37 mK [51]. Together, these observations suggest that the measured two-mode squeezing is a signature of entangled SAW modes.

Calibration and multimode state reconstruction
For the multimode entanglement experiment we perform a calibration of the gain and added noise in the amplification chain. We substitute a resistor at the mixing chamber for the device and measure the noise power We apply the PPT criterion to data presented in Fig. 3 in the main text. We vary the assumed phonon temperature, estimate G and n, and evaluate the PPT criteria at each step. The results suggest entanglement is present up to a phonon temperature of roughly 63 mK. Each line corresponds to different detunings from the pump.
as a function of temperature. Fits to the expression give the gain and noise parameters. The calibration is performed without the resonator connected at each frequency used in the measurement and the heating raises the temperature of the entire mixing chamber stage of the cryostat. Figure 11 shows the noise power at the frequency of mode f 2 (cf. Fig. 4). As the amplification chain makes use of a parametric amplifier, care must be taken to account for the idler noise in the analysis [52,53]. To accurately account for the variation in gain with frequency, the gain and noise terms appearing in Eq. F1 are extended to T = n i=1 √ G i I and N = n i=1 (G i − 1)(2n + 1)I + (G I,i − 1)(2n I,i + 1)I. The index i denotes the frequency modes and the subscript I denotes the idler contribution. To obtain reasonable fit parameters, we restrict the signal-idler gain to G I,i = G i and assume an idler noise to originate from a thermal state at temperature T I = 30 mK. From our analysis we obtain a lower than expected added noise temperature of T n ≤ 300 mK. In case our calibration underestimates the real added noise in the amplification chain, that should imply more entanglement in the reconstructed multimode state as the added noise obscures the correlations.
The fits yield uncertainties for the estimated gain and noise, which influence the entanglement significance by error propagation. The error accounting for the calibration as well as measurement error can be written as The contribution related to uncertainty in the signal gain is given by Additional contributions arise from uncertainty in the added noise σ ii,B = 2σ ni as well as the measured fluctuations in the covariance matrix Here, σṼ ij is given by the standard error of the mean of each covariance matrix element as calculated from the measured data. Furthermore, we cannot assume that errors in gain and noise are uncorrelated, which leads to the diagonal error term The covariance matrix for a physical state must satisfy the Heisenberg uncertainty relations, which may be expressed as where Ω is the symplectic matrix (cf. Eq. F5). Due to measurement noise and drift, this is not guaranteed to hold for the covariance matrix V obtained by inverting Eq. F1. To ensure a physical state before applying entanglement tests, we apply a reconstruction to find the most probable physical state V given a noisy measured state V [38]. This V is obtained by solving the optimization problem The pump configuration used in our measurement would suggest including more probe modes at higher frequency in the analysis. Including these modes renders the covariance matrix obtained unphysical by multiple standard deviations, presumably due to error in our calibration.
Appendix G: Two-mode quadrature histograms Figure 3b shows the two-mode quadrature histograms in the I + − I − plane. Here, the squeezing axis corresponds to the diagonal. In Fig. 12 we plot all two-mode quadrature histograms for the mode pair closest to the pump. The squeezing is manifest also in the Q + − Q − histogram, while the other quadratures show amplified noise.
FIG. 12. Single pump two-mode quadrature histograms. Measured for the mode pair adjacent to the pump mode (∆ = F SR). The pump off data has been subtracted. The squeezing appears in the planes combining the I or Q quadratures of both modes. Plotting the I and Q quadratures of different modes shows amplified noise, resulting in donut-shaped subtracted histograms.
The I + − I − histogram is shown without subtraction in Fig. 13a. In Fig. 13b we plot the reference histogram obtained with the pump off that is subtracted to generate the data shown in Fig. 12.

Appendix H: Scattering measurements
To verify the mode couplings induced by parametric pumping we perform scattering measurements in sample B. A signal is injected into one mode via the IDT and the scattering into other modes is measured. The scattering matrix for a single pump tone (four tones) is shown in Fig. 14 (Fig. 15). The scattering measurements verify that the parametric couplings relied on to generate entanglement are present. For an evenly-spaced pump comb the couplings are not all-to-all due to the deviation from uniform SAW mode spacing. The phase and amplitude of the scattered signal is sensitive to the pump configuration. Figure 16 shows the I and Q quadratures of a single scattering matrix element as a function of the pump spacing. The I and Q quadrature amplitudes of the S67 scattering matrix element as a function of pump tone spacing. The four-tone pump comb is uniform in frequency spacing and amplitude.