Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Generation and Structuring of Multipartite Entanglement in a Josephson Parametric System

Generation and Structuring of Multipartite Entanglement in a Josephson Parametric System IntroductionEntanglement is a crucial resource in advanced information processing based on quantum mechanical concepts.[1,2] To exceed the computational power of classical devices, quantum devices typically need to employ highly entangled states. Controlled generation of entanglement is the key resource not only in quantum computing,[3–8] but also in sensing[9–12] and secure communications.[13–15]One of the most easily accessible and, at the same time, reliable sources for entanglement generation is the vacuum state of a quantum field.[16] Squeezing, which is the fundamental operation for the continuous variable (CV) states production, allows generation of coherence and entanglement from vacuum fluctuations.[17] While two‐mode squeezing produces bipartite quantum states, multipartite states can be generated by applying similar operations.[18] Intriguingly, multipartite CV states are shown to enable various promising phenomena such as quantum state sharing[19] and secret sharing,[20] dense coding,[21] error correction,[22] and quantum teleportation.[23] Alongside with phase sensing[24] and quantum sensor networks,[25] multipartite entangled states have significant potential in multiparameter quantum metrology applications.[25,26] Furthermore, CV cluster states show potential as a universal quantum computing platform,[27–29] which has been under active development for the past 20 years. Cluster state calculus, foremost utilizing optical resources,[28,30–33] realize measurement‐based quantum computing.[3]While optical‐mode schemes for generation of multipartite states lack versatility, in‐situ tunability and are limited to optical frequencies, the microwave platform allows for full control of operations via input RF‐signals and integration with the existing silicon‐based circuitry. During the past few years, significant progress in processing of CV multipartite states at microwaves has been achieved; for instance, squeezed states produced by microwave cavities have been shown to exhibit correlations between photons in separate frequency bands[18] and strong entanglement between different modes.[34]In this work, we experimentally generate genuinely entangled tripartite and quadripartite states using a superconducting parametric cavity, operated under steady‐state conditions. Using the Gaussian‐mode formalism,[35–37] we characterize the generated states and verify entanglement from the covariance matrix. We develop an analytical description, which allows us to determine the entanglement structure of the generated state and establish a connection to H(amiltonian)‐graph representation.[7,29,38,39] All of the experimental results are in good agreement with the theoretical predictions, in which all circuit parameters were set in accordance with the measured characteristics of the device.The paper is organized as follows. Section 2 describes the quantum dynamics of a Josephson parametric system and demonstrates the generation protocol for CV multipartite states, such as fully inseparable and genuinely entangled states. Using analytical methods, we provide the entanglement structure, which is described by graphs and their corresponding adjacency matrices. In Section 3 we present the experimental setup and explain our data analysis methods. In Section 4 we present the experimental and theoretical results on the generation of multipartite entanglement using a Josephson parametric system in both tripartite and quadripartite cases. In Appendices A.1–A.5 we present details of our analysis, experimental techniques, and parameter extraction. Appendix A.6 deals with analytical solution of the equations of motion for tri‐/quadripartite states using simplified linear model of parametric amplifier and establish correlations (graph connections) between spectral modes, which allows us to classify and control the entanglement structure. Furthermore, analytical forms of relevant covariance matrices are given.Theoretical FoundationsQuantum Dynamics of the DeviceOur work employs a parametric system comprised of a superconducting λ/4$\lambda /4$ resonator terminated in a SQUID loop. Such a setup forms the archetype of a narrow‐band superconducting Josephson parametric amplifier (JPA).[40,41] In our setting, we pump the JPA using multi‐tone external RF magnetic flux through the SQUID at frequencies that are approximately twice the frequency of the resonator ωd∼2ωr$\omega _d \sim 2 \omega _\text{r}$ (three‐wave mixing).[18,42] In the rotating frame, the Hamiltonian of the system, as derived in refs. [40, 43], is given by:1Hsys,rwa(t)=ℏΔra∼†a∼+ℏ2∑d=1p(αd∗eiΔdta∼2+αde−iΔdta∼†2)+6ℏKa∼†a∼†a∼a∼$$\begin{eqnarray} {H}_{\text{sys,rwa}}(t) &=& {\hbar }\Delta _r \tilde{a}^\dagger \tilde{a} +\nonumber \\ &&\frac{\hbar }{2}\sum _{d=1}^{p}(\alpha ^{*}_d e^{i\Delta _d t}\tilde{a}^{2} + \alpha _d e^{-i\Delta _d t} \tilde{a}^{\dagger 2}) + 6{\hbar }{K} \tilde{a}^{\dagger }\tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$where a∼$\tilde{a}$ (a∼†$\tilde{a}^{\dagger }$) is the annihilation (creation) operator for cavity photons in the rotating frame at angular frequency ωΣ/2$\omega _{\Sigma }/2$, αd$\alpha _d$ is the complex amplitude of the d‐th pump tone and Δd=ωd−ωΣ$\Delta _d=\omega _d-\omega _{\Sigma }$ is the angular frequency detuning of the corresponding tone. Possible extra phase factors in different pump tones are included in the complex pump amplitude αd=|α|eiφd$\alpha _d=|\alpha | e^{i\varphi _d}$. Here Δr$\Delta _\text{r}$ denotes the detuning between half of the average pump angular frequency ωΣ$\omega _\Sigma$ and the resonator angular frequency ωr$\omega _\text{r}$: Δr=ωr−ωΣ/2$\Delta _\text{r}=\omega _r-\omega _{\Sigma }/2$, with ωΣ$\omega _{\Sigma }$ representing the average angular frequency in multi‐tone driven case: ωΣ=(1/p)∑d=1pωd$\omega _{\Sigma }=(1/p)\sum _{d=1}^{p}\omega _d$ with d={1,⋯,p}$d=\lbrace 1,\dots ,p \rbrace$ as the pump tone index.Strongly driven SQUIDs are notoriously nonlinear. Therefore, we also include the nonlinear Kerr term with strength K to the description of our parametric system. The Kerr constant controls the parametric behavior close to and above the critical pumping threshold α≥αcrit$\alpha \ge \alpha _{\text{crit}}$. Several effects are accounted for by the Kerr nonlinearity, such as limited maximum gain, compression, observed at α≲αcrit$\alpha \lesssim \alpha _{\text{crit}}$, broadening and shifting of the resonance curve, and parametric oscillation above the critical point. In our experiment, we employ the pump‐power‐dependent gain coefficient to extract the Kerr constant (see Appendix A.4).In order to describe the coupling of the cavity resonator to the incoming transmission line and to an intrinsic thermal bath, we include two additional terms in the full Hamiltonian:2H(t)=Hsys,rwa(t)+Hsig+Hloss$$\begin{equation} {H}(t)=H_{\text{sys,rwa}}(t)+H_{\text{sig}}+H_{\text{loss}} \end{equation}$$where Hsig$H_{\text{sig}}$ includes the coupling to the signal port transmission line with dissipation rate κ, while Hloss$H_{\text{loss}}$ includes the coupling to the internal loss port with linear dissipation rate γ.Using the quantum Langevin equation (QLE), we obtain the output modes in our parametric system. We employ the standard input/output formalism in the rotating frame, which yields3a∼̇(t)=(−iΔr−κ+γ2)a∼−i∑d=1pαdeiΔdta∼†+κb∼in+γc∼in−12iKa∼†a∼a∼$$\begin{eqnarray} \dot{\tilde{a}}(t)=(-i\Delta _r &-& \frac{\kappa +\gamma }{2})\tilde{a} - i \sum _{d=1}^{p}\alpha _d e^{i\Delta _d t} \tilde{a}^\dagger \nonumber \\ && +\ \sqrt {\kappa }\tilde{b}_{\text{in}}+\sqrt {\gamma }\tilde{c}_{\text{in}} - 12i{K} \tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$where b∼in$\tilde{b}_{\text{in}}$ and c∼in$\tilde{c}_{\text{in}}$ are the ladder (annihilation) operators for the signal and linear dissipation ports, respectively.The output mode b∼out$\tilde{b}_{\text{out}}$ is obtained using the following relation between the incoming and outgoing modes:4b∼out(t)=b∼in(t)−κa∼(t)$$\begin{equation} \tilde{b}_{out}(t)=\tilde{b}_{\text{in}}(t)-\sqrt {\kappa }\tilde{a}(t) \end{equation}$$We are interested in the correlations embedded in the output mode given by Equation (4) in the time domain. The correlations can be revealed in full after Fourier transformation to the frequency domain. By defining finite‐band spectral modes (see Section 2.2) in the frequency domain and examining correlations between these spectral ranges, we can verify the presence of entanglement in the band‐limited microwave signals.Spectral Modes DefinitionParametric downconversion processes and the definition of employed spectral modes within the fundamental cavity resonance in a multi‐pump JPA are illustrated in Figure 1. The spacing and width of the spectral modes are selected in such a manner that that the modes are generated within the linewidth of the JPA resonance. Each pump that acts on the JPA triggers spontaneous parametric downconversion of a pump photon (vertical red arrows) into two photons, with their energies summing up to the energy of the pump photon (blue arrows). This process is stimulated by vacuum fluctuations, whose existence is a fundamental feature of quantum electrodynamics. One might expect that the downconversion processes would be random, occurring independently for each pump in the multi‐tone pumping situation. In such a case, the result would simply be a sum of the downconversion processes, but this turns out not to be the case. Instead, the photons are fundamentally correlated, even if they originate from different pump tones, because they were “born into existence” by the same quantum fluctuation. In other words, one spectral mode contains photons correlated with the quanta in the other spectral modes and, consequently, we expect multipartite correlations to appear (depicted schematically via the zigzag lines).1FigureDefinition of spectral modes and their correlations in a multi‐tone pump setting. Spectral modes in a multi‐pump JPA, where the pump tones (red arrows) trigger parametric downconversion process (PDC) leading to the appearance of multipartite correlation between microwaves (blue arrows), extracted from vacuum fluctuations. Numbered spectral modes depicted in green are also correlated due to the continuous pumping of the JPA resulting in multipartite entanglement between microwaves in the spectral modes. The bandwidth of each spectral mode Δ in Fourier analysis is chosen to be much narrower than the cavity resonance width.We consider the fundamental resonance of a transmission line JPA centered at ωr$\omega _r$ frequency with a bandwidth 2δω$2\delta \omega$, within which ([ωr−δω,ωr+δω]$[\omega _\text{r}-\delta \omega , \omega _\text{r}+\delta \omega ]$) we define N spectral modes as depicted in Figure 1. Let us define a∼$\tilde{a}$ as a vector of spectral modes:5a∼={a∼1,⋯,a∼N,a∼1†,⋯,a∼N†}T$$\begin{equation} \tilde{a}=\lbrace \tilde{a}_1,\dots ,\tilde{a}_N,\tilde{a}^\dagger _1,\dots ,\tilde{a}^\dagger _N\rbrace ^T \end{equation}$$where N is a total mode number and the creation a∼i†=a∼i†(t)$\tilde{a}_i^{\dagger }=\tilde{a}_i^{\dagger }(t)$ and annihilation a∼i=a∼i(t)$\tilde{a}_i=\tilde{a}_i(t)$ operators are time‐dependent. In general, the frequency difference between half of pth and (p+1)$(p+1)$th pump tones defines the bandwidth of the spectral mode a∼i$\tilde{a}_i$. We employ an equidistant pump scheme where bandwidth of each mode a∼i$\tilde{a}_i$ is defined as Δ such that the spectrum of a full set of modes a∼i$\tilde{a}_i$ covers the bandwidth 2δω$2\delta \omega$ of the cavity mode a∼$\tilde{a}$. In the experiment, we collect the emitted power over the whole [ωr−δω,ωr+δω]$[\omega _\text{r}-\delta \omega , \omega _\text{r}+\delta \omega ]$ frequency range and separate signals in the N modes using numerical postprocessing. The same operation could be implemented by using accurate bandpass filters with bandwidth Δ.Within the scope of this work, we consider only tripartite (N=3$N=3$, p=2$p=2$) and quadripartite (N=4$N=4$, p=3$p=3$) quantum states. In the following, we elucidate the internal structure of the generated states through a graph representation based on the quantum Langevin equation.Graphical Description of Quantum StatesIn order to construct a comprehensive graphical representation, we consider interactions between cavity and input vacuum modes by solving the QLE given in Equation (3) for N cavity modes defined in Equation (5); for details see Appendix A.6. Assuming the strong coupling regime with κ≫Δ$\kappa \gg \Delta$ while neglecting dissipation losses γ and the Kerr‐nonlinearity K, Fourier transformation yields a system of linear equations which can be cast in a matrix form:6Ma∼(ω)=κb∼in(ω)$$\begin{equation} \mathbf {M}\tilde{a}(\omega )=\sqrt {\kappa }\tilde{b}_{\text{in}}(\omega ) \end{equation}$$Here the interaction matrix M$\mathbf {M}$ contains diagonal entries provided by cavity‐related part of the Hamiltonian, whereas the parametric terms appear in the off‐diagonal entries. For example, one obtains for the tripartite case with Δr=0$\Delta _r=0$ and having different phases for the pump tones α1=αeiφ1,α2=αeiφ2$\alpha _1=\alpha e^{i\varphi _1}, \alpha _2 =\alpha e^{i\varphi _2}$:7M=c000iαeiφ100c0iαeiφ10iαeiφ200c0iαeiφ200−iα†e−iφ10c00−iα†e−iφ10−iα†e−iφ20c00−iα†e−iφ2000c$$\begin{eqnarray} \mathbf {M}= \def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & i\alpha e^{i\varphi _1} & 0\\ 0 & c & 0 & i\alpha e^{i\varphi _1} & 0 & i\alpha e^{i\varphi _2} \\ 0 & 0 & c & 0 & i\alpha e^{i\varphi _2} & 0\\ 0 & -i\alpha ^\dagger e^{-i\varphi _1} & 0 & c & 0 & 0\\ -i\alpha ^\dagger e^{-i\varphi _1} & 0 & -i\alpha ^\dagger e^{-i\varphi _2} & 0 & c & 0\\ 0 & -i\alpha ^\dagger e^{-i\varphi _2} & 0 & 0 & 0 & c \end{bmatrix} \nonumber\hspace*{-10pt}\\ \end{eqnarray}$$where the frequency dependency enters through the c=−iω+κ/2$c=-i\omega +\kappa /2$ coefficient. To express the intracavity modes through the vacuum input, we use the inverse matrix M−1$\mathbf {M}^{-1}$:8a∼(ω)=κM−1b∼in(ω)$$\begin{equation} \tilde{a}(\omega )=\sqrt {\kappa }\mathbf {M}^{-1}\tilde{b}_{\text{in}}(\omega ) \end{equation}$$For the tripartite case, such a matrix is written in the following way9M−1=1c2−2α2·c−α2c0α2eiΔφc0−iαeiφ100c0−iαeiφ10−iαeiφ2α2eiΔφc0c−α2c0−iαeiφ200iαe−iφ10c−α2c0α2eiΔφciαe−iφ10iαe−iφ20c00iαe−iφ20α2eiΔφc0c−α2c$$\begin{eqnarray} \mathbf {M}^{-1} &=& \frac{1}{c^2-2\alpha ^2} \nonumber \cdot \\ &&\def\eqcellsep{&}\begin{bmatrix} c-\frac{\alpha ^2}{c} & 0 & \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}} & 0 & -i\alpha e^{i\varphi _1}& 0\\ 0 & c & 0 & -i\alpha e^{i\varphi _1} & 0 & -i\alpha e^{i\varphi _2}\\ \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}} & 0 & c-\frac{\alpha ^2}{c} & 0 & -i\alpha e^{i\varphi _2} & 0\\ 0 & i\alpha e^{-i\varphi _1} & 0 & c-\frac{\alpha ^2}{c} & 0 & \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}}\\ i\alpha e^{-i\varphi _1} & 0 & i\alpha e^{-i\varphi _2} & 0 & c & 0\\ 0 & i\alpha e^{-i\varphi _2} & 0 & \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}} & 0 & c-\frac{\alpha ^2}{c} \end{bmatrix}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$where Δφ=φ1−φ2$\Delta \varphi =\varphi _1-\varphi _2$. Besides two mode squeezing (TMS) correlations proportional to α, the matrix contains also beamsplitter correlations (BS) ∝α2$\propto \alpha ^2$ (denoted in bold). Note that the α2‐terms are absent in the matrix M$\mathbf {M}$ introduced in Equation (6). In this tripartite example, the phase difference Δφ$\Delta \varphi$ contributes only to the BS connection a∼i(ω)↔a∼j(ω)$\tilde{a}_i(\omega ) \leftrightarrow \tilde{a}_j(\omega )$ or a∼i†(ω)↔a∼j†(ω)$\tilde{a}^\dagger _i(\omega ) \leftrightarrow \tilde{a}^\dagger _j(\omega )$. TMS connections are defined by entries a∼i(ω)↔a∼j†(ω)$\tilde{a}_i(\omega ) \leftrightarrow \tilde{a}^\dagger _j(\omega )$ in M−1$\mathbf {M}^{-1}$. In Appendix A.6 we discuss how the phase shifts between pumps influence the structure of subspaces within the covariance matrix. In the quadripartite case, it turns out that those products of M−1$\mathbf {M}^{-1}$ responsible for beamsplitter interaction can be suppressed fully by choosing pump phases properly in certain pump tone configurations.Interestingly, the matrix κM−1$\sqrt {\kappa }\mathbf {M}^{-1}$ can also be interpreted as an adjacency matrix, which is used in graph theory for characterization of the connections, the graph edges. In regular H(amiltonian)‐graph,[7,29,38,39] each vertex represents a vacuum mode and the adjacency matrix describes correlations produced by two mode squeezing (TMS) between the vacuum modes. However, our scheme produces additional correlations, BS correlations, that can no longer be described purely by the TMS correlations and, therefore, the standard H‐graph theory needs a more generalized approach.In our approach, we introduce generalized H∼$\mathcal {\tilde{H}}$‐graphs formed by both two mode squeezing and beamsplitter correlations (see Appendix A.6). Examples of such graphs are presented in Figure 2a,b for the tripartite and quadripartite case, respectively. Intriguingly, the famous Greenberger–Horne–Zeilinger (GHZ) states[7,38,44] have different structure since they are devoid of BS correlations. However, by applying additional pump tones and adjusting phase difference between them, one can generate tripartite and quadripartite GHZ states consisting of only SQ correlations as shown in Figure 2c,d. In general, the considered multipumping scheme allows us to control SQ and BS bonds providing access to more complex structures of CV quantum states beyond GHZ‐like states.2FigureGraph representation for the entangled tripartite and quadripartite systems. Generalized H‐graph (H∼$\mathcal {\tilde{H}}$‐graph) representation of a) bisqueezed tripartite CV entangled state and b) quadripartite CV entangled H∼$\mathcal {\tilde{H}}$‐graph state obtained in our experiments. Vacuum modes (red circles) are connected via two‐mode squeezing (TMS, solid lines on graph) and beamsplitter (BS, dashed lines) correlation. Graphs (c,d) represent tri‐ and quadripartite GHZ states, which can be obtained via introducing an additional pump tone with Δd=0$\Delta _d=0$ in the tripartite case and two additional tones at −Δ$-\Delta$ and Δ in the quadripartite setting. The additional pumps supply missing TMS connections to the entangled states. For details, see Appendix A.6.From the experimental point of view, the observer is interested in the adjacency matrix for out‐coming modes b∼out$\tilde{b}_\text{out}$, which are obtained from a∼$\tilde{a}$ using the input–output relationship in Equation (4). This equation yields10b∼out(ω)=(I−κM−1)b∼in(ω)$$\begin{equation} \tilde{b}_{\text{out}}(\omega )=(\mathbf {I}-\kappa \mathbf {M}^{-1})\tilde{b}_{\text{in}}(\omega ) \end{equation}$$on the basis of which we may define the adjacency matrix M∼=I−κM−1$\mathbf {\tilde{M}}= \mathbf {I}-\kappa \mathbf {M}^{-1}$ for input–output mode graphs. Due to the linear nature of the equation, the unit matrix does not change the intrinsic form of interactions between the vacuum modes, but the correlation structure of intracavity and output spectral modes is equivalent. Consequently, analyzing a graph defined by the matrix M−1$\mathbf {M}^{-1}$ is sufficient to characterize the BS and TMS connections between vacuum spectral modes.Connection to Hamiltonian GraphCV cluster states with square‐lattice graph structure provides a foundation for measurement‐based continuous variable quantum computation (CVQC).[38] Cluster states can be asymptotically reached from H‐graph states in the case of infinite squeezing.[39] The H‐graph structure is defined by its adjacency matrix G$\mathbf {G}$, whose entries Gij$G_{ij}$ specify the multimode squeezing Hamiltonian as:11HS=ℏα∑i,jGij(a∼i†a∼j†+a∼ia∼j)$$\begin{equation} H_\text{S}= \hbar \alpha \sum _{i,j} G_{ij}(\tilde{a}_i^\dagger \tilde{a}_j^\dagger + \tilde{a}_i \tilde{a}_j) \end{equation}$$Here, the pump tone amplitudes are considered to have equal strength α. The matrix G$\mathbf {G}$ involves TMS correlations between modes a∼i↔a∼j†$\tilde{a}_i\leftrightarrow \tilde{a}^\dagger _j$, but as was pointed out before, the BS correlations do not show up in the H‐graph representation. The equations of motion for the operators are given by ia∼̇k†=α∑jGjka∼j$i\dot{\tilde{a}}_k^\dag =\alpha \sum \nolimits _j {{G_{jk}}} {\tilde{a}_j}$ and ia∼̇k=−α∑jGjk∗a∼j†$i{\dot{\tilde{a}}_k} = -\alpha \sum \nolimits _j {G_{jk}^*} \tilde{a}_j^\dag$. Taking the Fourier transform, the left hand side equals ω×a∼(ω)$ \omega \times \tilde{a}(\omega )$, and the combination ωa∼k(ω)=α∑jGjka∼j†(−ω)$\omega \tilde{a}_k(\omega )=\alpha \sum \nolimits _j {{G_{jk}}} {\tilde{a}_j}^\dag (-\omega )$ provides the connection to the QLE treatment in Equation (6): a∼(ω)$\tilde{a}(\omega )$ here is the cavity signal defined by the graph connections given by Gjk$G_{jk}$. Consequently, the basic graph structures are the same, but the form of M−1$\mathbf {M}^{-1}$ in the QLE analysis yields higher order correlations which are experimentally relevant.A standard description for graphs is based on the complex symmetric matrix Z=ie−G$\mathbf {Z}=ie^{-\mathbf {G}}$, which is interpreted as the adjacency matrix for an undirected Gaussian graph with complex‐valued edge weights.[38,39] Decomposing such a matrix up to quadratic terms Z=iI−iG+i2G2+O(h3)$\mathbf {Z}=i\mathbf {I}-i\mathbf {G}+\frac{i}{2}\mathbf {G}^2+\mathcal {O}(h^3)$, we obtain corrections to the adjacency matrix, which correspond to additional correlations, the BS correlations, obtained in our QLE analysis. Indeed, BS transformations embody interactions to second‐order, which provides classical correlations between corresponding nodes.[45]Let us now show the origin of BS correlations using the multimode squeezing Hamiltonian of Equation (11) in R=e−iℏHτ$R = e^{-\frac{i}{\hbar } H \tau }$ where we have considered that the system is pumped for a finite time τ. The multimode squeezing operator R can be decomposed to a combination of TMS operators, containing Bij=a∼i†a∼j†+a∼ia∼j$B_{ij}=\tilde{a}_i^\dagger \tilde{a}_j^\dagger +\tilde{a}_i\tilde{a}_j$, and BS transformations based on Tij=a∼ia∼j†−a∼i†a∼j$T_{ij}=\tilde{a}_i \tilde{a}_j^\dagger -\tilde{a}_i^\dagger \tilde{a}_j$. By utilizing the Zassenhaus expansion (up to first order of commutation relationship),[45] we obtain for the tripartite squeezing operator12R=e−iατ∑i,j=13Gij(a∼i†a∼j†+a∼ia∼j)=e−iατ(B12+B23)=e−iατB12e−iατB23eα2τ22[B12,B23]=e−iατB12e−iατB23eθ13T13$$\begin{eqnarray} R &=& e^{-i\alpha \tau \sum ^3_{i,j=1} G_{ij}(\tilde{a}_i^\dagger \tilde{a}_j^\dagger + \tilde{a}_i \tilde{a}_j)}=e^{-i\alpha \tau (B_{12}+ B_{23})}=e^{-i\alpha \tau B_{12}}\nonumber \\ && e^{-i\alpha \tau B_{23}}e^{\frac{\alpha ^2 \tau ^2 }{2}[B_{12},B_{23}]}= e^{-i\alpha \tau B_{12}}e^{-i\alpha \tau B_{23}}e^{\theta _{13}T_{13}} \end{eqnarray}$$Here, θ13 specifies the relative phase shift between the two pump tones. For detailed information on the expansion coefficients for a bisqueezed state we refer the reader to ref. [45].The total multimode squeezing operation can be considered as a combination of TMS operators, acting on the respective bipartitions, and BS transformations between the other modes. The beamsplitter correlations are phase dependent, and the strength of the BS contribution can be tuned down to zero in certain cases via proper choice of the phase difference between the pumps.The general decomposition of the multimode squeezing operator can be expressed as13R=∏1NTMSe−iατB12e−iατB23⋯∏1NBSeθ13T13eθ24T24⋯$$\begin{equation} R= \prod _{1}^{N_\text{TMS}}e^{-i\alpha \tau B_{12}}e^{-i\alpha \tau B_{23}}\dots \prod _{1}^{N_\text{BS}}e^{\theta _{13} T_{13}}e^{\theta _{24}T_{24}}\dots \end{equation}$$where the total number of TMS and BS operators in the decomposition is NBS=∑nN−2n$ N_{\text{BS}}=\sum _{n} N-2n$ and NTMS=∑nN−2n+1$N_{\text{TMS}}=\sum _{n} N-2n+1$; where n=1,2,⋯,⌊N2⌋$n=1,2,\dots ,\lfloor \frac{N}{2}\rfloor$ for the configuration introduced in Figure 1 with N−1$N-1$ pump tones and N spectral modes. Collecting all of the entries of Bij$B_{ij}$ and Tij$T_{ij}$, we obtain a generalized adjacency matrix G∼$\mathbf {\tilde{G}}$ with entries G∼ij=∑1NTMSBij+∑1NBSTij$\tilde{G}_{ij}=\sum _1^{N_{\text{TMS}}}B_{ij}+\sum _1^{N_{\text{BS}}}T_{ij}$. Thus, the beamsplitter correlations in adjacency matrix arise naturally from the squeezing operator formalism when the Hamiltonian is supplied with the second order terms in pump amplitude. The structure of the matrix G∼$\mathbf {\tilde{G}}$ defines the edge connections in the generalized H∼$\mathcal {\tilde{H}}$‐graph.Verification of the Multipartite EntanglementThe generalized graph analysis allows us to visualize the structure of entanglement in the quantum state generated by simultaneous multiple pump tones. However, in order to estimate the amount of quantum resources embedded in the state, we have to investigate and quantify the classical and quantum correlations and determine how they reflect the genuine multipartite entanglement of the state.Within the framework of parametric amplifiers, all microwave fields produced by a JPA below the critical threshold are Gaussian.[17,46] Therefore, the output states of a N‐mode JPA can be fully characterized by its covariance matrix of 2N‐length column vector with quadratures r∼=(x∼1,p∼1,⋯x∼N,p∼N)T$\tilde{r}\nobreakspace =\nobreakspace (\tilde{x}_1, \tilde{p}_1, \dots \tilde{x}_N, \tilde{p}_N )^{T}$, where x∼i=(a∼i+a∼i†)/2$ \tilde{x}_i = (\tilde{a}_i+\tilde{a}^{\dagger }_i)/2$ and p∼i=(a∼i−a∼i†)/2i$\tilde{p}_i=(\tilde{a}_i-\tilde{a}^{\dagger }_i)/2i$. The covariance matrix V$\mathcal {V}$, whose elements are given by14Vij=2Δr∼iΔr∼j+Δr∼jΔr∼i−4Δr∼iΔr∼j$$\begin{equation} \mathcal {V}_{ij}=2\left&lt;\Delta \tilde{r}_i \Delta \tilde{r}_j + \Delta \tilde{r}_j \Delta \tilde{r}_i \right&gt; - 4\left&lt;\Delta \tilde{r}_i \right&gt; \left&lt;\Delta \tilde{r}_j\right&gt; \end{equation}$$is sufficient for detection of the entanglement, eliminating the need for analysis of the full density matrix. The last term can be ignored as we take Δr∼i=r∼i−<r∼i>$\Delta \tilde{r}_i=\tilde{r}_i-&lt;\tilde{r}_i&gt;$. Obtaining the covariance matrix, we can analyze the entanglement and examine the structure of the quantum state. In this work, we consider fully inseparable states and genuinely entangled states, for which the covariance‐based detection is, in general, more robust than detection via complete determination of the state.[47] While the covariance matrix is sufficient for evaluating entanglement of Gaussian states, it is necessary to include higher‐order correlations in the evaluation of non‐Gaussian states.[48,49]To examine inseparability properties of the quantum state,[47,50–55] we apply symplectic transformations to the covariance matrix and calculate its symplectic eigenvalues—the positive partial transpose (PPT) criteria.[56,57] Such transformations are equivalent to a phase space reflection of a single party in the N‐partite state.[50] All minimum symplectic eigenvalues {νi}i=1N$\lbrace \nu _i\rbrace _{i=1}^N$ would be less than one, which indicates that this partially time‐reversed state is unphysical; in other words, the original state is fully inseparable. As has been pointed out in ref. [54], if the purity of states cannot be guaranteed in an experimental setting, verification of full inseparability in a multimode system does not imply genuine multipartite entanglement (GME).The entanglement structure becomes more involved with increasing number of parties. While the symplectic transform approach indicates that any one partite were inseparable from the whole, a state that is a mixture of separable states would show full inseparability based on this PPT criterion. The states that cannot be written in such a way are called genuinely entangled[17] and the verification of such states differs from full inseparability. Using generalized position and momentum observables, an entanglement criterion has been derived and applied to confirm tripartite energy‐time entanglement of three spatially separated photons.[58] In particular, there is an universal GME criterion derived in ref. [47] and further refined in ref. [54]. This GME criterion utilizes only variances of quadrature operators and it can be used for entanglement verification without any additional measurements. This general criterion was recently employed for verification of genuine tripartite entanglement of microwaves in a double superconducting cavity setting.[34]The GME criterion is based on the weighted variance of the quadratures, u=∑ihix∼i$u = \sum _i h_i \tilde{x}_i$ and v=∑kgkp∼k$v = \sum _k g_k \tilde{p}_k$; i,k={1,2,3,⋯,N}$i,k=\lbrace 1,2,3,\dots , N \rbrace$. Violation of the inequality15S≡Δu2+Δv2f3(hi,gi)≥1$$\begin{eqnarray} S \equiv \frac{\left&lt; \Delta u^2 \right&gt; + \left&lt; \Delta v^2 \right&gt;}{f_3(h_i,g_i)} \ge 1 \end{eqnarray}$$where16f3(hi,gi)=12min{|h1g1+h2g2|+|h3g3|,|h3g3+h2g2|+|h1g1|,|h1g1+h3g3|+|h2g2|,}$$\begin{eqnarray} f_3(h_i,g_i) = \frac{1}{2} \min \lbrace & |h_1 g_1 + h_2 g_2| + |h_3 g_3|, \nonumber \\ & |h_3 g_3 + h_2 g_2| + |h_1 g_1|, \nonumber \\ & |h_1 g_1 + h_3 g_3| + |h_2 g_2|, \rbrace \end{eqnarray}$$is sufficient to confirm genuine tripartite entanglement (N=3$N=3$) and violation of17S≡Δu2+Δv2f4(hi,gi)≥1$$\begin{eqnarray} S \equiv \frac{\left&lt; \Delta u^2 \right&gt; + \left&lt; \Delta v^2 \right&gt;}{f_4(h_i,g_i)} \ge 1 \end{eqnarray}$$where18f4(hi,gi)=min{|h1g1+h2g2+h3g3|+|h4g4|,|h4g4+h2g2+h3g3|+|h1g1|,|h4g4+h1g1+h3g3|+|h2g2|,|h4g4+h1g1+h2g2|+|h3g3|,|h1g1+h2g2|+|h3g3+h4g4|,|h1g1+h3g3|+|h2g2+h4g4|,|h2g2+h3g3|+|h1g1+h4g4|}$$\begin{eqnarray} f_4(h_i,g_i) = \min \lbrace & |h_1 g_1 + h_2 g_2 + h_3 g_3| + |h_4 g_4|, \nonumber \\ & |h_4 g_4 + h_2 g_2 + h_3 g_3| + |h_1 g_1|, \nonumber \\ & |h_4 g_4 + h_1 g_1 + h_3 g_3| + |h_2 g_2|, \nonumber \\ & |h_4 g_4 + h_1 g_1 + h_2 g_2| + |h_3 g_3|, \nonumber \\ & |h_1 g_1 + h_2 g_2| + |h_3 g_3 + h_4 g_4|, \nonumber \\ & |h_1 g_1 + h_3 g_3| + |h_2 g_2 + h_4 g_4|, \nonumber \\ & |h_2 g_2 + h_3 g_3| + |h_1 g_1 + h_4 g_4|\rbrace \end{eqnarray}$$is sufficient to confirm genuine quadripartite entanglement (N=4$N=4$) with weights hi,gk$h_i,g_k$ being in range [ − 1, 1]. To simplify the search domain, we set h1=g1=1$h_1=g_1=1$ and hi=h,gi=g$h_i = h, g_i = g$, for i={2,3}$i=\lbrace 2,3 \rbrace$ or i={2,3,4}$i=\lbrace 2,3,4 \rbrace$ with respect to the number of parties N.[54]In the double and triple pumping schemes, we generate generalized tripartite and quadripartite H$\mathcal {H}$‐graph states, which have a different structure compared with Greenberger‐Horne‐Zeilinger (GHZ) states.[7,38,44] The generalization deals with addition of BS correlations in the H$\mathcal {H}$‐graph structure. However, by applying additional pump tones and adjusting the phase difference between them, we can obtain regular GHZ type of entangled states. Thus, our scheme facilitates control of TMS and BS correlations and, thereby, allows tuning of the structure of the entangled state.Typically, the experimental weights of the graph edge connections are slightly non‐symmetric due to imperfections in the measurement settings. This results in a difference in the optimal weight values in the GME criterion. In order to find the full violation of the criterion in our analysis, we swap over all possible “base” modes (with weights hi=gi=1$h_i=g_i=1$) in order to detect the minimum value for S.ExperimentExperimental MethodsIn our microwave experiments, we employ a niobium λ/4$\lambda /4$ coplanar 50Ω$50\,\Omega$ transmission line terminated into a SQUID loop (QWJPA), forming a quarter‐wave Josephson parametric amplifier. The SQUID's junctions are formed by Nb/Al‐Al2O3/Nb 1× 1 μ$\umu$m tunnel barriers with Ic≃4$I_\text{c}\simeq 4$ μ$\umu$A critical current. The JPA, operating in three‐wave‐mixing mode around the cavity frequency ωr$\omega _\text{r}$, is pumped by an external RF magnetic flux through the SQUID loop using a single turn pump coil at frequency 2ωr$2\omega _\text{r}$ (marked as 2ωLO in Figure 3a).[40,59] We chose this operation regime, because for four‐wave mixing (typically using a current pump near ωr$\omega _\text{r}$), the large amplitude pump is within the amplification bandwidth, whereas the three‐wave mixing process separates the pump tone from the amplified signals, thus simplifying the practical use of the JPA. The loaded quality factor at the operation frequency is ≈900, while the internal Q is by a factor of three larger. The basic (zero‐flux) resonator frequency is 6.115 GHz, it can be tuned below 5.5 GHz by imposing external DC magnetic flux through the SQUID.3FigureExperimental scheme and device characterization. a) Principle of the experimental setup for tripartite entanglement measurements. The device is connected to the test ports via circulators. The DC bias current and AC pumping of flux are combined and reseparated in bias‐T components. Depending on the measurement type, input is connected either to a vector network analyzer (VNA) or a 50 Ω termination, whereas the output is directed either to a VNA, a signal analyzer, or a analog/digital converter. The frequency span of spectral modes and their separation is given for tripartite and quadripartite case in frames (b,c), respectively.Our measurement setup is illustrated in Figure 3a. The experiments were conducted at 20 mK using a BlueFors LD400 dry dilution cryostat. The JPA was protected from external magnetic fields using a Cryoperm shield. The DC flux bias and the RF pump shared a common on‐chip flux line, and the signals were combined in an external bias‐tee. Since our basic microwave setting is for reflection measurements, the sample is connected to the input and output ports via a circulator having a frequency band of 4 − 12 GHz. A vector network analyzer (VNA) was used to characterize the sample, whereas during the entanglement generation measurements, the signal port was kept terminated. By applying a multitone pump to JPA, correlated microwaves are generated from vacuum fluctuations. In the tripartite setting illustrated in Figure 3a, we had control over the relative phases of the pumps directly whereas, in the quadripartite case, phase rotation was possible in the data analysis only. Basic experimental data in the tripartite case as well as determination of the cavity parameters κ, γ, and K are discussed in Appendix A.4.Tripartite case: In the tripartite case, phase‐controlled pump signals from the RF waveform generator are mixed with the frequency‐doubled local oscillator (LO) frequency, filtered by a pair of home‐made, tunable bandpass cavity filters. In order to avoid spurious pumping at ω2LO, a band rejection filter tuned to 2ωLO is employed. The filtering ensures passage only for the desired two pump signal at angular frequencies ω1=ωr−Δ/2$\omega _1=\omega _\text{r}-\Delta /2$ and ω2=ωr+Δ/2$\omega _2=\omega _\text{r}+\Delta /2$. Sufficient noise thermalization is ensured by the 46 dB attenuation because the pump coil is only weakly coupled to the SQUID.We apply a DC magnetic flux of ΦDC=0.383Φ0$\Phi _{\text{DC}}=0.383\Phi _0$ through the SQUID loop, resulting in ωr2π=6.024$\frac{\omega _\text{r}}{2\pi }=6.024$ GHz for the cavity frequency. In the three‐mode experiment, we apply two pump tones at (2×ωr2π−2)$(2\times \frac{\omega _\text{r}}{2\pi }-2)$ MHz and (2×ωr2π+2)$(2\times \frac{\omega _\text{r}}{2\pi }+2)$ MHz, with the half‐frequencies positioned as depicted in Figure 3b and the correlated spectral modes are defined symmetrically with respect to the pump half‐frequencies. Each mode has a bandwidth of 1.9 MHz and is separated from the other modes by 0.1 MHz. The phase control in the measurement is facilitated by phase‐locking of microwave generators to a 10 MHz rubidium reference clock and by using a joint external trigger.To collect data for correlation analysis, we mix down the output signal using a synchronized LO signal and record the output quadratures using two channels of a Teledyne SP Devices ADQ14 digitizer with sample rate of 50 MSa s−1 per channel covering the bandwidth of 25 MHz. Furthermore, we employ an overall detuning of 14 MHz, that is, a heterodyne detection scheme, in order to avoid 1/f$1/f$ noise from the measurement devices and the IQ mixer in the frequency conversion part of the setup. Using digital postprocessing, we can easily shift the center frequency of the heterodyned MHz signal to zero, ready for final correlation analysis of the modes.Our three‐mode measurement scheme provides the remarkable advantage of physical control of the phase difference between the two pump tones, which is essential for the analysis of phase dependence in the entangled states. The phase difference between 2 LO signal used as the carrier of the pump tones and LO readout signal remains fixed. Therefore, a change of the initial phase of the IF‐signal in one of the pump generator channels (Agilent 33600B in Figure 3a) relative to the second IF channel creates an effective phase difference Δφ$\Delta \varphi$ between two pump tones. Importantly, this difference Δφ$\Delta \varphi$ is preserved after mixing with 2LO (using simplified notation): e−i(ω1t+Δφ)e−iω2LOt=e−i((ω1+ω2LO)t+Δφ);e−i(ω2t)e−iω2LOt=e−i((ω2+ω2LO)t)$e^{-i(\omega _1 t + \Delta \varphi )}e^{-i\omega _{\text{2LO}} t}=e^{-i((\omega _1 +\omega _{\text{2LO}})t+\Delta \varphi )}; e^{-i(\omega _2 t)}e^{-i\omega _{\text{2LO}} t}=e^{-i((\omega _2 +\omega _{\text{2LO}})t)}$.Since our fully phase‐locked scheme preserves the phases of the received, demodulated output quadratures, the measurement of covariance matrix components can be averaged for reducing noise in the elements. Finally, the reference phase of a single mode (defining the basis for I and Q) can be adjusted in postprocessing step in such a way that the corresponding subspace of the covariance matrix becomes a diagonal 2 × 2 matrix.In the tripartite case, indeed, we find that the hardware‐controlled relative phase rotation (in addition to the reference phase value to both channels of the pump generator) is equivalent to a proper phase rotation in the postprocessing step. The postprocessing will be discussed in more detail in Section 4.Quadripartite case: In the four spectral mode case, we simplify the experimental setup by eliminating the physical phase control, and replaced it by postprocessing of the received signals. This simplification possibility highlights the scalability of our entanglement generation method. The employed digital postprocessing is equivalent to hardware‐level selective separation of spectral modes into four channels, for example, using bandpass filters in conjunction with power splitters, and additional tunable delay lines for each selected spectral mode frequency.We apply a DC magnetic flux of ΦDC=0.417Φ0$\Phi _{\text{DC}}=0.417\,\Phi _0$ through a SQUID resulting in ωr=5.978$\omega _\text{r}=5.978$ GHz cavity frequency that slightly differs from the tripartite case, see Figure 3b,c. In order to generate quadripartite correlations, we apply three pump tones using Anapico APMS 12G generator, using strong high‐pass filtering (2 of Mini‐Circuits VHF‐8400+) to avoid subharmonic transmission to the circuitry. In this scheme, we avoid any external mixers for the input pump microwaves. By applying three phase‐locked pump tones at frequencies 2×ωr2π$2\times \frac{\omega _\text{r}}{2\pi }$ MHz, (2×ωr2π+1)$(2\times \frac{\omega _\text{r}}{2\pi }+1)$ MHz, and (2×ωr2π−1)$(2\times \frac{\omega _\text{r}}{2\pi }-1)$ MHz, we generated four correlated spectral modes out from the ground state of the microwave cavity. Each mode has a bandwidth of 0.4 MHz and is separated from adjacent modes by 0.2 MHz. The output microwaves are captured, mixed down, and digitized by Anritsu MS2830A Signal Analyzer with a bandwidth of 2 MHz. Again, averaging is needed to lower the noise in the covariance matrix elements, and in this scheme, digital postprocessing is necessary to unify the phase settings in the covariance matrices before summation.The experimental detection of H∼$\mathcal {\tilde{H}}$‐graph states and their genuine multipartite entanglement depends on relations among the covariance matrix elements as discussed in Section 2.5. The degree of violation in the GME condition S<1$S&lt;1$ depends strongly on the magnitude ratio of the diagonal covariance elements to the off‐diagonal ones. Therefore, calibration of the detected signal powers is decisive, which is discussed in Appendix A.3.Scaled Covariance MatrixThe system gain <GΣ>$&lt;G_{\Sigma }&gt;$ determined in Appendix A.3 refers to measured power per unit bandwidth. Since the measured spectral mode quadratures Ii$I_i$ and Qi$Q_i$ are determined over the band Δfi$\Delta f_i$, the scaled quadrature xi$x_i$, equivalent to the amplitude of the quantum mechanical operators x$\mathbf {x}$, is given by the formula19xi=IiGiZ0hfiΔfi$$\begin{equation} {x}_i=\frac{I_i}{\sqrt {G_i Z_0 h f_i \Delta f_i}} \end{equation}$$where Gi=<GΣ,i>$G_i=&lt;G_{\Sigma ,i}&gt;$ is the system gain for ith spectral mode, Z0=50$Z_0 = 50$ Ω is transmission line impedance and Δfi$\Delta f_i$ is the bandwidth of the spectral mode: Δfi=2$\Delta f_i = 2$ MHz or Δfi=0.4$\Delta f_i = 0.4$ MHz for tripartite and quadripartite cases, respectively. Similar scaling is applied also to the quadrature component Qi$Q_i$.Similar to our earlier work,[43] the noise added by the preamplifier is subtracted from the diagonal elements of the covariance matrix (see Equation (14))20V=4Von−Voff+Icothhfi2kbTi$$\begin{equation} \mathcal {V}=4{\left(\mathbf {V}_{\text{on}}-\mathbf {V}_{\text{off}} \right)} + \mathbf {I} \coth \frac{hf_i}{2 k_\text{b} T_i} \end{equation}$$where Voff$\mathbf {V}_{\text{off}}$ denotes the covariance matrix measured in the absence of the pump. Due to scaling of the covariance matrix V$\mathcal {V}$, this equation yields a unity diagonal matrix in the absence of pumping at T→0$T \rightarrow 0$. The average physical temperature in our experiments is Ti=20$T_i=20$ mK resulting in cothhfi2kbTi=1.000$\coth \frac{hf_i}{2 k_\text{b} T_i}=1.000$.Multipartite EntanglementTo characterize the structure of the entanglement in output states, we analyze the resulting covariance matrices using PPT formalism and GME criteria discussed in Section 2.5 for tripartite and quadripartite cases.Tripartite caseLeveraging the amplitude and phase control of the pump signals, we experimentally evaluate the PPT and GME criteria values at different pump parameters. For comparison, we also conducted detailed numerical simulations based on the QLE in Equation (A8) using experimentally determined JPA parameters in the measurements. In general, we find good agreement between simulations and the experimental data, which is reassuring concerning the validity of the results.Figure 4 depicts our experimental results on genuine tripartite entanglement and their comparison with simulations. Figure 4a illustrates results of numerical simulations on GME in terms of S defined in Equation (15). At weak pumping, the condition for genuine entanglement S<1$S &lt; 1$ is fulfilled almost independent of the pump phases, but with increasing A, the simulations reveal an even smaller range of Δφ$\Delta \varphi$ yielding S<1$S &lt; 1$ (see the inset in Figure 4a). The strongest genuine tripartite entanglement is reached at Δφ≃−120∘$\Delta \varphi \simeq -120^{\circ }$ under normalized pumping amplitude A≃0.22$A \simeq 0.22$, at which the simulations reach S=0.70$S=0.70$. At the minimum of S, the corresponding weights are hi={1,−0.65,−0.65}$h_i=\lbrace 1,-0.65,-0.65\rbrace$ and gi={1,0.65,0.65}$g_i=\lbrace 1,0.65,0.65\rbrace$. It is noteworthy that the phase setting Δφ=+90∘$\Delta \varphi = + 90^{\circ }$ yields clearly worse entanglement than Δφ=−90∘$\Delta \varphi = -90^{\circ }$. This asymmetry in GME between Δφ=±90∘$\Delta \varphi = \pm 90^{\circ }$ arises from differences in the covariance matrices which is illustrated in Figure 5.4FigurePhase‐dependent genuine entanglement of tripartite bisqueezed state. a) Simulation results for GME criterion as a function of normalized pump amplitude for three different values of the phase difference Δφ$\Delta \varphi$ between pump signals indicated in the figure; the simulation parameters were set to match the experiment (see Appendix A.4). The inset illustrates S(A,Δφ)$S(A,\Delta \varphi )$ up to critical amplitude A=0.5$A=0.5$. The weights hi,gi$h_i, g_i$ were optimized in the calculation of S as discussed in the text. b) Experimental values for S as a function of A at the same phase difference values Δφ$\Delta \varphi$ between pump signals as in frame (a). The inset illustrates measured S(A,Δφ)$S(A,\Delta \varphi )$ up to A=0.5$A=0.5$. In general, the measured S(A,Δφ)$S(A,\Delta \varphi )$ corresponds quite well to the inset in frame (a). Due to noise, the measured GME nearly vanishes around Δφ≃+90∘$\Delta \varphi \simeq +90^\circ$ where even the simulated S is only slightly below 1. The best genuine multipartite entanglement is reached at Δφ=−90∘⋯−120∘$\Delta \varphi = -90^\circ \dots -120^\circ$ owing to phase shifts introduced by the cavity (see Figure A3 in Appendix A.5). The parametric drive changes the phase response of the cavity which leads to a shift in the optimum conditions for GME as a function of A.5FigureCovariance matrix of the genuinely entangled tripartite bisqueezed state. Experimentally obtained tripartite covariance matrices for genuinely entangled bisqueezed states. The phase difference between the two displayed cases is 180°. Via this control of pump phases we demonstrate the rotation of the desired subspace elements, corresponding to TMS type of correlations (1 − 2 or 2 − 3). Subspace elements related to BS correlations (1 − 3), shows always dependency on the distribution of the elements in the corresponding driven TMS subspaces. However, the eigenvalues of the BS subspace remain constant.Our experimental data on S in Figure 4b displays similar features as Figure 4a. The measured GME criterion S as function of normalized pump amplitude for three phase differences is shown in Figure 4b. In the experimental data, nearly no GME is observed at positive phase differences, whereas Δφ=−90∘$\Delta \varphi = - 90^{\circ }$ and Δφ=−120∘$\Delta \varphi = - 120^{\circ }$ yield suppression down to S=0.75±0.05$S=0.75 \pm 0.05$. The measured result at Δφ=−120∘$\Delta \varphi = - 120^{\circ }$ follows quite well the simulated behavior as a function of the pump amplitude, and genuine entanglement is observed in the normalized amplitude range A∈[∼0.01,0.4]$A \in [\sim 0.01, 0.4]$. Overall, the pattern of S(φ,A)$S(\varphi , A)$ in the inset of Figure 4b coincides with the simulated pattern in Figure 4a. The agreement strongly supports the presence of genuinely entangled bisqueezed state in our experiment. We emphasize that the optimum entanglement at Δφ=−120∘$\Delta \varphi = - 120^{\circ }$, observed both in our simulations and in the experiment, cannot be obtained from a simple analytical calculations for the lossless, strongly coupled model. The reason is the frequency‐dependent phase response of the cavity due to finite coupling and dissipation rates (see Section A.5), which, when included in the simulations, result in very good matching with the experiment.Covariance matrices measured at the pump phase difference Δφ=+90∘$\Delta \varphi = + 90^{\circ }$ and Δφ=−90∘$\Delta \varphi = - 90^{\circ }$ are illustrated in Figure 5. Technically, by rotating the phase of a pump signal, we selectively control certain subspace of the covariance matrix, which can be seen in Figure 5. An applied phase shift to the first or the second pump rotates directly the subspace corresponding to two mode squeezing correlations, modes 1 − 2 or 2 − 3, respectively. If no phase shift to the selected pump tone is applied, the corresponding TMS subspace preserves its distribution of covariances. The subspace spanned by modes 1 and 3, corresponding to the beamsplitter type of correlations, has a structure according to products of the involved TMS subspaces. Distinct control of the BS subspace alone (leaving the TMS subspaces fixed) using a rotation of the pump phases is not possible.Comparison of Figures 5a,b, reveals how the subspaces transform with the phase difference from Δφ=−90∘$\Delta \varphi = - 90^{\circ }$ to Δφ=+90∘$\Delta \varphi = + 90^{\circ }$. Subspace 1 − 3 corresponding to BS correlations shows a sign inversion in its elements. Subspace 2 − 3 corresponds to the pump, the phase of which has not been changed and, thus, its elements remain fixed. The phase of the first pump has been changed by π which inverts the TMS correlation in subspace 1 − 2. These subspaces which are controlled by the pump phase settings, can expressed in symmetric <I1I2>≃<I2I3>$&lt;I_1I_2&gt; \simeq &lt;I_2I_3&gt;$ or antisymmetric form <I1I2>≃−<I2I3>$&lt;I_1I_2&gt; \simeq -&lt;I_2I_3&gt;$, see Figures 5a,b, respectively.The inseparability of the covariance matrix was investigated using the PPT criteria (see Section 2). Symplectic eigenvalues obtained for mode partitions 1 − 23, 2 − 13, and 3 − 12 are depicted in Figure 6a for simulations, while experimental data are displayed in Figure 6b; here the first index specifies the mode in which the sign of the momentum has been reversed. The results are plotted as a function of normalized pump amplitude since the phase difference between the pump tones does not play a role. Indeed, while the genuine entanglement is sensitive to both pump amplitude and phase difference, the PPT criterion is phase independent—the minimum symplectic eigenvalues remain constant when the phase difference is varied at fixed pump amplitude. Therefore, by exercising phase control over each pump, we gain the ability to switch from fully inseparable state to genuinely entangled state, without making any changes to the type of interaction between the modes. According to experimental results in Figure 6b, the middle frequency acted on by both pumps is the most inseparable part of the covariance matrix.6FigurePhase‐independent full inseparability of tripartite bisqueezed state. a) PPT criteria in terms of the minimum symplectic eigenvalues simulated for our double‐pump QWJPA using experimentally determined parameters. Eigenvalues min{νi}$\min \lbrace \nu _i\rbrace$ are traces over normalized pump amplitude A; permutations 1 − 23, 2 − 13, and 3 − 12 have been considered. The symplectic eigenvalues are the smallest for time‐reversed second mode (ν2−13$\nu _{2-13}$), which participates to both TMS processes. b) Experimentally determined symplectic eigenvalues for the same permutations ( in the experiment. Gray dashed line displays the full inseparability threshold. The difference in the simulated behavior of ν1−23$\nu _{1-23}$ and ν3−12$\nu _{3-12}$ is caused by asymmetry due to finite value of resonance detuning Δr$\Delta _\text{r}$.Quadripartite caseFor the quadripartite case, we apply three pump drives with identical amplitude α1=α2=α3=α$\alpha _1=\alpha _2=\alpha _3=\alpha$. While our goal is to demonstrate genuine entanglement generation of cluster states (mode structure depicted on Figure 2), we reject direct, physical phase control, and use digital postprocessing to transform the covariance matrix to the desired form on which we then verify its entanglement properties. However, we do preserve the coherence between pump tones by phase locking so that the relative phases do not fluctuate over time. By applying a postprocessing phase rotation for each mode separately, we bring the covariance matrix into the symmetric form (see Section 4).For the analysis of full inseparability of the covariance matrix according to the PPT criterion, we evaluate the minimum symplectic eigenvalues min{νi}$\min \lbrace \nu _i\rbrace$ of the following mode permutations: 1 − 234, 2 − 134, 3 − 124, 4 − 123. The experimentally obtained symplectic eigenvalues as function of normalized pump amplitude A are displayed in Figure 7 alongside with the corresponding predictions given by our numerical simulations. The minimum symplectic eigenvalue min{νi}=0.79±0.018$\min \lbrace \nu _i\rbrace =0.79\pm 0.018$ is reached, while all of the eigenvalues in the normalized pump amplitude range 0.01≲A≲0.15$0.01\lesssim A \lesssim 0.15$ are less than 1. Compared with the minimum symplectic eigenvalues in Figure 6, we may conclude that the influence of BS correlations on min{νi}$\min \lbrace \nu _i\rbrace$ value is less in the quadripartite state than for the tripartite case.7FigureFull inseparability and genuine entanglement of quadripartite H∼$\mathcal {\tilde{H}}$‐graph state (generalized H‐graph). a) Results on PPT criterion for four permutations of a four‐mode Gaussian state. Minimum eigenvalues min{νi}$\min \lbrace \nu _i\rbrace$, indicated as open circles, are traced over normalized pump amplitude A. Results of our QLE simulations, plotted as solid curves, exhibit good correspondence with the experimentally obtained values. The full inseparability condition {ν1−234<1}⋀{ν2−134<1}⋀{ν3−124<1}⋀{ν4−123<1}$\lbrace \nu _{1-234}&lt;1\rbrace \bigwedge \lbrace \nu _{2-134}&lt;1\rbrace \bigwedge \lbrace \nu _{3-124}&lt;1\rbrace \bigwedge \lbrace \nu _{4-123}&lt;1\rbrace$ is fulfilled in the range of 0.01≲A≲0.15$0.01\lesssim A\lesssim 0.15$. The green dashed line displays the entanglement threshold.The GME criterion for four modes as a function of the normalized pump amplitude is depicted in Figure 7b; the symbols display data while the simulation result is indicated by the solid curve. As was discussed in Section 2, the optimized weights in GME inequality Equation (17) are chosen in the same manner as in the tripartite case: h1=g1=1$h_1=g_1=1$ and hi=h,gi=g$h_i = h, g_i = g$, i={2,3,4}$i=\lbrace {2,3,4\rbrace }$. The strongest genuine entanglement S=0.84±0.02$S=0.84 \pm 0.02$ is observed at A≂0.08$A\eqsim 0.08$ pump amplitude using the weights hi={1,−0.51,−0.51,−0.51}$h_i=\lbrace 1,-0.51,-0.51,-0.51\rbrace$ and gi={1,0.69,0.69,0.69}$g_i=\lbrace 1,0.69,0.69,0.69\rbrace$. The numerical simulation provides h and g coefficients that coincide with the experimental values with 1% error, which strongly establishes that the states produced in the experiment coincide with the ones that were obtained and analyzed in the numerical model.The covariance matrices obtained in the experiment and using numerical simulation are presented in Figure 8a,b, respectively. They are determined at the strongest entanglement point reached at A≂0.08$A\eqsim 0.08$. TMS type of correlations are seen in the mode combinations 1 − 2, 2 − 3, 3 − 4, and 1 − 4. Subspaces corresponding to BS correlations are visible in the plot as product distributions in 1 − 3 and 2 − 4 subpartitions. The covariance matrix structures illustrated in Figure 8 correspond directly to the H∼$\mathcal {\tilde{H}}$‐graph structures shown on Figure 2b. In general, we conclude that for the employed pump configuration, the genuine quadripartite entanglement appears in the amplitude range 0.01≲A<0.13$0.01\lesssim A&lt;0.13$.8FigureCovariance matrix of genuinely entangled quadripartite H∼$\mathcal {\tilde{H}}$‐graph state. a) Experimental covariance where the rotation of the TMS subspaces 1−2, 2−3, and 3−4 have been made in such a way that the structure coincides with the matrix in Equation (A36) of Appendix A.6 (each pump has phase π2$\frac{\pi }{2}$). The employed pump amplitude A≂0.08$A\eqsim 0.08$ yields the smallest value for S. b) Simulated covariance matrix using equal pump phases π2$\frac{\pi }{2}$ at A≂0.08$A\eqsim 0.08$. The difference from the matrix in Equation (A36) is due to the cavity response that induces extra phase shifts.DiscussionThe control of bisqueezed tripartite and generalized H‐graph (H∼$\mathcal {\tilde{H}}$‐graph) quadripartite states by relative positioning of the pump frequencies and their phases is indicative of the strong potential of these methods for CV quantum state processing. The basic parametric microwave setting allows for enhancement in the number of spectral modes by additional pump tones, which leads to generation of more complex, entangled H∼$\mathcal {\tilde{H}}$‐graph states. Enhanced number of modes requires larger bandwidth, which calls for broadband parametric devices such as TWPAs[60–64] or broadband JPAs[59,65] in order to avoid problems with spectral mode crowding.Our approach based on QLE puts in evidence additional correlations, which are captured by the definition of H∼$\mathcal {\tilde{H}}$‐graph states. The correlations arise naturally from the connection between intracavity modes and the input vacuum modes, due to which the same vacuum fluctuations may act in the downconversion of more than one quanta. In the literature on cluster and H‐graph states, the adjacency matrix for H‐graphs is defined via the matrix specified in the multimode squeezing Hamiltonian. The QLE analysis corresponds to the expansion of the multimode squeezing operator up to second order, which leads to the appearance of beamsplitter correlations in the adjacency matrix. In Appendix A.6 (Equations (A37)–(A39)) we show how to use well‐chosen relative pump phase values in the quadripartite case to prepare an entangled square lattice state—that is, a state without BS correlations. For the case of very large squeezing, H∼$\mathcal {\tilde{H}}$‐graph state can be regarded as an approximation of a 4‐node cluster state, minimizing errors in gate operations of measurement‐based CV quantum computing.Cluster states form a promising platform for scalable quantum information processing. In one‐way quantum computing,[3] the entire computational resource is provided by the entanglement of the cluster state. The processing is based on quantum measurements which facilitate gate operations as well as the read‐out of the final result. However, cluster states can be obtained from graph states only in the mathematical limit of large squeezing parameter.[27–29,32] For quantum information processing steps, it is sufficient to perform sub‐cluster measurements in specified order using a suitable computational basis. In refs. [28, 66, 67] different computation scenarios based on resources provided by squeezing generators and beamsplitters are described. Encoding, gate and measurement operations have been so far considered in optical circuits for continuous variable quantum data and can be efficiently extended to the microwave realm. In this work, we have utilized this correspondence between optics and microwaves and demonstrated H∼$\mathcal {\tilde{H}}$‐graph state encoding.In contrast to computational models for graph states[38] considered as ideal clusters, hardware based on finite squeezing with noise and decoherence requires error correction procedures[6,68,69] to provide reliable CV computation. Using the presented scheme one can implement error correction codes based on the idea of repetitions of selective measurements and new encoding of H∼$\mathcal {\tilde{H}}$‐graph states before each gate operation. In ref.[ [70]] a multidimensional platform for scalable quantum computing has been proposed, based on cluster states created using microring resonators; also multiple frequency combs[67] created by optical parametric amplifiers and beamsplitters can serve as an excellent platform for quantum computation. Our work shows that the methods of generation of highly‐entangled CV states are not restricted to just optical parametric amplifiers, but the methods can be carried over into the microwave domain by employing parametric Josephson junction devices for creation of topologically involved and structurally versatile H∼$\mathcal {\tilde{H}}$‐graph states.An implementation of the universal quantum computer based on bosonic modes with the possibility of hardware‐efficient quantum error correction[71] requires efficient generation of continuous‐variable quantum resources. The genuine entanglement between several bosonic modes could potentially be employed for error‐correctable codeword states.[72] Besides potential in error correction, the introduction of entanglement into quantum measurement implementations leads to a quantum advantage in the detection process when detection is performed in the presence of high level of noise and loss.[73]Increase of the number of entangled spectral modes is essential for future technological application of these CV quantum state generation methods. The limiting factors are the requirements of high precision for the pump frequency and its phase, the stability of the biasing flux, and possible crowding of modes within a narrow‐band JPA resonance. However, recently it has been demonstrated that entanglement can be generated in low‐loss traveling wave parametric amplifiers.[] This opens a way to significant increase in the number of entangled modes.ConclusionIn this work, we presented a practical scheme for generation of controllable multipartite entanglement from vacuum fluctuations, based on multitone pumping scheme of a JPA, which facilitates pivotal resources for quantum technologies at microwave frequencies. While optical schemes for multipartite entanglement generation operate on even larger clusters, they lack versatility and are limited to optical frequencies as such. On the other hand, our scheme allows for a flexible increase in the number of modes and control of the entanglement configuration among the modes by adjusting pumping on the same device, whereas optical setups call for massive hardware reconfiguration when the entanglement structure is altered. Through phase and amplitude variation of the microwave pump tones, we reach a comprehensive control over the entanglement structure within the spectral modes of a single JPA cavity mode, which we experimentally verify in detail for the tripartite case.Using the developed scheme, we made the first successful demonstration of an on‐demand tunable, fully inseparable, and phase controllable genuinely entangled tripartite and quadripartite states in a superconducting system. The presence of multipartite quantum correlations was verified using the covariance matrix formalism and genuine entanglement criteria constructed from the measured quadratures. Experimental results were accurately reproduced by calculating symplectic eigenvalues of a partially‐transposed covariance matrix for full inseparability detection as well as computing GME criteria in normalized pump amplitude in range 0<A<0.5$0&lt;A&lt;0.5$ (0<A<0.25$0&lt;A &lt;0.25$) and verified genuine entanglement in the range of 0.01≲A<0.4$0.01\lesssim A&lt;0.4$ (0.01≲A<0.13$0.01\lesssim A&lt;0.13$) for the tripartite (quadripartite) state.We provided results of phase‐dependent GME criterion for bisqueezed state. With optimal phase shift between two pumping tones Δφ=−120∘$\Delta \varphi =-120^\circ$ minimum value of criterion S=0.75±0.05$S=0.75 \pm 0.05$ was obtained. This result was also faithfully reproduced by numerical simulations.In our analytical derivations, we demonstrated additional control possibilities over the BS correlations in the covariance matrix of quadripartite H∼$\mathcal {\tilde{H}}$‐graph state. To visualize the formed entanglement structure, we provided an extension for the known H‐graph adjacency matrix: besides TMS, it includes BS correlations between the vacuum modes. The QLE approach was used to introduce such an adjacency matrix and to connect it to the general approach starting from multimode squeezing operator and the TMS Hamiltonian for the multi‐mode case with multiple pumps. As shown in Appendix A.6, BS correlations can be fully suppressed by implementing a 180° phase shift of one pump. Such a phase combination creates a distinct square‐lattice H‐graph state which, for the limit of infinite squeezing parameter, transforms to a square‐lattice cluster state.Additional TMS correlations can be introduced by inserting new pump tones, which can change the nature of the entangled states drastically. For example, using two additional pump tones with half frequencies at {−Δ2;Δ2}$\lbrace -\frac{\Delta }{2};\frac{\Delta }{2}\rbrace$, we are able to connect all 4 modes with TMS correlations and thereby achieve a GHZ‐like state. Furthermore, by tuning the phases of the pumps, the state can be converted into square‐lattice H‐graph state. With the bandwidth improvements provided by the state‐of‐the‐art superconducting parametric devices, such as the broadband, low‐loss travelling wave parametric amplifier,[62–64] we expect a substantial increase in the number of entangled modes, which facilitates generation of highly‐squeezed square‐lattice H‐graph states for CV quantum computation at microwave frequencies.AAppendixAnalytical Methods, Experimental Procedures, and Entanglement AnalysisA.1Details of Theoretical DescriptionThe Hamiltonian of JPA system is given byA1Ĥsys=ℏωrâ†â+ℏ2∑d=1pαd∗eiωdt+αde−iωdt(â2+â†2)+ℏK(â+â†)4$$\begin{eqnarray} {\hat{H}}_{\text{sys}}=\hbar \omega _r \hat{a}^\dagger \hat{a} & + \frac{\hbar }{2}\sum _{d=1}^{p}{\left[\alpha ^{*}_d e^{i\omega _d t} + \alpha _d e^{-i\omega _d t}\right]}(\hat{a}^{2} + \hat{a}^{\dagger 2}) \nonumber \\ & + \hbar K (\hat{a}+ \hat{a}^{\dagger })^{4} \end{eqnarray}$$where â(â†$\hat{a}(\hat{a}^{\dagger }$) is the annihilation (creation) operators for cavity photons, αd$\alpha _d$ is the complex amplitude for pump tone d, and K denotes the strength of the Kerr nonlinearity term. Using the average of p pump tones ωd$\omega _d$, d={1,⋯,p}$d=\lbrace 1,\dots ,p \rbrace$, we define the detuning between the half pump frequency and the resonator frequency: Δr=ωr−ωΣ2$\Delta _\text{r}=\omega _\text{r}-\frac{\omega _{\Sigma }}{2}$, ωΣ=∑d=1pωdp$\omega _{\Sigma }=\frac{\sum _{d=1}^{p}\omega _d}{p}$.For each of the p pump tones, we define the detuning from the average frequency Δd=ωd−ωΣ,d={1,⋯,p}$\Delta _d=\omega _d-\omega _{\Sigma }, d=\lbrace 1,\dots ,p \rbrace$. By applying the rotating wave approximation in the frame ωΣ/2$\omega _\Sigma /2$ (a∼(t)=â(t)eiωΣt/2$\tilde{a}(t) = \hat{a}(t)e^{i\omega _\Sigma t/2}$) and leaving only the effective high‐order terms, we obtain for the nonlinear part of the HamiltonianA2Hsys,rwa(t)=ℏΔra∼†a∼+ℏ2∑d=1p(αd∗eiΔdta∼2+αde−iΔdta∼†2)+6ℏKa∼†a∼†a∼a∼$$\begin{eqnarray} {H}_{\text{sys,rwa}}(t)={\hbar }\Delta _r \tilde{a}^\dagger \tilde{a} & + \frac{\hbar }{2}\sum _{d=1}^{p}(\alpha ^{*}_d e^{i\Delta _d t}\tilde{a}^{2} + \alpha _d e^{-i\Delta _d t} \tilde{a}^{\dagger 2}) \nonumber \\ & + 6{\hbar }{K} \tilde{a}^{\dagger }\tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$As usual, the bosonic commutation relationships are valid for the cavity modes [a∼,a∼†]=1$[ \tilde{a},\tilde{a}^{\dagger }]=1$.The parametric resonator is coupled to a transmission line via the signal port and to the thermal bath via a linear dissipation port. The coupling Hamiltonian associated with the signal port is given byA3Hsig=ℏ∫dωb∼†b∼+κa∼†b∼−κ∗b∼†a∼$$\begin{equation} H_{\text{sig}}=\hbar \int d\omega {\left(\tilde{b}^{\dagger }\tilde{b} + \kappa \tilde{a}^{\dagger } \tilde{b} - \kappa ^{*} \tilde{b}^{\dagger } \tilde{a}\right)} \end{equation}$$where creation and annihilation operators b∼†$\tilde{b}^{\dagger }$ and b∼$\tilde{b}$ refer to modes in the transmission line, and κ denotes the coupling rate. The Hamiltonian related to the linear dissipation portA4Hloss=ℏ∫dωc∼†c∼+γa∼†c∼−γ∗c∼†a∼$$\begin{equation} H_{\text{loss}}=\hbar \int d\omega {\left(\tilde{c}^{\dagger }\tilde{c} + \gamma \tilde{a}^{\dagger } \tilde{c} - \gamma ^{*} \tilde{c}^{\dagger }\tilde{a}\right)} \end{equation}$$where c∼†$\tilde{c}^{\dagger }$ and c∼$\tilde{c}$ describe creation and annihilation of thermal bath modes and the rate γ represent the coupling of cavity modes to the linear dissipation port. The transmission line and bath operators obey the commutation relationsA5b∼(ω),b∼†(ω′)=c∼(ω),c∼†(ω′)=δ(ω−ω′)$$\begin{equation} {\left[ \tilde{b}(\omega ),\tilde{b}^{\dagger }(\omega ^{\prime })\right]}={\left[ \tilde{c}(\omega ),\tilde{c}^{\dagger }(\omega ^{\prime })\right]}=\delta (\omega -\omega ^{\prime }) \end{equation}$$andA6b∼(ω),c∼†(ω′)=c∼(ω),b∼†(ω′)=b∼(ω),b∼(ω′)=c∼(ω),c∼(ω′)=0$$\begin{eqnarray} {\left[ \tilde{b}(\omega ),\tilde{c}^{\dagger }(\omega ^{\prime })\right]} = & {\left[ \tilde{c}(\omega ),\tilde{b}^{\dagger }(\omega ^{\prime })\right]} = {\left[ \tilde{b}(\omega ),\tilde{b}(\omega ^{\prime })\right]} \nonumber \\ & ={\left[ \tilde{c}(\omega ),\tilde{c}(\omega ^{\prime })\right]}=0 \end{eqnarray}$$The total Hamiltonian can be conveniently written as a sum of the separate parts given aboveA7H(t)=Hsys,rwa(t)+Hsig+Hloss$$\begin{equation} {H}(t)=H_{\text{sys,rwa}}(t)+H_{\text{sig}}+H_{\text{loss}} \end{equation}$$For further analysis and for our simulations, we use the Quantum Langevin Equation (QLE) for the cavity operator a∼(t)$\tilde{a}(t)$A8a∼̇(t)=(−iΔr−κ+γ2)a∼−i∑d=1pαdeiΔdta∼†+κb∼in+γc∼in−12iKa∼†a∼a∼$$\begin{eqnarray} \dot{\tilde{a}}(t)=(-i\Delta _r &- \frac{\kappa +\gamma }{2})\tilde{a} - i \sum _{d=1}^{p}\alpha _d e^{i\Delta _d t} \tilde{a}^\dagger \nonumber \\ & + \sqrt {\kappa }\tilde{b}_{in}+\sqrt {\gamma }\tilde{c}_{in} - 12i{K} \tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$where the presence of the Kerr term allows us to consider dynamics of the parametric resonator above the critical oscillation threshold. To obtain the modes coming out from the cavity, we employ the standard input–output formalism which yields the relationshipA9b∼out(t)=b∼in(t)−κa∼(t)$$\begin{equation} \tilde{b}_{\text{out}}(t)=\tilde{b}_{\text{in}}(t)-\sqrt {\kappa }\tilde{a}(t) \end{equation}$$Equations (A8) and (A9) are used in our numerical simulations with Matlab ODE45 solver.A.2Full InseparabilityAssuming that the microwave fields produced by the JPA below the critical pumping threshold are Gaussian,[34] the states with multiple spectral modes can be fully characterized by measuring the covariance matrix of corresponding in‐phase I and quadrature Q voltages. For the measurement of tripartite correlations, we collect quadrature data for 0.8 s at every phase difference and pump amplitude value, without averaging. For the quadripartite case, we repeat the experiment 20 times at each pump power value and every quadrature sequence has a duration of 1.3 s.The quantum quadratures x∼i=a∼i+a∼i†2$ \tilde{x}_i = \frac{\tilde{a}_i+\tilde{a}_i^{\dagger }}{2}$ and p∼i=a∼i−a∼i†2i$\tilde{p}_i=\frac{\tilde{a}_i-\tilde{a}_i^{\dagger }}{2i}$ can be combined into a 2N‐long column vector operator for the N‐mode state r∼=(x∼1,p∼1,⋯x∼N,p∼N)T$\tilde{r} = (\tilde{x}_1, \tilde{p}_1, \dots \tilde{x}_N, \tilde{p}_N )^{T}$. The commutation relations can be written down in a skew‐symmetric, block‐diagonal matrix form[50]A10r∼i,r∼j=i2ΩijandΩ=⨂i=1N01−10$$\begin{equation} {\left[\tilde{r}_i, \tilde{r}_j\right]}=\frac{i}{2}\Omega _{ij} \text{ and } \mathbf {\Omega } = \bigotimes _{i=1}^N \def\eqcellsep{&}\begin{bmatrix} 0 & 1\\ -1 & 0 \end{bmatrix} \end{equation}$$The covariance matrix V$\mathbf {V}$ is given by elements Vij=12<Δr∼iΔr∼j+Δr∼jΔr∼i>−<Δr∼i><Δr∼j>$V_{ij}=\frac{1}{2}&lt;\Delta \tilde{r}_i \Delta \tilde{r}_j + \Delta \tilde{r}_j \Delta \tilde{r}_i &gt; - &lt;\Delta \tilde{r}_i &gt; &lt;\Delta \tilde{r}_j&gt;$ where we have defined standard error Δr∼i=r∼i−<r∼i>$\Delta \tilde{r}_i=\tilde{r}_i-&lt;\tilde{r}_i&gt;$ and <r∼i>=tr(r∼iρ̂)$&lt;\tilde{r}_i&gt;= \text{tr}(\tilde{r}_i \hat{\rho })$. The uncertainty principle requires thatA11V+i4Ω≥0$$\begin{equation} \mathbf {V}+\frac{i}{4}\mathbf {\Omega }\ge 0 \end{equation}$$applies for a physical covariance matrix.For verification of entanglement, we may investigate a modified equationA12V′+i4Ω≥0$$\begin{equation} \mathbf {V}^\prime +\frac{i}{4}\mathbf {\Omega }\ge 0 \end{equation}$$where Vk′=λkVλk$\mathbf {V}^\prime _k=\pmb {\lambda }_k \mathbf {V} \pmb {\lambda }_k$, λk$\pmb {\lambda }_{k}$ is diagonal matrix with ones entries, except of that related to kth mode, with value of −1. For example, transformation with λk≡N=diag(1,1,…,1,−1)$\pmb {\lambda }_{k\equiv N}=\text{diag}(1,1,\ldots,1,-1)$ means a partial transposition of the covariance matrix with respect to the last mode. The positive partial transpose (PPT) criterion for multipartite case requires that there is a violation of Equation (A12) when applying a partial transposition with respect to each from full set of modes: Vk′≥i4Ω$\mathbf {V}^\prime _k \ge \frac{i}{4}\mathbf {\Omega }$. In ref. [47], the entangled states are classified in accordance to the number of modes for which the condition Equation (A12) is broken. We follow this approach to demonstrate the highest class—full inseparability—in four mode case.Unitary operations which retain the Gaussian character of the states, for example, squeezing, are of particular importance. Such operations on the Hilbert space correspond to a linear transformation P in the phase‐space which preserve the symplectic form, that is,A13Ω=PTΩP$$\begin{equation} \mathbf {\Omega } = \mathbf {P}^T \mathbf {\Omega } \mathbf {P} \end{equation}$$Symplectic transformations on a 2N‐dimensional phase‐space form the real symplectic group denoted as Sp(2N;R)$ Sp(2N; {\bf R})$, which is a proper subgroup of the special linear group of 2N × 2N matrices.[74] By utilizing Williamson's theorem,[75] any covariance matrix can be expressed in the Williamson normal form:A14V∼k=PTVk′P$$\begin{equation} \mathbf {\tilde{V}}_k= \mathbf {P}^T \mathbf {V}^\prime _k \mathbf {P} \end{equation}$$where V∼k$\mathbf {\tilde{V}}_k$ is a 2N‐dimensional diagonal matrix consisting of the symplectic eigenvalues, ν∼k$\tilde{\nu }_k$, of the covariance matrix. The symplectic eigenvalues are called the symplectic spectrum which provides a practical means to verify physicality and various entanglement criteria. Separability is in force, when condition ν∼k≥1/4$\tilde{\nu }_k \ge 1/4$ fulfilled for V∼k$\mathbf {\tilde{V}}_k$.For convenience, we insert an additional factor of 4 to the covariance matrix and work with fluctuations with zero mean values: Vij∗=2<Δr∼iΔr∼j+Δr∼jΔr∼i>${V}^*_{ij}=2&lt;\Delta \tilde{r}_i \Delta \tilde{r}_j + \Delta \tilde{r}_j \Delta \tilde{r}_i &gt;$. Consequently, for evidence of “fully inseparable” states, we need to find minimum symplectic eigenvalues with ν∼k∗<1$\tilde{\nu }_k^* &lt; 1$ for each partial transposition k.A.3System Gain CalibrationOur system gain calibration procedure consists of a measurement of Johnson–Nyquist noise spectral density emitted by a 50 Ω termination at different temperatures. Assuming perfect matching of the source and load impedances, the received power per unit of bandwidth can be written by applying the Friis formula: the measured noise is given by the noise temperature of the source Ts$T_\text{s}$, the contribution of the cooled amplifier THEMT$T_{\text{HEMT}}$, and the noise of the room‐temperature amplifiers TRT$T_{\text{RT}}$ multiplied by the system gain GΣ,i=GHEMTiGRTi$G_{\Sigma , i}=G_{\text{HEMT}_i}G_{\text{RT}_i}$:A15Ii2+Qi2Z0Δfi=kbGΣ,iTi+THEMT+TRTGHEMT$$\begin{equation} \frac{\left&lt; I_i^2 + Q_i^2 \right&gt;}{ Z_0 \Delta f_i } = k_b G_{\Sigma ,i} {\left(T_i+T_{\text{HEMT}}+\frac{T_{\text{RT}}}{G_{\text{HEMT}}}\right)} \end{equation}$$Here i refers to the frequency of the spectral mode and Δfi$\Delta f_i$ refers to the bandwidth of the detection of quadratures Ii$I_i$ and Qi$Q_i$. The total gain GΣ,i$G_{\Sigma , i}$ was separately determined for different spectral modes.Figure A1 displays the measured noise power per unit band as a function of sample temperature Ts$T_\text{s}$, averaged over frequencies covering the resonance curve. By fitting a line to the data, we obtain <GΣ,i>=94.4±0.2$&lt;G_{\Sigma , i}&gt;=94.4 \pm 0.2$ dB for the average total gain. The linear fit in Figure A1 is performed at T>0.2$T&gt;0.2$ K, which allows us to neglect the corrections from the coth(ℏω/2kbTs)$\coth (\hbar \omega /2k_\text{b} T_\text{s})$.A1FigureGain calibration using linear temperature dependence of the measured thermal noise spectral density of a 50 Ω terminator measured as a function of Ts$T_\text{s}$, the source temperature. The average total gain <GΣ,i>=94.4±0.2$&lt;G_{\Sigma , i}&gt;=94.4 \pm 0.2$ dB over the cavity resonance is obtained from the linear fit (in red) to the data. This gain value <GΣ,i>$&lt;G_{\Sigma , i}&gt;$ also includes frequency mixing losses/amplification in the signal analyzer circuit part. The term Tpreamp=THEMT+TRTGHEMT=5.2±0.25$T_{\text{preamp}}=T_{\text{HEMT}}+\frac{T_\text{RT}}{G_\text{HEMT}}=5.2\pm 0.25$K characterizes the equivalent noise temperature of the amplifiers; the largest contribution originates from the cooled HEMT amplifier at 4 K. The value coth(hfi2kbTpreamp)$\coth (\frac{hf_i}{2 k_b T_\text{preamp}})$ sets the background for the diagonal elements in the covariance matrix 4Voff$4\mathbf {V}_\text{off}$.The error in the system gain calibration results in uncertainty in the symplectic eigenvalues on the order of 2%, that is, the eigenvalues fall in the range of min{ν∼k∗}$\min \lbrace \tilde{\nu }_k^*\rbrace$=min{ν∼k∗}±0.018$\min \lbrace \tilde{\nu }_k^*\rbrace \pm 0.018$ for each partial bipartition. Random variations of the system parameters were reduced by averaging the outcome by ten to twenty times.A.4System Parameter FittingIn order to determine coupling rates γ and κ introduced in Section 2 A, we characterized our nonlinear resonator as a two‐port device using a vector network analyzer. For the characterization, we chose the optimal DC operating point ΦDC=0.383Φ0$\Phi _{\text{DC}} = 0.383\,\Phi _0$ depicted in Figure 3b. At this DC‐flux, we measured the resonance curve in the absence of the pump in order to estimate the external and internal loss rates κ and γ, respectively. By fitting the measured resonance curve to the analytical solution of the QLE (b∼out(ω)/b∼in(ω)$\tilde{b}_{\text{out}}(\omega )/\tilde{b}_{\text{in}}(\omega )$), derived for the linear case without any pump drive, we obtain the coupling coefficients κ2π=4.44$\frac{\kappa }{2\pi } = 4.44$ MHz and γ2π=2.30$\frac{\gamma }{2\pi }=2.30$ MHz. The employed analytical solution, displayed in Equation (A16), was derived from the full QLE in Equation (A8) without taking the nonlinear part −iKa∼†a∼a∼$-iK \tilde{a}^{\dagger } \tilde{a} \tilde{a}$ into accountA16b∼out(ω)b∼in(ω)=1−κ(−i(ω−ωr)+κ+γ2)$$\begin{equation} \frac{\tilde{b}_{\text{out}}(\omega )}{\tilde{b}_{\text{in}}(\omega )}=1-\frac{\kappa }{(-i(\omega -\omega _\text{r}) +\frac{\kappa +\gamma }{2})} \end{equation}$$For fitting of the Kerr constant K, we employed the whole form of the QLE in the rotating wave approximation Equation (A8). By comparing the measured and simulated gain coefficients G(ωprobe−ωr,A)$G(\omega _\text{probe}-\omega _\text{r},A)$ (Figure A2) in the cavity at large pump amplitudes, we obtain an estimate K=6.5ωr$K = 6.5\omega _\text{r}$ for the Kerr constant.A2Figurea,b) Gain coefficient measurement and QLE's simulation results (Equation (A8)) which are used for fitting coupling and Kerr constants. Vertical axis represents normalized pump amplitude A=ακ+γ$A=\frac{\alpha }{\kappa +\gamma }$. Horizontal axis represents detuning between probe signal frequency and resonance frequency ωprobe−ωr2π$\frac{\omega _\text{probe} - \omega _\text{r}}{2\pi }$. Pumping carried out in degenerate mode, ωd=2ωr$\omega _d=2\omega _\text{r}$. Fano resonance picture, given in experimental plot, explained by phase shift between the cavity and input modes and described by complex rate value of κ.A.5Cavity Phase ResponseOptimal value of GME criterion, which governed by “symmetric” covariance matrix view in tripartite mode case, can be obtained with {π2;π2}$\lbrace \frac{\pi }{2};\frac{\pi }{2}\rbrace$ only if modes reshuffling suffers no additional phase rotations (see next subsection). However, in experiments we deal with finite values of coupling and dissipation loss rates. Cavity phase response becomes crucial figure in pump tones phase shift adjustments. Cavity phase response illustrated on Figure A3.A3FigureCavity phase response given by fitted experimental parameters κ2π=4.44$\frac{\kappa }{2\pi } = 4.44$ MHz and γ2π=2.30$\frac{\gamma }{2\pi }=2.30$ MHz. Vertical dash lines show center frequencies of first and last modes. Corresponding phase shifts applied to pump tones to reach “symmetric” covariance matrix view on double frequencies are Δφ1=π2−π4=π4$\Delta \varphi _1=\frac{\pi }{2}-\frac{\pi }{4}=\frac{\pi }{4}$ and Δφ2=π2+π4=3π4$\Delta \varphi _2=\frac{\pi }{2}+\frac{\pi }{4}=\frac{3\pi }{4}$ with corresponding phase shift between pump tones π2$\frac{\pi }{2}$, given in results (see Figure 4). Half of applied phase shift described by pump tones on double resonance frequencies. Additional phase shift with increase of A relates to modification of phase response curve during pumping.A.6Multifrequency Correlations in Terms of QLE with 3 and 4 Spectral ModesAs discussed in Section 2 C, our measurement setting probes outgoing waves from the parametric resonator, which brings about slight differences with standard quantum optics schemes where the entanglement analysis is based on the Hamiltonian of the system. In our case, the QLE provides a good description, and here we derive the relevant matrix equations describing the coupling of the different outgoing spectral modes under two and three pump tones (3 and 4 spectral modes, respectively).3 Mode Case: Let us define a∼$\tilde{a}$ as a vector of spectral modes:A17a∼={a∼1,a∼2,a∼3,a∼1†,a∼2†,a∼3†}T;$$\begin{equation} \tilde{a}=\lbrace {\tilde{a}_1,\tilde{a}_2,\tilde{a}_3,\tilde{a}^\dagger _1,\tilde{a}^\dagger _2,\tilde{a}^\dagger _3\rbrace }^{T}; \end{equation}$$where the creation a∼i†=a∼i†(t)$\tilde{a}_i^{\dagger }=\tilde{a}_i^{\dagger }(t)$ and annihilation a∼i=a∼i(t)$\tilde{a}_i=\tilde{a}_i(t)$ operators are time‐dependent. After Fourier transform,A18a∼(ω)={a∼1(ω),a∼2(ω),a∼3(ω),a∼1†(−ω),a∼2†(−ω),a∼3†(−ω)}T.$$\begin{equation} \tilde{a}(\omega )=\lbrace {\tilde{a}_1(\omega ),\tilde{a}_2(\omega ),\tilde{a}_3(\omega ),\tilde{a}^\dagger _1(-\omega ),\tilde{a}^\dagger _2(-\omega ),\tilde{a}^\dagger _3(-\omega )\rbrace }^{T}. \end{equation}$$We define our spectral modes ai$a_i$ as {(−3Δ2,−Δ2);(−Δ2,Δ2);(Δ2,3Δ2)}$\lbrace (-\frac{3\Delta }{2},-\frac{\Delta }{2});(-\frac{\Delta }{2},\frac{\Delta }{2});(\frac{\Delta }{2},\frac{3\Delta }{2})\rbrace$ according pump tone positions {−Δ,Δ}$\lbrace -\Delta ,\Delta \rbrace$ (see Figure 1). Similarly, we define for the input and output modes b̂in/out$\hat{b}_{\text{in}/\text{out}}$:A19b∼in/out={b∼in1/out1,b∼in2/out2,b∼in3/out3,b∼in2/out2†,b∼in1/out1†,b∼in3/out3†T.$$\begin{eqnarray} \tilde{b}_{\text{in/out}}=& \lbrace {\tilde{b}_{\text{in}1/\text{out}1},\tilde{b}_{\text{in}2/\text{out}2},\tilde{b}_{\text{in}3/\text{out}3}}, \nonumber \\ &{\left. \tilde{b}^\dagger _{\text{in}2/\text{out}2},\tilde{b}^\dagger _{\text{in}1/\text{out}1},\tilde{b}^\dagger _{\text{in}3/\text{out}3}\right\rbrace} ^{T}. \end{eqnarray}$$The commutation relationships for the case of N modes can be conveniently expressed in matrix form. We use the common convention for [a∼i,a∼j]$[\tilde{a}_i,\tilde{a}_j ]$ from ref. [74].The effect of Kerr nonlinearity is significant only at large pump amplitudes. Hence, we may take the QLE Equation (A8) without the nonlinear part for our treatment. In theoretical analysis we assume, that spectral modes lay down deep in cavity mode, such that Δ≪κ$\Delta \ll \kappa$; we also neglect internal dissipation expressed by loss rate γ. For that case phase shift between modes, provided by phase response of the cavity, can be neglected. Guided by standard Fourier transform technique for solving linear QLE,[40] we denote a∼i(ω)=∫a∼i(t)eiωtdt$\tilde{a}_i(\omega )=\int \tilde{a}_i(t) e^{i\omega t}\text{d}t$ and Fourier transform the QLE term by term. Owing to detuning of the pump tones in the rotating frame, there will be coupling of spectral modes and we have mode index exchange. For example, for a∼1,2†$\tilde{a}^\dagger _{1,2}$: ∫a∼1†(t)eiωteiΔ1tdt=a∼2†(−ω)$\int \tilde{a}^\dagger _1(t)e^{i\omega t}e^{i\Delta _1 t}\text{d}t=\tilde{a}^{\dagger }_2(-\omega )$; ∫a∼2†(t)eiωteiΔ2tdt=a∼1†(−ω)$\int \tilde{a}^\dagger _2(t)e^{i\omega t}e^{i\Delta _2 t}\text{d}t=\tilde{a}^{\dagger }_1(-\omega )$, while for a∼2,3†$\tilde{a}^\dagger _{2,3}$: ∫a∼2†(t)eiωte−iΔ2tdt=a∼3†(−ω)$\int \tilde{a}^\dagger _2(t)e^{i\omega t}e^{-i\Delta _2 t}\text{d}t=\tilde{a}^{\dagger }_3(-\omega )$; ∫a∼3†(t)eiωte−iΔ3tdt=a∼2†(−ω)$\int \tilde{a}^\dagger _3(t)e^{i\omega t}e^{-i\Delta _3 t}\text{d}t=\tilde{a}^{\dagger }_2(-\omega )$. Thus, it is seen that each pump creates a two‐mode squeezed state (TMS) between two neighboring spectral modes independently from RWA's zero‐frequency position.After Fourier transforming,A20(−i(ω−Δr)+κ2)a∼(ω)+iα(∫a∼†(t)eiωte−iΔdtdt+∫a∼†(t)eiωteiΔdtdt)=κb∼in(ω)$$\begin{eqnarray} (-i(\omega -\Delta _r)+&\frac{\kappa }{2})\tilde{a}(\omega )+i\alpha (\int \tilde{a}^\dagger (t)e^{i\omega t}e^{-i\Delta _d t}\text{d}t+ \nonumber \\ &\int \tilde{a}^\dagger (t)e^{i\omega t}e^{i\Delta _d t}\text{d}t)=\sqrt {\kappa }\tilde{b}_{\text{in}}(\omega ) \end{eqnarray}$$the QLE yields the following system of linear equations:A21κb∼in1(ω)=(−i(ω−Δr)+κ2)a∼1(ω)+iαa∼2†(−ω)κb∼in2(ω)=(−i(ω−Δr)+κ2)a∼2(ω)+iα(a∼1†(−ω)+a∼3†(−ω))κb∼in3(ω)=(−i(ω−Δr)+κ2)a∼3(ω)+iαa∼2†(−ω)κb∼in1†(−ω)=(−i(ω+Δr)+κ2)a∼1†(−ω)−iα†a∼2(ω)κb∼in2†(−ω)=(−i(ω+Δr)+κ2)a∼2†(−ω)−iα†(a∼1(ω)+a∼3(ω))κb∼in3†(−ω)=(−i(ω+Δr)+κ2)a∼3†(−ω)−iα†a∼2(ω).$$\begin{eqnarray} \sqrt {\kappa }\tilde{b}_{\text{in}1}(\omega )=(-i(\omega -\Delta _r)+\frac{\kappa }{2})\tilde{a}_1(\omega )+i\alpha \tilde{a}^\dagger _2(-\omega ) \nonumber \\ \sqrt {\kappa }\tilde{b}_{\text{in}2}(\omega )=(-i(\omega -\Delta _r)+\frac{\kappa }{2})\tilde{a}_2(\omega )+i\alpha (\tilde{a}^\dagger _1(-\omega )+\tilde{a}^\dagger _3(-\omega )) \nonumber \\ \sqrt {\kappa }\tilde{b}_{\text{in}3}(\omega )=(-i(\omega -\Delta _r)+\frac{\kappa }{2})\tilde{a}_3(\omega )+i\alpha \tilde{a}^\dagger _2(-\omega ) \nonumber \\ \sqrt {\kappa }\tilde{b}^\dagger _{\text{in}1}(-\omega )=(-i(\omega +\Delta _r)+\frac{\kappa }{2})\tilde{a}^\dagger _1(-\omega )-i\alpha ^\dagger \tilde{a}_2(\omega ) \nonumber \\ \sqrt {\kappa }\tilde{b}^\dagger _{\text{in}2}(-\omega )=(-i(\omega +\Delta _r)+\frac{\kappa }{2})\tilde{a}^\dagger _2(-\omega )-i\alpha ^\dagger (\tilde{a}_1(\omega )+\tilde{a}_3(\omega )) \nonumber \\ \sqrt {\kappa }\tilde{b}^\dagger _{\text{in}3}(-\omega )=(-i(\omega +\Delta _r)+\frac{\kappa }{2})\tilde{a}^\dagger _3(-\omega )-i\alpha ^\dagger \tilde{a}_2(\omega ). \end{eqnarray}$$We cast Equation (A21) into matrix form:A22Ma∼(ω)=κb∼in(ω)$$\begin{equation} \mathbf {M}\tilde{a}(\omega )=\sqrt {\kappa }\tilde{b}_{\text{in}}(\omega ) \end{equation}$$A23M=c1000iα00c10iα0iα00c10iα00−iα†0c200−iα†0−iα†0c200−iα†000c2$$\begin{equation} \mathbf {M}=\def\eqcellsep{&}\begin{bmatrix} c_1 & 0 & 0 & 0 & i\alpha & 0\\ 0 & c_1 & 0 & i\alpha & 0 & i\alpha \\ 0 & 0 & c_1 & 0 & i\alpha & 0\\ 0 & -i\alpha ^\dagger & 0 & c_2 & 0 & 0\\ -i\alpha ^\dagger & 0 & -i\alpha ^\dagger & 0 & c_2 & 0\\ 0 & -i\alpha ^\dagger & 0 & 0 & 0 & c_2 \end{bmatrix} \end{equation}$$where c1=−i(ω−Δr)+κ2$c_1=-i(\omega -\Delta _r)+\frac{\kappa }{2}$ and c2=−i(ω+Δr)+κ2$c_2=-i(\omega +\Delta _r)+\frac{\kappa }{2}$ Solving for the inverse of M̂$\hat{M}$ and using Equation (A9), we obtainA24b∼out(ω)=(I−κM−1)b∼in(ω)$$\begin{equation} \tilde{b}_{\text{out}}(\omega )=(\mathbf {I}-\kappa \mathbf {M}^{-1})\tilde{b}_{\text{in}}(\omega ) \end{equation}$$for the outgoing radiation b∼out(ω)$\tilde{b}_{\text{out}}(\omega )$ in terms of incoming waves b∼in(ω)$\tilde{b}_{\text{in}}(\omega )$.Because our goal is to determine the structure of the experimental covariance matrix, it is unsatisfactory to consider cavity modes a∼$\tilde{a}$ with equation a∼(ω)=κM−1b∼in(ω)$\tilde{a}(\omega )=\sqrt {\kappa }\mathbf {M}^{-1}\tilde{b}_{\text{in}}(\omega )$ though it has a more compact final form. However, the presence of the identity matrix I$\mathbf {I}$ and the multiplication factor κ do not change the final structure.Assuming that the pump amplitude α is a real number and c1=c2=c$c_1=c_2=c$ (zero detuning case), we haveA25M−1=1c2−2α2·c−α2c0α2c0−iα00c0−iα0−iαα2c0c−α2c0−iα00iα0c−α2c0α2ciα0iα0c00iα0α2c0c−α2c$$\begin{eqnarray} \mathbf {M}^{-1}= \frac{1}{c^2-2\alpha ^2} \nonumber \cdot \\ \def\eqcellsep{&}\begin{bmatrix} c-\frac{\alpha ^2}{c} & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & -i\alpha & 0\\ 0 & c & 0 & -i\alpha & 0 & -i\alpha \\ \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} & 0 & -i\alpha & 0\\ 0 & i\alpha & 0 & c-\frac{\alpha ^2}{c} & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}}\\ i\alpha & 0 & i\alpha & 0 & c & 0\\ 0 & i\alpha & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} \end{bmatrix} \end{eqnarray}$$This allows us to draw the generalized H∼$\mathcal {\tilde{H}}$‐graph for describing the parametric interaction between the spectral modes, Figure 2. The off‐diagonal beamsplitter elements proportional to α2 are set in bold in Equation (A25).Still, we want to construct the parametric interaction matrix S−1$\mathbf {S}^{-1}$ for quadrature vector operator r∼$\tilde{r}$. Using a linear operator matrix K$\mathbf {K}$ to implement a change of basisA26K=12100100−i00i000100100−i00i000100100−i00i$$\begin{eqnarray} \mathbf {K}= \frac{1}{2} \def\eqcellsep{&}\begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0\\ -i & 0 & 0 & i & 0 & 0\\ 0 & 1 & 0 & 0 & 1 & 0\\ 0 & -i & 0 & 0 & i & 0 \\ 0 & 0 & 1 & 0 & 0 & 1\\ 0 & 0 & -i & 0 & 0 & i \end{bmatrix} \end{eqnarray}$$we obtain by a canonical transformation S−1=κKM−1K−1$\mathbf {S}^{-1}=\sqrt {\kappa }\mathbf {K}\mathbf {M}^{-1}\mathbf {K}^{-1}$:A27S−1=κc2−2α2·c+α2c00−αα2c00c+α2c−α00α2c0−αc00−α−α00c−α0α2c00−αc+α2c00α2c−α00c+α2c$$\begin{eqnarray} &&\mathbf {S}^{-1}= \frac{\sqrt {\kappa }}{c^2-2\alpha ^2} \nonumber \cdot \\ &&\def\eqcellsep{&}\begin{bmatrix} c+\frac{\alpha ^2}{c} & 0 & 0 & -\alpha & \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 \\ 0 & c+\frac{\alpha ^2}{c} & -\alpha & 0 & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}}\\ 0 & -\alpha & c & 0 & 0 & -\alpha \\ -\alpha & 0 & 0 & c & -\alpha & 0 \\ \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & 0 & -\alpha & c+\frac{\alpha ^2}{c} & 0\\ 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}} & -\alpha & 0 & 0 & c+\frac{\alpha ^2}{c} \end{bmatrix} \end{eqnarray}$$Note, that here the overall structure of the matrix has changed because of the basis change from ladder to quadrature operators. This is seen, for example, in the distribution of the off‐diagonal beamsplitter correlations (shown in bold).Since the environment of the cavity is in the ground state, b∼in$\tilde{b}_\text{in}$ has a Gaussian covariance matrix of the form Vin=14I$\mathbf {V}_\text{in}=\frac{1}{4}\mathbf {I}$. Consequently, the covariance matrix of the cavity spectral modes a∼i$\tilde{a}_i$ can be represented asA28Va=S−1Vin(S−1)T$$\begin{equation} \mathbf {V}_\text{a}=\mathbf {S}^{-1}\mathbf {V}_\text{in}(\mathbf {S}^{-1})^{T} \end{equation}$$or, equivalently, for output modes b∼out$\tilde{b}_\text{out}$: Vout=(I−κS−1)Vin(I−κS−1)T$\mathbf {V}_\text{out}= (\mathbf {I}-\sqrt {\kappa }\mathbf {S}^{-1})\mathbf {V}_\text{in}(\mathbf {I}-\sqrt {\kappa }\mathbf {S}^{-1})^{T}$. Both forms Va$\mathbf {V}_\text{a}$ and Vout$\mathbf {V}_\text{out}$ can be employed for studying the structure of parametric interactions between the quadratures, because input–output relationship doesnot change the general structure of the couplings between the quadratures (see below).As shown in Section 4 experimentally, phase shift between pumps changes the appearance of the covariance matrix (see Figure 5) as well as the strength of genuine multipartite entanglement. A change in the matrix M$\mathbf {M}$ due to a phase shift is illustrated in Equation (A29), in which the phase of the first pump has been rotated by eiπ2$e^{\frac{i\pi }{2}}$.A29M̂=c000−α00c0−α0iα00c0iα00−α0c00−α0−iα0c00−iα000c$$\begin{equation} \hat{M}=\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & -\pmb {\alpha } & 0\\ 0 & c & 0 & -\pmb {\alpha } & 0 & i\alpha \\ 0 & 0 & c & 0 & i\alpha & 0\\ 0 & -\pmb {\alpha } & 0 & c & 0 & 0\\ -\pmb {\alpha } & 0 & -i\alpha & 0 & c & 0\\ 0 & -i\alpha & 0 & 0 & 0 & c \end{bmatrix} \end{equation}$$The elements affected by the rotation are indicated in bold in the matrix. The elements in bold face indicate coupling between modes a∼1(ω)↔a∼2(ω)$\tilde{a}_1(\omega )\leftrightarrow \tilde{a}_2(\omega )$ while the other off‐diagonal elements indicate squeezing across a∼2(ω)↔a∼3(ω)$\tilde{a}_2(\omega )\leftrightarrow \tilde{a}_3(\omega )$. Note, that phase rotation operates in opposite direction on rows related to b∼in(ω)$\tilde{b}_{\text{in}}(\omega )$ and b∼in†(ω)$\tilde{b}_{\text{in}}^{\dagger }(\omega )$.The inversion of the rotated matrix M$\mathbf {M}$ yields for the parametric interaction matrix, where all the beamsplitter elements (in bold) have acquired a π/2$\pi /2$ phase shift. This phase shift can be unwound by a phase shift on the second pump, which indicates different phase dependence of the beamsplitter correlations compared with the TMS correlations. The structure of matrix M−1$\mathbf {M}^{-1}$ in Equation (A30) shows that phase rotation of specified pump tones does not change parametric interaction form between modes, preserving structure of a bisqueezed tripartite state. However, as shown in the main text, the criterion describing the strength of GME (see Equation (15)) depends on the difference of pump phases and strong genuine entanglement is reached only at specific phase settings.A30M−1=1c2−2α2·c−α2c0iα2c0α00c0α0−iα−iα2c0c−α2c0−iα00α0c−α2c0−iα2cα0iα0c00iα0iα2c0c−α2c$$\begin{eqnarray} \mathbf {M}^{-1}= \frac{1}{c^2-2\alpha ^2}\cdot \nonumber \\ \def\eqcellsep{&}\begin{bmatrix} c-\frac{\alpha ^2}{c} & 0 & \mathbf {\frac{\pmb {i\alpha }^2}{c}} & 0 & \alpha & 0\\ 0 & c & 0 & \alpha & 0 & -i\alpha \\ -\mathbf {\frac{\pmb {i\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} & 0 & -i\alpha & 0\\ 0 & \alpha & 0 & c-\frac{\alpha ^2}{c} & 0 & -\mathbf {\frac{\pmb {i\alpha }^2}{c}}\\ \alpha & 0 & i\alpha & 0 & c & 0\\ 0 & i\alpha & 0 & \mathbf {\frac{\pmb {i\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} \end{bmatrix} \end{eqnarray}$$The covariance matrix Va$\mathbf {V}_\text{a}$ obtained from matrix M$\mathbf {M}$ in Equation (A25) with zero pump phase shifts is given in Equation (A31). The corresponding covariance matrix for π/2$\pi /2$ phase rotation in the first pump is displayed in Equation (A32). The matrix Va$\mathbf {V}_\text{a}$ in Equation (A32) has one rotated subspace, corresponding to two quadrature pairs; these rotated components are indicated by bold face. Based on these analytical relationships we conclude that control over desired covariance matrix TMS‐subspace can be provided by phase rotation of corresponding pump tone. Finally, we introduce the same phase rotation eiπ2$e^{\frac{i\pi }{2}}$ to the second pump. This brings the covariance matrix for the spectral cavity modes to the “standard‐symmetric” form displayed in Equation (A33). By comparing Equation (A31) and Equation (A33) we note that the beamsplitter elements (in bold) in the covariance matrix are unchanged (the phase difference between the pumps is the same) while the TMS elements are different.A31Va=κ4(c2−2α2)2·−α2+2α4c2+c200−2αc3α2−2α4c200−α2+2α4c2+c2−2αc003α2−2α4c20−2αc2α2+c200−2αc−2αc002α2+c2−2αc03α2−2α4c200−2αc−α2+2α4c2+c2003α2−2α4c2−2αc00−α2+2α4c2+c2$$\begin{equation} \mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \def\eqcellsep{&}\begin{bmatrix} -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & 0 & -2 \alpha c & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 \\ 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & -2 \alpha c & 0 & 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} \\ 0 & -2 \alpha c & 2 \alpha ^2+c^2 & 0 & 0 & -2 \alpha c \\ -2 \alpha c & 0 & 0 & 2 \alpha ^2+c^2 & -2 \alpha c & 0 \\ \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & 0 & -2 \alpha c & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 \\ 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & -2 \alpha c & 0 & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 \\ \end{bmatrix} \end{equation}$$A32Va=κ4(c2−2α2)2·−α2+2α4c2+c202αc00−3α2+2α4c20−α2+2α4c2+c20−2αc3α2−2α4c202αc02α2+c200−2αc0−2αc02α2+c2−2αc003α2−2α4c20−2αc−α2+2α4c2+c20−3α2+2α4c20−2αc00−α2+2α4c2+c2$$\begin{equation} \mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \def\eqcellsep{&}\begin{bmatrix} -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & \pmb {2 \alpha c} & 0 & 0 & \pmb {-3 \alpha ^2+\frac{2 \alpha ^4}{c^2}} \\ 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & \pmb {-2 \alpha c} & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 \\ \pmb {2 \alpha c} & 0& 2 \alpha ^2+c^2 & 0 & 0 & -2 \alpha c \\ 0 & \pmb {-2 \alpha c} & 0 & 2 \alpha ^2+c^2 & -2 \alpha c & 0 \\ 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & -2 \alpha c & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 \\ \pmb {-3 \alpha ^2+\frac{2 \alpha ^4}{c^2}} & 0 &-2 \alpha c & 0 & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 \\ \end{bmatrix} \end{equation}$$A33Va=κ4(c2−2α2)2·−α2+2α4c2+c202αc03α2−2α4c200−α2+2α4c2+c20−2αc03α2−2α4c22αc02α2+c202αc00−2αc02α2+c20−2αc3α2−2α4c202αc0−α2+2α4c2+c2003α2−2α4c20−2αc0−α2+2α4c2+c2.$$\begin{equation} \mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \def\eqcellsep{&}\begin{bmatrix} -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 \\ 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c & 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} \\ 2 \alpha c & 0 & 2 \alpha ^2+c^2 & 0 & 2 \alpha c & 0 \\ 0 & -2 \alpha c & 0 & 2 \alpha ^2+c^2 & 0 & -2 \alpha c \\ \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & 2 \alpha c & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 \\ 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & -2 \alpha c & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 \\ \end{bmatrix}. \end{equation}$$4 Mode Case: We treat the case with three pump tones (p=3$p=3$) in the same way as we did with two pumps above. By taking ωΣ=∑d=1pωdp=2ωr$\omega _{\Sigma }=\frac{\sum _{d=1}^{p}\omega _d}{p} = 2 \omega _\text{r}$, we have a configuration of pump tones {−2Δ,0,2Δ}$\lbrace {-2\Delta ,0,2\Delta \rbrace }$ with respect to 2ωr$2 \omega _\text{r}$. We define our spectral modes around ωr$\omega _\text{r}$ as {(−2Δ,−Δ);(−Δ,0);(0,Δ);(Δ,2Δ)}$\lbrace (-2\Delta ,-\Delta );(-\Delta ,0);(0,\Delta );(\Delta ,2\Delta )\rbrace$ =̂{a∼1(ω);a∼2(ω);a∼3(ω);a∼4(ω)}$\widehat{=} \lbrace \tilde{a}_1(\omega ); \tilde{a}_2(\omega ); \tilde{a}_3(\omega ); \tilde{a}_4(\omega ) \rbrace$. Individual pump tones at ω1 and ω3 create TMS states between neighboring spectral modes as in the two pump case above. However, the middle pump creates TMS correlations between two pairs of spectral modes: a∼1(ω)↔a∼4†(−ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}^\dagger _4(-\omega )$ and a∼2(ω)↔a∼3†(−ω)$\tilde{a}_2(\omega ) \leftrightarrow \tilde{a}^\dagger _3(-\omega )$. Rotation of the middle pump phase has an effect on both corresponding subspaces of the covariance matrix. Consequently, this pump configuration is quite suitable for producing square H∼$\mathcal {\tilde{H}}$‐graph states (see Figure 2).A34M=c0000−α0−α0c00−α0−α000c00−α0−α000c−α0−α00−α0−αc000−α0−α00c000−α0−α00c0−α0−α0000c$$\begin{align} \mathbf {M}=\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & 0 & -\alpha & 0 & \pmb {-\alpha }\\ 0 & c & 0 & 0 & -\alpha & 0 & \pmb {-\alpha } & 0\\ 0 & 0 & c & 0 & 0 & \pmb {-\alpha }& 0 & -\alpha \\ 0 & 0 & 0 & c & \pmb {-\alpha }& 0 & -\alpha & 0\\ 0 & -\alpha & 0 & \pmb {-\alpha }& c & 0 & 0 & 0\\ -\alpha & 0 & \pmb {-\alpha } & 0 & 0 & c & 0 & 0\\ 0 & \pmb {-\alpha }& 0 & -\alpha & 0 & 0 & c & 0\\ \pmb {-\alpha }& 0 & -\alpha & 0 & 0 & 0 & 0 & c\\ \end{bmatrix} \end{align}$$The parametric interactions in the covariance matrix can be analyzed in the same way as above, but now the number of phase differences influencing the beamsplitter correlations has increased. The system of linear equations ban be written as for three modes in Equation (A21), but we skip it and write down the interaction matrix M$\mathbf {M}$ (Equation (A34)), where all c coefficient are equal since we have assumed Δr=ωr−ωΣ2=0$\Delta _\text{r}=\omega _\text{r}-\frac{\omega _{\Sigma }}{2}=0$. The signs of α's are governed by the choice of pump phases as {αeiπ2;αeiπ2;αeiπ2}$ \lbrace { \alpha e^{\frac{i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}\rbrace }$. The correlations produced by the pump at ω2=2ωr$\omega _2=2\omega _\text{r}$ are indicated in bold. The special role of the central pump is seen because its correlations cover the whole ascending diagonal.The inverse matrix M−1$\mathbf {M}^{-1}$ reveals the beamsplitter correlations between a∼1(ω)↔a∼3(ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}_3(\omega )$ and a∼2(ω)↔a∼4(ω)$\tilde{a}_2(\omega ) \leftrightarrow \tilde{a}_4(\omega )$:A35M−1=1(c2−4α2)·c−2α2c02α2c00α0α0c−2α2c02α2cα0α02α2c0c−2α2c00α0α02α2c0c−2α2cα0α00α0αc−2α2c02α2c0α0α00c−2α2c02α2c0α0α2α2c0c−2α2c0α0α002α2c0c−2α2c$$\begin{eqnarray} &&\mathbf {M}^{-1}=\frac{1}{(c^2-4\alpha ^2)}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} {{ \def\eqcellsep{&}\begin{array}{cccccccc}c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} & 0 & 0 & \alpha & 0 & \alpha \\[3pt] 0 & c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} & \alpha & 0 & \alpha & 0 \\[3pt] \pmb {\frac{2 \alpha ^2}{c}} & 0 & c-\frac{2 \alpha ^2}{c} & 0 & 0 & \alpha & 0 & \alpha \\[3pt] 0 & \pmb {\frac{2 \alpha ^2}{c}}& 0 & c-\frac{2 \alpha ^2}{c} & \alpha & 0 & \alpha & 0 \\[3pt] 0 & \alpha & 0 & \alpha & c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} & 0 \\[3pt] \alpha & 0 & \alpha & 0 & 0 & c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} \\[3pt] 0 & \alpha & 0 & \alpha & \pmb {\frac{2 \alpha ^2}{c}} & 0 & c-\frac{2 \alpha ^2}{c} & 0 \\[3pt] \alpha & 0 & \alpha & 0 & 0 & \pmb {\frac{2 \alpha ^2}{c}}& 0 & c-\frac{2 \alpha ^2}{c} \\[3pt] \end{array} }} \end{bmatrix}\nonumber\\ \end{eqnarray}$$The beamsplitter correlations are indicated in bold in this matrix M−1$\mathbf {M}^{-1}$. We see that there are two sequences of pump transformations that yield BS correlations between modes a∼1(ω)↔a∼2(ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}_2(\omega )$ and a∼3(ω)↔a∼4(ω)$\tilde{a}_3(\omega ) \leftrightarrow \tilde{a}_4(\omega )$. This agrees with the simple argument that indicates BS correlations to exist when two spectral bands are connected across squeezing action by two pumps with one joint frequency. Higher order correlations via three pumps exist also, but these are neglected in our analysis. Note that the number of beamsplitter correlations also coincides with the independent number of phase differences between the pumps. Connection of the cavity spectral mode correlations to H∼$\mathcal {\tilde{H}}$‐graphs is illustrated in Figure 2.The beamsplitter correlations are prominent also in the covariance matrix Va$\mathbf {V}_\text{a}$ (see Equation (A28)):A36Va=κ4(c2−4α2)2·−2α2+8α4c2+c202αc06α2−8α4c202αc00−2α2+8α4c2+c20−2αc06α2−8α4c20−2αc2αc0−2α2+8α4c2+c202αc06α2−8α4c200−2αc0−2α2+8α4c2+c20−2αc06α2−8α4c26α2−8α4c202αc0−2α2+8α4c2+c202αc006α2−8α4c20−2αc0−2α2+8α4c2+c20−2αc2αc06α2−8α4c202αc0−2α2+8α4c2+c200−2αc06α2−8α4c20−2αc0−2α2+8α4c2+c2.$$\begin{eqnarray} &&\mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-4\alpha ^2)^2}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} {{ \def\eqcellsep{&}\begin{array}{cccccccc}-2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} }& 0 & 2 \alpha c & 0 \\[3pt] 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & -2 \alpha c \\[3pt] 2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 \\[3pt] 0 & -2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } \\[3pt] \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & 2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 \\[3pt] 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & -2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c \\[3pt] 2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & 2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 \\[3pt] 0 & -2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & -2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 \end{array} }} \end{bmatrix}.\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$Bold font marks the beamsplitter correlations which display a different structure in comparison to Equation (A35) owing to the base change to quadratures ordered as (x∼1,p∼1,⋯x∼N,p∼N)T$(\tilde{x}_1, \tilde{p}_1, \dots \tilde{x}_N, \tilde{p}_N )^{T}$. So the BS correlations are between quadratures of a∼1(ω)↔a∼3(ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}_3(\omega )$ and a∼2(ω)↔a∼4(ω)$\tilde{a}_2(\omega ) \leftrightarrow \tilde{a}_4(\omega )$.By choosing the phase of the first pump to be opposite to that of the second and the third {αe−iπ2;αeiπ2;αeiπ2}$ \lbrace { \alpha e^{\frac{-i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}\rbrace }$ we are able to flip the sign of one minor diagonal indicated by bold font in Equation (A37).A37M=c0000α0−α0c00α0−α000c00−α0−α000c−α0−α00α0−αc000α0−α00c000−α0−α00c0−α0−α0000c$$\begin{equation} \mathbf {M}=\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & 0 & \pmb {\alpha } & 0 & -\alpha \\ 0 & c & 0 & 0 & \pmb {\alpha } & 0 & -\alpha & 0\\ 0 & 0 & c & 0 & 0 & -\alpha & 0 & -\alpha \\ 0 & 0 & 0 & c & -\alpha & 0 & -\alpha & 0\\ 0 & \pmb {\alpha } & 0 & -\alpha & c & 0 & 0 & 0\\ \pmb {\alpha } & 0 & -\alpha & 0 & 0 & c & 0 & 0\\ 0 & -\alpha & 0 & -\alpha & 0 & 0 & c & 0\\ -\alpha & 0 & -\alpha & 0 & 0 & 0 & 0 & c\\ \end{bmatrix} \end{equation}$$Interestingly, this choice of phases leads to full cancellation of the beamsplitter correlation terms. This is seen in the structure of the inverse matrix:A38M−1=1(c2−2α2)·c0000−α0α0c00−α0α000c00α0α000cα0α00−α0αc000−α0α00c000α0α00c0α0α0000c$$\begin{eqnarray} &&\mathbf {M}^{-1}=\frac{1}{(c^2-2\alpha ^2)}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & 0 & -\alpha & 0 & \alpha \\ 0 & c & 0 & 0 & -\alpha & 0 & \alpha & 0\\ 0 & 0 & c & 0 & 0 & \alpha & 0 & \alpha \\ 0 & 0 & 0 & c & \alpha & 0 & \alpha & 0\\ 0 & -\alpha & 0 & \alpha & c & 0 & 0 & 0\\ -\alpha & 0 & \alpha & 0 & 0 & c & 0 & 0\\ 0 & \alpha & 0 & \alpha & 0 & 0 & c & 0\\ \alpha & 0 & \alpha & 0 & 0 & 0 & 0 & c\\ \end{bmatrix} \end{eqnarray}$$which does not have any elements proportional to α2. Without any beamsplitter correlations, Equation (A38) indicates a clear connection to a square‐lattice H∼$\mathcal {\tilde{H}}$‐graph.Finally, the corresponding covariance matrix Va$\mathbf {V}_\text{a}$ for cavity spectral modes is given by Equation (A39). This structure for the covariance matrix is obtained when all the pump signals have an additional phase shift of eiπ2$e^{\frac{i\pi }{2}}$. Such a choice of phases will result in a covariance matrix with “symmetric” structure as shown in experimental data in Figure 8a,b. By controlling the phases of the pump tones separately, we can rotate and adjust certain subspaces of the 8 × 8 covariance matrix. In particular, the influence of the beamsplitter correlations can be eliminated from Va$\mathbf {V}_\text{a}$ in the four pump case.Regarding the quadripartite covariance matrix structures, the relative phase shift between the pump tones are not influenced by the cavity response in the limit of vanishing band widths or with the assumption of huge coupling rate and tiny internal dissipation loss rate. However, additional phase shifts will appear if these conditions are not met, which has to be taken into account in the generation of the desired entangled states.In principle, it would be possible to evaluate the criteria for GME from the analytical expressions derived in this Appendix (see e.g. Equations (A33) and (A39)). However, we leave the conclusions about genuine entanglement both in the tripartite and quadripartite case for analysis based on numerical simulations where even the nonlinear terms can be taken into account. The nonlinear terms are of central importance when increasing the pump drive past the critical pumping amplitude.A39Va=κ4(c2−2α2)2·c2+2α20−2αc0002αc00c2+2α202αc000−2αc−2αc0c2+2α202αc00002αc0c2+2α20−2αc00002αc0c2+2α202αc0000−2αc0c2+2α20−2αc2αc0002αc0c2+2α200−2αc000−2αc0c2+2α2.$$\begin{eqnarray} &&\mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} {{ \def\eqcellsep{&}\begin{array}{cccccccc}c^2+2 \alpha ^2 & 0 & -2 \alpha c & 0 & 0 & 0 & 2 \alpha c & 0 \\[3pt] 0 & c^2+2 \alpha ^2 & 0 & 2 \alpha c & 0 & 0 & 0 & -2 \alpha c \\[3pt] -2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & 2 \alpha c & 0 & 0 & 0 \\[3pt] 0 & 2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & -2 \alpha c & 0 & 0 \\[3pt] 0 & 0 & 2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & 2 \alpha c & 0 \\[3pt] 0 & 0 & 0 & -2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & -2 \alpha c \\[3pt] 2 \alpha c & 0 & 0 & 0 & 2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 \\[3pt] 0 & -2 \alpha c & 0 & 0 & 0 & -2 \alpha c & 0 & c^2+2 \alpha ^2 \end{array} }} \end{bmatrix}.\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$AcknowledgementsThe authors thank Mark Dykman, Stefano Pirandola, Olivier Pfister, Aidar Sultanov, Shruti Dogra, and Marko Kuzmanovié for fruitful discussions and useful comments. This work was supported by the Academy of Finland through the Finnish Center of Excellence in Quantum Technology QTF (projects 312295, 312296, 336810, 336813). This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 862644 (FET‐Open project QUARTET) and grant agreement no. 824109 (European Microkelvin Platform), as well as ERC grant agreement no. 670743 (QuDeT). G.S.P and K.V.P. also would like to thank Saab for support, under a research agreement between Saab and Aalto University. The work of PJH and IL was supported by MATINE grant VN/13774/2020‐PLM‐49.Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.R. Jozsa, N. Linden, Proc. R. Soc. London, Ser. A 2003, 459, 2011.R. Horodecki, P. Horodecki, M. Horodecki, K. Horodecki, Rev. Mod. Phys. 2009, 81, 865.R. Raussendorf, H. J. Briegel, Phys. Rev. Lett. 2001, 86, 5188.N. C. Menicucci, S. T. Flammia, O. Pfister, Phys. Rev. Lett. 2008, 101, 130501.T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, J. L. O'Brien, Nature 2010, 464, 45.N. C. Menicucci, Phys. Rev. Lett. 2014, 112, 120504.O. Pfister, J. Phys. B: At. Mol. Opt. Phys. 2019, 53, 012001.40 years of quantum computing, Nat. Rev. Phys. 2022, 4, https://doi.org/10.1038/s42254-021-00410-6.V. Giovannetti, S. Lloyd, L. Maccone, Nat. Photonics 2011, 5, 222.Z. Zhang, S. Mouradian, F. N. C. Wong, J. H. Shapiro, Phys. Rev. Lett. 2015, 114, 110506.C. L. Degen, F. Reinhard, P. Cappellaro, Rev. Mod. Phys. 2017, 89, 035002.L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, P. Treutlein, Rev. Mod. Phys. 2018, 90, 035005.X. Su, J. Jing, Q. Pan, C. Xie, Phys. Rev. A 2006, 74, 062305.N. Gisin, R. Thew, Nat. Photonics 2007, 1, 165.J. Yin, Y.‐H. Li, S.‐K. Liao, M. Yang, Y. Cao, L. Zhang, J.‐G. Ren, W.‐Q. Cai, W.‐Y. Liu, S.‐L. Li, R. Shu, Y.‐M. Huang, L. Deng, L. Li, Q. Zhang, N.‐L. Liu, Y.‐A. Chen, C.‐Y. Lu, X.‐B. Wang, F. Xu, J.‐Y. Wang, Ch.‐Z. Peng, A. K. Ekert, J.‐W. Pan, Nature 2020, 582, 501.S. J. Summers, R. Werner, Phys. Lett. A 1985, 110, 257.S. L. Braunstein, P. van Loock, Rev. Mod. Phys. 2005, 77, 513.P. Lähteenmäki, G. S. Paraoanu, J. Hassel, P. J. Hakonen, Nat. Commun. 2016, 7, 12548.A. M. Lance, T. Symul, W. P. Bowen, B. C. Sanders, T. Tyc, T. C. Ralph, P. K. Lam, Phys. Rev. A 2005, 71, 033814.M. Hillery, V. Bužek, A. Berthiaume, Phys. Rev. A 1999, 59, 1829.S. L. Braunstein, H. J. Kimble, Phys. Rev. A 2000, 61, 042302.J. Dias, T. C. Ralph, Phys. Rev. A 2018, 97, 032335.G. J. Milburn, S. L. Braunstein, Phys. Rev. A 1999, 60, 937.X. Guo, C. R. Breum, J. Borregaard, S. Izumi, M. V. Larsen, T. Gehring, M. Christandl, J. S. Neergaard‐Nielsen, U. L. Andersen, Nat. Phys. 2020, 16, 281.R. Nichols, P. Liuzzo‐Scorpo, P. A. Knott, G. Adesso, Phys. Rev. A 2018, 98, 012114.M. Gessner, A. Smerzi, L. Pezzè, Nat. Commun. 2020, 11, 3817.J. Zhang, S. L. Braunstein, Phys. Rev. A 2006, 73, 032318.N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, M. A. Nielsen, Phys. Rev. Lett. 2006, 97, 110501.M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, P. van Loock, Phys. Rev. A 2009, 79, 062318.A. Mantri, T. F. Demarie, J. F. Fitzsimons, Sci. Rep. 2017, 7, 42861.P. Walther, K. J. Resch, T. Rudolph, E. Schenck, H. Weinfurter, V. Vedral, M. Aspelmeyer, A. Zeilinger, Nature 2005, 434, 169.S. Zippilli, D. Vitali, Phys. Rev. A 2020, 102, 052424.S. Zippilli, D. Vitali, Phys. Rev. Lett. 2021, 126, 020402.C. W. Sandbo Chang, M. Simoen, J. Aumentado, C. Sabín, P. Forn‐Díaz, A. M. Vadiraj, F. Quijandría, G. Johansson, I. Fuentes, C. M. Wilson, Phys. Rev. Appl. 2018, 10, 044019.G. Adesso, F. Illuminati, J. Phys. A: Math. Theor. 2007, 40, 7821.C. Weedbrook, S. Pirandola, R. García‐Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, S. Lloyd, Rev. Mod. Phys. 2012, 84, 621.C. Fabre, N. Treps, Rev. Mod. Phys. 2020, 92, 035005.N. C. Menicucci, S. T. Flammia, P. van Loock, Phys. Rev. A 2011, 83, 042335.D. E. Bruschi, C. Sabín, P. Kok, G. Johansson, P. Delsing, I. Fuentes, Sci. Rep. 2016, 6, 18349.T. Yamamoto, K. Koshino, Y. Nakamura, in Parametric Amplifier and Oscillator Based on Josephson Junction Circuitry, Springer, Japan, Tokyo 2016, pp. 495–513.Z. Wang, M. Pechal, E. A. Wollack, P. Arrangoiz‐Arriola, M. Gao, N. R. Lee, A. H. Safavi‐Naeini, Phys. Rev. X 2019, 9, 021049.T. Korkalainen, I. Lilja, M. R. Perelshtein, K. V. Petrovnin, G. S. Paraoanu, P. J. Hakonen, AIP Conf. Proc. 2021, 2362, 030001.P. Lähteenmäki, G. Paraoanu, J. Hassel, P. J. Hakonen, PNAS 2013, 110, 4234.H. J. Briegel, R. Raussendorf, Phys. Rev. Lett. 2001, 86, 910.D. E. Bruschi, C. Sabín, G. S. Paraoanu, Phys. Rev. A 2017, 95, 062324.G. Adesso, S. Ragy, A. R. Lee, Open Syst. Inf. Dyn. 2014, 21, 1440001.P. van Loock, A. Furusawa, Phys. Rev. A 2003, 67, 052315.M. Hillery, H. T. Dung, H. Zheng, Phys. Rev. A 2010, 81, 062322.A. Agustí, C. W. S. Chang, F. Quijandría, G. Johansson, C. M. Wilson, C. Sabín, Phys. Rev. Lett. 2020, 125, 020502.R. Simon, Phys. Rev. Lett. 2000, 84, 2726.L.‐M. Duan, G. Giedke, J. I. Cirac, P. Zoller, Phys. Rev. Lett. 2000, 84, 2722.S. Mancini, V. Giovannetti, D. Vitali, P. Tombesi, Phys. Rev. Lett. 2002, 88, 120401.P. Hyllus, J. Eisert, New J. Phys. 2006, 8, 51.R. Y. Teh, M. D. Reid, Phys. Rev. A 2014, 90, 062337.E. Shchukin, P. van Loock, Phys. Rev. A 2015, 92, 042328.A. Peres, Phys. Rev. Lett. 1996, 77, 1413.M. Horodecki, P. Horodecki, R. Horodecki, Phys. Lett. A 1996, 223, 1.L. K. Shalm, D. R. Hamel, Z. Yan, C. Simon, K. J. Resch, T. Jennewein, Nat. Phys. 2013, 9, 19.T. Elo, T. S. Abhilash, M. R. Perelshtein, I. Lilja, E. V. Korostylev, P. J. Hakonen, Appl. Phys. Lett. 2019, 114, 152601.C. Macklin, K. O'Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, I. Siddiqi, Science 2015, 350, 307.A. Zorin, Phys. Rev. Appl. 2019, 12, 044051.M. Perelshtein, K. Petrovnin, V. Vesterinen, S. Hamedani Raja, I. Lilja, M. Will, A. Savin, S. Simbierowicz, R. Jabdaraghi, J. Lehtinen, L. Grönberg, J. Hassel, M. Prunnila, J. Govenius, G. Paraoanu, P. Hakonen, Phys. Rev. Appl. 2022, 18, 024063.M. Esposito, A. Ranadive, L. Planat, S. Leger, D. Fraudet, V. Jouanny, O. Buisson, W. Guichard, C. Naud, J. Aumentado, F. Lecocq, N. Roch, Phys. Rev. Lett. 2022, 128, 153603.J. Y. Qiu, A. Grimsmo, K. Peng, B. Kannan, B. Lienhard, Y. Sung, P. Krantz, V. Bolkhovsky, G. Calusine, D. Kim, A. Melville, B. M. Niedzielski, J. Yoder, M. E. Schwartz, T. P. Orlando, I. Siddiqi, S. Gustavsson, K. P. O'Brien, W. D. Oliver, arXiv:2201.11261, 2022.T. Roy, S. Kundu, M. Chand, A. M. Vadiraj, A. Ranadive, N. Nehra, M. P. Patankar, J. Aumentado, A. A. Clerk, R. Vijay, Appl. Phys. Lett. 2015, 107, 262601.R. N. Alexander, N. C. Menicucci, Phys. Rev. A 2016, 93, 062326.R. Pooser, J. Jing, Phys. Rev. A 2014, 90, 043841.M. V. Larsen, J. S. Neergaard‐Nielsen, U. L. Andersen, Phys. Rev. A 2020, 102, 042608.J. Wu, Q. Zhuang, Phys. Rev. Appl. 2021, 15, 034073.B.‐H. Wu, R. N. Alexander, S. Liu, Z. Zhang, Phys. Rev. Res. 2020, 2, 023138.T. Hillmann, F. Quijandría, G. Johansson, A. Ferraro, S. Gasparinetti, G. Ferrini, Phys. Rev. Lett. 2020, 125, 160501.Y. Y. Gao, B. J. Lester, K. S. Chou, L. Frunzio, M. H. Devoret, L. Jiang, S. Girvin, R. J. Schoelkopf, Nature 2019, 566, 509.S. Lloyd, Science 2008, 321, 1463.R. Simon, N. Mukunda, B. Dutta, Phys. Rev. A 1994, 49, 1567.J. Williamson, Am. J. Math. 1936, 58, 141. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Advanced Quantum Technologies Wiley

Generation and Structuring of Multipartite Entanglement in a Josephson Parametric System

Download PDF
19 pages

Loading next page...
 
/lp/wiley/generation-and-structuring-of-multipartite-entanglement-in-a-josephson-wxe2qVENFq
Publisher
Wiley
Copyright
© 2023 Wiley‐VCH GmbH
eISSN
2511-9044
DOI
10.1002/qute.202200031
Publisher site
See Article on Publisher Site

Abstract

IntroductionEntanglement is a crucial resource in advanced information processing based on quantum mechanical concepts.[1,2] To exceed the computational power of classical devices, quantum devices typically need to employ highly entangled states. Controlled generation of entanglement is the key resource not only in quantum computing,[3–8] but also in sensing[9–12] and secure communications.[13–15]One of the most easily accessible and, at the same time, reliable sources for entanglement generation is the vacuum state of a quantum field.[16] Squeezing, which is the fundamental operation for the continuous variable (CV) states production, allows generation of coherence and entanglement from vacuum fluctuations.[17] While two‐mode squeezing produces bipartite quantum states, multipartite states can be generated by applying similar operations.[18] Intriguingly, multipartite CV states are shown to enable various promising phenomena such as quantum state sharing[19] and secret sharing,[20] dense coding,[21] error correction,[22] and quantum teleportation.[23] Alongside with phase sensing[24] and quantum sensor networks,[25] multipartite entangled states have significant potential in multiparameter quantum metrology applications.[25,26] Furthermore, CV cluster states show potential as a universal quantum computing platform,[27–29] which has been under active development for the past 20 years. Cluster state calculus, foremost utilizing optical resources,[28,30–33] realize measurement‐based quantum computing.[3]While optical‐mode schemes for generation of multipartite states lack versatility, in‐situ tunability and are limited to optical frequencies, the microwave platform allows for full control of operations via input RF‐signals and integration with the existing silicon‐based circuitry. During the past few years, significant progress in processing of CV multipartite states at microwaves has been achieved; for instance, squeezed states produced by microwave cavities have been shown to exhibit correlations between photons in separate frequency bands[18] and strong entanglement between different modes.[34]In this work, we experimentally generate genuinely entangled tripartite and quadripartite states using a superconducting parametric cavity, operated under steady‐state conditions. Using the Gaussian‐mode formalism,[35–37] we characterize the generated states and verify entanglement from the covariance matrix. We develop an analytical description, which allows us to determine the entanglement structure of the generated state and establish a connection to H(amiltonian)‐graph representation.[7,29,38,39] All of the experimental results are in good agreement with the theoretical predictions, in which all circuit parameters were set in accordance with the measured characteristics of the device.The paper is organized as follows. Section 2 describes the quantum dynamics of a Josephson parametric system and demonstrates the generation protocol for CV multipartite states, such as fully inseparable and genuinely entangled states. Using analytical methods, we provide the entanglement structure, which is described by graphs and their corresponding adjacency matrices. In Section 3 we present the experimental setup and explain our data analysis methods. In Section 4 we present the experimental and theoretical results on the generation of multipartite entanglement using a Josephson parametric system in both tripartite and quadripartite cases. In Appendices A.1–A.5 we present details of our analysis, experimental techniques, and parameter extraction. Appendix A.6 deals with analytical solution of the equations of motion for tri‐/quadripartite states using simplified linear model of parametric amplifier and establish correlations (graph connections) between spectral modes, which allows us to classify and control the entanglement structure. Furthermore, analytical forms of relevant covariance matrices are given.Theoretical FoundationsQuantum Dynamics of the DeviceOur work employs a parametric system comprised of a superconducting λ/4$\lambda /4$ resonator terminated in a SQUID loop. Such a setup forms the archetype of a narrow‐band superconducting Josephson parametric amplifier (JPA).[40,41] In our setting, we pump the JPA using multi‐tone external RF magnetic flux through the SQUID at frequencies that are approximately twice the frequency of the resonator ωd∼2ωr$\omega _d \sim 2 \omega _\text{r}$ (three‐wave mixing).[18,42] In the rotating frame, the Hamiltonian of the system, as derived in refs. [40, 43], is given by:1Hsys,rwa(t)=ℏΔra∼†a∼+ℏ2∑d=1p(αd∗eiΔdta∼2+αde−iΔdta∼†2)+6ℏKa∼†a∼†a∼a∼$$\begin{eqnarray} {H}_{\text{sys,rwa}}(t) &=& {\hbar }\Delta _r \tilde{a}^\dagger \tilde{a} +\nonumber \\ &&\frac{\hbar }{2}\sum _{d=1}^{p}(\alpha ^{*}_d e^{i\Delta _d t}\tilde{a}^{2} + \alpha _d e^{-i\Delta _d t} \tilde{a}^{\dagger 2}) + 6{\hbar }{K} \tilde{a}^{\dagger }\tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$where a∼$\tilde{a}$ (a∼†$\tilde{a}^{\dagger }$) is the annihilation (creation) operator for cavity photons in the rotating frame at angular frequency ωΣ/2$\omega _{\Sigma }/2$, αd$\alpha _d$ is the complex amplitude of the d‐th pump tone and Δd=ωd−ωΣ$\Delta _d=\omega _d-\omega _{\Sigma }$ is the angular frequency detuning of the corresponding tone. Possible extra phase factors in different pump tones are included in the complex pump amplitude αd=|α|eiφd$\alpha _d=|\alpha | e^{i\varphi _d}$. Here Δr$\Delta _\text{r}$ denotes the detuning between half of the average pump angular frequency ωΣ$\omega _\Sigma$ and the resonator angular frequency ωr$\omega _\text{r}$: Δr=ωr−ωΣ/2$\Delta _\text{r}=\omega _r-\omega _{\Sigma }/2$, with ωΣ$\omega _{\Sigma }$ representing the average angular frequency in multi‐tone driven case: ωΣ=(1/p)∑d=1pωd$\omega _{\Sigma }=(1/p)\sum _{d=1}^{p}\omega _d$ with d={1,⋯,p}$d=\lbrace 1,\dots ,p \rbrace$ as the pump tone index.Strongly driven SQUIDs are notoriously nonlinear. Therefore, we also include the nonlinear Kerr term with strength K to the description of our parametric system. The Kerr constant controls the parametric behavior close to and above the critical pumping threshold α≥αcrit$\alpha \ge \alpha _{\text{crit}}$. Several effects are accounted for by the Kerr nonlinearity, such as limited maximum gain, compression, observed at α≲αcrit$\alpha \lesssim \alpha _{\text{crit}}$, broadening and shifting of the resonance curve, and parametric oscillation above the critical point. In our experiment, we employ the pump‐power‐dependent gain coefficient to extract the Kerr constant (see Appendix A.4).In order to describe the coupling of the cavity resonator to the incoming transmission line and to an intrinsic thermal bath, we include two additional terms in the full Hamiltonian:2H(t)=Hsys,rwa(t)+Hsig+Hloss$$\begin{equation} {H}(t)=H_{\text{sys,rwa}}(t)+H_{\text{sig}}+H_{\text{loss}} \end{equation}$$where Hsig$H_{\text{sig}}$ includes the coupling to the signal port transmission line with dissipation rate κ, while Hloss$H_{\text{loss}}$ includes the coupling to the internal loss port with linear dissipation rate γ.Using the quantum Langevin equation (QLE), we obtain the output modes in our parametric system. We employ the standard input/output formalism in the rotating frame, which yields3a∼̇(t)=(−iΔr−κ+γ2)a∼−i∑d=1pαdeiΔdta∼†+κb∼in+γc∼in−12iKa∼†a∼a∼$$\begin{eqnarray} \dot{\tilde{a}}(t)=(-i\Delta _r &-& \frac{\kappa +\gamma }{2})\tilde{a} - i \sum _{d=1}^{p}\alpha _d e^{i\Delta _d t} \tilde{a}^\dagger \nonumber \\ && +\ \sqrt {\kappa }\tilde{b}_{\text{in}}+\sqrt {\gamma }\tilde{c}_{\text{in}} - 12i{K} \tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$where b∼in$\tilde{b}_{\text{in}}$ and c∼in$\tilde{c}_{\text{in}}$ are the ladder (annihilation) operators for the signal and linear dissipation ports, respectively.The output mode b∼out$\tilde{b}_{\text{out}}$ is obtained using the following relation between the incoming and outgoing modes:4b∼out(t)=b∼in(t)−κa∼(t)$$\begin{equation} \tilde{b}_{out}(t)=\tilde{b}_{\text{in}}(t)-\sqrt {\kappa }\tilde{a}(t) \end{equation}$$We are interested in the correlations embedded in the output mode given by Equation (4) in the time domain. The correlations can be revealed in full after Fourier transformation to the frequency domain. By defining finite‐band spectral modes (see Section 2.2) in the frequency domain and examining correlations between these spectral ranges, we can verify the presence of entanglement in the band‐limited microwave signals.Spectral Modes DefinitionParametric downconversion processes and the definition of employed spectral modes within the fundamental cavity resonance in a multi‐pump JPA are illustrated in Figure 1. The spacing and width of the spectral modes are selected in such a manner that that the modes are generated within the linewidth of the JPA resonance. Each pump that acts on the JPA triggers spontaneous parametric downconversion of a pump photon (vertical red arrows) into two photons, with their energies summing up to the energy of the pump photon (blue arrows). This process is stimulated by vacuum fluctuations, whose existence is a fundamental feature of quantum electrodynamics. One might expect that the downconversion processes would be random, occurring independently for each pump in the multi‐tone pumping situation. In such a case, the result would simply be a sum of the downconversion processes, but this turns out not to be the case. Instead, the photons are fundamentally correlated, even if they originate from different pump tones, because they were “born into existence” by the same quantum fluctuation. In other words, one spectral mode contains photons correlated with the quanta in the other spectral modes and, consequently, we expect multipartite correlations to appear (depicted schematically via the zigzag lines).1FigureDefinition of spectral modes and their correlations in a multi‐tone pump setting. Spectral modes in a multi‐pump JPA, where the pump tones (red arrows) trigger parametric downconversion process (PDC) leading to the appearance of multipartite correlation between microwaves (blue arrows), extracted from vacuum fluctuations. Numbered spectral modes depicted in green are also correlated due to the continuous pumping of the JPA resulting in multipartite entanglement between microwaves in the spectral modes. The bandwidth of each spectral mode Δ in Fourier analysis is chosen to be much narrower than the cavity resonance width.We consider the fundamental resonance of a transmission line JPA centered at ωr$\omega _r$ frequency with a bandwidth 2δω$2\delta \omega$, within which ([ωr−δω,ωr+δω]$[\omega _\text{r}-\delta \omega , \omega _\text{r}+\delta \omega ]$) we define N spectral modes as depicted in Figure 1. Let us define a∼$\tilde{a}$ as a vector of spectral modes:5a∼={a∼1,⋯,a∼N,a∼1†,⋯,a∼N†}T$$\begin{equation} \tilde{a}=\lbrace \tilde{a}_1,\dots ,\tilde{a}_N,\tilde{a}^\dagger _1,\dots ,\tilde{a}^\dagger _N\rbrace ^T \end{equation}$$where N is a total mode number and the creation a∼i†=a∼i†(t)$\tilde{a}_i^{\dagger }=\tilde{a}_i^{\dagger }(t)$ and annihilation a∼i=a∼i(t)$\tilde{a}_i=\tilde{a}_i(t)$ operators are time‐dependent. In general, the frequency difference between half of pth and (p+1)$(p+1)$th pump tones defines the bandwidth of the spectral mode a∼i$\tilde{a}_i$. We employ an equidistant pump scheme where bandwidth of each mode a∼i$\tilde{a}_i$ is defined as Δ such that the spectrum of a full set of modes a∼i$\tilde{a}_i$ covers the bandwidth 2δω$2\delta \omega$ of the cavity mode a∼$\tilde{a}$. In the experiment, we collect the emitted power over the whole [ωr−δω,ωr+δω]$[\omega _\text{r}-\delta \omega , \omega _\text{r}+\delta \omega ]$ frequency range and separate signals in the N modes using numerical postprocessing. The same operation could be implemented by using accurate bandpass filters with bandwidth Δ.Within the scope of this work, we consider only tripartite (N=3$N=3$, p=2$p=2$) and quadripartite (N=4$N=4$, p=3$p=3$) quantum states. In the following, we elucidate the internal structure of the generated states through a graph representation based on the quantum Langevin equation.Graphical Description of Quantum StatesIn order to construct a comprehensive graphical representation, we consider interactions between cavity and input vacuum modes by solving the QLE given in Equation (3) for N cavity modes defined in Equation (5); for details see Appendix A.6. Assuming the strong coupling regime with κ≫Δ$\kappa \gg \Delta$ while neglecting dissipation losses γ and the Kerr‐nonlinearity K, Fourier transformation yields a system of linear equations which can be cast in a matrix form:6Ma∼(ω)=κb∼in(ω)$$\begin{equation} \mathbf {M}\tilde{a}(\omega )=\sqrt {\kappa }\tilde{b}_{\text{in}}(\omega ) \end{equation}$$Here the interaction matrix M$\mathbf {M}$ contains diagonal entries provided by cavity‐related part of the Hamiltonian, whereas the parametric terms appear in the off‐diagonal entries. For example, one obtains for the tripartite case with Δr=0$\Delta _r=0$ and having different phases for the pump tones α1=αeiφ1,α2=αeiφ2$\alpha _1=\alpha e^{i\varphi _1}, \alpha _2 =\alpha e^{i\varphi _2}$:7M=c000iαeiφ100c0iαeiφ10iαeiφ200c0iαeiφ200−iα†e−iφ10c00−iα†e−iφ10−iα†e−iφ20c00−iα†e−iφ2000c$$\begin{eqnarray} \mathbf {M}= \def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & i\alpha e^{i\varphi _1} & 0\\ 0 & c & 0 & i\alpha e^{i\varphi _1} & 0 & i\alpha e^{i\varphi _2} \\ 0 & 0 & c & 0 & i\alpha e^{i\varphi _2} & 0\\ 0 & -i\alpha ^\dagger e^{-i\varphi _1} & 0 & c & 0 & 0\\ -i\alpha ^\dagger e^{-i\varphi _1} & 0 & -i\alpha ^\dagger e^{-i\varphi _2} & 0 & c & 0\\ 0 & -i\alpha ^\dagger e^{-i\varphi _2} & 0 & 0 & 0 & c \end{bmatrix} \nonumber\hspace*{-10pt}\\ \end{eqnarray}$$where the frequency dependency enters through the c=−iω+κ/2$c=-i\omega +\kappa /2$ coefficient. To express the intracavity modes through the vacuum input, we use the inverse matrix M−1$\mathbf {M}^{-1}$:8a∼(ω)=κM−1b∼in(ω)$$\begin{equation} \tilde{a}(\omega )=\sqrt {\kappa }\mathbf {M}^{-1}\tilde{b}_{\text{in}}(\omega ) \end{equation}$$For the tripartite case, such a matrix is written in the following way9M−1=1c2−2α2·c−α2c0α2eiΔφc0−iαeiφ100c0−iαeiφ10−iαeiφ2α2eiΔφc0c−α2c0−iαeiφ200iαe−iφ10c−α2c0α2eiΔφciαe−iφ10iαe−iφ20c00iαe−iφ20α2eiΔφc0c−α2c$$\begin{eqnarray} \mathbf {M}^{-1} &=& \frac{1}{c^2-2\alpha ^2} \nonumber \cdot \\ &&\def\eqcellsep{&}\begin{bmatrix} c-\frac{\alpha ^2}{c} & 0 & \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}} & 0 & -i\alpha e^{i\varphi _1}& 0\\ 0 & c & 0 & -i\alpha e^{i\varphi _1} & 0 & -i\alpha e^{i\varphi _2}\\ \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}} & 0 & c-\frac{\alpha ^2}{c} & 0 & -i\alpha e^{i\varphi _2} & 0\\ 0 & i\alpha e^{-i\varphi _1} & 0 & c-\frac{\alpha ^2}{c} & 0 & \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}}\\ i\alpha e^{-i\varphi _1} & 0 & i\alpha e^{-i\varphi _2} & 0 & c & 0\\ 0 & i\alpha e^{-i\varphi _2} & 0 & \pmb {\frac{\alpha ^2 e^{i\Delta \varphi }}{c}} & 0 & c-\frac{\alpha ^2}{c} \end{bmatrix}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$where Δφ=φ1−φ2$\Delta \varphi =\varphi _1-\varphi _2$. Besides two mode squeezing (TMS) correlations proportional to α, the matrix contains also beamsplitter correlations (BS) ∝α2$\propto \alpha ^2$ (denoted in bold). Note that the α2‐terms are absent in the matrix M$\mathbf {M}$ introduced in Equation (6). In this tripartite example, the phase difference Δφ$\Delta \varphi$ contributes only to the BS connection a∼i(ω)↔a∼j(ω)$\tilde{a}_i(\omega ) \leftrightarrow \tilde{a}_j(\omega )$ or a∼i†(ω)↔a∼j†(ω)$\tilde{a}^\dagger _i(\omega ) \leftrightarrow \tilde{a}^\dagger _j(\omega )$. TMS connections are defined by entries a∼i(ω)↔a∼j†(ω)$\tilde{a}_i(\omega ) \leftrightarrow \tilde{a}^\dagger _j(\omega )$ in M−1$\mathbf {M}^{-1}$. In Appendix A.6 we discuss how the phase shifts between pumps influence the structure of subspaces within the covariance matrix. In the quadripartite case, it turns out that those products of M−1$\mathbf {M}^{-1}$ responsible for beamsplitter interaction can be suppressed fully by choosing pump phases properly in certain pump tone configurations.Interestingly, the matrix κM−1$\sqrt {\kappa }\mathbf {M}^{-1}$ can also be interpreted as an adjacency matrix, which is used in graph theory for characterization of the connections, the graph edges. In regular H(amiltonian)‐graph,[7,29,38,39] each vertex represents a vacuum mode and the adjacency matrix describes correlations produced by two mode squeezing (TMS) between the vacuum modes. However, our scheme produces additional correlations, BS correlations, that can no longer be described purely by the TMS correlations and, therefore, the standard H‐graph theory needs a more generalized approach.In our approach, we introduce generalized H∼$\mathcal {\tilde{H}}$‐graphs formed by both two mode squeezing and beamsplitter correlations (see Appendix A.6). Examples of such graphs are presented in Figure 2a,b for the tripartite and quadripartite case, respectively. Intriguingly, the famous Greenberger–Horne–Zeilinger (GHZ) states[7,38,44] have different structure since they are devoid of BS correlations. However, by applying additional pump tones and adjusting phase difference between them, one can generate tripartite and quadripartite GHZ states consisting of only SQ correlations as shown in Figure 2c,d. In general, the considered multipumping scheme allows us to control SQ and BS bonds providing access to more complex structures of CV quantum states beyond GHZ‐like states.2FigureGraph representation for the entangled tripartite and quadripartite systems. Generalized H‐graph (H∼$\mathcal {\tilde{H}}$‐graph) representation of a) bisqueezed tripartite CV entangled state and b) quadripartite CV entangled H∼$\mathcal {\tilde{H}}$‐graph state obtained in our experiments. Vacuum modes (red circles) are connected via two‐mode squeezing (TMS, solid lines on graph) and beamsplitter (BS, dashed lines) correlation. Graphs (c,d) represent tri‐ and quadripartite GHZ states, which can be obtained via introducing an additional pump tone with Δd=0$\Delta _d=0$ in the tripartite case and two additional tones at −Δ$-\Delta$ and Δ in the quadripartite setting. The additional pumps supply missing TMS connections to the entangled states. For details, see Appendix A.6.From the experimental point of view, the observer is interested in the adjacency matrix for out‐coming modes b∼out$\tilde{b}_\text{out}$, which are obtained from a∼$\tilde{a}$ using the input–output relationship in Equation (4). This equation yields10b∼out(ω)=(I−κM−1)b∼in(ω)$$\begin{equation} \tilde{b}_{\text{out}}(\omega )=(\mathbf {I}-\kappa \mathbf {M}^{-1})\tilde{b}_{\text{in}}(\omega ) \end{equation}$$on the basis of which we may define the adjacency matrix M∼=I−κM−1$\mathbf {\tilde{M}}= \mathbf {I}-\kappa \mathbf {M}^{-1}$ for input–output mode graphs. Due to the linear nature of the equation, the unit matrix does not change the intrinsic form of interactions between the vacuum modes, but the correlation structure of intracavity and output spectral modes is equivalent. Consequently, analyzing a graph defined by the matrix M−1$\mathbf {M}^{-1}$ is sufficient to characterize the BS and TMS connections between vacuum spectral modes.Connection to Hamiltonian GraphCV cluster states with square‐lattice graph structure provides a foundation for measurement‐based continuous variable quantum computation (CVQC).[38] Cluster states can be asymptotically reached from H‐graph states in the case of infinite squeezing.[39] The H‐graph structure is defined by its adjacency matrix G$\mathbf {G}$, whose entries Gij$G_{ij}$ specify the multimode squeezing Hamiltonian as:11HS=ℏα∑i,jGij(a∼i†a∼j†+a∼ia∼j)$$\begin{equation} H_\text{S}= \hbar \alpha \sum _{i,j} G_{ij}(\tilde{a}_i^\dagger \tilde{a}_j^\dagger + \tilde{a}_i \tilde{a}_j) \end{equation}$$Here, the pump tone amplitudes are considered to have equal strength α. The matrix G$\mathbf {G}$ involves TMS correlations between modes a∼i↔a∼j†$\tilde{a}_i\leftrightarrow \tilde{a}^\dagger _j$, but as was pointed out before, the BS correlations do not show up in the H‐graph representation. The equations of motion for the operators are given by ia∼̇k†=α∑jGjka∼j$i\dot{\tilde{a}}_k^\dag =\alpha \sum \nolimits _j {{G_{jk}}} {\tilde{a}_j}$ and ia∼̇k=−α∑jGjk∗a∼j†$i{\dot{\tilde{a}}_k} = -\alpha \sum \nolimits _j {G_{jk}^*} \tilde{a}_j^\dag$. Taking the Fourier transform, the left hand side equals ω×a∼(ω)$ \omega \times \tilde{a}(\omega )$, and the combination ωa∼k(ω)=α∑jGjka∼j†(−ω)$\omega \tilde{a}_k(\omega )=\alpha \sum \nolimits _j {{G_{jk}}} {\tilde{a}_j}^\dag (-\omega )$ provides the connection to the QLE treatment in Equation (6): a∼(ω)$\tilde{a}(\omega )$ here is the cavity signal defined by the graph connections given by Gjk$G_{jk}$. Consequently, the basic graph structures are the same, but the form of M−1$\mathbf {M}^{-1}$ in the QLE analysis yields higher order correlations which are experimentally relevant.A standard description for graphs is based on the complex symmetric matrix Z=ie−G$\mathbf {Z}=ie^{-\mathbf {G}}$, which is interpreted as the adjacency matrix for an undirected Gaussian graph with complex‐valued edge weights.[38,39] Decomposing such a matrix up to quadratic terms Z=iI−iG+i2G2+O(h3)$\mathbf {Z}=i\mathbf {I}-i\mathbf {G}+\frac{i}{2}\mathbf {G}^2+\mathcal {O}(h^3)$, we obtain corrections to the adjacency matrix, which correspond to additional correlations, the BS correlations, obtained in our QLE analysis. Indeed, BS transformations embody interactions to second‐order, which provides classical correlations between corresponding nodes.[45]Let us now show the origin of BS correlations using the multimode squeezing Hamiltonian of Equation (11) in R=e−iℏHτ$R = e^{-\frac{i}{\hbar } H \tau }$ where we have considered that the system is pumped for a finite time τ. The multimode squeezing operator R can be decomposed to a combination of TMS operators, containing Bij=a∼i†a∼j†+a∼ia∼j$B_{ij}=\tilde{a}_i^\dagger \tilde{a}_j^\dagger +\tilde{a}_i\tilde{a}_j$, and BS transformations based on Tij=a∼ia∼j†−a∼i†a∼j$T_{ij}=\tilde{a}_i \tilde{a}_j^\dagger -\tilde{a}_i^\dagger \tilde{a}_j$. By utilizing the Zassenhaus expansion (up to first order of commutation relationship),[45] we obtain for the tripartite squeezing operator12R=e−iατ∑i,j=13Gij(a∼i†a∼j†+a∼ia∼j)=e−iατ(B12+B23)=e−iατB12e−iατB23eα2τ22[B12,B23]=e−iατB12e−iατB23eθ13T13$$\begin{eqnarray} R &=& e^{-i\alpha \tau \sum ^3_{i,j=1} G_{ij}(\tilde{a}_i^\dagger \tilde{a}_j^\dagger + \tilde{a}_i \tilde{a}_j)}=e^{-i\alpha \tau (B_{12}+ B_{23})}=e^{-i\alpha \tau B_{12}}\nonumber \\ && e^{-i\alpha \tau B_{23}}e^{\frac{\alpha ^2 \tau ^2 }{2}[B_{12},B_{23}]}= e^{-i\alpha \tau B_{12}}e^{-i\alpha \tau B_{23}}e^{\theta _{13}T_{13}} \end{eqnarray}$$Here, θ13 specifies the relative phase shift between the two pump tones. For detailed information on the expansion coefficients for a bisqueezed state we refer the reader to ref. [45].The total multimode squeezing operation can be considered as a combination of TMS operators, acting on the respective bipartitions, and BS transformations between the other modes. The beamsplitter correlations are phase dependent, and the strength of the BS contribution can be tuned down to zero in certain cases via proper choice of the phase difference between the pumps.The general decomposition of the multimode squeezing operator can be expressed as13R=∏1NTMSe−iατB12e−iατB23⋯∏1NBSeθ13T13eθ24T24⋯$$\begin{equation} R= \prod _{1}^{N_\text{TMS}}e^{-i\alpha \tau B_{12}}e^{-i\alpha \tau B_{23}}\dots \prod _{1}^{N_\text{BS}}e^{\theta _{13} T_{13}}e^{\theta _{24}T_{24}}\dots \end{equation}$$where the total number of TMS and BS operators in the decomposition is NBS=∑nN−2n$ N_{\text{BS}}=\sum _{n} N-2n$ and NTMS=∑nN−2n+1$N_{\text{TMS}}=\sum _{n} N-2n+1$; where n=1,2,⋯,⌊N2⌋$n=1,2,\dots ,\lfloor \frac{N}{2}\rfloor$ for the configuration introduced in Figure 1 with N−1$N-1$ pump tones and N spectral modes. Collecting all of the entries of Bij$B_{ij}$ and Tij$T_{ij}$, we obtain a generalized adjacency matrix G∼$\mathbf {\tilde{G}}$ with entries G∼ij=∑1NTMSBij+∑1NBSTij$\tilde{G}_{ij}=\sum _1^{N_{\text{TMS}}}B_{ij}+\sum _1^{N_{\text{BS}}}T_{ij}$. Thus, the beamsplitter correlations in adjacency matrix arise naturally from the squeezing operator formalism when the Hamiltonian is supplied with the second order terms in pump amplitude. The structure of the matrix G∼$\mathbf {\tilde{G}}$ defines the edge connections in the generalized H∼$\mathcal {\tilde{H}}$‐graph.Verification of the Multipartite EntanglementThe generalized graph analysis allows us to visualize the structure of entanglement in the quantum state generated by simultaneous multiple pump tones. However, in order to estimate the amount of quantum resources embedded in the state, we have to investigate and quantify the classical and quantum correlations and determine how they reflect the genuine multipartite entanglement of the state.Within the framework of parametric amplifiers, all microwave fields produced by a JPA below the critical threshold are Gaussian.[17,46] Therefore, the output states of a N‐mode JPA can be fully characterized by its covariance matrix of 2N‐length column vector with quadratures r∼=(x∼1,p∼1,⋯x∼N,p∼N)T$\tilde{r}\nobreakspace =\nobreakspace (\tilde{x}_1, \tilde{p}_1, \dots \tilde{x}_N, \tilde{p}_N )^{T}$, where x∼i=(a∼i+a∼i†)/2$ \tilde{x}_i = (\tilde{a}_i+\tilde{a}^{\dagger }_i)/2$ and p∼i=(a∼i−a∼i†)/2i$\tilde{p}_i=(\tilde{a}_i-\tilde{a}^{\dagger }_i)/2i$. The covariance matrix V$\mathcal {V}$, whose elements are given by14Vij=2Δr∼iΔr∼j+Δr∼jΔr∼i−4Δr∼iΔr∼j$$\begin{equation} \mathcal {V}_{ij}=2\left&lt;\Delta \tilde{r}_i \Delta \tilde{r}_j + \Delta \tilde{r}_j \Delta \tilde{r}_i \right&gt; - 4\left&lt;\Delta \tilde{r}_i \right&gt; \left&lt;\Delta \tilde{r}_j\right&gt; \end{equation}$$is sufficient for detection of the entanglement, eliminating the need for analysis of the full density matrix. The last term can be ignored as we take Δr∼i=r∼i−<r∼i>$\Delta \tilde{r}_i=\tilde{r}_i-&lt;\tilde{r}_i&gt;$. Obtaining the covariance matrix, we can analyze the entanglement and examine the structure of the quantum state. In this work, we consider fully inseparable states and genuinely entangled states, for which the covariance‐based detection is, in general, more robust than detection via complete determination of the state.[47] While the covariance matrix is sufficient for evaluating entanglement of Gaussian states, it is necessary to include higher‐order correlations in the evaluation of non‐Gaussian states.[48,49]To examine inseparability properties of the quantum state,[47,50–55] we apply symplectic transformations to the covariance matrix and calculate its symplectic eigenvalues—the positive partial transpose (PPT) criteria.[56,57] Such transformations are equivalent to a phase space reflection of a single party in the N‐partite state.[50] All minimum symplectic eigenvalues {νi}i=1N$\lbrace \nu _i\rbrace _{i=1}^N$ would be less than one, which indicates that this partially time‐reversed state is unphysical; in other words, the original state is fully inseparable. As has been pointed out in ref. [54], if the purity of states cannot be guaranteed in an experimental setting, verification of full inseparability in a multimode system does not imply genuine multipartite entanglement (GME).The entanglement structure becomes more involved with increasing number of parties. While the symplectic transform approach indicates that any one partite were inseparable from the whole, a state that is a mixture of separable states would show full inseparability based on this PPT criterion. The states that cannot be written in such a way are called genuinely entangled[17] and the verification of such states differs from full inseparability. Using generalized position and momentum observables, an entanglement criterion has been derived and applied to confirm tripartite energy‐time entanglement of three spatially separated photons.[58] In particular, there is an universal GME criterion derived in ref. [47] and further refined in ref. [54]. This GME criterion utilizes only variances of quadrature operators and it can be used for entanglement verification without any additional measurements. This general criterion was recently employed for verification of genuine tripartite entanglement of microwaves in a double superconducting cavity setting.[34]The GME criterion is based on the weighted variance of the quadratures, u=∑ihix∼i$u = \sum _i h_i \tilde{x}_i$ and v=∑kgkp∼k$v = \sum _k g_k \tilde{p}_k$; i,k={1,2,3,⋯,N}$i,k=\lbrace 1,2,3,\dots , N \rbrace$. Violation of the inequality15S≡Δu2+Δv2f3(hi,gi)≥1$$\begin{eqnarray} S \equiv \frac{\left&lt; \Delta u^2 \right&gt; + \left&lt; \Delta v^2 \right&gt;}{f_3(h_i,g_i)} \ge 1 \end{eqnarray}$$where16f3(hi,gi)=12min{|h1g1+h2g2|+|h3g3|,|h3g3+h2g2|+|h1g1|,|h1g1+h3g3|+|h2g2|,}$$\begin{eqnarray} f_3(h_i,g_i) = \frac{1}{2} \min \lbrace & |h_1 g_1 + h_2 g_2| + |h_3 g_3|, \nonumber \\ & |h_3 g_3 + h_2 g_2| + |h_1 g_1|, \nonumber \\ & |h_1 g_1 + h_3 g_3| + |h_2 g_2|, \rbrace \end{eqnarray}$$is sufficient to confirm genuine tripartite entanglement (N=3$N=3$) and violation of17S≡Δu2+Δv2f4(hi,gi)≥1$$\begin{eqnarray} S \equiv \frac{\left&lt; \Delta u^2 \right&gt; + \left&lt; \Delta v^2 \right&gt;}{f_4(h_i,g_i)} \ge 1 \end{eqnarray}$$where18f4(hi,gi)=min{|h1g1+h2g2+h3g3|+|h4g4|,|h4g4+h2g2+h3g3|+|h1g1|,|h4g4+h1g1+h3g3|+|h2g2|,|h4g4+h1g1+h2g2|+|h3g3|,|h1g1+h2g2|+|h3g3+h4g4|,|h1g1+h3g3|+|h2g2+h4g4|,|h2g2+h3g3|+|h1g1+h4g4|}$$\begin{eqnarray} f_4(h_i,g_i) = \min \lbrace & |h_1 g_1 + h_2 g_2 + h_3 g_3| + |h_4 g_4|, \nonumber \\ & |h_4 g_4 + h_2 g_2 + h_3 g_3| + |h_1 g_1|, \nonumber \\ & |h_4 g_4 + h_1 g_1 + h_3 g_3| + |h_2 g_2|, \nonumber \\ & |h_4 g_4 + h_1 g_1 + h_2 g_2| + |h_3 g_3|, \nonumber \\ & |h_1 g_1 + h_2 g_2| + |h_3 g_3 + h_4 g_4|, \nonumber \\ & |h_1 g_1 + h_3 g_3| + |h_2 g_2 + h_4 g_4|, \nonumber \\ & |h_2 g_2 + h_3 g_3| + |h_1 g_1 + h_4 g_4|\rbrace \end{eqnarray}$$is sufficient to confirm genuine quadripartite entanglement (N=4$N=4$) with weights hi,gk$h_i,g_k$ being in range [ − 1, 1]. To simplify the search domain, we set h1=g1=1$h_1=g_1=1$ and hi=h,gi=g$h_i = h, g_i = g$, for i={2,3}$i=\lbrace 2,3 \rbrace$ or i={2,3,4}$i=\lbrace 2,3,4 \rbrace$ with respect to the number of parties N.[54]In the double and triple pumping schemes, we generate generalized tripartite and quadripartite H$\mathcal {H}$‐graph states, which have a different structure compared with Greenberger‐Horne‐Zeilinger (GHZ) states.[7,38,44] The generalization deals with addition of BS correlations in the H$\mathcal {H}$‐graph structure. However, by applying additional pump tones and adjusting the phase difference between them, we can obtain regular GHZ type of entangled states. Thus, our scheme facilitates control of TMS and BS correlations and, thereby, allows tuning of the structure of the entangled state.Typically, the experimental weights of the graph edge connections are slightly non‐symmetric due to imperfections in the measurement settings. This results in a difference in the optimal weight values in the GME criterion. In order to find the full violation of the criterion in our analysis, we swap over all possible “base” modes (with weights hi=gi=1$h_i=g_i=1$) in order to detect the minimum value for S.ExperimentExperimental MethodsIn our microwave experiments, we employ a niobium λ/4$\lambda /4$ coplanar 50Ω$50\,\Omega$ transmission line terminated into a SQUID loop (QWJPA), forming a quarter‐wave Josephson parametric amplifier. The SQUID's junctions are formed by Nb/Al‐Al2O3/Nb 1× 1 μ$\umu$m tunnel barriers with Ic≃4$I_\text{c}\simeq 4$ μ$\umu$A critical current. The JPA, operating in three‐wave‐mixing mode around the cavity frequency ωr$\omega _\text{r}$, is pumped by an external RF magnetic flux through the SQUID loop using a single turn pump coil at frequency 2ωr$2\omega _\text{r}$ (marked as 2ωLO in Figure 3a).[40,59] We chose this operation regime, because for four‐wave mixing (typically using a current pump near ωr$\omega _\text{r}$), the large amplitude pump is within the amplification bandwidth, whereas the three‐wave mixing process separates the pump tone from the amplified signals, thus simplifying the practical use of the JPA. The loaded quality factor at the operation frequency is ≈900, while the internal Q is by a factor of three larger. The basic (zero‐flux) resonator frequency is 6.115 GHz, it can be tuned below 5.5 GHz by imposing external DC magnetic flux through the SQUID.3FigureExperimental scheme and device characterization. a) Principle of the experimental setup for tripartite entanglement measurements. The device is connected to the test ports via circulators. The DC bias current and AC pumping of flux are combined and reseparated in bias‐T components. Depending on the measurement type, input is connected either to a vector network analyzer (VNA) or a 50 Ω termination, whereas the output is directed either to a VNA, a signal analyzer, or a analog/digital converter. The frequency span of spectral modes and their separation is given for tripartite and quadripartite case in frames (b,c), respectively.Our measurement setup is illustrated in Figure 3a. The experiments were conducted at 20 mK using a BlueFors LD400 dry dilution cryostat. The JPA was protected from external magnetic fields using a Cryoperm shield. The DC flux bias and the RF pump shared a common on‐chip flux line, and the signals were combined in an external bias‐tee. Since our basic microwave setting is for reflection measurements, the sample is connected to the input and output ports via a circulator having a frequency band of 4 − 12 GHz. A vector network analyzer (VNA) was used to characterize the sample, whereas during the entanglement generation measurements, the signal port was kept terminated. By applying a multitone pump to JPA, correlated microwaves are generated from vacuum fluctuations. In the tripartite setting illustrated in Figure 3a, we had control over the relative phases of the pumps directly whereas, in the quadripartite case, phase rotation was possible in the data analysis only. Basic experimental data in the tripartite case as well as determination of the cavity parameters κ, γ, and K are discussed in Appendix A.4.Tripartite case: In the tripartite case, phase‐controlled pump signals from the RF waveform generator are mixed with the frequency‐doubled local oscillator (LO) frequency, filtered by a pair of home‐made, tunable bandpass cavity filters. In order to avoid spurious pumping at ω2LO, a band rejection filter tuned to 2ωLO is employed. The filtering ensures passage only for the desired two pump signal at angular frequencies ω1=ωr−Δ/2$\omega _1=\omega _\text{r}-\Delta /2$ and ω2=ωr+Δ/2$\omega _2=\omega _\text{r}+\Delta /2$. Sufficient noise thermalization is ensured by the 46 dB attenuation because the pump coil is only weakly coupled to the SQUID.We apply a DC magnetic flux of ΦDC=0.383Φ0$\Phi _{\text{DC}}=0.383\Phi _0$ through the SQUID loop, resulting in ωr2π=6.024$\frac{\omega _\text{r}}{2\pi }=6.024$ GHz for the cavity frequency. In the three‐mode experiment, we apply two pump tones at (2×ωr2π−2)$(2\times \frac{\omega _\text{r}}{2\pi }-2)$ MHz and (2×ωr2π+2)$(2\times \frac{\omega _\text{r}}{2\pi }+2)$ MHz, with the half‐frequencies positioned as depicted in Figure 3b and the correlated spectral modes are defined symmetrically with respect to the pump half‐frequencies. Each mode has a bandwidth of 1.9 MHz and is separated from the other modes by 0.1 MHz. The phase control in the measurement is facilitated by phase‐locking of microwave generators to a 10 MHz rubidium reference clock and by using a joint external trigger.To collect data for correlation analysis, we mix down the output signal using a synchronized LO signal and record the output quadratures using two channels of a Teledyne SP Devices ADQ14 digitizer with sample rate of 50 MSa s−1 per channel covering the bandwidth of 25 MHz. Furthermore, we employ an overall detuning of 14 MHz, that is, a heterodyne detection scheme, in order to avoid 1/f$1/f$ noise from the measurement devices and the IQ mixer in the frequency conversion part of the setup. Using digital postprocessing, we can easily shift the center frequency of the heterodyned MHz signal to zero, ready for final correlation analysis of the modes.Our three‐mode measurement scheme provides the remarkable advantage of physical control of the phase difference between the two pump tones, which is essential for the analysis of phase dependence in the entangled states. The phase difference between 2 LO signal used as the carrier of the pump tones and LO readout signal remains fixed. Therefore, a change of the initial phase of the IF‐signal in one of the pump generator channels (Agilent 33600B in Figure 3a) relative to the second IF channel creates an effective phase difference Δφ$\Delta \varphi$ between two pump tones. Importantly, this difference Δφ$\Delta \varphi$ is preserved after mixing with 2LO (using simplified notation): e−i(ω1t+Δφ)e−iω2LOt=e−i((ω1+ω2LO)t+Δφ);e−i(ω2t)e−iω2LOt=e−i((ω2+ω2LO)t)$e^{-i(\omega _1 t + \Delta \varphi )}e^{-i\omega _{\text{2LO}} t}=e^{-i((\omega _1 +\omega _{\text{2LO}})t+\Delta \varphi )}; e^{-i(\omega _2 t)}e^{-i\omega _{\text{2LO}} t}=e^{-i((\omega _2 +\omega _{\text{2LO}})t)}$.Since our fully phase‐locked scheme preserves the phases of the received, demodulated output quadratures, the measurement of covariance matrix components can be averaged for reducing noise in the elements. Finally, the reference phase of a single mode (defining the basis for I and Q) can be adjusted in postprocessing step in such a way that the corresponding subspace of the covariance matrix becomes a diagonal 2 × 2 matrix.In the tripartite case, indeed, we find that the hardware‐controlled relative phase rotation (in addition to the reference phase value to both channels of the pump generator) is equivalent to a proper phase rotation in the postprocessing step. The postprocessing will be discussed in more detail in Section 4.Quadripartite case: In the four spectral mode case, we simplify the experimental setup by eliminating the physical phase control, and replaced it by postprocessing of the received signals. This simplification possibility highlights the scalability of our entanglement generation method. The employed digital postprocessing is equivalent to hardware‐level selective separation of spectral modes into four channels, for example, using bandpass filters in conjunction with power splitters, and additional tunable delay lines for each selected spectral mode frequency.We apply a DC magnetic flux of ΦDC=0.417Φ0$\Phi _{\text{DC}}=0.417\,\Phi _0$ through a SQUID resulting in ωr=5.978$\omega _\text{r}=5.978$ GHz cavity frequency that slightly differs from the tripartite case, see Figure 3b,c. In order to generate quadripartite correlations, we apply three pump tones using Anapico APMS 12G generator, using strong high‐pass filtering (2 of Mini‐Circuits VHF‐8400+) to avoid subharmonic transmission to the circuitry. In this scheme, we avoid any external mixers for the input pump microwaves. By applying three phase‐locked pump tones at frequencies 2×ωr2π$2\times \frac{\omega _\text{r}}{2\pi }$ MHz, (2×ωr2π+1)$(2\times \frac{\omega _\text{r}}{2\pi }+1)$ MHz, and (2×ωr2π−1)$(2\times \frac{\omega _\text{r}}{2\pi }-1)$ MHz, we generated four correlated spectral modes out from the ground state of the microwave cavity. Each mode has a bandwidth of 0.4 MHz and is separated from adjacent modes by 0.2 MHz. The output microwaves are captured, mixed down, and digitized by Anritsu MS2830A Signal Analyzer with a bandwidth of 2 MHz. Again, averaging is needed to lower the noise in the covariance matrix elements, and in this scheme, digital postprocessing is necessary to unify the phase settings in the covariance matrices before summation.The experimental detection of H∼$\mathcal {\tilde{H}}$‐graph states and their genuine multipartite entanglement depends on relations among the covariance matrix elements as discussed in Section 2.5. The degree of violation in the GME condition S<1$S&lt;1$ depends strongly on the magnitude ratio of the diagonal covariance elements to the off‐diagonal ones. Therefore, calibration of the detected signal powers is decisive, which is discussed in Appendix A.3.Scaled Covariance MatrixThe system gain <GΣ>$&lt;G_{\Sigma }&gt;$ determined in Appendix A.3 refers to measured power per unit bandwidth. Since the measured spectral mode quadratures Ii$I_i$ and Qi$Q_i$ are determined over the band Δfi$\Delta f_i$, the scaled quadrature xi$x_i$, equivalent to the amplitude of the quantum mechanical operators x$\mathbf {x}$, is given by the formula19xi=IiGiZ0hfiΔfi$$\begin{equation} {x}_i=\frac{I_i}{\sqrt {G_i Z_0 h f_i \Delta f_i}} \end{equation}$$where Gi=<GΣ,i>$G_i=&lt;G_{\Sigma ,i}&gt;$ is the system gain for ith spectral mode, Z0=50$Z_0 = 50$ Ω is transmission line impedance and Δfi$\Delta f_i$ is the bandwidth of the spectral mode: Δfi=2$\Delta f_i = 2$ MHz or Δfi=0.4$\Delta f_i = 0.4$ MHz for tripartite and quadripartite cases, respectively. Similar scaling is applied also to the quadrature component Qi$Q_i$.Similar to our earlier work,[43] the noise added by the preamplifier is subtracted from the diagonal elements of the covariance matrix (see Equation (14))20V=4Von−Voff+Icothhfi2kbTi$$\begin{equation} \mathcal {V}=4{\left(\mathbf {V}_{\text{on}}-\mathbf {V}_{\text{off}} \right)} + \mathbf {I} \coth \frac{hf_i}{2 k_\text{b} T_i} \end{equation}$$where Voff$\mathbf {V}_{\text{off}}$ denotes the covariance matrix measured in the absence of the pump. Due to scaling of the covariance matrix V$\mathcal {V}$, this equation yields a unity diagonal matrix in the absence of pumping at T→0$T \rightarrow 0$. The average physical temperature in our experiments is Ti=20$T_i=20$ mK resulting in cothhfi2kbTi=1.000$\coth \frac{hf_i}{2 k_\text{b} T_i}=1.000$.Multipartite EntanglementTo characterize the structure of the entanglement in output states, we analyze the resulting covariance matrices using PPT formalism and GME criteria discussed in Section 2.5 for tripartite and quadripartite cases.Tripartite caseLeveraging the amplitude and phase control of the pump signals, we experimentally evaluate the PPT and GME criteria values at different pump parameters. For comparison, we also conducted detailed numerical simulations based on the QLE in Equation (A8) using experimentally determined JPA parameters in the measurements. In general, we find good agreement between simulations and the experimental data, which is reassuring concerning the validity of the results.Figure 4 depicts our experimental results on genuine tripartite entanglement and their comparison with simulations. Figure 4a illustrates results of numerical simulations on GME in terms of S defined in Equation (15). At weak pumping, the condition for genuine entanglement S<1$S &lt; 1$ is fulfilled almost independent of the pump phases, but with increasing A, the simulations reveal an even smaller range of Δφ$\Delta \varphi$ yielding S<1$S &lt; 1$ (see the inset in Figure 4a). The strongest genuine tripartite entanglement is reached at Δφ≃−120∘$\Delta \varphi \simeq -120^{\circ }$ under normalized pumping amplitude A≃0.22$A \simeq 0.22$, at which the simulations reach S=0.70$S=0.70$. At the minimum of S, the corresponding weights are hi={1,−0.65,−0.65}$h_i=\lbrace 1,-0.65,-0.65\rbrace$ and gi={1,0.65,0.65}$g_i=\lbrace 1,0.65,0.65\rbrace$. It is noteworthy that the phase setting Δφ=+90∘$\Delta \varphi = + 90^{\circ }$ yields clearly worse entanglement than Δφ=−90∘$\Delta \varphi = -90^{\circ }$. This asymmetry in GME between Δφ=±90∘$\Delta \varphi = \pm 90^{\circ }$ arises from differences in the covariance matrices which is illustrated in Figure 5.4FigurePhase‐dependent genuine entanglement of tripartite bisqueezed state. a) Simulation results for GME criterion as a function of normalized pump amplitude for three different values of the phase difference Δφ$\Delta \varphi$ between pump signals indicated in the figure; the simulation parameters were set to match the experiment (see Appendix A.4). The inset illustrates S(A,Δφ)$S(A,\Delta \varphi )$ up to critical amplitude A=0.5$A=0.5$. The weights hi,gi$h_i, g_i$ were optimized in the calculation of S as discussed in the text. b) Experimental values for S as a function of A at the same phase difference values Δφ$\Delta \varphi$ between pump signals as in frame (a). The inset illustrates measured S(A,Δφ)$S(A,\Delta \varphi )$ up to A=0.5$A=0.5$. In general, the measured S(A,Δφ)$S(A,\Delta \varphi )$ corresponds quite well to the inset in frame (a). Due to noise, the measured GME nearly vanishes around Δφ≃+90∘$\Delta \varphi \simeq +90^\circ$ where even the simulated S is only slightly below 1. The best genuine multipartite entanglement is reached at Δφ=−90∘⋯−120∘$\Delta \varphi = -90^\circ \dots -120^\circ$ owing to phase shifts introduced by the cavity (see Figure A3 in Appendix A.5). The parametric drive changes the phase response of the cavity which leads to a shift in the optimum conditions for GME as a function of A.5FigureCovariance matrix of the genuinely entangled tripartite bisqueezed state. Experimentally obtained tripartite covariance matrices for genuinely entangled bisqueezed states. The phase difference between the two displayed cases is 180°. Via this control of pump phases we demonstrate the rotation of the desired subspace elements, corresponding to TMS type of correlations (1 − 2 or 2 − 3). Subspace elements related to BS correlations (1 − 3), shows always dependency on the distribution of the elements in the corresponding driven TMS subspaces. However, the eigenvalues of the BS subspace remain constant.Our experimental data on S in Figure 4b displays similar features as Figure 4a. The measured GME criterion S as function of normalized pump amplitude for three phase differences is shown in Figure 4b. In the experimental data, nearly no GME is observed at positive phase differences, whereas Δφ=−90∘$\Delta \varphi = - 90^{\circ }$ and Δφ=−120∘$\Delta \varphi = - 120^{\circ }$ yield suppression down to S=0.75±0.05$S=0.75 \pm 0.05$. The measured result at Δφ=−120∘$\Delta \varphi = - 120^{\circ }$ follows quite well the simulated behavior as a function of the pump amplitude, and genuine entanglement is observed in the normalized amplitude range A∈[∼0.01,0.4]$A \in [\sim 0.01, 0.4]$. Overall, the pattern of S(φ,A)$S(\varphi , A)$ in the inset of Figure 4b coincides with the simulated pattern in Figure 4a. The agreement strongly supports the presence of genuinely entangled bisqueezed state in our experiment. We emphasize that the optimum entanglement at Δφ=−120∘$\Delta \varphi = - 120^{\circ }$, observed both in our simulations and in the experiment, cannot be obtained from a simple analytical calculations for the lossless, strongly coupled model. The reason is the frequency‐dependent phase response of the cavity due to finite coupling and dissipation rates (see Section A.5), which, when included in the simulations, result in very good matching with the experiment.Covariance matrices measured at the pump phase difference Δφ=+90∘$\Delta \varphi = + 90^{\circ }$ and Δφ=−90∘$\Delta \varphi = - 90^{\circ }$ are illustrated in Figure 5. Technically, by rotating the phase of a pump signal, we selectively control certain subspace of the covariance matrix, which can be seen in Figure 5. An applied phase shift to the first or the second pump rotates directly the subspace corresponding to two mode squeezing correlations, modes 1 − 2 or 2 − 3, respectively. If no phase shift to the selected pump tone is applied, the corresponding TMS subspace preserves its distribution of covariances. The subspace spanned by modes 1 and 3, corresponding to the beamsplitter type of correlations, has a structure according to products of the involved TMS subspaces. Distinct control of the BS subspace alone (leaving the TMS subspaces fixed) using a rotation of the pump phases is not possible.Comparison of Figures 5a,b, reveals how the subspaces transform with the phase difference from Δφ=−90∘$\Delta \varphi = - 90^{\circ }$ to Δφ=+90∘$\Delta \varphi = + 90^{\circ }$. Subspace 1 − 3 corresponding to BS correlations shows a sign inversion in its elements. Subspace 2 − 3 corresponds to the pump, the phase of which has not been changed and, thus, its elements remain fixed. The phase of the first pump has been changed by π which inverts the TMS correlation in subspace 1 − 2. These subspaces which are controlled by the pump phase settings, can expressed in symmetric <I1I2>≃<I2I3>$&lt;I_1I_2&gt; \simeq &lt;I_2I_3&gt;$ or antisymmetric form <I1I2>≃−<I2I3>$&lt;I_1I_2&gt; \simeq -&lt;I_2I_3&gt;$, see Figures 5a,b, respectively.The inseparability of the covariance matrix was investigated using the PPT criteria (see Section 2). Symplectic eigenvalues obtained for mode partitions 1 − 23, 2 − 13, and 3 − 12 are depicted in Figure 6a for simulations, while experimental data are displayed in Figure 6b; here the first index specifies the mode in which the sign of the momentum has been reversed. The results are plotted as a function of normalized pump amplitude since the phase difference between the pump tones does not play a role. Indeed, while the genuine entanglement is sensitive to both pump amplitude and phase difference, the PPT criterion is phase independent—the minimum symplectic eigenvalues remain constant when the phase difference is varied at fixed pump amplitude. Therefore, by exercising phase control over each pump, we gain the ability to switch from fully inseparable state to genuinely entangled state, without making any changes to the type of interaction between the modes. According to experimental results in Figure 6b, the middle frequency acted on by both pumps is the most inseparable part of the covariance matrix.6FigurePhase‐independent full inseparability of tripartite bisqueezed state. a) PPT criteria in terms of the minimum symplectic eigenvalues simulated for our double‐pump QWJPA using experimentally determined parameters. Eigenvalues min{νi}$\min \lbrace \nu _i\rbrace$ are traces over normalized pump amplitude A; permutations 1 − 23, 2 − 13, and 3 − 12 have been considered. The symplectic eigenvalues are the smallest for time‐reversed second mode (ν2−13$\nu _{2-13}$), which participates to both TMS processes. b) Experimentally determined symplectic eigenvalues for the same permutations ( in the experiment. Gray dashed line displays the full inseparability threshold. The difference in the simulated behavior of ν1−23$\nu _{1-23}$ and ν3−12$\nu _{3-12}$ is caused by asymmetry due to finite value of resonance detuning Δr$\Delta _\text{r}$.Quadripartite caseFor the quadripartite case, we apply three pump drives with identical amplitude α1=α2=α3=α$\alpha _1=\alpha _2=\alpha _3=\alpha$. While our goal is to demonstrate genuine entanglement generation of cluster states (mode structure depicted on Figure 2), we reject direct, physical phase control, and use digital postprocessing to transform the covariance matrix to the desired form on which we then verify its entanglement properties. However, we do preserve the coherence between pump tones by phase locking so that the relative phases do not fluctuate over time. By applying a postprocessing phase rotation for each mode separately, we bring the covariance matrix into the symmetric form (see Section 4).For the analysis of full inseparability of the covariance matrix according to the PPT criterion, we evaluate the minimum symplectic eigenvalues min{νi}$\min \lbrace \nu _i\rbrace$ of the following mode permutations: 1 − 234, 2 − 134, 3 − 124, 4 − 123. The experimentally obtained symplectic eigenvalues as function of normalized pump amplitude A are displayed in Figure 7 alongside with the corresponding predictions given by our numerical simulations. The minimum symplectic eigenvalue min{νi}=0.79±0.018$\min \lbrace \nu _i\rbrace =0.79\pm 0.018$ is reached, while all of the eigenvalues in the normalized pump amplitude range 0.01≲A≲0.15$0.01\lesssim A \lesssim 0.15$ are less than 1. Compared with the minimum symplectic eigenvalues in Figure 6, we may conclude that the influence of BS correlations on min{νi}$\min \lbrace \nu _i\rbrace$ value is less in the quadripartite state than for the tripartite case.7FigureFull inseparability and genuine entanglement of quadripartite H∼$\mathcal {\tilde{H}}$‐graph state (generalized H‐graph). a) Results on PPT criterion for four permutations of a four‐mode Gaussian state. Minimum eigenvalues min{νi}$\min \lbrace \nu _i\rbrace$, indicated as open circles, are traced over normalized pump amplitude A. Results of our QLE simulations, plotted as solid curves, exhibit good correspondence with the experimentally obtained values. The full inseparability condition {ν1−234<1}⋀{ν2−134<1}⋀{ν3−124<1}⋀{ν4−123<1}$\lbrace \nu _{1-234}&lt;1\rbrace \bigwedge \lbrace \nu _{2-134}&lt;1\rbrace \bigwedge \lbrace \nu _{3-124}&lt;1\rbrace \bigwedge \lbrace \nu _{4-123}&lt;1\rbrace$ is fulfilled in the range of 0.01≲A≲0.15$0.01\lesssim A\lesssim 0.15$. The green dashed line displays the entanglement threshold.The GME criterion for four modes as a function of the normalized pump amplitude is depicted in Figure 7b; the symbols display data while the simulation result is indicated by the solid curve. As was discussed in Section 2, the optimized weights in GME inequality Equation (17) are chosen in the same manner as in the tripartite case: h1=g1=1$h_1=g_1=1$ and hi=h,gi=g$h_i = h, g_i = g$, i={2,3,4}$i=\lbrace {2,3,4\rbrace }$. The strongest genuine entanglement S=0.84±0.02$S=0.84 \pm 0.02$ is observed at A≂0.08$A\eqsim 0.08$ pump amplitude using the weights hi={1,−0.51,−0.51,−0.51}$h_i=\lbrace 1,-0.51,-0.51,-0.51\rbrace$ and gi={1,0.69,0.69,0.69}$g_i=\lbrace 1,0.69,0.69,0.69\rbrace$. The numerical simulation provides h and g coefficients that coincide with the experimental values with 1% error, which strongly establishes that the states produced in the experiment coincide with the ones that were obtained and analyzed in the numerical model.The covariance matrices obtained in the experiment and using numerical simulation are presented in Figure 8a,b, respectively. They are determined at the strongest entanglement point reached at A≂0.08$A\eqsim 0.08$. TMS type of correlations are seen in the mode combinations 1 − 2, 2 − 3, 3 − 4, and 1 − 4. Subspaces corresponding to BS correlations are visible in the plot as product distributions in 1 − 3 and 2 − 4 subpartitions. The covariance matrix structures illustrated in Figure 8 correspond directly to the H∼$\mathcal {\tilde{H}}$‐graph structures shown on Figure 2b. In general, we conclude that for the employed pump configuration, the genuine quadripartite entanglement appears in the amplitude range 0.01≲A<0.13$0.01\lesssim A&lt;0.13$.8FigureCovariance matrix of genuinely entangled quadripartite H∼$\mathcal {\tilde{H}}$‐graph state. a) Experimental covariance where the rotation of the TMS subspaces 1−2, 2−3, and 3−4 have been made in such a way that the structure coincides with the matrix in Equation (A36) of Appendix A.6 (each pump has phase π2$\frac{\pi }{2}$). The employed pump amplitude A≂0.08$A\eqsim 0.08$ yields the smallest value for S. b) Simulated covariance matrix using equal pump phases π2$\frac{\pi }{2}$ at A≂0.08$A\eqsim 0.08$. The difference from the matrix in Equation (A36) is due to the cavity response that induces extra phase shifts.DiscussionThe control of bisqueezed tripartite and generalized H‐graph (H∼$\mathcal {\tilde{H}}$‐graph) quadripartite states by relative positioning of the pump frequencies and their phases is indicative of the strong potential of these methods for CV quantum state processing. The basic parametric microwave setting allows for enhancement in the number of spectral modes by additional pump tones, which leads to generation of more complex, entangled H∼$\mathcal {\tilde{H}}$‐graph states. Enhanced number of modes requires larger bandwidth, which calls for broadband parametric devices such as TWPAs[60–64] or broadband JPAs[59,65] in order to avoid problems with spectral mode crowding.Our approach based on QLE puts in evidence additional correlations, which are captured by the definition of H∼$\mathcal {\tilde{H}}$‐graph states. The correlations arise naturally from the connection between intracavity modes and the input vacuum modes, due to which the same vacuum fluctuations may act in the downconversion of more than one quanta. In the literature on cluster and H‐graph states, the adjacency matrix for H‐graphs is defined via the matrix specified in the multimode squeezing Hamiltonian. The QLE analysis corresponds to the expansion of the multimode squeezing operator up to second order, which leads to the appearance of beamsplitter correlations in the adjacency matrix. In Appendix A.6 (Equations (A37)–(A39)) we show how to use well‐chosen relative pump phase values in the quadripartite case to prepare an entangled square lattice state—that is, a state without BS correlations. For the case of very large squeezing, H∼$\mathcal {\tilde{H}}$‐graph state can be regarded as an approximation of a 4‐node cluster state, minimizing errors in gate operations of measurement‐based CV quantum computing.Cluster states form a promising platform for scalable quantum information processing. In one‐way quantum computing,[3] the entire computational resource is provided by the entanglement of the cluster state. The processing is based on quantum measurements which facilitate gate operations as well as the read‐out of the final result. However, cluster states can be obtained from graph states only in the mathematical limit of large squeezing parameter.[27–29,32] For quantum information processing steps, it is sufficient to perform sub‐cluster measurements in specified order using a suitable computational basis. In refs. [28, 66, 67] different computation scenarios based on resources provided by squeezing generators and beamsplitters are described. Encoding, gate and measurement operations have been so far considered in optical circuits for continuous variable quantum data and can be efficiently extended to the microwave realm. In this work, we have utilized this correspondence between optics and microwaves and demonstrated H∼$\mathcal {\tilde{H}}$‐graph state encoding.In contrast to computational models for graph states[38] considered as ideal clusters, hardware based on finite squeezing with noise and decoherence requires error correction procedures[6,68,69] to provide reliable CV computation. Using the presented scheme one can implement error correction codes based on the idea of repetitions of selective measurements and new encoding of H∼$\mathcal {\tilde{H}}$‐graph states before each gate operation. In ref.[ [70]] a multidimensional platform for scalable quantum computing has been proposed, based on cluster states created using microring resonators; also multiple frequency combs[67] created by optical parametric amplifiers and beamsplitters can serve as an excellent platform for quantum computation. Our work shows that the methods of generation of highly‐entangled CV states are not restricted to just optical parametric amplifiers, but the methods can be carried over into the microwave domain by employing parametric Josephson junction devices for creation of topologically involved and structurally versatile H∼$\mathcal {\tilde{H}}$‐graph states.An implementation of the universal quantum computer based on bosonic modes with the possibility of hardware‐efficient quantum error correction[71] requires efficient generation of continuous‐variable quantum resources. The genuine entanglement between several bosonic modes could potentially be employed for error‐correctable codeword states.[72] Besides potential in error correction, the introduction of entanglement into quantum measurement implementations leads to a quantum advantage in the detection process when detection is performed in the presence of high level of noise and loss.[73]Increase of the number of entangled spectral modes is essential for future technological application of these CV quantum state generation methods. The limiting factors are the requirements of high precision for the pump frequency and its phase, the stability of the biasing flux, and possible crowding of modes within a narrow‐band JPA resonance. However, recently it has been demonstrated that entanglement can be generated in low‐loss traveling wave parametric amplifiers.[] This opens a way to significant increase in the number of entangled modes.ConclusionIn this work, we presented a practical scheme for generation of controllable multipartite entanglement from vacuum fluctuations, based on multitone pumping scheme of a JPA, which facilitates pivotal resources for quantum technologies at microwave frequencies. While optical schemes for multipartite entanglement generation operate on even larger clusters, they lack versatility and are limited to optical frequencies as such. On the other hand, our scheme allows for a flexible increase in the number of modes and control of the entanglement configuration among the modes by adjusting pumping on the same device, whereas optical setups call for massive hardware reconfiguration when the entanglement structure is altered. Through phase and amplitude variation of the microwave pump tones, we reach a comprehensive control over the entanglement structure within the spectral modes of a single JPA cavity mode, which we experimentally verify in detail for the tripartite case.Using the developed scheme, we made the first successful demonstration of an on‐demand tunable, fully inseparable, and phase controllable genuinely entangled tripartite and quadripartite states in a superconducting system. The presence of multipartite quantum correlations was verified using the covariance matrix formalism and genuine entanglement criteria constructed from the measured quadratures. Experimental results were accurately reproduced by calculating symplectic eigenvalues of a partially‐transposed covariance matrix for full inseparability detection as well as computing GME criteria in normalized pump amplitude in range 0<A<0.5$0&lt;A&lt;0.5$ (0<A<0.25$0&lt;A &lt;0.25$) and verified genuine entanglement in the range of 0.01≲A<0.4$0.01\lesssim A&lt;0.4$ (0.01≲A<0.13$0.01\lesssim A&lt;0.13$) for the tripartite (quadripartite) state.We provided results of phase‐dependent GME criterion for bisqueezed state. With optimal phase shift between two pumping tones Δφ=−120∘$\Delta \varphi =-120^\circ$ minimum value of criterion S=0.75±0.05$S=0.75 \pm 0.05$ was obtained. This result was also faithfully reproduced by numerical simulations.In our analytical derivations, we demonstrated additional control possibilities over the BS correlations in the covariance matrix of quadripartite H∼$\mathcal {\tilde{H}}$‐graph state. To visualize the formed entanglement structure, we provided an extension for the known H‐graph adjacency matrix: besides TMS, it includes BS correlations between the vacuum modes. The QLE approach was used to introduce such an adjacency matrix and to connect it to the general approach starting from multimode squeezing operator and the TMS Hamiltonian for the multi‐mode case with multiple pumps. As shown in Appendix A.6, BS correlations can be fully suppressed by implementing a 180° phase shift of one pump. Such a phase combination creates a distinct square‐lattice H‐graph state which, for the limit of infinite squeezing parameter, transforms to a square‐lattice cluster state.Additional TMS correlations can be introduced by inserting new pump tones, which can change the nature of the entangled states drastically. For example, using two additional pump tones with half frequencies at {−Δ2;Δ2}$\lbrace -\frac{\Delta }{2};\frac{\Delta }{2}\rbrace$, we are able to connect all 4 modes with TMS correlations and thereby achieve a GHZ‐like state. Furthermore, by tuning the phases of the pumps, the state can be converted into square‐lattice H‐graph state. With the bandwidth improvements provided by the state‐of‐the‐art superconducting parametric devices, such as the broadband, low‐loss travelling wave parametric amplifier,[62–64] we expect a substantial increase in the number of entangled modes, which facilitates generation of highly‐squeezed square‐lattice H‐graph states for CV quantum computation at microwave frequencies.AAppendixAnalytical Methods, Experimental Procedures, and Entanglement AnalysisA.1Details of Theoretical DescriptionThe Hamiltonian of JPA system is given byA1Ĥsys=ℏωrâ†â+ℏ2∑d=1pαd∗eiωdt+αde−iωdt(â2+â†2)+ℏK(â+â†)4$$\begin{eqnarray} {\hat{H}}_{\text{sys}}=\hbar \omega _r \hat{a}^\dagger \hat{a} & + \frac{\hbar }{2}\sum _{d=1}^{p}{\left[\alpha ^{*}_d e^{i\omega _d t} + \alpha _d e^{-i\omega _d t}\right]}(\hat{a}^{2} + \hat{a}^{\dagger 2}) \nonumber \\ & + \hbar K (\hat{a}+ \hat{a}^{\dagger })^{4} \end{eqnarray}$$where â(â†$\hat{a}(\hat{a}^{\dagger }$) is the annihilation (creation) operators for cavity photons, αd$\alpha _d$ is the complex amplitude for pump tone d, and K denotes the strength of the Kerr nonlinearity term. Using the average of p pump tones ωd$\omega _d$, d={1,⋯,p}$d=\lbrace 1,\dots ,p \rbrace$, we define the detuning between the half pump frequency and the resonator frequency: Δr=ωr−ωΣ2$\Delta _\text{r}=\omega _\text{r}-\frac{\omega _{\Sigma }}{2}$, ωΣ=∑d=1pωdp$\omega _{\Sigma }=\frac{\sum _{d=1}^{p}\omega _d}{p}$.For each of the p pump tones, we define the detuning from the average frequency Δd=ωd−ωΣ,d={1,⋯,p}$\Delta _d=\omega _d-\omega _{\Sigma }, d=\lbrace 1,\dots ,p \rbrace$. By applying the rotating wave approximation in the frame ωΣ/2$\omega _\Sigma /2$ (a∼(t)=â(t)eiωΣt/2$\tilde{a}(t) = \hat{a}(t)e^{i\omega _\Sigma t/2}$) and leaving only the effective high‐order terms, we obtain for the nonlinear part of the HamiltonianA2Hsys,rwa(t)=ℏΔra∼†a∼+ℏ2∑d=1p(αd∗eiΔdta∼2+αde−iΔdta∼†2)+6ℏKa∼†a∼†a∼a∼$$\begin{eqnarray} {H}_{\text{sys,rwa}}(t)={\hbar }\Delta _r \tilde{a}^\dagger \tilde{a} & + \frac{\hbar }{2}\sum _{d=1}^{p}(\alpha ^{*}_d e^{i\Delta _d t}\tilde{a}^{2} + \alpha _d e^{-i\Delta _d t} \tilde{a}^{\dagger 2}) \nonumber \\ & + 6{\hbar }{K} \tilde{a}^{\dagger }\tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$As usual, the bosonic commutation relationships are valid for the cavity modes [a∼,a∼†]=1$[ \tilde{a},\tilde{a}^{\dagger }]=1$.The parametric resonator is coupled to a transmission line via the signal port and to the thermal bath via a linear dissipation port. The coupling Hamiltonian associated with the signal port is given byA3Hsig=ℏ∫dωb∼†b∼+κa∼†b∼−κ∗b∼†a∼$$\begin{equation} H_{\text{sig}}=\hbar \int d\omega {\left(\tilde{b}^{\dagger }\tilde{b} + \kappa \tilde{a}^{\dagger } \tilde{b} - \kappa ^{*} \tilde{b}^{\dagger } \tilde{a}\right)} \end{equation}$$where creation and annihilation operators b∼†$\tilde{b}^{\dagger }$ and b∼$\tilde{b}$ refer to modes in the transmission line, and κ denotes the coupling rate. The Hamiltonian related to the linear dissipation portA4Hloss=ℏ∫dωc∼†c∼+γa∼†c∼−γ∗c∼†a∼$$\begin{equation} H_{\text{loss}}=\hbar \int d\omega {\left(\tilde{c}^{\dagger }\tilde{c} + \gamma \tilde{a}^{\dagger } \tilde{c} - \gamma ^{*} \tilde{c}^{\dagger }\tilde{a}\right)} \end{equation}$$where c∼†$\tilde{c}^{\dagger }$ and c∼$\tilde{c}$ describe creation and annihilation of thermal bath modes and the rate γ represent the coupling of cavity modes to the linear dissipation port. The transmission line and bath operators obey the commutation relationsA5b∼(ω),b∼†(ω′)=c∼(ω),c∼†(ω′)=δ(ω−ω′)$$\begin{equation} {\left[ \tilde{b}(\omega ),\tilde{b}^{\dagger }(\omega ^{\prime })\right]}={\left[ \tilde{c}(\omega ),\tilde{c}^{\dagger }(\omega ^{\prime })\right]}=\delta (\omega -\omega ^{\prime }) \end{equation}$$andA6b∼(ω),c∼†(ω′)=c∼(ω),b∼†(ω′)=b∼(ω),b∼(ω′)=c∼(ω),c∼(ω′)=0$$\begin{eqnarray} {\left[ \tilde{b}(\omega ),\tilde{c}^{\dagger }(\omega ^{\prime })\right]} = & {\left[ \tilde{c}(\omega ),\tilde{b}^{\dagger }(\omega ^{\prime })\right]} = {\left[ \tilde{b}(\omega ),\tilde{b}(\omega ^{\prime })\right]} \nonumber \\ & ={\left[ \tilde{c}(\omega ),\tilde{c}(\omega ^{\prime })\right]}=0 \end{eqnarray}$$The total Hamiltonian can be conveniently written as a sum of the separate parts given aboveA7H(t)=Hsys,rwa(t)+Hsig+Hloss$$\begin{equation} {H}(t)=H_{\text{sys,rwa}}(t)+H_{\text{sig}}+H_{\text{loss}} \end{equation}$$For further analysis and for our simulations, we use the Quantum Langevin Equation (QLE) for the cavity operator a∼(t)$\tilde{a}(t)$A8a∼̇(t)=(−iΔr−κ+γ2)a∼−i∑d=1pαdeiΔdta∼†+κb∼in+γc∼in−12iKa∼†a∼a∼$$\begin{eqnarray} \dot{\tilde{a}}(t)=(-i\Delta _r &- \frac{\kappa +\gamma }{2})\tilde{a} - i \sum _{d=1}^{p}\alpha _d e^{i\Delta _d t} \tilde{a}^\dagger \nonumber \\ & + \sqrt {\kappa }\tilde{b}_{in}+\sqrt {\gamma }\tilde{c}_{in} - 12i{K} \tilde{a}^{\dagger } \tilde{a} \tilde{a} \end{eqnarray}$$where the presence of the Kerr term allows us to consider dynamics of the parametric resonator above the critical oscillation threshold. To obtain the modes coming out from the cavity, we employ the standard input–output formalism which yields the relationshipA9b∼out(t)=b∼in(t)−κa∼(t)$$\begin{equation} \tilde{b}_{\text{out}}(t)=\tilde{b}_{\text{in}}(t)-\sqrt {\kappa }\tilde{a}(t) \end{equation}$$Equations (A8) and (A9) are used in our numerical simulations with Matlab ODE45 solver.A.2Full InseparabilityAssuming that the microwave fields produced by the JPA below the critical pumping threshold are Gaussian,[34] the states with multiple spectral modes can be fully characterized by measuring the covariance matrix of corresponding in‐phase I and quadrature Q voltages. For the measurement of tripartite correlations, we collect quadrature data for 0.8 s at every phase difference and pump amplitude value, without averaging. For the quadripartite case, we repeat the experiment 20 times at each pump power value and every quadrature sequence has a duration of 1.3 s.The quantum quadratures x∼i=a∼i+a∼i†2$ \tilde{x}_i = \frac{\tilde{a}_i+\tilde{a}_i^{\dagger }}{2}$ and p∼i=a∼i−a∼i†2i$\tilde{p}_i=\frac{\tilde{a}_i-\tilde{a}_i^{\dagger }}{2i}$ can be combined into a 2N‐long column vector operator for the N‐mode state r∼=(x∼1,p∼1,⋯x∼N,p∼N)T$\tilde{r} = (\tilde{x}_1, \tilde{p}_1, \dots \tilde{x}_N, \tilde{p}_N )^{T}$. The commutation relations can be written down in a skew‐symmetric, block‐diagonal matrix form[50]A10r∼i,r∼j=i2ΩijandΩ=⨂i=1N01−10$$\begin{equation} {\left[\tilde{r}_i, \tilde{r}_j\right]}=\frac{i}{2}\Omega _{ij} \text{ and } \mathbf {\Omega } = \bigotimes _{i=1}^N \def\eqcellsep{&}\begin{bmatrix} 0 & 1\\ -1 & 0 \end{bmatrix} \end{equation}$$The covariance matrix V$\mathbf {V}$ is given by elements Vij=12<Δr∼iΔr∼j+Δr∼jΔr∼i>−<Δr∼i><Δr∼j>$V_{ij}=\frac{1}{2}&lt;\Delta \tilde{r}_i \Delta \tilde{r}_j + \Delta \tilde{r}_j \Delta \tilde{r}_i &gt; - &lt;\Delta \tilde{r}_i &gt; &lt;\Delta \tilde{r}_j&gt;$ where we have defined standard error Δr∼i=r∼i−<r∼i>$\Delta \tilde{r}_i=\tilde{r}_i-&lt;\tilde{r}_i&gt;$ and <r∼i>=tr(r∼iρ̂)$&lt;\tilde{r}_i&gt;= \text{tr}(\tilde{r}_i \hat{\rho })$. The uncertainty principle requires thatA11V+i4Ω≥0$$\begin{equation} \mathbf {V}+\frac{i}{4}\mathbf {\Omega }\ge 0 \end{equation}$$applies for a physical covariance matrix.For verification of entanglement, we may investigate a modified equationA12V′+i4Ω≥0$$\begin{equation} \mathbf {V}^\prime +\frac{i}{4}\mathbf {\Omega }\ge 0 \end{equation}$$where Vk′=λkVλk$\mathbf {V}^\prime _k=\pmb {\lambda }_k \mathbf {V} \pmb {\lambda }_k$, λk$\pmb {\lambda }_{k}$ is diagonal matrix with ones entries, except of that related to kth mode, with value of −1. For example, transformation with λk≡N=diag(1,1,…,1,−1)$\pmb {\lambda }_{k\equiv N}=\text{diag}(1,1,\ldots,1,-1)$ means a partial transposition of the covariance matrix with respect to the last mode. The positive partial transpose (PPT) criterion for multipartite case requires that there is a violation of Equation (A12) when applying a partial transposition with respect to each from full set of modes: Vk′≥i4Ω$\mathbf {V}^\prime _k \ge \frac{i}{4}\mathbf {\Omega }$. In ref. [47], the entangled states are classified in accordance to the number of modes for which the condition Equation (A12) is broken. We follow this approach to demonstrate the highest class—full inseparability—in four mode case.Unitary operations which retain the Gaussian character of the states, for example, squeezing, are of particular importance. Such operations on the Hilbert space correspond to a linear transformation P in the phase‐space which preserve the symplectic form, that is,A13Ω=PTΩP$$\begin{equation} \mathbf {\Omega } = \mathbf {P}^T \mathbf {\Omega } \mathbf {P} \end{equation}$$Symplectic transformations on a 2N‐dimensional phase‐space form the real symplectic group denoted as Sp(2N;R)$ Sp(2N; {\bf R})$, which is a proper subgroup of the special linear group of 2N × 2N matrices.[74] By utilizing Williamson's theorem,[75] any covariance matrix can be expressed in the Williamson normal form:A14V∼k=PTVk′P$$\begin{equation} \mathbf {\tilde{V}}_k= \mathbf {P}^T \mathbf {V}^\prime _k \mathbf {P} \end{equation}$$where V∼k$\mathbf {\tilde{V}}_k$ is a 2N‐dimensional diagonal matrix consisting of the symplectic eigenvalues, ν∼k$\tilde{\nu }_k$, of the covariance matrix. The symplectic eigenvalues are called the symplectic spectrum which provides a practical means to verify physicality and various entanglement criteria. Separability is in force, when condition ν∼k≥1/4$\tilde{\nu }_k \ge 1/4$ fulfilled for V∼k$\mathbf {\tilde{V}}_k$.For convenience, we insert an additional factor of 4 to the covariance matrix and work with fluctuations with zero mean values: Vij∗=2<Δr∼iΔr∼j+Δr∼jΔr∼i>${V}^*_{ij}=2&lt;\Delta \tilde{r}_i \Delta \tilde{r}_j + \Delta \tilde{r}_j \Delta \tilde{r}_i &gt;$. Consequently, for evidence of “fully inseparable” states, we need to find minimum symplectic eigenvalues with ν∼k∗<1$\tilde{\nu }_k^* &lt; 1$ for each partial transposition k.A.3System Gain CalibrationOur system gain calibration procedure consists of a measurement of Johnson–Nyquist noise spectral density emitted by a 50 Ω termination at different temperatures. Assuming perfect matching of the source and load impedances, the received power per unit of bandwidth can be written by applying the Friis formula: the measured noise is given by the noise temperature of the source Ts$T_\text{s}$, the contribution of the cooled amplifier THEMT$T_{\text{HEMT}}$, and the noise of the room‐temperature amplifiers TRT$T_{\text{RT}}$ multiplied by the system gain GΣ,i=GHEMTiGRTi$G_{\Sigma , i}=G_{\text{HEMT}_i}G_{\text{RT}_i}$:A15Ii2+Qi2Z0Δfi=kbGΣ,iTi+THEMT+TRTGHEMT$$\begin{equation} \frac{\left&lt; I_i^2 + Q_i^2 \right&gt;}{ Z_0 \Delta f_i } = k_b G_{\Sigma ,i} {\left(T_i+T_{\text{HEMT}}+\frac{T_{\text{RT}}}{G_{\text{HEMT}}}\right)} \end{equation}$$Here i refers to the frequency of the spectral mode and Δfi$\Delta f_i$ refers to the bandwidth of the detection of quadratures Ii$I_i$ and Qi$Q_i$. The total gain GΣ,i$G_{\Sigma , i}$ was separately determined for different spectral modes.Figure A1 displays the measured noise power per unit band as a function of sample temperature Ts$T_\text{s}$, averaged over frequencies covering the resonance curve. By fitting a line to the data, we obtain <GΣ,i>=94.4±0.2$&lt;G_{\Sigma , i}&gt;=94.4 \pm 0.2$ dB for the average total gain. The linear fit in Figure A1 is performed at T>0.2$T&gt;0.2$ K, which allows us to neglect the corrections from the coth(ℏω/2kbTs)$\coth (\hbar \omega /2k_\text{b} T_\text{s})$.A1FigureGain calibration using linear temperature dependence of the measured thermal noise spectral density of a 50 Ω terminator measured as a function of Ts$T_\text{s}$, the source temperature. The average total gain <GΣ,i>=94.4±0.2$&lt;G_{\Sigma , i}&gt;=94.4 \pm 0.2$ dB over the cavity resonance is obtained from the linear fit (in red) to the data. This gain value <GΣ,i>$&lt;G_{\Sigma , i}&gt;$ also includes frequency mixing losses/amplification in the signal analyzer circuit part. The term Tpreamp=THEMT+TRTGHEMT=5.2±0.25$T_{\text{preamp}}=T_{\text{HEMT}}+\frac{T_\text{RT}}{G_\text{HEMT}}=5.2\pm 0.25$K characterizes the equivalent noise temperature of the amplifiers; the largest contribution originates from the cooled HEMT amplifier at 4 K. The value coth(hfi2kbTpreamp)$\coth (\frac{hf_i}{2 k_b T_\text{preamp}})$ sets the background for the diagonal elements in the covariance matrix 4Voff$4\mathbf {V}_\text{off}$.The error in the system gain calibration results in uncertainty in the symplectic eigenvalues on the order of 2%, that is, the eigenvalues fall in the range of min{ν∼k∗}$\min \lbrace \tilde{\nu }_k^*\rbrace$=min{ν∼k∗}±0.018$\min \lbrace \tilde{\nu }_k^*\rbrace \pm 0.018$ for each partial bipartition. Random variations of the system parameters were reduced by averaging the outcome by ten to twenty times.A.4System Parameter FittingIn order to determine coupling rates γ and κ introduced in Section 2 A, we characterized our nonlinear resonator as a two‐port device using a vector network analyzer. For the characterization, we chose the optimal DC operating point ΦDC=0.383Φ0$\Phi _{\text{DC}} = 0.383\,\Phi _0$ depicted in Figure 3b. At this DC‐flux, we measured the resonance curve in the absence of the pump in order to estimate the external and internal loss rates κ and γ, respectively. By fitting the measured resonance curve to the analytical solution of the QLE (b∼out(ω)/b∼in(ω)$\tilde{b}_{\text{out}}(\omega )/\tilde{b}_{\text{in}}(\omega )$), derived for the linear case without any pump drive, we obtain the coupling coefficients κ2π=4.44$\frac{\kappa }{2\pi } = 4.44$ MHz and γ2π=2.30$\frac{\gamma }{2\pi }=2.30$ MHz. The employed analytical solution, displayed in Equation (A16), was derived from the full QLE in Equation (A8) without taking the nonlinear part −iKa∼†a∼a∼$-iK \tilde{a}^{\dagger } \tilde{a} \tilde{a}$ into accountA16b∼out(ω)b∼in(ω)=1−κ(−i(ω−ωr)+κ+γ2)$$\begin{equation} \frac{\tilde{b}_{\text{out}}(\omega )}{\tilde{b}_{\text{in}}(\omega )}=1-\frac{\kappa }{(-i(\omega -\omega _\text{r}) +\frac{\kappa +\gamma }{2})} \end{equation}$$For fitting of the Kerr constant K, we employed the whole form of the QLE in the rotating wave approximation Equation (A8). By comparing the measured and simulated gain coefficients G(ωprobe−ωr,A)$G(\omega _\text{probe}-\omega _\text{r},A)$ (Figure A2) in the cavity at large pump amplitudes, we obtain an estimate K=6.5ωr$K = 6.5\omega _\text{r}$ for the Kerr constant.A2Figurea,b) Gain coefficient measurement and QLE's simulation results (Equation (A8)) which are used for fitting coupling and Kerr constants. Vertical axis represents normalized pump amplitude A=ακ+γ$A=\frac{\alpha }{\kappa +\gamma }$. Horizontal axis represents detuning between probe signal frequency and resonance frequency ωprobe−ωr2π$\frac{\omega _\text{probe} - \omega _\text{r}}{2\pi }$. Pumping carried out in degenerate mode, ωd=2ωr$\omega _d=2\omega _\text{r}$. Fano resonance picture, given in experimental plot, explained by phase shift between the cavity and input modes and described by complex rate value of κ.A.5Cavity Phase ResponseOptimal value of GME criterion, which governed by “symmetric” covariance matrix view in tripartite mode case, can be obtained with {π2;π2}$\lbrace \frac{\pi }{2};\frac{\pi }{2}\rbrace$ only if modes reshuffling suffers no additional phase rotations (see next subsection). However, in experiments we deal with finite values of coupling and dissipation loss rates. Cavity phase response becomes crucial figure in pump tones phase shift adjustments. Cavity phase response illustrated on Figure A3.A3FigureCavity phase response given by fitted experimental parameters κ2π=4.44$\frac{\kappa }{2\pi } = 4.44$ MHz and γ2π=2.30$\frac{\gamma }{2\pi }=2.30$ MHz. Vertical dash lines show center frequencies of first and last modes. Corresponding phase shifts applied to pump tones to reach “symmetric” covariance matrix view on double frequencies are Δφ1=π2−π4=π4$\Delta \varphi _1=\frac{\pi }{2}-\frac{\pi }{4}=\frac{\pi }{4}$ and Δφ2=π2+π4=3π4$\Delta \varphi _2=\frac{\pi }{2}+\frac{\pi }{4}=\frac{3\pi }{4}$ with corresponding phase shift between pump tones π2$\frac{\pi }{2}$, given in results (see Figure 4). Half of applied phase shift described by pump tones on double resonance frequencies. Additional phase shift with increase of A relates to modification of phase response curve during pumping.A.6Multifrequency Correlations in Terms of QLE with 3 and 4 Spectral ModesAs discussed in Section 2 C, our measurement setting probes outgoing waves from the parametric resonator, which brings about slight differences with standard quantum optics schemes where the entanglement analysis is based on the Hamiltonian of the system. In our case, the QLE provides a good description, and here we derive the relevant matrix equations describing the coupling of the different outgoing spectral modes under two and three pump tones (3 and 4 spectral modes, respectively).3 Mode Case: Let us define a∼$\tilde{a}$ as a vector of spectral modes:A17a∼={a∼1,a∼2,a∼3,a∼1†,a∼2†,a∼3†}T;$$\begin{equation} \tilde{a}=\lbrace {\tilde{a}_1,\tilde{a}_2,\tilde{a}_3,\tilde{a}^\dagger _1,\tilde{a}^\dagger _2,\tilde{a}^\dagger _3\rbrace }^{T}; \end{equation}$$where the creation a∼i†=a∼i†(t)$\tilde{a}_i^{\dagger }=\tilde{a}_i^{\dagger }(t)$ and annihilation a∼i=a∼i(t)$\tilde{a}_i=\tilde{a}_i(t)$ operators are time‐dependent. After Fourier transform,A18a∼(ω)={a∼1(ω),a∼2(ω),a∼3(ω),a∼1†(−ω),a∼2†(−ω),a∼3†(−ω)}T.$$\begin{equation} \tilde{a}(\omega )=\lbrace {\tilde{a}_1(\omega ),\tilde{a}_2(\omega ),\tilde{a}_3(\omega ),\tilde{a}^\dagger _1(-\omega ),\tilde{a}^\dagger _2(-\omega ),\tilde{a}^\dagger _3(-\omega )\rbrace }^{T}. \end{equation}$$We define our spectral modes ai$a_i$ as {(−3Δ2,−Δ2);(−Δ2,Δ2);(Δ2,3Δ2)}$\lbrace (-\frac{3\Delta }{2},-\frac{\Delta }{2});(-\frac{\Delta }{2},\frac{\Delta }{2});(\frac{\Delta }{2},\frac{3\Delta }{2})\rbrace$ according pump tone positions {−Δ,Δ}$\lbrace -\Delta ,\Delta \rbrace$ (see Figure 1). Similarly, we define for the input and output modes b̂in/out$\hat{b}_{\text{in}/\text{out}}$:A19b∼in/out={b∼in1/out1,b∼in2/out2,b∼in3/out3,b∼in2/out2†,b∼in1/out1†,b∼in3/out3†T.$$\begin{eqnarray} \tilde{b}_{\text{in/out}}=& \lbrace {\tilde{b}_{\text{in}1/\text{out}1},\tilde{b}_{\text{in}2/\text{out}2},\tilde{b}_{\text{in}3/\text{out}3}}, \nonumber \\ &{\left. \tilde{b}^\dagger _{\text{in}2/\text{out}2},\tilde{b}^\dagger _{\text{in}1/\text{out}1},\tilde{b}^\dagger _{\text{in}3/\text{out}3}\right\rbrace} ^{T}. \end{eqnarray}$$The commutation relationships for the case of N modes can be conveniently expressed in matrix form. We use the common convention for [a∼i,a∼j]$[\tilde{a}_i,\tilde{a}_j ]$ from ref. [74].The effect of Kerr nonlinearity is significant only at large pump amplitudes. Hence, we may take the QLE Equation (A8) without the nonlinear part for our treatment. In theoretical analysis we assume, that spectral modes lay down deep in cavity mode, such that Δ≪κ$\Delta \ll \kappa$; we also neglect internal dissipation expressed by loss rate γ. For that case phase shift between modes, provided by phase response of the cavity, can be neglected. Guided by standard Fourier transform technique for solving linear QLE,[40] we denote a∼i(ω)=∫a∼i(t)eiωtdt$\tilde{a}_i(\omega )=\int \tilde{a}_i(t) e^{i\omega t}\text{d}t$ and Fourier transform the QLE term by term. Owing to detuning of the pump tones in the rotating frame, there will be coupling of spectral modes and we have mode index exchange. For example, for a∼1,2†$\tilde{a}^\dagger _{1,2}$: ∫a∼1†(t)eiωteiΔ1tdt=a∼2†(−ω)$\int \tilde{a}^\dagger _1(t)e^{i\omega t}e^{i\Delta _1 t}\text{d}t=\tilde{a}^{\dagger }_2(-\omega )$; ∫a∼2†(t)eiωteiΔ2tdt=a∼1†(−ω)$\int \tilde{a}^\dagger _2(t)e^{i\omega t}e^{i\Delta _2 t}\text{d}t=\tilde{a}^{\dagger }_1(-\omega )$, while for a∼2,3†$\tilde{a}^\dagger _{2,3}$: ∫a∼2†(t)eiωte−iΔ2tdt=a∼3†(−ω)$\int \tilde{a}^\dagger _2(t)e^{i\omega t}e^{-i\Delta _2 t}\text{d}t=\tilde{a}^{\dagger }_3(-\omega )$; ∫a∼3†(t)eiωte−iΔ3tdt=a∼2†(−ω)$\int \tilde{a}^\dagger _3(t)e^{i\omega t}e^{-i\Delta _3 t}\text{d}t=\tilde{a}^{\dagger }_2(-\omega )$. Thus, it is seen that each pump creates a two‐mode squeezed state (TMS) between two neighboring spectral modes independently from RWA's zero‐frequency position.After Fourier transforming,A20(−i(ω−Δr)+κ2)a∼(ω)+iα(∫a∼†(t)eiωte−iΔdtdt+∫a∼†(t)eiωteiΔdtdt)=κb∼in(ω)$$\begin{eqnarray} (-i(\omega -\Delta _r)+&\frac{\kappa }{2})\tilde{a}(\omega )+i\alpha (\int \tilde{a}^\dagger (t)e^{i\omega t}e^{-i\Delta _d t}\text{d}t+ \nonumber \\ &\int \tilde{a}^\dagger (t)e^{i\omega t}e^{i\Delta _d t}\text{d}t)=\sqrt {\kappa }\tilde{b}_{\text{in}}(\omega ) \end{eqnarray}$$the QLE yields the following system of linear equations:A21κb∼in1(ω)=(−i(ω−Δr)+κ2)a∼1(ω)+iαa∼2†(−ω)κb∼in2(ω)=(−i(ω−Δr)+κ2)a∼2(ω)+iα(a∼1†(−ω)+a∼3†(−ω))κb∼in3(ω)=(−i(ω−Δr)+κ2)a∼3(ω)+iαa∼2†(−ω)κb∼in1†(−ω)=(−i(ω+Δr)+κ2)a∼1†(−ω)−iα†a∼2(ω)κb∼in2†(−ω)=(−i(ω+Δr)+κ2)a∼2†(−ω)−iα†(a∼1(ω)+a∼3(ω))κb∼in3†(−ω)=(−i(ω+Δr)+κ2)a∼3†(−ω)−iα†a∼2(ω).$$\begin{eqnarray} \sqrt {\kappa }\tilde{b}_{\text{in}1}(\omega )=(-i(\omega -\Delta _r)+\frac{\kappa }{2})\tilde{a}_1(\omega )+i\alpha \tilde{a}^\dagger _2(-\omega ) \nonumber \\ \sqrt {\kappa }\tilde{b}_{\text{in}2}(\omega )=(-i(\omega -\Delta _r)+\frac{\kappa }{2})\tilde{a}_2(\omega )+i\alpha (\tilde{a}^\dagger _1(-\omega )+\tilde{a}^\dagger _3(-\omega )) \nonumber \\ \sqrt {\kappa }\tilde{b}_{\text{in}3}(\omega )=(-i(\omega -\Delta _r)+\frac{\kappa }{2})\tilde{a}_3(\omega )+i\alpha \tilde{a}^\dagger _2(-\omega ) \nonumber \\ \sqrt {\kappa }\tilde{b}^\dagger _{\text{in}1}(-\omega )=(-i(\omega +\Delta _r)+\frac{\kappa }{2})\tilde{a}^\dagger _1(-\omega )-i\alpha ^\dagger \tilde{a}_2(\omega ) \nonumber \\ \sqrt {\kappa }\tilde{b}^\dagger _{\text{in}2}(-\omega )=(-i(\omega +\Delta _r)+\frac{\kappa }{2})\tilde{a}^\dagger _2(-\omega )-i\alpha ^\dagger (\tilde{a}_1(\omega )+\tilde{a}_3(\omega )) \nonumber \\ \sqrt {\kappa }\tilde{b}^\dagger _{\text{in}3}(-\omega )=(-i(\omega +\Delta _r)+\frac{\kappa }{2})\tilde{a}^\dagger _3(-\omega )-i\alpha ^\dagger \tilde{a}_2(\omega ). \end{eqnarray}$$We cast Equation (A21) into matrix form:A22Ma∼(ω)=κb∼in(ω)$$\begin{equation} \mathbf {M}\tilde{a}(\omega )=\sqrt {\kappa }\tilde{b}_{\text{in}}(\omega ) \end{equation}$$A23M=c1000iα00c10iα0iα00c10iα00−iα†0c200−iα†0−iα†0c200−iα†000c2$$\begin{equation} \mathbf {M}=\def\eqcellsep{&}\begin{bmatrix} c_1 & 0 & 0 & 0 & i\alpha & 0\\ 0 & c_1 & 0 & i\alpha & 0 & i\alpha \\ 0 & 0 & c_1 & 0 & i\alpha & 0\\ 0 & -i\alpha ^\dagger & 0 & c_2 & 0 & 0\\ -i\alpha ^\dagger & 0 & -i\alpha ^\dagger & 0 & c_2 & 0\\ 0 & -i\alpha ^\dagger & 0 & 0 & 0 & c_2 \end{bmatrix} \end{equation}$$where c1=−i(ω−Δr)+κ2$c_1=-i(\omega -\Delta _r)+\frac{\kappa }{2}$ and c2=−i(ω+Δr)+κ2$c_2=-i(\omega +\Delta _r)+\frac{\kappa }{2}$ Solving for the inverse of M̂$\hat{M}$ and using Equation (A9), we obtainA24b∼out(ω)=(I−κM−1)b∼in(ω)$$\begin{equation} \tilde{b}_{\text{out}}(\omega )=(\mathbf {I}-\kappa \mathbf {M}^{-1})\tilde{b}_{\text{in}}(\omega ) \end{equation}$$for the outgoing radiation b∼out(ω)$\tilde{b}_{\text{out}}(\omega )$ in terms of incoming waves b∼in(ω)$\tilde{b}_{\text{in}}(\omega )$.Because our goal is to determine the structure of the experimental covariance matrix, it is unsatisfactory to consider cavity modes a∼$\tilde{a}$ with equation a∼(ω)=κM−1b∼in(ω)$\tilde{a}(\omega )=\sqrt {\kappa }\mathbf {M}^{-1}\tilde{b}_{\text{in}}(\omega )$ though it has a more compact final form. However, the presence of the identity matrix I$\mathbf {I}$ and the multiplication factor κ do not change the final structure.Assuming that the pump amplitude α is a real number and c1=c2=c$c_1=c_2=c$ (zero detuning case), we haveA25M−1=1c2−2α2·c−α2c0α2c0−iα00c0−iα0−iαα2c0c−α2c0−iα00iα0c−α2c0α2ciα0iα0c00iα0α2c0c−α2c$$\begin{eqnarray} \mathbf {M}^{-1}= \frac{1}{c^2-2\alpha ^2} \nonumber \cdot \\ \def\eqcellsep{&}\begin{bmatrix} c-\frac{\alpha ^2}{c} & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & -i\alpha & 0\\ 0 & c & 0 & -i\alpha & 0 & -i\alpha \\ \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} & 0 & -i\alpha & 0\\ 0 & i\alpha & 0 & c-\frac{\alpha ^2}{c} & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}}\\ i\alpha & 0 & i\alpha & 0 & c & 0\\ 0 & i\alpha & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} \end{bmatrix} \end{eqnarray}$$This allows us to draw the generalized H∼$\mathcal {\tilde{H}}$‐graph for describing the parametric interaction between the spectral modes, Figure 2. The off‐diagonal beamsplitter elements proportional to α2 are set in bold in Equation (A25).Still, we want to construct the parametric interaction matrix S−1$\mathbf {S}^{-1}$ for quadrature vector operator r∼$\tilde{r}$. Using a linear operator matrix K$\mathbf {K}$ to implement a change of basisA26K=12100100−i00i000100100−i00i000100100−i00i$$\begin{eqnarray} \mathbf {K}= \frac{1}{2} \def\eqcellsep{&}\begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0\\ -i & 0 & 0 & i & 0 & 0\\ 0 & 1 & 0 & 0 & 1 & 0\\ 0 & -i & 0 & 0 & i & 0 \\ 0 & 0 & 1 & 0 & 0 & 1\\ 0 & 0 & -i & 0 & 0 & i \end{bmatrix} \end{eqnarray}$$we obtain by a canonical transformation S−1=κKM−1K−1$\mathbf {S}^{-1}=\sqrt {\kappa }\mathbf {K}\mathbf {M}^{-1}\mathbf {K}^{-1}$:A27S−1=κc2−2α2·c+α2c00−αα2c00c+α2c−α00α2c0−αc00−α−α00c−α0α2c00−αc+α2c00α2c−α00c+α2c$$\begin{eqnarray} &&\mathbf {S}^{-1}= \frac{\sqrt {\kappa }}{c^2-2\alpha ^2} \nonumber \cdot \\ &&\def\eqcellsep{&}\begin{bmatrix} c+\frac{\alpha ^2}{c} & 0 & 0 & -\alpha & \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 \\ 0 & c+\frac{\alpha ^2}{c} & -\alpha & 0 & 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}}\\ 0 & -\alpha & c & 0 & 0 & -\alpha \\ -\alpha & 0 & 0 & c & -\alpha & 0 \\ \mathbf {\frac{\pmb {\alpha }^2}{c}} & 0 & 0 & -\alpha & c+\frac{\alpha ^2}{c} & 0\\ 0 & \mathbf {\frac{\pmb {\alpha }^2}{c}} & -\alpha & 0 & 0 & c+\frac{\alpha ^2}{c} \end{bmatrix} \end{eqnarray}$$Note, that here the overall structure of the matrix has changed because of the basis change from ladder to quadrature operators. This is seen, for example, in the distribution of the off‐diagonal beamsplitter correlations (shown in bold).Since the environment of the cavity is in the ground state, b∼in$\tilde{b}_\text{in}$ has a Gaussian covariance matrix of the form Vin=14I$\mathbf {V}_\text{in}=\frac{1}{4}\mathbf {I}$. Consequently, the covariance matrix of the cavity spectral modes a∼i$\tilde{a}_i$ can be represented asA28Va=S−1Vin(S−1)T$$\begin{equation} \mathbf {V}_\text{a}=\mathbf {S}^{-1}\mathbf {V}_\text{in}(\mathbf {S}^{-1})^{T} \end{equation}$$or, equivalently, for output modes b∼out$\tilde{b}_\text{out}$: Vout=(I−κS−1)Vin(I−κS−1)T$\mathbf {V}_\text{out}= (\mathbf {I}-\sqrt {\kappa }\mathbf {S}^{-1})\mathbf {V}_\text{in}(\mathbf {I}-\sqrt {\kappa }\mathbf {S}^{-1})^{T}$. Both forms Va$\mathbf {V}_\text{a}$ and Vout$\mathbf {V}_\text{out}$ can be employed for studying the structure of parametric interactions between the quadratures, because input–output relationship doesnot change the general structure of the couplings between the quadratures (see below).As shown in Section 4 experimentally, phase shift between pumps changes the appearance of the covariance matrix (see Figure 5) as well as the strength of genuine multipartite entanglement. A change in the matrix M$\mathbf {M}$ due to a phase shift is illustrated in Equation (A29), in which the phase of the first pump has been rotated by eiπ2$e^{\frac{i\pi }{2}}$.A29M̂=c000−α00c0−α0iα00c0iα00−α0c00−α0−iα0c00−iα000c$$\begin{equation} \hat{M}=\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & -\pmb {\alpha } & 0\\ 0 & c & 0 & -\pmb {\alpha } & 0 & i\alpha \\ 0 & 0 & c & 0 & i\alpha & 0\\ 0 & -\pmb {\alpha } & 0 & c & 0 & 0\\ -\pmb {\alpha } & 0 & -i\alpha & 0 & c & 0\\ 0 & -i\alpha & 0 & 0 & 0 & c \end{bmatrix} \end{equation}$$The elements affected by the rotation are indicated in bold in the matrix. The elements in bold face indicate coupling between modes a∼1(ω)↔a∼2(ω)$\tilde{a}_1(\omega )\leftrightarrow \tilde{a}_2(\omega )$ while the other off‐diagonal elements indicate squeezing across a∼2(ω)↔a∼3(ω)$\tilde{a}_2(\omega )\leftrightarrow \tilde{a}_3(\omega )$. Note, that phase rotation operates in opposite direction on rows related to b∼in(ω)$\tilde{b}_{\text{in}}(\omega )$ and b∼in†(ω)$\tilde{b}_{\text{in}}^{\dagger }(\omega )$.The inversion of the rotated matrix M$\mathbf {M}$ yields for the parametric interaction matrix, where all the beamsplitter elements (in bold) have acquired a π/2$\pi /2$ phase shift. This phase shift can be unwound by a phase shift on the second pump, which indicates different phase dependence of the beamsplitter correlations compared with the TMS correlations. The structure of matrix M−1$\mathbf {M}^{-1}$ in Equation (A30) shows that phase rotation of specified pump tones does not change parametric interaction form between modes, preserving structure of a bisqueezed tripartite state. However, as shown in the main text, the criterion describing the strength of GME (see Equation (15)) depends on the difference of pump phases and strong genuine entanglement is reached only at specific phase settings.A30M−1=1c2−2α2·c−α2c0iα2c0α00c0α0−iα−iα2c0c−α2c0−iα00α0c−α2c0−iα2cα0iα0c00iα0iα2c0c−α2c$$\begin{eqnarray} \mathbf {M}^{-1}= \frac{1}{c^2-2\alpha ^2}\cdot \nonumber \\ \def\eqcellsep{&}\begin{bmatrix} c-\frac{\alpha ^2}{c} & 0 & \mathbf {\frac{\pmb {i\alpha }^2}{c}} & 0 & \alpha & 0\\ 0 & c & 0 & \alpha & 0 & -i\alpha \\ -\mathbf {\frac{\pmb {i\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} & 0 & -i\alpha & 0\\ 0 & \alpha & 0 & c-\frac{\alpha ^2}{c} & 0 & -\mathbf {\frac{\pmb {i\alpha }^2}{c}}\\ \alpha & 0 & i\alpha & 0 & c & 0\\ 0 & i\alpha & 0 & \mathbf {\frac{\pmb {i\alpha }^2}{c}} & 0 & c-\frac{\alpha ^2}{c} \end{bmatrix} \end{eqnarray}$$The covariance matrix Va$\mathbf {V}_\text{a}$ obtained from matrix M$\mathbf {M}$ in Equation (A25) with zero pump phase shifts is given in Equation (A31). The corresponding covariance matrix for π/2$\pi /2$ phase rotation in the first pump is displayed in Equation (A32). The matrix Va$\mathbf {V}_\text{a}$ in Equation (A32) has one rotated subspace, corresponding to two quadrature pairs; these rotated components are indicated by bold face. Based on these analytical relationships we conclude that control over desired covariance matrix TMS‐subspace can be provided by phase rotation of corresponding pump tone. Finally, we introduce the same phase rotation eiπ2$e^{\frac{i\pi }{2}}$ to the second pump. This brings the covariance matrix for the spectral cavity modes to the “standard‐symmetric” form displayed in Equation (A33). By comparing Equation (A31) and Equation (A33) we note that the beamsplitter elements (in bold) in the covariance matrix are unchanged (the phase difference between the pumps is the same) while the TMS elements are different.A31Va=κ4(c2−2α2)2·−α2+2α4c2+c200−2αc3α2−2α4c200−α2+2α4c2+c2−2αc003α2−2α4c20−2αc2α2+c200−2αc−2αc002α2+c2−2αc03α2−2α4c200−2αc−α2+2α4c2+c2003α2−2α4c2−2αc00−α2+2α4c2+c2$$\begin{equation} \mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \def\eqcellsep{&}\begin{bmatrix} -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & 0 & -2 \alpha c & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 \\ 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & -2 \alpha c & 0 & 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} \\ 0 & -2 \alpha c & 2 \alpha ^2+c^2 & 0 & 0 & -2 \alpha c \\ -2 \alpha c & 0 & 0 & 2 \alpha ^2+c^2 & -2 \alpha c & 0 \\ \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & 0 & -2 \alpha c & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 \\ 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & -2 \alpha c & 0 & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 \\ \end{bmatrix} \end{equation}$$A32Va=κ4(c2−2α2)2·−α2+2α4c2+c202αc00−3α2+2α4c20−α2+2α4c2+c20−2αc3α2−2α4c202αc02α2+c200−2αc0−2αc02α2+c2−2αc003α2−2α4c20−2αc−α2+2α4c2+c20−3α2+2α4c20−2αc00−α2+2α4c2+c2$$\begin{equation} \mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \def\eqcellsep{&}\begin{bmatrix} -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & \pmb {2 \alpha c} & 0 & 0 & \pmb {-3 \alpha ^2+\frac{2 \alpha ^4}{c^2}} \\ 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & \pmb {-2 \alpha c} & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 \\ \pmb {2 \alpha c} & 0& 2 \alpha ^2+c^2 & 0 & 0 & -2 \alpha c \\ 0 & \pmb {-2 \alpha c} & 0 & 2 \alpha ^2+c^2 & -2 \alpha c & 0 \\ 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & -2 \alpha c & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 \\ \pmb {-3 \alpha ^2+\frac{2 \alpha ^4}{c^2}} & 0 &-2 \alpha c & 0 & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 \\ \end{bmatrix} \end{equation}$$A33Va=κ4(c2−2α2)2·−α2+2α4c2+c202αc03α2−2α4c200−α2+2α4c2+c20−2αc03α2−2α4c22αc02α2+c202αc00−2αc02α2+c20−2αc3α2−2α4c202αc0−α2+2α4c2+c2003α2−2α4c20−2αc0−α2+2α4c2+c2.$$\begin{equation} \mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \def\eqcellsep{&}\begin{bmatrix} -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 \\ 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c & 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} \\ 2 \alpha c & 0 & 2 \alpha ^2+c^2 & 0 & 2 \alpha c & 0 \\ 0 & -2 \alpha c & 0 & 2 \alpha ^2+c^2 & 0 & -2 \alpha c \\ \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & 2 \alpha c & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 & 0 \\ 0 & \pmb {3 \alpha ^2-\frac{2 \alpha ^4}{c^2}} & 0 & -2 \alpha c & 0 & -\alpha ^2+\frac{2 \alpha ^4}{c^2}+c^2 \\ \end{bmatrix}. \end{equation}$$4 Mode Case: We treat the case with three pump tones (p=3$p=3$) in the same way as we did with two pumps above. By taking ωΣ=∑d=1pωdp=2ωr$\omega _{\Sigma }=\frac{\sum _{d=1}^{p}\omega _d}{p} = 2 \omega _\text{r}$, we have a configuration of pump tones {−2Δ,0,2Δ}$\lbrace {-2\Delta ,0,2\Delta \rbrace }$ with respect to 2ωr$2 \omega _\text{r}$. We define our spectral modes around ωr$\omega _\text{r}$ as {(−2Δ,−Δ);(−Δ,0);(0,Δ);(Δ,2Δ)}$\lbrace (-2\Delta ,-\Delta );(-\Delta ,0);(0,\Delta );(\Delta ,2\Delta )\rbrace$ =̂{a∼1(ω);a∼2(ω);a∼3(ω);a∼4(ω)}$\widehat{=} \lbrace \tilde{a}_1(\omega ); \tilde{a}_2(\omega ); \tilde{a}_3(\omega ); \tilde{a}_4(\omega ) \rbrace$. Individual pump tones at ω1 and ω3 create TMS states between neighboring spectral modes as in the two pump case above. However, the middle pump creates TMS correlations between two pairs of spectral modes: a∼1(ω)↔a∼4†(−ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}^\dagger _4(-\omega )$ and a∼2(ω)↔a∼3†(−ω)$\tilde{a}_2(\omega ) \leftrightarrow \tilde{a}^\dagger _3(-\omega )$. Rotation of the middle pump phase has an effect on both corresponding subspaces of the covariance matrix. Consequently, this pump configuration is quite suitable for producing square H∼$\mathcal {\tilde{H}}$‐graph states (see Figure 2).A34M=c0000−α0−α0c00−α0−α000c00−α0−α000c−α0−α00−α0−αc000−α0−α00c000−α0−α00c0−α0−α0000c$$\begin{align} \mathbf {M}=\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & 0 & -\alpha & 0 & \pmb {-\alpha }\\ 0 & c & 0 & 0 & -\alpha & 0 & \pmb {-\alpha } & 0\\ 0 & 0 & c & 0 & 0 & \pmb {-\alpha }& 0 & -\alpha \\ 0 & 0 & 0 & c & \pmb {-\alpha }& 0 & -\alpha & 0\\ 0 & -\alpha & 0 & \pmb {-\alpha }& c & 0 & 0 & 0\\ -\alpha & 0 & \pmb {-\alpha } & 0 & 0 & c & 0 & 0\\ 0 & \pmb {-\alpha }& 0 & -\alpha & 0 & 0 & c & 0\\ \pmb {-\alpha }& 0 & -\alpha & 0 & 0 & 0 & 0 & c\\ \end{bmatrix} \end{align}$$The parametric interactions in the covariance matrix can be analyzed in the same way as above, but now the number of phase differences influencing the beamsplitter correlations has increased. The system of linear equations ban be written as for three modes in Equation (A21), but we skip it and write down the interaction matrix M$\mathbf {M}$ (Equation (A34)), where all c coefficient are equal since we have assumed Δr=ωr−ωΣ2=0$\Delta _\text{r}=\omega _\text{r}-\frac{\omega _{\Sigma }}{2}=0$. The signs of α's are governed by the choice of pump phases as {αeiπ2;αeiπ2;αeiπ2}$ \lbrace { \alpha e^{\frac{i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}\rbrace }$. The correlations produced by the pump at ω2=2ωr$\omega _2=2\omega _\text{r}$ are indicated in bold. The special role of the central pump is seen because its correlations cover the whole ascending diagonal.The inverse matrix M−1$\mathbf {M}^{-1}$ reveals the beamsplitter correlations between a∼1(ω)↔a∼3(ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}_3(\omega )$ and a∼2(ω)↔a∼4(ω)$\tilde{a}_2(\omega ) \leftrightarrow \tilde{a}_4(\omega )$:A35M−1=1(c2−4α2)·c−2α2c02α2c00α0α0c−2α2c02α2cα0α02α2c0c−2α2c00α0α02α2c0c−2α2cα0α00α0αc−2α2c02α2c0α0α00c−2α2c02α2c0α0α2α2c0c−2α2c0α0α002α2c0c−2α2c$$\begin{eqnarray} &&\mathbf {M}^{-1}=\frac{1}{(c^2-4\alpha ^2)}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} {{ \def\eqcellsep{&}\begin{array}{cccccccc}c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} & 0 & 0 & \alpha & 0 & \alpha \\[3pt] 0 & c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} & \alpha & 0 & \alpha & 0 \\[3pt] \pmb {\frac{2 \alpha ^2}{c}} & 0 & c-\frac{2 \alpha ^2}{c} & 0 & 0 & \alpha & 0 & \alpha \\[3pt] 0 & \pmb {\frac{2 \alpha ^2}{c}}& 0 & c-\frac{2 \alpha ^2}{c} & \alpha & 0 & \alpha & 0 \\[3pt] 0 & \alpha & 0 & \alpha & c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} & 0 \\[3pt] \alpha & 0 & \alpha & 0 & 0 & c-\frac{2 \alpha ^2}{c} & 0 & \pmb {\frac{2 \alpha ^2}{c}} \\[3pt] 0 & \alpha & 0 & \alpha & \pmb {\frac{2 \alpha ^2}{c}} & 0 & c-\frac{2 \alpha ^2}{c} & 0 \\[3pt] \alpha & 0 & \alpha & 0 & 0 & \pmb {\frac{2 \alpha ^2}{c}}& 0 & c-\frac{2 \alpha ^2}{c} \\[3pt] \end{array} }} \end{bmatrix}\nonumber\\ \end{eqnarray}$$The beamsplitter correlations are indicated in bold in this matrix M−1$\mathbf {M}^{-1}$. We see that there are two sequences of pump transformations that yield BS correlations between modes a∼1(ω)↔a∼2(ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}_2(\omega )$ and a∼3(ω)↔a∼4(ω)$\tilde{a}_3(\omega ) \leftrightarrow \tilde{a}_4(\omega )$. This agrees with the simple argument that indicates BS correlations to exist when two spectral bands are connected across squeezing action by two pumps with one joint frequency. Higher order correlations via three pumps exist also, but these are neglected in our analysis. Note that the number of beamsplitter correlations also coincides with the independent number of phase differences between the pumps. Connection of the cavity spectral mode correlations to H∼$\mathcal {\tilde{H}}$‐graphs is illustrated in Figure 2.The beamsplitter correlations are prominent also in the covariance matrix Va$\mathbf {V}_\text{a}$ (see Equation (A28)):A36Va=κ4(c2−4α2)2·−2α2+8α4c2+c202αc06α2−8α4c202αc00−2α2+8α4c2+c20−2αc06α2−8α4c20−2αc2αc0−2α2+8α4c2+c202αc06α2−8α4c200−2αc0−2α2+8α4c2+c20−2αc06α2−8α4c26α2−8α4c202αc0−2α2+8α4c2+c202αc006α2−8α4c20−2αc0−2α2+8α4c2+c20−2αc2αc06α2−8α4c202αc0−2α2+8α4c2+c200−2αc06α2−8α4c20−2αc0−2α2+8α4c2+c2.$$\begin{eqnarray} &&\mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-4\alpha ^2)^2}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} {{ \def\eqcellsep{&}\begin{array}{cccccccc}-2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} }& 0 & 2 \alpha c & 0 \\[3pt] 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & -2 \alpha c \\[3pt] 2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 \\[3pt] 0 & -2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } \\[3pt] \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & 2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & 2 \alpha c & 0 \\[3pt] 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & -2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 & -2 \alpha c \\[3pt] 2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & 2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 & 0 \\[3pt] 0 & -2 \alpha c & 0 & \pmb { 6 \alpha ^2-\frac{8 \alpha ^4}{c^2} } & 0 & -2 \alpha c & 0 & -2 \alpha ^2+\frac{8 \alpha ^4}{c^2}+c^2 \end{array} }} \end{bmatrix}.\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$Bold font marks the beamsplitter correlations which display a different structure in comparison to Equation (A35) owing to the base change to quadratures ordered as (x∼1,p∼1,⋯x∼N,p∼N)T$(\tilde{x}_1, \tilde{p}_1, \dots \tilde{x}_N, \tilde{p}_N )^{T}$. So the BS correlations are between quadratures of a∼1(ω)↔a∼3(ω)$\tilde{a}_1(\omega ) \leftrightarrow \tilde{a}_3(\omega )$ and a∼2(ω)↔a∼4(ω)$\tilde{a}_2(\omega ) \leftrightarrow \tilde{a}_4(\omega )$.By choosing the phase of the first pump to be opposite to that of the second and the third {αe−iπ2;αeiπ2;αeiπ2}$ \lbrace { \alpha e^{\frac{-i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}; \alpha e^{\frac{i\pi }{2}}\rbrace }$ we are able to flip the sign of one minor diagonal indicated by bold font in Equation (A37).A37M=c0000α0−α0c00α0−α000c00−α0−α000c−α0−α00α0−αc000α0−α00c000−α0−α00c0−α0−α0000c$$\begin{equation} \mathbf {M}=\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & 0 & \pmb {\alpha } & 0 & -\alpha \\ 0 & c & 0 & 0 & \pmb {\alpha } & 0 & -\alpha & 0\\ 0 & 0 & c & 0 & 0 & -\alpha & 0 & -\alpha \\ 0 & 0 & 0 & c & -\alpha & 0 & -\alpha & 0\\ 0 & \pmb {\alpha } & 0 & -\alpha & c & 0 & 0 & 0\\ \pmb {\alpha } & 0 & -\alpha & 0 & 0 & c & 0 & 0\\ 0 & -\alpha & 0 & -\alpha & 0 & 0 & c & 0\\ -\alpha & 0 & -\alpha & 0 & 0 & 0 & 0 & c\\ \end{bmatrix} \end{equation}$$Interestingly, this choice of phases leads to full cancellation of the beamsplitter correlation terms. This is seen in the structure of the inverse matrix:A38M−1=1(c2−2α2)·c0000−α0α0c00−α0α000c00α0α000cα0α00−α0αc000−α0α00c000α0α00c0α0α0000c$$\begin{eqnarray} &&\mathbf {M}^{-1}=\frac{1}{(c^2-2\alpha ^2)}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} c & 0 & 0 & 0 & 0 & -\alpha & 0 & \alpha \\ 0 & c & 0 & 0 & -\alpha & 0 & \alpha & 0\\ 0 & 0 & c & 0 & 0 & \alpha & 0 & \alpha \\ 0 & 0 & 0 & c & \alpha & 0 & \alpha & 0\\ 0 & -\alpha & 0 & \alpha & c & 0 & 0 & 0\\ -\alpha & 0 & \alpha & 0 & 0 & c & 0 & 0\\ 0 & \alpha & 0 & \alpha & 0 & 0 & c & 0\\ \alpha & 0 & \alpha & 0 & 0 & 0 & 0 & c\\ \end{bmatrix} \end{eqnarray}$$which does not have any elements proportional to α2. Without any beamsplitter correlations, Equation (A38) indicates a clear connection to a square‐lattice H∼$\mathcal {\tilde{H}}$‐graph.Finally, the corresponding covariance matrix Va$\mathbf {V}_\text{a}$ for cavity spectral modes is given by Equation (A39). This structure for the covariance matrix is obtained when all the pump signals have an additional phase shift of eiπ2$e^{\frac{i\pi }{2}}$. Such a choice of phases will result in a covariance matrix with “symmetric” structure as shown in experimental data in Figure 8a,b. By controlling the phases of the pump tones separately, we can rotate and adjust certain subspaces of the 8 × 8 covariance matrix. In particular, the influence of the beamsplitter correlations can be eliminated from Va$\mathbf {V}_\text{a}$ in the four pump case.Regarding the quadripartite covariance matrix structures, the relative phase shift between the pump tones are not influenced by the cavity response in the limit of vanishing band widths or with the assumption of huge coupling rate and tiny internal dissipation loss rate. However, additional phase shifts will appear if these conditions are not met, which has to be taken into account in the generation of the desired entangled states.In principle, it would be possible to evaluate the criteria for GME from the analytical expressions derived in this Appendix (see e.g. Equations (A33) and (A39)). However, we leave the conclusions about genuine entanglement both in the tripartite and quadripartite case for analysis based on numerical simulations where even the nonlinear terms can be taken into account. The nonlinear terms are of central importance when increasing the pump drive past the critical pumping amplitude.A39Va=κ4(c2−2α2)2·c2+2α20−2αc0002αc00c2+2α202αc000−2αc−2αc0c2+2α202αc00002αc0c2+2α20−2αc00002αc0c2+2α202αc0000−2αc0c2+2α20−2αc2αc0002αc0c2+2α200−2αc000−2αc0c2+2α2.$$\begin{eqnarray} &&\mathbf {V}_\text{a}= \frac{\kappa }{4(c^2-2\alpha ^2)^2}\cdot \nonumber \\ &&\def\eqcellsep{&}\begin{bmatrix} {{ \def\eqcellsep{&}\begin{array}{cccccccc}c^2+2 \alpha ^2 & 0 & -2 \alpha c & 0 & 0 & 0 & 2 \alpha c & 0 \\[3pt] 0 & c^2+2 \alpha ^2 & 0 & 2 \alpha c & 0 & 0 & 0 & -2 \alpha c \\[3pt] -2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & 2 \alpha c & 0 & 0 & 0 \\[3pt] 0 & 2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & -2 \alpha c & 0 & 0 \\[3pt] 0 & 0 & 2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & 2 \alpha c & 0 \\[3pt] 0 & 0 & 0 & -2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 & -2 \alpha c \\[3pt] 2 \alpha c & 0 & 0 & 0 & 2 \alpha c & 0 & c^2+2 \alpha ^2 & 0 \\[3pt] 0 & -2 \alpha c & 0 & 0 & 0 & -2 \alpha c & 0 & c^2+2 \alpha ^2 \end{array} }} \end{bmatrix}.\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$AcknowledgementsThe authors thank Mark Dykman, Stefano Pirandola, Olivier Pfister, Aidar Sultanov, Shruti Dogra, and Marko Kuzmanovié for fruitful discussions and useful comments. This work was supported by the Academy of Finland through the Finnish Center of Excellence in Quantum Technology QTF (projects 312295, 312296, 336810, 336813). This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 862644 (FET‐Open project QUARTET) and grant agreement no. 824109 (European Microkelvin Platform), as well as ERC grant agreement no. 670743 (QuDeT). G.S.P and K.V.P. also would like to thank Saab for support, under a research agreement between Saab and Aalto University. The work of PJH and IL was supported by MATINE grant VN/13774/2020‐PLM‐49.Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.R. Jozsa, N. Linden, Proc. R. Soc. London, Ser. A 2003, 459, 2011.R. Horodecki, P. Horodecki, M. Horodecki, K. Horodecki, Rev. Mod. Phys. 2009, 81, 865.R. Raussendorf, H. J. Briegel, Phys. Rev. Lett. 2001, 86, 5188.N. C. Menicucci, S. T. Flammia, O. Pfister, Phys. Rev. Lett. 2008, 101, 130501.T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, J. L. O'Brien, Nature 2010, 464, 45.N. C. Menicucci, Phys. Rev. Lett. 2014, 112, 120504.O. Pfister, J. Phys. B: At. Mol. Opt. Phys. 2019, 53, 012001.40 years of quantum computing, Nat. Rev. Phys. 2022, 4, https://doi.org/10.1038/s42254-021-00410-6.V. Giovannetti, S. Lloyd, L. Maccone, Nat. Photonics 2011, 5, 222.Z. Zhang, S. Mouradian, F. N. C. Wong, J. H. Shapiro, Phys. Rev. Lett. 2015, 114, 110506.C. L. Degen, F. Reinhard, P. Cappellaro, Rev. Mod. Phys. 2017, 89, 035002.L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, P. Treutlein, Rev. Mod. Phys. 2018, 90, 035005.X. Su, J. Jing, Q. Pan, C. Xie, Phys. Rev. A 2006, 74, 062305.N. Gisin, R. Thew, Nat. Photonics 2007, 1, 165.J. Yin, Y.‐H. Li, S.‐K. Liao, M. Yang, Y. Cao, L. Zhang, J.‐G. Ren, W.‐Q. Cai, W.‐Y. Liu, S.‐L. Li, R. Shu, Y.‐M. Huang, L. Deng, L. Li, Q. Zhang, N.‐L. Liu, Y.‐A. Chen, C.‐Y. Lu, X.‐B. Wang, F. Xu, J.‐Y. Wang, Ch.‐Z. Peng, A. K. Ekert, J.‐W. Pan, Nature 2020, 582, 501.S. J. Summers, R. Werner, Phys. Lett. A 1985, 110, 257.S. L. Braunstein, P. van Loock, Rev. Mod. Phys. 2005, 77, 513.P. Lähteenmäki, G. S. Paraoanu, J. Hassel, P. J. Hakonen, Nat. Commun. 2016, 7, 12548.A. M. Lance, T. Symul, W. P. Bowen, B. C. Sanders, T. Tyc, T. C. Ralph, P. K. Lam, Phys. Rev. A 2005, 71, 033814.M. Hillery, V. Bužek, A. Berthiaume, Phys. Rev. A 1999, 59, 1829.S. L. Braunstein, H. J. Kimble, Phys. Rev. A 2000, 61, 042302.J. Dias, T. C. Ralph, Phys. Rev. A 2018, 97, 032335.G. J. Milburn, S. L. Braunstein, Phys. Rev. A 1999, 60, 937.X. Guo, C. R. Breum, J. Borregaard, S. Izumi, M. V. Larsen, T. Gehring, M. Christandl, J. S. Neergaard‐Nielsen, U. L. Andersen, Nat. Phys. 2020, 16, 281.R. Nichols, P. Liuzzo‐Scorpo, P. A. Knott, G. Adesso, Phys. Rev. A 2018, 98, 012114.M. Gessner, A. Smerzi, L. Pezzè, Nat. Commun. 2020, 11, 3817.J. Zhang, S. L. Braunstein, Phys. Rev. A 2006, 73, 032318.N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, M. A. Nielsen, Phys. Rev. Lett. 2006, 97, 110501.M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, P. van Loock, Phys. Rev. A 2009, 79, 062318.A. Mantri, T. F. Demarie, J. F. Fitzsimons, Sci. Rep. 2017, 7, 42861.P. Walther, K. J. Resch, T. Rudolph, E. Schenck, H. Weinfurter, V. Vedral, M. Aspelmeyer, A. Zeilinger, Nature 2005, 434, 169.S. Zippilli, D. Vitali, Phys. Rev. A 2020, 102, 052424.S. Zippilli, D. Vitali, Phys. Rev. Lett. 2021, 126, 020402.C. W. Sandbo Chang, M. Simoen, J. Aumentado, C. Sabín, P. Forn‐Díaz, A. M. Vadiraj, F. Quijandría, G. Johansson, I. Fuentes, C. M. Wilson, Phys. Rev. Appl. 2018, 10, 044019.G. Adesso, F. Illuminati, J. Phys. A: Math. Theor. 2007, 40, 7821.C. Weedbrook, S. Pirandola, R. García‐Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, S. Lloyd, Rev. Mod. Phys. 2012, 84, 621.C. Fabre, N. Treps, Rev. Mod. Phys. 2020, 92, 035005.N. C. Menicucci, S. T. Flammia, P. van Loock, Phys. Rev. A 2011, 83, 042335.D. E. Bruschi, C. Sabín, P. Kok, G. Johansson, P. Delsing, I. Fuentes, Sci. Rep. 2016, 6, 18349.T. Yamamoto, K. Koshino, Y. Nakamura, in Parametric Amplifier and Oscillator Based on Josephson Junction Circuitry, Springer, Japan, Tokyo 2016, pp. 495–513.Z. Wang, M. Pechal, E. A. Wollack, P. Arrangoiz‐Arriola, M. Gao, N. R. Lee, A. H. Safavi‐Naeini, Phys. Rev. X 2019, 9, 021049.T. Korkalainen, I. Lilja, M. R. Perelshtein, K. V. Petrovnin, G. S. Paraoanu, P. J. Hakonen, AIP Conf. Proc. 2021, 2362, 030001.P. Lähteenmäki, G. Paraoanu, J. Hassel, P. J. Hakonen, PNAS 2013, 110, 4234.H. J. Briegel, R. Raussendorf, Phys. Rev. Lett. 2001, 86, 910.D. E. Bruschi, C. Sabín, G. S. Paraoanu, Phys. Rev. A 2017, 95, 062324.G. Adesso, S. Ragy, A. R. Lee, Open Syst. Inf. Dyn. 2014, 21, 1440001.P. van Loock, A. Furusawa, Phys. Rev. A 2003, 67, 052315.M. Hillery, H. T. Dung, H. Zheng, Phys. Rev. A 2010, 81, 062322.A. Agustí, C. W. S. Chang, F. Quijandría, G. Johansson, C. M. Wilson, C. Sabín, Phys. Rev. Lett. 2020, 125, 020502.R. Simon, Phys. Rev. Lett. 2000, 84, 2726.L.‐M. Duan, G. Giedke, J. I. Cirac, P. Zoller, Phys. Rev. Lett. 2000, 84, 2722.S. Mancini, V. Giovannetti, D. Vitali, P. Tombesi, Phys. Rev. Lett. 2002, 88, 120401.P. Hyllus, J. Eisert, New J. Phys. 2006, 8, 51.R. Y. Teh, M. D. Reid, Phys. Rev. A 2014, 90, 062337.E. Shchukin, P. van Loock, Phys. Rev. A 2015, 92, 042328.A. Peres, Phys. Rev. Lett. 1996, 77, 1413.M. Horodecki, P. Horodecki, R. Horodecki, Phys. Lett. A 1996, 223, 1.L. K. Shalm, D. R. Hamel, Z. Yan, C. Simon, K. J. Resch, T. Jennewein, Nat. Phys. 2013, 9, 19.T. Elo, T. S. Abhilash, M. R. Perelshtein, I. Lilja, E. V. Korostylev, P. J. Hakonen, Appl. Phys. Lett. 2019, 114, 152601.C. Macklin, K. O'Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, I. Siddiqi, Science 2015, 350, 307.A. Zorin, Phys. Rev. Appl. 2019, 12, 044051.M. Perelshtein, K. Petrovnin, V. Vesterinen, S. Hamedani Raja, I. Lilja, M. Will, A. Savin, S. Simbierowicz, R. Jabdaraghi, J. Lehtinen, L. Grönberg, J. Hassel, M. Prunnila, J. Govenius, G. Paraoanu, P. Hakonen, Phys. Rev. Appl. 2022, 18, 024063.M. Esposito, A. Ranadive, L. Planat, S. Leger, D. Fraudet, V. Jouanny, O. Buisson, W. Guichard, C. Naud, J. Aumentado, F. Lecocq, N. Roch, Phys. Rev. Lett. 2022, 128, 153603.J. Y. Qiu, A. Grimsmo, K. Peng, B. Kannan, B. Lienhard, Y. Sung, P. Krantz, V. Bolkhovsky, G. Calusine, D. Kim, A. Melville, B. M. Niedzielski, J. Yoder, M. E. Schwartz, T. P. Orlando, I. Siddiqi, S. Gustavsson, K. P. O'Brien, W. D. Oliver, arXiv:2201.11261, 2022.T. Roy, S. Kundu, M. Chand, A. M. Vadiraj, A. Ranadive, N. Nehra, M. P. Patankar, J. Aumentado, A. A. Clerk, R. Vijay, Appl. Phys. Lett. 2015, 107, 262601.R. N. Alexander, N. C. Menicucci, Phys. Rev. A 2016, 93, 062326.R. Pooser, J. Jing, Phys. Rev. A 2014, 90, 043841.M. V. Larsen, J. S. Neergaard‐Nielsen, U. L. Andersen, Phys. Rev. A 2020, 102, 042608.J. Wu, Q. Zhuang, Phys. Rev. Appl. 2021, 15, 034073.B.‐H. Wu, R. N. Alexander, S. Liu, Z. Zhang, Phys. Rev. Res. 2020, 2, 023138.T. Hillmann, F. Quijandría, G. Johansson, A. Ferraro, S. Gasparinetti, G. Ferrini, Phys. Rev. Lett. 2020, 125, 160501.Y. Y. Gao, B. J. Lester, K. S. Chou, L. Frunzio, M. H. Devoret, L. Jiang, S. Girvin, R. J. Schoelkopf, Nature 2019, 566, 509.S. Lloyd, Science 2008, 321, 1463.R. Simon, N. Mukunda, B. Dutta, Phys. Rev. A 1994, 49, 1567.J. Williamson, Am. J. Math. 1936, 58, 141.

Journal

Advanced Quantum Technologies – Wiley

Published: Jan 1, 2023

Keywords: continuous variables; Josephson parametric amplifier; microwaves; multipartite entanglement; quantum correlations; quantum entanglement

References

  • Parametric Amplifier and Oscillator Based on Josephson Junction Circuitry
    Yamamoto, T.; Koshino, K.; Nakamura, Y.