Generation and Structuring of Multipartite Entanglement in a Josephson Parametric System

Quantum correlations are a vital resource in advanced information processing based on quantum phenomena. Remarkably, the vacuum state of a quantum field may act as a key element for the generation of multipartite quantum entanglement. In this work, generation of genuine tripartite entangled state and its control is achieved by the use of the phase difference between two continuous pump tones. Control of the subspaces of the covariance matrix for tripartite bisqueezed state is demonstrated. Furthermore, by optimizing the phase relationships in a three‐tone pumping scheme genuine quadripartite entanglement of a generalized H‐graph state ( H∼$\mathcal {\tilde{H}}$ ‐graph) is explored. This scheme provides a comprehensive control toolbox for the entanglement structure and allows to demonstrate, for first time to the authors' knowledge, genuine quadripartite entanglement of microwave modes. All experimental results are verified with numerical simulations of the nonlinear quantum Langevin equation. It is envisioned that quantum resources facilitated by multi‐pump configurations offer enhanced prospects for quantum data processing using parametric microwave cavities.

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][28][29], which has been under active development for the past 20 years.Cluster state calculus, foremost utilizing optical resources [28,[30][31][32][33], realize measurementbased 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 steadystate conditions.Using the Gaussian-mode formalism [35][36][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 II 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 III we present the experimental setup and explain our data analysis methods.In Section IV we present the experimental and theoretical results on the generation of multipartite entanglement using a Josephson paramet-arXiv:2203.09247v2[quant-ph] 18 Apr 2023 ric system in both tripartite and quadripartite cases.In Appendices VIIA-E we present details of our analysis, experimental techniques, and parameter extraction.Appendix VIIF 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.

II. THEORETICAL FOUNDATIONS
A. Quantum dynamics of the device Our work employs a parametric system comprised of a superconducting λ/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 multitone external RF magnetic flux through the SQUID at frequencies that are approximately twice the frequency of the resonator ω d ∼ 2ω r (three-wave mixing) [18,42].In the rotating frame, the Hamiltonian of the system, as derived in Refs.40 and 43, is given by: where ã (ã † ) is the annihilation (creation) operator for cavity photons in the rotating frame at angular frequency ω Σ /2, α d is the complex amplitude of the d-th pump tone and ∆ d = ω d − ω Σ 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 = |α|e iϕ d .Here ∆ r denotes the detuning between half of the average pump angular frequency ω Σ and the resonator angular frequency ω r : ∆ r = ω r −ω Σ /2, with ω Σ representing the average angular frequency in multi-tone driven case: ω Σ = (1/p) p d=1 ω d with d = {1, . . ., p} 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 .Several effects are accounted for by the Kerr nonlinearity, such as limited maximum gain, compression, observed at α α 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 VII D).
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: where H sig includes the coupling to the signal port transmission line with dissipation rate κ, while H 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 yields where bin and cin are the ladder (annihilation) operators for the signal and linear dissipation ports, respectively.The output mode bout is obtained using the following relation between the incoming and outgoing modes: We are interested in the correlations embedded in the output mode given by Eq. ( 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 II B) in the frequency domain and examining correlations between these spectral ranges, we can verify the presence of entanglement in the band-limited microwave signals.

B. Spectral modes definition
Parametric downconversion processes and the definition of employed spectral modes within the fundamental cavity resonance in a multi-pump JPA are illustrated in Fig. 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).

JPA resonance
ω1/2 FIG. 1. Definition 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 frequency with a bandwidth 2δω, within which ([ω r − δω, ω r + δω]) we define N spectral modes as depicted in Fig. 1.Let us define ã as a vector of spectral modes: where N is a total mode number and the creation ã † i = ã † i (t) and annihilation ãi = ãi (t) operators are timedependent.In general, the frequency difference between half of p-th and (p + 1)-th pump tones defines the bandwidth of the spectral mode ãi .We employ an equidistant pump scheme where bandwidth of each mode ãi is defined as ∆ such that the spectrum of a full set of modes ãi covers the bandwidth 2δω of the cavity mode ã.In the experiment, we collect the emitted power over the whole [ω r − δω, ω r + δω] 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, p = 2) and quadripartite (N = 4, 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.

C. Graphical description of quantum states
In order to construct a comprehensive graphical representation, we consider interactions between cavity and input vacuum modes by solving the QLE given in Eq. (3) for N cavity modes defined in Eq. ( 5); for details see Appendix VII F. Assuming the strong coupling regime with κ ∆ while neglecting dissipation losses γ and the Kerrnonlinearity K, Fourier transformation yields a system of linear equations which can be cast in a matrix form: Here the interaction matrix 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 and having different phases for the pump tones α 1 = αe iϕ1 , α 2 = αe iϕ2 : where the frequency dependency enters through the c = −iω + κ/2 coefficient.To express the intracavity modes through the vacuum input, we use the inverse matrix For the tripartite case, such a matrix is written in the following way where ∆ϕ = ϕ 1 −ϕ 2 .Besides two mode squeezing (TMS) correlations proportional to α, the matrix contains also beamsplitter correlations (BS) ∝ α 2 (denoted in bold).Note that the α 2 -terms are absent in the matrix M introduced in Eq. ( 6).In this tripartite example, the phase difference ∆ϕ contributes only to the BS connection ãi (ω) ↔ ãj (ω) or ã † i (ω) ↔ ã † j (ω).TMS connections are defined by entries ãi (ω) ↔ ã † j (ω) in M −1 .In Appendix VII F 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 responsible for beamsplitter interaction can be suppressed fully by choosing pump phases properly in certain pump tone configurations.
Interestingly, the matrix √ κ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-graphs formed by both two mode squeezing and beamsplitter correlations (see Appendix VII F).Examples of such graphs are presented in Fig. 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 Fig. 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.
From the experimental point of view, the observer is interested in the adjacency matrix for out-coming modes bout , which are obtained from ã using the input-output relationship in Eq. ( 4).This equation yields bout (ω) = (I − κM −1 ) bin (ω), (10) on the basis of which we may define the adjacency matrix M = I − κ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 is sufficient to characterize the BS and TMS connections between vacuum spectral modes.

D. Connection to Hamiltonian graph
CV cluster states with square-lattice graph structure provide 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, whose entries G ij specify the multimode squeezing Hamiltonian as: Here, the pump tone amplitudes are considered to have equal strength α.The matrix G involves TMS correlations between modes ãi ↔ ã † 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 i ȧ † k = α j G jk ãj and i ȧk = −α j G * jk ã † j .Taking the Fourier transform, the left hand side equals ω × ã(ω), and the combination ωã k (ω) = α j G jk ã † j (−ω) provides the connection to the QLE treatment in Eq. ( 6): ã(ω) here is the cavity signal defined by the graph connections given by G jk .Consequently, the basic graph structures are the same, but the form of 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 , 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 + i 2 G 2 + 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 Eq. 11 in R = e − i Hτ 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 B ij = ã † i ã † j + ãi ãj , and BS transformations based on T ij = ãi ã † j − ã † i ãj .By utilizing the Zassenhaus expansion (up to first order of commutation relationship) [45], we obtain for the tripartite squeezing operator 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 as e θ13T13 e θ24T24 . . ., (13) where the total number of TMS and BS operators in the decomposition is N BS = n N − 2n and N TMS = n N − 2n + 1; where n = 1, 2, . . ., N   2   for the configuration introduced in Fig. 1 with N − 1 pump tones and N spectral modes.Collecting all of the entries of B ij and T ij , we obtain a generalized adjacency matrix G with entries Gij = NTMS 1 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 defines the edge connections in the generalized H-graph.

E. Verification of the multipartite entanglement
The generalized graph analysis allows us to visualise 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 , p1 , . . .xN , pN ) T , where xi = (ã i + ã † i )/2 and pi = (ã i − ã † i )/2i.The covariance matrix V, whose elements are given by 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 = ri − ri .
Obtaining the covariance matrix, we can analyse 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][51][52][53][54][55], we apply symplectic transformations to the covariance matrix and calculate its symplectic eigenvalues -the 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 } N i=1 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, where is sufficient to confirm genuine tripartite entanglement (N = 3) and violation of where is sufficient to confirm genuine quadripartite entanglement (N = 4) with weights h i , g k being in range [−1, 1].
To simplify the search domain, we set h 1 = g 1 = 1 and In the double and triple pumping schemes, we generate generalized tripartite and quadripartite 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-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 h i = g i = 1) in order to detect the minimum value for S.

III. EXPERIMENT A. Experimental methods
In our microwave experiments, we employ a niobium λ/4 coplanar 50 Ω transmission line terminated into a SQUID loop (QWJPA), forming a quarter-wave Josephson parametric amplifier.The SQUID's junctions are formed by Nb/Al-Al 2 O 3 /Nb 1×1 µm tunnel barriers with I c 4 µA critical current.The JPA, operating in threewave-mixing mode around the cavity frequency ω r , is pumped by an external RF magnetic flux through the SQUID loop using a single turn pump coil at frequency 2ω r (marked as 2ω LO in Fig. 3a) [40,59].We chose this operation regime, because for four-wave mixing (typically using a current pump near ω r ), the large amplitude pump is within the amplification bandwidth, whereas the threewave 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 FIG. 3. Experimental 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) and (c), respectively.
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.
Our measurement setup is illustrated in Fig. 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 Fig. 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 VII D.
a. Tripartite case.In the tripartite case, phasecontrolled 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 and ω 2 = ω r + ∆/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 through the SQUID loop, resulting in ωr 2π = 6.024GHz for the cavity frequency.In the three-mode experiment, we apply two pump tones at (2 × ωr 2π − 2) MHz and (2 × ωr 2π + 2) MHz, with the half-frequencies positioned as depicted in Fig. 3b and the correlated spectral modes are defined symmetrically with respect to the pump halffrequencies.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 per channel covering the bandwidth of 25 MHz.Furthermore, we employ an overall detuning of 14 MHz, i.e. a heterodyne detection scheme, in order to avoid 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.
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 IV.
b. 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, e.g. 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 through a SQUID resulting in ω r = 5.978 GHz cavity frequency that slightly differs from the tripartite case, see Fig. 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× ωr 2π MHz, (2 × ωr 2π + 1) MHz, and (2 × ωr 2π − 1) MHz, we generate 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-graph states and their genuine multipartite entanglement depends on relations among the covariance matrix elements as discussed in Section II E. The degree of violation in the GME con-dition S < 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 VII C.

B. Scaled covariance matrix
The system gain G Σ determined in Appendix VII C refers to measured power per unit band width.Since the measured spectral mode quadratures I i and Q i are determined over the band ∆f i , the scaled quadrature x i , equivalent to the amplitude of the quantum mechanical operators x, is given by the formula where G i = G Σ,i is the system gain for i th spectral mode, Z 0 = 50 Ω is transmission line impedance and ∆f i is the bandwidth of the spectral mode: ∆f i = 2 MHz or ∆f i = 0.4 MHz for tripartite and quadripartite case, respectively.Similar scaling is applied also to the quadrature component 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 Eq. 14): where V off denotes the covariance matrix measured in the absence of the pump.Due to scaling of the covariance matrix V, this equation yields a unity diagonal matrix in the absence of pumping at T → 0. The average physical temperature in our experiments is T i = 20 mK resulting in coth hfi 2k b Ti = 1.000.

IV. MULTIPARTITE ENTANGLEMENT
To characterize the structure of the entanglement in output states, we analyze the resulting covariance matrices using positive partial transpose (PPT) formalism and GME criteria discussed in Section II E for tripartite and quadripartite cases.a. Tripartite case.Leveraging 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 Eq. 26 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.
Fig. 4 depicts our experimental results on genuine tripartite entanglement and their comparison with simulations.Fig. 4a illustrates results of numerical simulations on GME in terms of S defined in Eq. 15.At weak pumping, the condition for genuine entanglement S < 1 is fulfilled almost independent of the pump phases, but with increasing A, the simulations reveal an even smaller range of ∆ϕ yielding S < 1 (see the inset in Fig. 4a).The strongest genuine tripartite entanglement is reached at ∆ϕ −120 • under normalized pumping amplitude A 0.22, at which the simulations reach S = 0.70.At the minimum of S, the corresponding weights are h i = {1, −0.65, −0.65} and g i = {1, 0.65, 0.65}.It is noteworthy that the phase setting ∆ϕ = +90 • yields clearly worse entanglement than ∆ϕ = −90 • .This asymmetry in GME between ∆ϕ = ±90 • arises from differences in the covariance matrices which is illustrated in Fig. 5.
Our experimental data on S in Fig. 4b displays similar features as Fig. 4a.The measured GME criterion S as function of normalized pump amplitude for three phase differences is shown in Fig. 4b.In the experimental data, nearly no GME is observed at positive phase differences, whereas ∆ϕ = −90 • and ∆ϕ = −120 • yield suppression down to S = 0.75 ± 0.05.The measured result at ∆ϕ = −120 • 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].Overall, the pattern of S(ϕ, A) in the inset of Fig. 4b coincides with the simulated pattern in Fig. 4a.The agreement strongly supports the presence of genuinely entangled bisqueezed state in our experiment.We emphasize that the optimum entanglement at ∆ϕ = −120 • , 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 VII E), which, when included in the simulations, result in very good matching with the experiment.
Covariance matrices measured at the pump phase difference ∆ϕ = +90 • and ∆ϕ = −90 • are illustrated in Fig. 5. Technically, by rotating the phase of a pump signal, we selectively control certain subspace of the covariance matrix, which can be seen in Fig. 5.An applied phase shift to the 1st or the 2nd 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 Figs.5a and 5b, reveals how the subspaces transform with the phase difference from ∆ϕ = −90 • to ∆ϕ = +90 • .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 re- The inseparability of the covariance matrix was investigated using the PPT criteria (see Sect.II).Symplectic eigenvalues obtained for mode partitions 1 − 23, 2 − 13, and 3 − 12 are depicted in Fig. 6a for simulations, while experimental data are displayed in Fig. 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 Fig. 6b, the middle frequency acted on by both pumps is the most inseparable part of the covariance matrix.
b. Quadripartite case For the quadripartite case, we apply three pump drives with identical amplitude α 1 = α 2 = α 3 = α.While our goal is to demonstrate genuine entanglement generation of cluster states (mode structure depicted on Fig. 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 Sect.IV 0 a).For the analysis of full inseparability of the covariance matrix according to the PPT criterion, we evaluate the minimum symplectic eigenvalues min{ν i } 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 Fig. 7 alongside with the corresponding predictions given by our numerical simulations.The minimum symplectic eigenvalue min{ν i } = 0.79 ± 0.018 is reached, while all of the eigenvalues in the normalized pump amplitude range 0.01 A 0.15 are less than 1.Compared with the minimum symplectic eigenvalues in Fig. 6, we may conclude that the influence of BS correlations on min{ν i } value is less in the quadripartite state than for the tripartite case.
The GME criterion for four modes as a function of the normalized pump amplitude is depicted in Fig. 7b; the symbols display data while the simulation result is indicated by the solid curve.As was discussed in Sect.II, the optimized weights in GME inequality Eq. ( 16) are chosen in the same manner as in the tripartite case: h 1 = g 1 = 1 and h i = h, g i = g, i = {2, 3, 4}.The strongest genuine entanglement S = 0.84 ± 0.02 is observed at A 0.08 pump amplitude using the weights h i = {1, −0.51, −0.51, −0.51} and g i = {1, 0.69, 0.69, 0.69}.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 analysed in the numerical model.
The covariance matrices obtained in the experiment and using numerical simulation are presented in Figs.8a  and 8b, respectively.They are determined at the strongest entanglement point reached at A 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 Fig. 8  Covariance matrix of genuinely entangled quadripartite H-graph state.a) Experimental covariance where the rotation of the TMS subspaces 1-2, 2-3, amd 3-4 have been made in such a way that the structure coincides with the matrix in Eq. ( 54) of Appendix VII F (each pump has phase π 2 ).The employed pump amplitude A 0.08 yields the smallest value for S. b) Simulated covariance matrix using equal pump phases π 2 at A 0.08.The difference from the matrix in Eq. ( 54) is due to the cavity response that induces extra phase shifts.
we conclude that for the employed pump configuration, the genuine quadripartite entanglement appears in the amplitude range 0.01 A < 0.13.

V. DISCUSSION
The control of bisqueezed tripartite and generalized H-graph ( 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-graph states.Enhanced number of modes requires larger bandwidth, which calls for broadband parametric devices such as TWPAs [60][61][62][63][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-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 VII F (Eqs. ( 55)-( 57)) 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-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][28][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, and 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 have utilized this correspondence between optics and microwaves and demonstrated 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-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-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 errorcorrectable 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 [62][63][64].This opens a way to significant increase in the number of entangled modes.

VI. CONCLUSION
In 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 < A < 0.25) and verified genuine entanglement in the range of 0.01 A < 0.4 (0.01 A < 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 • minimum value of criterion S = 0.75 ± 0.05 was obtained.This result were 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-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 VII F, BS correlations can be fully suppressed by implementing a 180 o 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 }, 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-ofthe-art superconducting parametric devices, such as the broadband, low-loss travelling wave parametric amplifier [62][63][64], we expect a substantial increase in the number of entangled modes, which facilitates generation of highlysqueezed square-lattice H-graph states for CV quantum computation at microwave frequencies.
The Hamiltonian of JPA system is given by where â(â † ) is the annihilation (creation) operators for cavity photons, α 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 , d = {1, . . ., p}, we define the detuning between the half pump frequency and the resonator frequency . For each of the p pump tones, we define the detuning from the average frequency ∆ d = ω d −ω Σ , d = {1, . . ., p}.By applying the rotating wave approximation in the frame ω Σ /2 (ã(t) = â(t)e iωΣt/2 ) and leaving only the effective high-order terms, we obtain for the nonlinear part of the Hamiltonian As usual, the bosonic commutation relationships are valid for the cavity modes ã, ã † = 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 by where creation and annihilation operators b † and b refer to modes in the transmission line, and κ denotes the coupling rate.The Hamiltonian related to the linear dissipation port where c † and 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 relations The total Hamiltonian can be conveniently written as a sum of the separate parts given above: For further analysis and for our simulations, we use the Quantum Langevin Equation (QLE) for the cavity operator ã(t): 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 relationship: Eqs. ( 26) and ( 27) are used in our numerical simulations with Matlab ODE45 solver.

B. Full inseparability
Assuming 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 seconds 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 seconds.
The quantum quadratures xi = ãi+ã † i 2 and pi = ãi−ã † i 2i can be combined into a 2N-long column vector operator for the N-mode state r = (x 1 , p1 , . . .xN , pN ) T .The commutation relations can be written down in a skewsymmetric, block-diagonal matrix form [50]: The covariance matrix V is given by elements V ij = 1 2 ∆r i ∆r j + ∆r j ∆r i − ∆r i ∆r j where we have defined standard error ∆r i = ri − ri and ri = tr (r i ρ).
The uncertainty principle requires that applies for a physical covariance matrix.For verification of entanglement, we may investigate a modified equation where V k = λ λ λ k Vλ λ λ k , λ λ λ k is diagonal matrix with ones entries, except of that related to k-th mode, with value of −1.For example, transformation with λ λ λ k≡N = diag(1, 1, ..., 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 Eq. ( 30) when applying a partial transposition with respect to each from full set of modes: V k ≥ i 4 Ω.In Ref. 47, the entangled states are classified in accordance to the number of modes for which the condition (30) 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, e.g.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, i.e., Symplectic transformations on a 2N-dimensional phasespace form the real symplectic group denoted as Sp(2N ; 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: where Ṽk is a 2N-dimensional diagonal matrix consisting of the symplectic eigenvalues, ν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 fulfilled for Ṽk .
For convenience, we insert an additional factor of 4 to the covariance matrix and work with fluctuations with zero mean values: V * ij = 2 ∆r i ∆r j + ∆r j ∆r i .Consequently, for evidence of 'fully inseparable' states, we need to find minimum symplectic eigenvalues with ν * k < 1 for each partial transposition k.

C. System gain calibration
Our 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 T s , the contribution of the cooled amplifier T HEMT , and the noise of the room-temperature amplifiers T RT multi- plied by the system gain G Σ,i = G HEMTi G RTi : Here i refers to the frequency of the spectral mode and ∆f i refers to the bandwidth of the detection of quadratures I i and Q i .The total gain G Σ,i was separately determined for different spectral modes.Fig. 9 displays the measured noise power per unit band as a function of sample temperature T s , averaged over frequencies covering the resonance curve.By fitting a line to the data, we obtain G Σ,i = 94.4 ± 0.2 dB for the average total gain.The linear fit in Fig. 9 is performed at T > 0.2 K, which allows us to neglect the corrections from the coth( ω/2k b T s ).
The error in the system gain calibration results in uncertainty in the symplectic eigenvalues on the order of 2%, i.e. the eigenvalues fall in the range of min{ν * k }=min{ν * k } ± 0.018 for each partial bipartition.Random variations of the system parameters were reduced by averaging the outcome by ten to twenty times.

D. System parameter fitting
In order to determine coupling rates γ and κ introduced in Section IIA, 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 depicted in Fig. 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 .Pumping carried out in degenerate mode, ω d = 2ω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 κ.
of the QLE ( bout (ω)/ bin (ω)), derived for the linear case without any pump drive, we obtain the coupling coefficients κ 2π = 4.44 MHz and γ 2π = 2.30 MHz.The employed analytical solution, displayed in Eq. ( 34), was derived from the full QLE in Eq. ( 26) without taking the nonlinear part −iKã † ãã into account: For fitting of the Kerr constant K, we employed the whole form of the QLE in the rotating wave approximation (26).By comparing the measured and simulated gain coefficients G(ω probe − ω r , A) (Fig. 10) in the cavity at large pump amplitudes, we obtain an estimate K = 6.5ω r for the Kerr constant.Corresponding phase shifts applied to pump tones to reach "symmetric" covariance matrix view on double frequencies are ∆ϕ1 = π 2 − π 4 = π 4 and ∆ϕ2 = π 2 + π 4 = 3π 4 with corresponding phase shift between pump tones π 2 , given in results (See.Fig. 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.

E. Cavity phase response
Optimal value of GME criterion, which governed by "symmetric" covariance matrix view in tripartite mode case, can be obtained with { π 2 ; π 2 } 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 Fig. 11.

F. Multifrequency correlations in terms of QLE with 3 and 4 spectral modes
As discussed in Section IIC, 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 ã as a vector of spectral modes: where the creation ã † i = ã † i (t) and annihilation ãi = ãi (t) operators are time-dependent.After Fourier transform, 2 )} according pump tone positions {−∆, ∆} (see Fig. 1).Similarly, we define for the input and output modes bin/out : The commutation relationships for the case of N modes can be conveniently expressed in matrix form.We use the common convention for [ã i , ãj ] from Ref. 74.
The effect of Kerr nonlinearity is significant only at large pump amplitudes.Hence, we may take the QLE (26) without the nonlinear part for our treatment.In theoretical analysis we assume, that spectral modes lay down deep in cavity mode, such that ∆ κ; 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 ãi (ω) = ãi (t)e iωt dt 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 ã †  (38) the QLE yields the following system of linear equations: We cast Eq. ( 39) into matrix form: where Solving for the inverse of M and using Eq. ( 27), we obtain bout (ω) = (I − κM −1 ) bin (ω) (42) for the outgoing radiation bout (ω) in terms of incoming waves bin (ω).
Because our goal is to determine the structure of the experimental covariance matrix, it is unsatisfactory to consider cavity modes ã with equation ã(ω) = √ κM −1b in (ω) though it has a more compact final form.However, the presence of the identity matrix I and the multiplication factor κ do not change the final structure.
Assuming that the pump amplitude α is a real number and c 1 = c 2 = c (zero detuning case), we have This allows us to draw the generalized H-graph for describing the parametric interaction between the spectral modes, Fig. 2. The off-diagonal beamsplitter elements proportional to α 2 are set in bold in Eq. (43).Still, we want to construct the parametric interaction matrix S −1 for quadrature vector operator r.Using a linear operator matrix K to implement a change of basis we obtain by a canonical transformation 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, bin has a Gaussian covariance matrix of the form V in = 1 4 I. Consequently, the covariance matrix of the cavity spectral modes ãi can be represented as or, equivalently, for output modes bout : Both forms V a and V out can be employed for studying the structure of parametric interactions between the quadratures, because inputoutput relationship doesn't change the general structure of the couplings between the quadratures (see below).
As shown in Section IV experimentally, phase shift between pumps changes the appearance of the covariance matrix (see Fig. 5) as well as the strength of genuine multipartite entanglement.A change in the matrix M due to a phase shift is illustrated in Eq. ( 47), in which the phase of the first pump has been rotated by e The elements affected by the rotation are indicated in bold in the matrix.The elements in bold face indicate coupling between modes ã1 (ω) ↔ ã2 (ω) while the other off-diagonal elements indicate squeezing across ã2 (ω) ↔ ã3 (ω).Note, that phase rotation operates in opposite direction on rows related to bin (ω) and b † in (ω).The inversion of the rotated matrix M yields for the parametric interaction matrix, where all the beamsplitter elements (in bold) have acquired a π/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 in Eq. (48) 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 Eq. ( 15)) depends on the difference of pump phases and strong genuine entanglement is reached only at specific phase settings.
The covariance matrix V a obtained from matrix M in Eq. ( 43) with zero pump phase shifts is given in Eq. ( 49).The corresponding covariance matrix for π/2 phase rotation in the first pump is displayed in Eq. (50).The matrix V a in Eq. ( 50) 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 e iπ 2 to the second pump.This brings the covariance matrix for the spectral cavity modes to the "standard-symmetric" form displayed in Eq. (51).By comparing Eq. ( 49) and Eq. ( 51) 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.
The parametric interactions in the covariance matrix can be analysed 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 Eq. ( 39), but we skip it and write down the interaction matrix M (Eq.( 52)), where all c coefficient are equal since we have assumed ∆ r = ω r − ωΣ 2 = 0.The signs of α's are governed by the choice of pump phases as {αe iπ 2 ; αe iπ 2 ; αe iπ 2 }.The correlations produced by the pump at ω 2 = 2ω 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 reveals the beamsplitter correlations between ã1 (ω) ↔ ã3 (ω) and ã2 (ω) ↔ ã4 (ω): The beamsplitter correlations are indicated in bold in this matrix M −1 .We see that there are two sequences of pump transformations that yield BS correlations between modes ã1 (ω) ↔ ã2 (ω) and ã3 (ω) ↔ ã4 (ω).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-graphs is illustrated in Fig. 2. The beamsplitter correlations are prominent also in the covariance matrix V a (see 46): Bold font marks the beamsplitter correlations which display a different structure in comparison to Eq. ( 53) owing to the base change to quadratures ordered as (x 1 , p1 , . . .xN , pN ) T .So the BS correlations are between quadratures of ã1 (ω) ↔ ã3 (ω) and ã2 (ω) ↔ ã4 (ω).
Finally, the corresponding covariance matrix V a for cavity spectral modes is given by Eq. ( 57).This structure for the covariance matrix is obtained when all the pump signals have an additional phase shift of e iπ 2 .Such a choice of phases will result in a covariance matrix with "symmetric" structure as shown in experimental data in Figs.8a and 8b.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 V 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.Eqs. ( 51) and ( 57)).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.

FIG. 2 .
FIG. 2. Graph representation for the entangled tripartite and quadripartite systems.Generalized H-graph ( H-graph) representation of bisqueezed tripartite CV entangled state (a) and quadripartite CV entangled H-graph state (b) 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) and (d) represent tri-and quadripartite GHZ states, which can be obtained via introducing an additional pump tone with ∆ d = 0 in the tripartite case and two additional tones at −∆ and ∆ in the quadripartite setting.The additional pumps supply missing TMS connections to the entangled states.For details, see Appendix VII F.

FIG. 4 .IG. 5 .
FIG. 4. Phase-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 ∆ϕ between pump signals indicated in the figure; the simulation parameters were set to match the experiment (see Appendix VII D).The inset illustrates S(A, ∆ϕ) up to critical amplitude A = 0.5.The weights hi, gi 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 ∆ϕ between pump signals as in frame (a).The inset illustrates measured S(A, ∆ϕ) up to A = 0.5.In general, the measured S(A, ∆ϕ) corresponds quite well to the inset in frame (a).Due to noise, the measured GME nearly vanishes around ∆ϕ +90 • where even the simulated S is only slightly below 1.The best genuine multipartite entanglement is reached at ∆ϕ = −90 • • • • − 120 • owing to phase shifts introduced by the cavity (see Fig.11in Appendix VII E).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.

FIG. 6 .
FIG.6.Phase-independent full inseparability of tripartite bisqueezed state.a) PPT criteria in terms of the minimum symplectic eigenvalues simulated for our doublepump QWJPA using experimentally determined parameters.Eigenvalues min{νi} 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), which participates to both TMS processes.b) Experimentally determined symplectic eigenvalues for the same permutations (•, •, •); the solid lines are just to guide the eyes.We find full inseparability at normalized pumping 0.05 A 0.3 in the experiment.Grey dashed line displays the full inseparability threshold.The difference in the simulated behavior of ν1−23 and ν3−12 is caused by asymmetry due to finite value of resonance detuning ∆r.
FIG. 8.Covariance matrix of genuinely entangled quadripartite H-graph state.a) Experimental covariance where the rotation of the TMS subspaces 1-2, 2-3, amd 3-4 have been made in such a way that the structure coincides with the matrix in Eq. (54) of Appendix VII F (each pump has phase π 2 ).The employed pump amplitude A 0.08 yields the smallest value for S. b) Simulated covariance matrix using equal pump phases π 2 at A 0.08.The difference from the matrix in Eq. (54) is due to the cavity response that induces extra phase shifts.

10 FIG. 9 .
FIG. 9. Gain calibration using linear temperature dependence of the measured thermal noise spectral density of a 50 Ohm terminator measured as a function of Ts, the source temperature.The average total gain GΣ,i = 94.4 ± 0.2 dB over the cavity resonance is obtained from the linear fit (in red) to the data.This gain value GΣ,i also includes frequency mixing losses/amplification in the signal analyzer circuit part.The term Tpreamp = THEMT + T RT G HEMT = 5.2 ± 0.25K characterizes the equivalent noise temperature of the amplifiers; the largest contribution originates from the cooled HEMT amplifier at 4 K.The value coth( hf i 2k b Tpreamp ) sets the background for the diagonal elements in the covariance matrix 4V off .

FIG. 11 .
FIG.11.Cavity phase response given by fitted experimental parameters κ 2π = 4.44 MHz and γ 2π = 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 and ∆ϕ2 = π 2 + π 4 = 3π 4 with corresponding phase shift between pump tones π 2 , given in results (See.Fig.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.

2 }
we are able to flip the sign of one minor diagonal indicated by bold font in Eq. (55).