# Generation and structuring of multipartite entanglement in
Josephson parametric system

K. V. Petrovnin<EMAIL_ADDRESS>QTF Centre of Excellence,
Department of Applied Physics, Aalto University, P.O. Box 15100, FI-00076
AALTO, Finland  M. R. Perelshtein  QTF Centre of Excellence, Department of
Applied Physics, Aalto University, P.O. Box 15100, FI-00076 AALTO, Finland  T.
Korkalainen  QTF Centre of Excellence, Department of Applied Physics, Aalto
University, P.O. Box 15100, FI-00076 AALTO, Finland  V. Vesterinen  QTF Centre
of Excellence, VTT Technical Research Centre of Finland Ltd,
P.O. Box 1000, FI-02044 VTT, Finland  I. Lilja  QTF Centre of Excellence,
Department of Applied Physics, Aalto University, P.O. Box 15100, FI-00076
AALTO, Finland  G. S. Paraoanu  QTF Centre of Excellence, Department of
Applied Physics, Aalto University, P.O. Box 15100, FI-00076 AALTO, Finland
InstituteQ – the Finnish Quantum Institute, Aalto University, Finland  P. J.
Hakonen<EMAIL_ADDRESS>QTF Centre of Excellence, Department of
Applied Physics, Aalto University, P.O. Box 15100, FI-00076 AALTO, Finland
InstituteQ – the Finnish Quantum Institute, Aalto University, Finland

###### Abstract

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, we achieve generation of genuine tripartite
entangled state and its control by the use of the phase difference between two
continuous pump tones. We demonstrate control of the subspaces of the
covariance matrix for tripartite bisqueezed state. Furthermore, by optimizing
the phase relationships in a three-tone pumping scheme we explore genuine
quadripartite entanglement of a generalized H-graph state
($\mathscr{\tilde{H}}$-graph). Our scheme provides a comprehensive control
toolbox for the entanglement structure and allows us to demonstrate, for first
time to our knowledge, genuine quadripartite entanglement of microwave modes.
All experimental results are verified with numerical simulations of the
nonlinear quantum Langevin equation. We envision that quantum resources
facilitated by multi-pump configurations offer enhanced prospects for quantum
data processing using parametric microwave cavities.

##  I Introduction

Entanglement is a crucial resource in advanced information processing based on
quantum mechanical concepts [1]. To exceed the computational power of
classical devices, quantum devices typically need to employ highly entangled
states. The controlled generation of entanglement is the key resource not only
in quantum computing, but also in sensing and secure communications.

One of the most easily accessible and, at the same time, strong sources for
entanglement generation is the vacuum state of a quantum field [2]. Squeezing,
which is the fundamental operation for the continuous variable (CV) states
production, allows generation of coherence and entanglement from vacuum
fluctuations [3]. While two-mode squeezing produces bipartite quantum states,
multipartite states can be generated by applying similar operations [4].
Intriguingly, multipartite CV states are shown to enable various promising
phenomena such as quantum state sharing [5] and secret sharing [6], dense
coding [7], error correction [8] and quantum teleportation [9]. Furthermore,
CV cluster states show potential as a universal quantum computing platform
[10, 11, 12], which has been under active development for the past 20 years.
Cluster state calculus, foremost utilizing optical resources [13, 11, 14],
realize measurement-based quantum computing [15].

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 [4] and strong entanglement between different modes
[16].

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 [17, 18, 19], 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 [12, 20, 21, 22]. 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
vacuum-induced 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 parametric 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

###  II.1 Quantum dynamics of the device

Our work employs a parametric system comprised of a superconducting
$\lambda/4$ resonator terminated in a SQUID loop. Such a setup forms the
archetype of a narrow-band superconducting Josephson parametric amplifier
(JPA) [23, 24]. 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 $\omega_{d}\sim 2\omega_{r}$ (three-wave
mixing) [4, 25]. In the rotating frame, the Hamiltonian of the system, as
derived in Refs. 26, 23, is given by:

$\displaystyle{H}_{\text{sys,rwa}}(t)={\hbar}\Delta_{r}\tilde{a}^{\dagger}\tilde{a}$
$\displaystyle+$
$\displaystyle\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},$ (1)

where $\tilde{a}$ ($\tilde{a}^{\dagger}$) is the annihilation (creation)
operator for cavity photons in the rotating frame at angular frequency
$\omega_{\Sigma}/2$, $\alpha_{d}$ is the complex amplitude of the $d$-th pump
tone and $\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
$\alpha_{d}=|\alpha|e^{i\varphi_{d}}$. Here $\Delta_{r}$ denotes the detuning
between half of the average pump angular frequency $\omega_{\Sigma}$ and the
resonator angular frequency $\omega_{r}$:
$\Delta_{r}=\omega_{r}-\omega_{\Sigma}/2$, with $\omega_{\Sigma}$ representing
the average angular frequency in multi-tone driven case:
$\omega_{\Sigma}=(1/p)\sum_{d=1}^{p}\omega_{d}$ with
$d=\left\\{1,\dots,p\right\\}$ 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 $\alpha\geq\alpha_{crit}$. Several effects are
accounted for by the Kerr nonlinearity, such as limited maximum gain,
compression, observed at $\alpha\lesssim\alpha_{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.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:

${H}(t)=H_{\text{sys,rwa}}(t)+H_{\text{sig}}+H_{\text{loss}},$ (2)

where $H_{\text{sig}}$ includes the coupling to the signal port transmission
line with dissipation rate $\kappa$, while $H_{\text{loss}}$ includes the
coupling to the internal loss port with linear dissipation rate $\gamma$.

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

$\displaystyle\dot{\tilde{a}}(t)=(-i\Delta_{r}$
$\displaystyle-\frac{\kappa+\gamma}{2})\tilde{a}-i\sum_{d=1}^{p}\alpha_{d}e^{i\Delta_{d}t}\tilde{a}^{\dagger}$
(3)
$\displaystyle+\sqrt{\kappa}\tilde{b}_{in}+\sqrt{\gamma}\tilde{c}_{in}-12i{K}\tilde{a}^{\dagger}\tilde{a}\tilde{a}$

where $\tilde{b}_{in}$ and $\tilde{c}_{in}$ are the ladder (annihilation)
operators for the signal and linear dissipation ports, respectively.

The output mode $\tilde{b}_{out}$ is obtained using the following relation
between the incoming and outgoing modes:

$\tilde{b}_{out}(t)=\tilde{b}_{in}(t)-\sqrt{\kappa}\tilde{a}(t).$ (4)

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.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.

Figure 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 vacuum-induced microwaves (blue
arrows). 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
$\Delta$ in Fourier analysis is chosen to be much narrower than the cavity
resonance width.

###  II.2 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).

We consider the fundamental resonance of a transmission line JPA centered at
$\omega_{r}$ frequency with a bandwidth $2\delta\omega$, within which
($[\omega_{r}-\delta\omega,\omega_{r}+\delta\omega]$) we define $N$ spectral
modes as depicted in Fig. 1. Let us define $\tilde{a}$ as a vector of spectral
modes:

$\tilde{a}=\\{\tilde{a}_{1},\dots,\tilde{a}_{N},\tilde{a}^{\dagger}_{1},\dots,\tilde{a}^{\dagger}_{N}\\}^{T},$
(5)

where $N$ is a total mode number and the creation
$\tilde{a}_{i}^{\dagger}=\tilde{a}_{i}^{\dagger}(t)$ and annihilation
$\tilde{a}_{i}=\tilde{a}_{i}(t)$ operators are time-dependent. In general, the
frequency difference between half of $p$-th and $(p+1)$-th pump tones defines
the bandwidth of the spectral mode $\tilde{a}_{i}$. We employ an equidistant
pump scheme where bandwidth of each mode $\tilde{a}_{i}$ is defined as
$\Delta$ such that the spectrum of a full set of modes $\tilde{a}_{i}$ covers
the bandwidth $2\delta\omega$ of the cavity mode $\tilde{a}$. In the
experiment, we collect the emitted power over the whole
$[\omega_{r}-\delta\omega,\omega_{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 $\Delta$.

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.

###  II.3 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.6. Assuming the strong coupling regime with $\kappa\gg\Delta$ while
neglecting dissipation losses $\gamma$ and the Kerr-nonlinearity $K$, Fourier
transformation yields a system of linear equations which can be cast in a
matrix form:

$\mathbf{M}\tilde{a}(\omega)=\sqrt{\kappa}\tilde{b}_{\text{in}}(\omega).$ (6)

Here the interaction matrix $\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 $\Delta_{r}=0$ and having different phases for the pump tones
$\alpha_{1}=\alpha e^{i\varphi_{1}},\alpha_{2}=\alpha e^{i\varphi_{2}}$:

$\mathbf{M}=\leavevmode\resizebox{355.57156pt}{}{{\hbox{ \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},}}
}$ (7)

where the frequency dependency enters through the $c=-i\omega+\kappa/2$
coefficient. To express the intracavity modes through the vacuum input, we use
the inverse matrix $\mathbf{M}^{-1}$:

$\tilde{a}(\omega)=\sqrt{\kappa}\mathbf{M}^{-1}\tilde{b}_{\text{in}}(\omega).$
(8)

For the tripartite case, such a matrix is written in the following way

$\displaystyle\mathbf{M}^{-1}=\frac{1}{c^{2}-2\alpha^{2}}\cdot$ (9)

where $\Delta\varphi=\varphi_{1}-\varphi_{2}$. Besides two mode squeezing
(TMS) correlations proportional to $\alpha$, the matrix contains also
beamsplitter correlations (BS) $\propto\alpha^{2}$ (denoted in bold). Note
that the $\alpha^{2}$-terms are absent in the matrix $\mathbf{M}$ introduced
in Eq. (6). In this tripartite example, the phase difference $\Delta\varphi$
contributes only to the BS connection
$\tilde{a}_{i}(\omega)\leftrightarrow\tilde{a}_{j}(\omega)$ or
$\tilde{a}^{\dagger}_{i}(\omega)\leftrightarrow\tilde{a}^{\dagger}_{j}(\omega)$.
TMS connections are defined by entries
$\tilde{a}_{i}(\omega)\leftrightarrow\tilde{a}^{\dagger}_{j}(\omega)$ in
$\mathbf{M}^{-1}$. In Appendix VII.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 $\mathbf{M}^{-1}$
responsible for beamsplitter interaction can be suppressed fully by choosing
pump phases properly in certain pump tone configurations.

Figure 2:  Graph representation for the entangled tripartite and quadripartite
systems. Generalized H-graph ($\mathscr{\tilde{H}}$-graph) representation of
bisqueezed tripartite CV entangled state (a) and quadripartite CV entangled
$\mathscr{\tilde{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 $\Delta_{d}=0$ in the tripartite case
and two additional tones at $-\Delta$ and $\Delta$ in the quadripartite
setting. The additional pumps supply missing TMS connections to the entangled
states. For details, see Appendix VII.6.

Interestingly, the matrix $\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 [12, 20, 21, 22], 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 $\mathscr{\tilde{H}}$-graphs formed
by both two mode squeezing and beamsplitter correlations (see Appendix VII.6).
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 [27, 20, 22] 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 $\tilde{b}_{\text{out}}$, which are
obtained from $\tilde{a}$ using the input-output relationship in Eq. (4). This
equation yields

$\tilde{b}_{\text{out}}(\omega)=(\mathbf{I}-\kappa\mathbf{M}^{-1})\tilde{b}_{\text{in}}(\omega),$
(10)

on the basis of which we may define the adjacency matrix
$\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 $\mathbf{M}^{-1}$ is
sufficient to characterize the BS and TMS connections between vacuum spectral
modes.

###  II.4 Connection to Hamiltonian graph

CV cluster states with square-lattice graph structure provide a foundation for
measurement-based continuous variable quantum computation (CVQC) [20]. Cluster
states can be asymptotically reached from H-graph states in the case of
infinite squeezing [21]. The H-graph structure is defined by its adjacency
matrix $\mathbf{G}$, whose entries $G_{ij}$ specify the multimode squeezing
Hamiltonian as:

$H_{\text{S}}=\hbar\alpha\sum_{i,j}G_{ij}(\tilde{a}_{i}^{\dagger}\tilde{a}_{j}^{\dagger}+\tilde{a}_{i}\tilde{a}_{j}).$
(11)

Here, the pump tone amplitudes are considered to have equal strength $\alpha$.
The matrix $\mathbf{G}$ involves TMS correlations between modes
$\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
$i\dot{\tilde{a}}_{k}^{\dagger}=\alpha\sum\nolimits_{j}{{G_{jk}}}{\tilde{a}_{j}}$
and
$i{\dot{\tilde{a}}_{k}}=-\alpha\sum\nolimits_{j}{G_{jk}^{*}}\tilde{a}_{j}^{\dagger}$.
Taking the Fourier transform, the left hand side equals
$\omega\times\tilde{a}(\omega)$, and the combination
$\omega\tilde{a}_{k}(\omega)=\alpha\sum\nolimits_{j}{{G_{jk}}}{\tilde{a}_{j}}^{\dagger}(-\omega)$
provides the connection to the QLE treatment in Eq. (6): $\tilde{a}(\omega)$
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
$\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
$\mathbf{Z}=ie^{-\mathbf{G}}$, which is interpreted as the adjacency matrix
for an undirected Gaussian graph with complex-valued edge weights [20, 21].
Decomposing such a matrix up to quadratic terms
$\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 [28].

Let us now show the origin of BS correlations using the multimode squeezing
Hamiltonian of Eq. 11 in $R=e^{-\frac{i}{\hbar}H\tau}$ where we have
considered that the system is pumped for a finite time $\tau$. The multimode
squeezing operator $R$ can be decomposed to a combination of TMS operators,
containing
$B_{ij}=\tilde{a}_{i}^{\dagger}\tilde{a}_{j}^{\dagger}+\tilde{a}_{i}\tilde{a}_{j}$,
and BS transformations based on
$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) [28], we obtain for the tripartite squeezing operator

$\displaystyle
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}}$ $\displaystyle 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}}.$ (12)

Here, $\theta_{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. 28.

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

$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,$
(13)

where the total number of TMS and BS operators in the decomposition is
$N_{\text{BS}}=\sum_{n}N-2n$ and $N_{\text{TMS}}=\sum_{n}N-2n+1$; where
$n=1,2,\dots,\lfloor\frac{N}{2}\rfloor$ 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
$\mathbf{\tilde{G}}$ with entries
$\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
$\mathbf{\tilde{G}}$ defines the edge connections in the generalized
$\mathscr{\tilde{H}}$-graph.

###  II.5 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 [3, 29]. Therefore, the
output states of a $N$-mode JPA can be fully characterized by its covariance
matrix of $2N$-length column vector with quadratures
$\tilde{r}~{}=~{}\left(\tilde{x}_{1},\tilde{p}_{1},\dots\tilde{x}_{N},\tilde{p}_{N}\right)^{T}$,
where $\tilde{x}_{i}=(\tilde{a}_{i}+\tilde{a}^{\dagger}_{i})/2$ and
$\tilde{p}_{i}=(\tilde{a}_{i}-\tilde{a}^{\dagger}_{i})/2i$. The covariance
matrix $\mathscr{V}$, whose elements are given by

$\mathscr{V}_{ij}=2\left<\Delta\tilde{r}_{i}\Delta\tilde{r}_{j}+\Delta\tilde{r}_{j}\Delta\tilde{r}_{i}\right>-4\left<\Delta\tilde{r}_{i}\right>\left<\Delta\tilde{r}_{j}\right>,$
(14)

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
$\Delta\tilde{r}_{i}=\tilde{r}_{i}-\left<\tilde{r}_{i}\right>$. 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
[30]. 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 [31, 32].

To examine inseparability properties of the quantum state [33, 34, 35, 30, 36,
37, 38], we apply symplectic transformations to the covariance matrix and
calculate its symplectic eigenvalues – the PPT criteria [39, 40]. Such
transformations are equivalent to a phase space reflection of a single party
in the $N$-partite state [33]. All minimum symplectic eigenvalues
$\\{\nu_{i}\\}_{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. 37, 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
[3] 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 [41]. In particular, there is an universal
GME criterion derived in Ref. [30] and further refined in Ref. [37]. 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 [16].

The GME criterion is based on the weighted variance of the quadratures,
$u=\sum_{i}h_{i}\tilde{x}_{i}$ and $v=\sum_{k}g_{k}\tilde{p}_{k}$;
$i,k=\left\\{1,2,3,\dots,N\right\\}$. Violation of the inequality

$\displaystyle S\equiv\frac{\left<\Delta u^{2}\right>+\left<\Delta
v^{2}\right>}{f_{3}(h_{i},g_{i})}\geq 1,$ (15)

where

$\displaystyle f_{3}(h_{i},g_{i})=\frac{1}{2}\min\\{$
$\displaystyle|h_{1}g_{1}+h_{2}g_{2}|+|h_{3}g_{3}|,$
$\displaystyle|h_{3}g_{3}+h_{2}g_{2}|+|h_{1}g_{1}|,$
$\displaystyle|h_{1}g_{1}+h_{3}g_{3}|+|h_{2}g_{2}|,\\}$

is sufficient to confirm genuine tripartite entanglement ($N=3$) and violation
of

$\displaystyle S\equiv\frac{\left<\Delta u^{2}\right>+\left<\Delta
v^{2}\right>}{f_{4}(h_{i},g_{i})}\geq 1,$ (16)

where

$\displaystyle f_{4}(h_{i},g_{i})=\min\\{$
$\displaystyle|h_{1}g_{1}+h_{2}g_{2}+h_{3}g_{3}|+|h_{4}g_{4}|,$
$\displaystyle|h_{4}g_{4}+h_{2}g_{2}+h_{3}g_{3}|+|h_{1}g_{1}|,$
$\displaystyle|h_{4}g_{4}+h_{1}g_{1}+h_{3}g_{3}|+|h_{2}g_{2}|,$
$\displaystyle|h_{4}g_{4}+h_{1}g_{1}+h_{2}g_{2}|+|h_{3}g_{3}|,$
$\displaystyle|h_{1}g_{1}+h_{2}g_{2}|+|h_{3}g_{3}+h_{4}g_{4}|,$
$\displaystyle|h_{1}g_{1}+h_{3}g_{3}|+|h_{2}g_{2}+h_{4}g_{4}|,$
$\displaystyle|h_{2}g_{2}+h_{3}g_{3}|+|h_{1}g_{1}+h_{4}g_{4}|\\}$

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 $h_{i}=h,g_{i}=g$, for $i=\left\\{2,3\right\\}$ or
$i=\left\\{2,3,4\right\\}$ with respect to the number of parties $N$ [37].

In the double and triple pumping schemes, we generate generalized tripartite
and quadripartite $\mathscr{H}$-graph states, which have a different structure
compared with Greenberger-Horne-Zeilinger (GHZ) states [27, 20, 22]. The
generalization deals with addition of BS correlations in the
$\mathscr{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

###  III.1 Experimental methods

Figure 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 $\Omega$ 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.

In our microwave experiments, we employ a niobium $\lambda/4$ coplanar
$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$\times$1 $\mu$m tunnel barriers with $I_{c}\simeq 4$
$\mu$A critical current. The JPA, operating in three-wave-mixing mode around
the cavity frequency $\omega_{r}$, is pumped by an external RF magnetic flux
through the SQUID loop using a single turn pump coil at frequency
$2\omega_{r}$ (marked as $2\omega_{\mathrm{LO}}$ in Fig. 3a) [23, 42]. We
chose this operation regime, because for four-wave mixing (typically using a
current pump near $\omega_{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 $\sim$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.

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 $\kappa$, $\gamma$ and $K$ are discussed in Appendix VII.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 $\omega_{\mathrm{2LO}}$, a band
rejection filter tuned to $2\omega_{\mathrm{LO}}$ is employed. The filtering
ensures passage only for the desired two pump signal at angular frequencies
$\omega_{1}=\omega_{r}-\Delta/2$ and $\omega_{2}=\omega_{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 $\Phi_{DC}=0.383\Phi_{0}$ through the SQUID
loop, resulting in $\frac{\omega_{r}}{2\pi}=6.024$ GHz for the cavity
frequency. In the three-mode experiment, we apply two pump tones at
$(2\times\frac{\omega_{r}}{2\pi}-2)$ MHz and
$(2\times\frac{\omega_{r}}{2\pi}+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 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 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. Therefore, a change of the initial phase of
the IF-signal in one of the pump generator channels (Agilent 33600B in Fig.
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(\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\times 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.

#### 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 $\Phi_{DC}=0.417\,\Phi_{0}$ through a SQUID
resulting in $\omega_{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\times\frac{\omega_{r}}{2\pi}$ MHz,
$(2\times\frac{\omega_{r}}{2\pi}+1)$ MHz, and
$(2\times\frac{\omega_{r}}{2\pi}-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 $\mathcal{\tilde{H}}$-graph states and their
genuine multipartite entanglement depends on relations among the covariance
matrix elements as discussed in Section II.5. The degree of violation in the
GME condition $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.3.

###  III.2 Scaled covariance matrix

The system gain $\left<G_{\Sigma}\right>$ determined in Appendix VII.3 refers
to measured power per unit band width. Since the measured spectral mode
quadratures $I_{i}$ and $Q_{i}$ are determined over the band $\Delta f_{i}$,
the scaled quadrature $x_{i}$, equivalent to the amplitude of the quantum
mechanical operators $\mathbf{x}$, is given by the formula

${x}_{i}=\frac{I_{i}}{\sqrt{G_{i}Z_{0}hf_{i}\Delta f_{i}}},$ (17)

where $G_{i}=\left<G_{\Sigma,i}\right>$ is the system gain for
$i^{\mathrm{th}}$ spectral mode, $Z_{0}=50$ $\Omega$ is transmission line
impedance and $\Delta f_{i}$ is the bandwidth of the spectral mode: $\Delta
f_{i}=2$ MHz or $\Delta 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 [26], the noise added by the preamplifier is
subtracted from the diagonal elements of the covariance matrix (see Eq. 14):

$\mathscr{V}=4\left(\mathbf{V}_{\text{on}}-\mathbf{V}_{\text{off}}\right)+\mathbf{I}\coth\frac{hf_{i}}{2k_{b}T_{i}}.$
(18)

where $\mathbf{V}_{\text{off}}$ denotes the covariance matrix measured in the
absence of the pump. Due to scaling of the covariance matrix $\mathscr{V}$,
this equation yields a unity diagonal matrix in the absence of pumping at
$T\rightarrow 0$. The average physical temperature in our experiments is
$T_{i}=20$ mK resulting in $\coth\frac{hf_{i}}{2k_{b}T_{i}}=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.5 for tripartite and
quadripartite cases.

#### 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.

Figure 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
$\Delta\varphi$ between pump signals indicated in the figure; the simulation
parameters were set to match the experiment (see Appendix VII.4). The inset
illustrates $S(A,\Delta\varphi)$ up to critical amplitude $A=0.5$. The weights
$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,\Delta\varphi)$ up to $A=0.5$. In general, the
measured $S(A,\Delta\varphi)$ corresponds quite well to the inset in frame
(a). Due to noise, the measured GME nearly vanishes around
$\Delta\varphi\simeq+90^{\circ}$ where even the simulated $S$ is only slightly
below 1. The best genuine multipartite entanglement is reached at
$\Delta\varphi=-90^{\circ}\dots-120^{\circ}$ owing to phase shifts introduced
by the cavity (see Fig. 11 in Appendix VII.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$.

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 $\Delta\varphi$ yielding $S<1$ (see the inset in Fig. 4a).
The strongest genuine tripartite entanglement is reached at
$\Delta\varphi\simeq-120^{\circ}$ under normalized pumping amplitude $A\simeq
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
$\Delta\varphi=+90^{\circ}$ yields clearly worse entanglement than
$\Delta\varphi=-90^{\circ}$. This asymmetry in GME between $\Delta\varphi=\pm
90^{\circ}$ 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
$\Delta\varphi=-90^{\circ}$ and $\Delta\varphi=-120^{\circ}$ yield suppression
down to $S=0.75\pm 0.05$. The measured result at $\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\in[\sim 0.01,0.4]$. Overall, the pattern of $S(\varphi,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
$\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 VII.5), which, when included in the simulations, result in very good
matching with the experiment.

Figure 5:  Covariance 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^{\circ}$. 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.

Covariance matrices measured at the pump phase difference
$\Delta\varphi=+90^{\circ}$ and $\Delta\varphi=-90^{\circ}$ 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 $\Delta\varphi=-90^{\circ}$ to
$\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 $\pi$ which inverts the
TMS correlation in subspace $1-2$. These subspaces which are controlled by the
pump phase settings, can expressed in symmetric
$\left<I_{1}I_{2}\right>\simeq\left<I_{2}I_{3}\right>$ or antisymmetric form
$\left<I_{1}I_{2}\right>\simeq-\left<I_{2}I_{3}\right>$, see Figs. 5a and 5b,
respectively.

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.

Figure 6:  Phase-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\\{\nu_{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 ($\nu_{2-13}$), which participates to both TMS processes. b)
Experimentally determined symplectic eigenvalues for the same permutations
($\circ$, $\circ$, $\circ$); the solid lines are just to guide the eyes. We
find full inseparability at normalized pumping $0.05\lesssim A\lesssim 0.3$ in
the experiment. Grey dashed line displays the full inseparability threshold.
The difference in the simulated behavior of $\nu_{1-23}$ and $\nu_{3-12}$ is
caused by asymmetry due to finite value of resonance detuning $\Delta_{r}$.
Figure 7:  Full inseparability and genuine entanglement of quadripartite
$\mathscr{\tilde{H}}$-graph state (generalized H-graph). a) Results on PPT
criterion for four permutations of a four-mode Gaussian state. Minimum
eigenvalues $\min\\{\nu_{i}\\}$, 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
$\\{\nu_{1-234}<1\\}\bigwedge\\{\nu_{2-134}<1\\}\bigwedge\\{\nu_{3-124}<1\\}\bigwedge\\{\nu_{4-123}<1\\}$
is fulfilled in the range of $0.01\lesssim A\lesssim 0.15$. The green dashed
line displays the entanglement threshold. Figure 8:  Covariance matrix of
genuinely entangled quadripartite $\mathscr{\tilde{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.6 (each pump has phase $\frac{\pi}{2}$). The
employed pump amplitude $A\eqsim 0.08$ yields the smallest value for $S$. b)
Simulated covariance matrix using equal pump phases $\frac{\pi}{2}$ at
$A\eqsim 0.08$. The difference from the matrix in Eq. (54) is due to the
cavity response that induces extra phase shifts.

#### Quadripartite case

For the quadripartite case, we apply three pump drives with identical
amplitude $\alpha_{1}=\alpha_{2}=\alpha_{3}=\alpha$. 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).

For the analysis of full inseparability of the covariance matrix according to
the PPT criterion, we evaluate the minimum symplectic eigenvalues
$\min\\{\nu_{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\\{\nu_{i}\\}=0.79\pm 0.018$ is reached,
while all of the eigenvalues in the normalized pump amplitude range
$0.01\lesssim A\lesssim 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\\{\nu_{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\pm 0.02$ is
observed at $A\eqsim 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\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 Fig. 8 correspond directly to the
$\mathscr{\tilde{H}}$-graph structures shown on Fig. 2b. In general, we
conclude that for the employed pump configuration, the genuine quadripartite
entanglement appears in the amplitude range $0.01\lesssim A<0.13$.

##  V Discussion

The control of bisqueezed tripartite and generalized H-graph
($\mathscr{\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
$\mathscr{\tilde{H}}$-graph states. Enhanced number of modes requires larger
bandwidth, which calls for broadband parametric devices such as TWPAs [43, 44,
45, 46, 47] or broadband JPAs [48, 42] 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 $\mathscr{\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 VII.6 (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, $\mathscr{\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 [15], 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 [10,
11, 12]. For quantum information processing steps, it is sufficient to perform
sub-cluster measurements in specified order using a suitable computational
basis. In Refs. 11, 49, 50 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
$\mathscr{\tilde{H}}$-graph state encoding.

In contrast to computational models for graph states [20] considered as ideal
clusters, hardware based on finite squeezing with noise and decoherence
requires error correction procedures [51, 52, 53] 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 $\mathscr{\tilde{H}}$-graph states before each gate operation. In
Ref. 54 a multidimensional platform for scalable quantum computing has been
proposed, based on cluster states created using microring resonators; also
multiple frequency combs [50] 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 $\mathscr{\tilde{H}}$-graph states.

An implementation of the universal quantum computer based on bosonic modes
with the possibility of hardware-efficient quantum error correction [55]
requires efficient generation of continuous-variable quantum resources. The
genuine entanglement between several bosonic modes could potentially be
employed for error-correctable codeword states [56]. 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 [57].

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 [45, 46, 47]. 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,
vacuum-induced multipartite entanglement 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\lesssim A<0.4$ ($0.01\lesssim 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
$\Delta\varphi=-120^{\circ}$ minimum value of criterion $S=0.75\pm 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 $\mathscr{\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 VII.6, BS correlations can be fully suppressed by implementing a
180o 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
$\\{-\frac{\Delta}{2};\frac{\Delta}{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-of-the-
art superconducting parametric devices, such as the broadband, low-loss
travelling wave parametric amplifier [45, 46, 47], 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.

## Acknowledgments

We thank Mark Dykman, Stefano Pirandola, Olivier Pfister, Aidar Sultanov,
Shruti Dogra and Marko Kuzmanović for fruitful discussions and useful
comments. The work of KVP has been funded from the European Union’s Horizon
2020 research and innovation programme under grant agreement no. 862644 (FET-
Open project: Quantum readout techniques and technologies QUARTET). PJH
acknowledges ERC grant agreement no. 670743 (QuDeT) and grant agreement no.
824109 (European Microkelvin Platform). The contribution of PJH and IL was
supported by MATINE grant VN/13774/2020-PLM-49. GSP and KVP would like to
thank Saab for scientific collaboration under a research agreement with Aalto
University. The work at VTT and Aalto was additionally supported by the
Academy of Finland through the Finnish Center of Excellence in Quantum
Technology QTF (projects 312295, 312296, 336810, 336813). This research
benefited of the facilities of Micronova Nanofabrication Centre and Low
Temperature Laboratory of the Aalto University OtaNano infrastructure.

## References

  * Jozsa and Linden [2003] R. Jozsa and N. Linden, On the role of entanglement in quantum-computational speed-up, Proc. R. Soc. Lond. A. 459, 2011 (2003).
  * Summers and Werner [1985] S. J. Summers and R. Werner, The vacuum violates bell’s inequalities, Phys. Lett. A 110, 257 (1985).
  * Braunstein and van Loock [2005] S. L. Braunstein and P. van Loock, Quantum information with continuous variables, Rev. Mod. Phys. 77, 513 (2005).
  * Lähteenmäki  _et al._ [2016] P. Lähteenmäki, G. S. Paraoanu, J. Hassel, and P. J. Hakonen, Coherence and multimode correlations from vacuum fluctuations in a microwave superconducting cavity, Nat. Commun. 7, 1 (2016).
  * Lance  _et al._ [2005] A. M. Lance, T. Symul, W. P. Bowen, B. C. Sanders, T. Tyc, T. C. Ralph, and P. K. Lam, Continuous-variable quantum-state sharing via quantum disentanglement, Phys. Rev. A 71, 033814 (2005).
  * Hillery  _et al._ [1999] M. Hillery, V. Bužek, and A. Berthiaume, Quantum secret sharing, Phys. Rev. A 59, 1829 (1999).
  * Braunstein and Kimble [2000] S. L. Braunstein and H. J. Kimble, Dense coding for continuous variables, Phys. Rev. A 61, 042302 (2000).
  * Dias and Ralph [2018] J. Dias and T. C. Ralph, Quantum error correction of continuous-variable states with realistic resources, Phys. Rev. A 97, 032335 (2018).
  * Milburn and Braunstein [1999] G. J. Milburn and S. L. Braunstein, Quantum teleportation with squeezed vacuum states, Phys. Rev. A 60, 937 (1999).
  * Zhang and Braunstein [2006] J. Zhang and S. L. Braunstein, Continuous-variable Gaussian analog of cluster states, Phys. Rev. A 73, 032318 (2006).
  * Menicucci  _et al._ [2006] N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, Universal quantum computation with continuous-variable cluster states, Phys. Rev. Lett. 97, 110501 (2006).
  * Gu  _et al._ [2009] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and P. van Loock, Quantum computing with continuous-variable clusters, Phys. Rev. A 79, 062318 (2009).
  * Mantri  _et al._ [2017] A. Mantri, T. F. Demarie, and J. F. Fitzsimons, Universality of quantum computation with cluster states and (x, y)-plane measurements, Sci. Rep. 7, 1 (2017).
  * Walther  _et al._ [2005] P. Walther, K. J. Resch, T. Rudolph, E. Schenck, H. Weinfurter, V. Vedral, M. Aspelmeyer, and A. Zeilinger, Experimental one-way quantum computing, Nature 434, 169 (2005).
  * Raussendorf and Briegel [2001] R. Raussendorf and H. J. Briegel, A one-way quantum computer, Phys. Rev. Lett. 86, 5188 (2001).
  * Sandbo Chang  _et al._ [2018] 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, and C. M. Wilson, Generating multimode entangled microwaves with a superconducting parametric cavity, Phys. Rev. Applied 10, 044019 (2018).
  * Adesso and Illuminati [2007] G. Adesso and F. Illuminati, Entanglement in continuous-variable systems: recent advances and current perspectives, J. Phys. A: Math. Theor. 40, 7821 (2007).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Rev. Mod. Phys. 84, 621 (2012).
  * Fabre and Treps [2020] C. Fabre and N. Treps, Modes and states in quantum optics, Rev. Mod. Phys. 92, 035005 (2020), 1912.09321 .
  * Menicucci  _et al._ [2011] N. C. Menicucci, S. T. Flammia, and P. van Loock, Graphical calculus for Gaussian pure states, Phys. Rev. A 83, 042335 (2011).
  * Bruschi  _et al._ [2016] D. E. Bruschi, C. Sabín, P. Kok, G. Johansson, P. Delsing, and I. Fuentes, Towards universal quantum computation through relativistic motion, Sci. Rep. 6, 1 (2016).
  * Pfister [2019] O. Pfister, Continuous-variable quantum computing in the quantum optical frequency comb, J. Phys. B: At. Mol. Opt. Phys. 53, 012001 (2019).
  * Yamamoto  _et al._ [2016] T. Yamamoto, K. Koshino, and Y. Nakamura, Parametric amplifier and oscillator based on josephson junction circuitry, in  _Principles and Methods of Quantum Information Technologies_, edited by Y. Yamamoto and K. Semba (Springer Japan, Tokyo, 2016) pp. 495–513.
  * Wang  _et al._ [2019] Z. Wang, M. Pechal, E. A. Wollack, P. Arrangoiz-Arriola, M. Gao, N. R. Lee, and A. H. Safavi-Naeini, Quantum dynamics of a few-photon parametric oscillator, Phys. Rev. X 9, 021049 (2019).
  * Korkalainen  _et al._ [2021] T. Korkalainen, I. Lilja, M. R. Perelshtein, K. V. Petrovnin, G. S. Paraoanu, and P. J. Hakonen, Vacuum-induced correlations in superconducting microwave cavity under multiple pump tones, AIP Conference Proceedings 2362, 030001 (2021).
  * Lähteenmäki  _et al._ [2013] P. Lähteenmäki, G. Paraoanu, J. Hassel, and P. J. Hakonen, Dynamical casimir effect in a Josephson metamaterial, PNAS 110, 4234 (2013).
  * Briegel and Raussendorf [2001] H. J. Briegel and R. Raussendorf, Persistent entanglement in arrays of interacting particles, Phys. Rev. Lett. 86, 910 (2001).
  * Bruschi  _et al._ [2017] D. E. Bruschi, C. Sabín, and G. S. Paraoanu, Entanglement, coherence, and redistribution of quantum resources in double spontaneous down-conversion processes, Phys. Rev. A 95, 062324 (2017).
  * Adesso  _et al._ [2014] G. Adesso, S. Ragy, and A. R. Lee, Continuous variable quantum information: Gaussian states and beyond, Open Systems & Information Dynamics 21, 1440001 (2014).
  * van Loock and Furusawa [2003] P. van Loock and A. Furusawa, Detecting genuine multipartite continuous-variable entanglement, Phys. Rev. A 67, 052315 (2003).
  * Hillery  _et al._ [2010] M. Hillery, H. T. Dung, and H. Zheng, Conditions for entanglement in multipartite systems, Phys. Rev. A 81, 062322 (2010).
  * Agustí  _et al._ [2020] A. Agustí, C. W. S. Chang, F. Quijandría, G. Johansson, C. M. Wilson, and C. Sabín, Tripartite genuine non-gaussian entanglement in three-mode spontaneous parametric down-conversion, Phys. Rev. Lett. 125, 020502 (2020).
  * Simon [2000] R. Simon, Peres-horodecki separability criterion for continuous variable systems, Phys. Rev. Lett. 84, 2726 (2000).
  * Duan  _et al._ [2000] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Inseparability criterion for continuous variable systems, Phys. Rev. Lett. 84, 2722 (2000).
  * Mancini  _et al._ [2002] S. Mancini, V. Giovannetti, D. Vitali, and P. Tombesi, Entangling macroscopic oscillators exploiting radiation pressure, Phys. Rev. Lett. 88, 120401 (2002).
  * Hyllus and Eisert [2006] P. Hyllus and J. Eisert, Optimal entanglement witnesses for continuous-variable systems, New Journal of Physics 8, 51 (2006).
  * Teh and Reid [2014] R. Y. Teh and M. D. Reid, Criteria for genuine $n$-partite continuous-variable entanglement and einstein-podolsky-rosen steering, Phys. Rev. A 90, 062337 (2014).
  * Shchukin and van Loock [2015] E. Shchukin and P. van Loock, Generalized conditions for genuine multipartite continuous-variable entanglement, Phys. Rev. A 92, 042328 (2015).
  * Peres [1996] A. Peres, Separability criterion for density matrices, Phys. Rev. Lett. 77, 1413 (1996).
  * Horodecki  _et al._ [1996] M. Horodecki, P. Horodecki, and R. Horodecki, Separability of mixed states: necessary and sufficient conditions, Phys. Lett. A 223, 1 (1996).
  * Shalm  _et al._ [2013] L. K. Shalm, D. R. Hamel, Z. Yan, C. Simon, K. J. Resch, and T. Jennewein, Three-photon energy–time entanglement, Nat. Phys. 9, 19 (2013).
  * Elo  _et al._ [2019] T. Elo, T. S. Abhilash, M. R. Perelshtein, I. Lilja, E. V. Korostylev, and P. J. Hakonen, Broadband lumped-element Josephson parametric amplifier with single-step lithography, Appl. Phys. Lett. 114, 152601 (2019).
  * Macklin  _et al._ [2015] C. Macklin, K. O’Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, and I. Siddiqi, A near–quantum-limited josephson traveling-wave parametric amplifier, Science 350, 307 (2015), https://www.science.org/doi/pdf/10.1126/science.aaa8525 .
  * Zorin [2019] A. Zorin, Flux-driven Josephson traveling-wave parametric amplifier, Phys. Rev. Appl. 12, 044051 (2019).
  * Perelshtein  _et al._ [2021] M. R. Perelshtein, K. V. Petrovnin, V. Vesterinen, S. H. Raja, I. Lilja, M. Will, A. Savin, S. Simbierowicz, R. Jabdaraghi, J. Lehtinen, L. Grönberg, J. Hassel, M. Prunnila, J. Govenius, S. Paraoanu, and P. Hakonen, Broadband continuous variable entanglement generationusing a Kerr-free Josephson metamaterial (2021), arXiv:2111.06145 [quant-ph] .
  * Esposito  _et al._ [2021] M. Esposito, A. Ranadive, L. Planat, S. Leger, D. Fraudet, V. Jouanny, O. Buisson, W. Guichard, C. Naud, J. Aumentado, F. Lecocq, and N. Roch, Observation of two-mode squeezing in a traveling wave parametric amplifier (2021), arXiv:2111.03696 [quant-ph] .
  * Qiu  _et al._ [2022] 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, and W. D. Oliver, Broadband squeezed microwaves and amplification with a Josephson traveling-wave parametric amplifier (2022), arXiv:2201.11261 [quant-ph] .
  * Roy  _et al._ [2015] T. Roy, S. Kundu, M. Chand, A. M. Vadiraj, A. Ranadive, N. Nehra, M. P. Patankar, J. Aumentado, A. A. Clerk, and R. Vijay, Broadband parametric amplification with impedance engineering: Beyond the gain-bandwidth product, Appl. Phys. Lett. 107, 262601 (2015).
  * Alexander and Menicucci [2016] R. N. Alexander and N. C. Menicucci, Flexible quantum circuits using scalable continuous-variable cluster states, Phys. Rev. A 93, 062326 (2016).
  * Pooser and Jing [2014] R. Pooser and J. Jing, Continuous-variable cluster-state generation over the optical spatial mode comb, Phys. Rev. A 90, 043841 (2014).
  * Menicucci [2014] N. C. Menicucci, Fault-tolerant measurement-based quantum computing with continuous-variable cluster states, Phys. Rev. Lett. 112, 120504 (2014).
  * Larsen  _et al._ [2020] M. V. Larsen, J. S. Neergaard-Nielsen, and U. L. Andersen, Architecture and noise analysis of continuous-variable quantum gates using two-dimensional cluster states, Phys. Rev. A 102, 042608 (2020).
  * Wu and Zhuang [2021] J. Wu and Q. Zhuang, Continuous-variable error correction for general Gaussian noises, Phys. Rev. Appl. 15, 034073 (2021).
  * Wu  _et al._ [2020] B.-H. Wu, R. N. Alexander, S. Liu, and Z. Zhang, Quantum computing with multidimensional continuous-variable cluster states in a scalable photonic platform, Phys. Rev. Res. 2, 023138 (2020).
  * Hillmann  _et al._ [2020] T. Hillmann, F. Quijandría, G. Johansson, A. Ferraro, S. Gasparinetti, and G. Ferrini, Universal gate set for continuous-variable quantum computation with microwave circuits, Phys. Rev. Lett. 125, 160501 (2020).
  * Gao  _et al._ [2019] Y. Y. Gao, B. J. Lester, K. S. Chou, L. Frunzio, M. H. Devoret, L. Jiang, S. Girvin, and R. J. Schoelkopf, Entanglement of bosonic modes through an engineered exchange interaction, Nature 566, 509 (2019).
  * Lloyd [2008] S. Lloyd, Enhanced sensitivity of photodetection via quantum illumination, Science 321, 1463 (2008).
  * Simon  _et al._ [1994] R. Simon, N. Mukunda, and B. Dutta, Quantum-noise matrix for multimode systems: U(n) invariance, squeezing, and normal forms, Phys. Rev. A 49, 1567 (1994).
  * Williamson [1936] J. Williamson, On the algebraic problem concerning the normal forms of linear dynamical systems, American Journal of Mathematics 58, 141 (1936).

##  VII Analytical methods, experimental procedures, and entanglement analysis

###  VII.1 Details of theoretical description

The Hamiltonian of JPA system is given by

$\displaystyle{\hat{H}}_{\text{sys}}=\hbar\omega_{r}\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\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})$ (19) $\displaystyle+\hbar K(\hat{a}+\hat{a}^{\dagger})^{4},$

where $\hat{a}$($\hat{a}^{\dagger}$) is the annihilation (creation) operators
for cavity photons, $\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 $\omega_{d}$, $d=\left\\{1,\dots,p\right\\}$, we define the
detuning between the half pump frequency and the resonator frequency:
$\Delta_{r}=\omega_{r}-\frac{\omega_{\Sigma}}{2}$,
$\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
$\Delta_{d}=\omega_{d}-\omega_{\Sigma},d=\left\\{1,\dots,p\right\\}$. By
applying the rotating wave approximation in the frame $\omega_{\Sigma}/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
Hamiltonian

$\displaystyle{H}_{\text{sys,rwa}}(t)={\hbar}\Delta_{r}\tilde{a}^{\dagger}\tilde{a}$
$\displaystyle+\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})$ (20)
$\displaystyle+6{\hbar}{K}\tilde{a}^{\dagger}\tilde{a}^{\dagger}\tilde{a}\tilde{a}.$

As usual, the bosonic commutation relationships are valid for the cavity modes
$\left[\tilde{a},\tilde{a}^{\dagger}\right]=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

$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),$
(21)

where creation and annihilation operators $\tilde{b}^{\dagger}$ and
$\tilde{b}$ refer to modes in the transmission line, and $\kappa$ denotes the
coupling rate. The Hamiltonian related to the linear dissipation port

$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),$
(22)

where $\tilde{c}^{\dagger}$ and $\tilde{c}$ describe creation and annihilation
of thermal bath modes and the rate $\gamma$ represent the coupling of cavity
modes to the linear dissipation port. The transmission line and bath operators
obey the commutation relations

$\left[\tilde{b}(\omega),\tilde{b}^{\dagger}(\omega^{\prime})\right]=\left[\tilde{c}(\omega),\tilde{c}^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime})$
(23)

and

$\displaystyle\left[\tilde{b}(\omega),\tilde{c}^{\dagger}(\omega^{\prime})\right]=$
$\displaystyle\left[\tilde{c}(\omega),\tilde{b}^{\dagger}(\omega^{\prime})\right]=\left[\tilde{b}(\omega),\tilde{b}(\omega^{\prime})\right]$
(24)
$\displaystyle=\left[\tilde{c}(\omega),\tilde{c}(\omega^{\prime})\right]=0.$

The total Hamiltonian can be conveniently written as a sum of the separate
parts given above:

${H}(t)=H_{\text{sys,rwa}}(t)+H_{\text{sig}}+H_{\text{loss}}$ (25)

For further analysis and for our simulations, we use the Quantum Langevin
Equation (QLE) for the cavity operator $\tilde{a}(t)$:

$\displaystyle\dot{\tilde{a}}(t)=(-i\Delta_{r}$
$\displaystyle-\frac{\kappa+\gamma}{2})\tilde{a}-i\sum_{d=1}^{p}\alpha_{d}e^{i\Delta_{d}t}\tilde{a}^{\dagger}$
(26)
$\displaystyle+\sqrt{\kappa}\tilde{b}_{in}+\sqrt{\gamma}\tilde{c}_{in}-12i{K}\tilde{a}^{\dagger}\tilde{a}\tilde{a},$

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:

$\tilde{b}_{out}(t)=\tilde{b}_{in}(t)-\sqrt{\kappa}\tilde{a}(t).$ (27)

Eqs. (26) and (27) are used in our numerical simulations with Matlab ODE45
solver.

###  VII.2 Full inseparability

Assuming that the microwave fields produced by the JPA below the critical
pumping threshold are Gaussian [16], 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
$\tilde{x}_{i}=\frac{\tilde{a}_{i}+\tilde{a}_{i}^{\dagger}}{2}$ and
$\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
$\tilde{r}=\left(\tilde{x}_{1},\tilde{p}_{1},\dots\tilde{x}_{N},\tilde{p}_{N}\right)^{T}$.
The commutation relations can be written down in a skew-symmetric, block-
diagonal matrix form [33]:

$\left[\tilde{r}_{i},\tilde{r}_{j}\right]=\frac{i}{2}\Omega_{ij}\text{ and
}\mathbf{\Omega}=\bigotimes_{i=1}^{N}\begin{bmatrix}0&1\\\ -1&0\end{bmatrix}.$
(28)

The covariance matrix $\mathbf{V}$ is given by elements
$V_{ij}=\frac{1}{2}\left<\Delta\tilde{r}_{i}\Delta\tilde{r}_{j}+\Delta\tilde{r}_{j}\Delta\tilde{r}_{i}\right>-\left<\Delta\tilde{r}_{i}\right>\left<\Delta\tilde{r}_{j}\right>$
where we have defined standard error
$\Delta\tilde{r}_{i}=\tilde{r}_{i}-\left<\tilde{r}_{i}\right>$ and
$\left<\tilde{r}_{i}\right>=\text{tr}\left(\tilde{r}_{i}\hat{\rho}\right)$.
The uncertainty principle requires that

$\mathbf{V}+\frac{i}{4}\mathbf{\Omega}\geq 0$ (29)

applies for a physical covariance matrix.

For verification of entanglement, we may investigate a modified equation

$\mathbf{V}^{\prime}+\frac{i}{4}\mathbf{\Omega}\geq 0,$ (30)

where
$\mathbf{V}^{\prime}_{k}=\boldsymbol{\lambda}_{k}\mathbf{V}\boldsymbol{\lambda}_{k}$,
$\boldsymbol{\lambda}_{k}$ is diagonal matrix with ones entries, except of
that related to $k$-th mode, with value of $-1$. For example, transformation
with $\boldsymbol{\lambda}_{k\equiv N}=\text{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:
$\mathbf{V}^{\prime}_{k}\geq\frac{i}{4}\mathbf{\Omega}$. In Ref. 30, 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.,

$\mathbf{\Omega}=\mathbf{P}^{T}\mathbf{\Omega}\mathbf{P}.$ (31)

Symplectic transformations on a 2N-dimensional phase-space form the real
symplectic group denoted as $Sp(2N;\textbf{R})$, which is a proper subgroup of
the special linear group of 2N $\times$ 2N matrices [58]. By utilizing
Williamson’s theorem [59], any covariance matrix can be expressed in the
Williamson normal form:

$\mathbf{\tilde{V}}_{k}=\mathbf{P}^{T}\mathbf{V}^{\prime}_{k}\mathbf{P}$ (32)

where $\mathbf{\tilde{V}}_{k}$ is a 2N-dimensional diagonal matrix consisting
of the symplectic eigenvalues, $\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 $\tilde{\nu}_{k}\geq 1/4$ fulfilled
for $\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:
${V}^{*}_{ij}=2\left<\Delta\tilde{r}_{i}\Delta\tilde{r}_{j}+\Delta\tilde{r}_{j}\Delta\tilde{r}_{i}\right>$.
Consequently, for evidence of ’fully inseparable’ states, we need to find
minimum symplectic eigenvalues with $\tilde{\nu}_{k}^{*}<1$ for each partial
transposition $k$.

###  VII.3 System gain calibration

Figure 9: Gain calibration using linear temperature dependence of the measured
thermal noise spectral density of a 50 Ohm terminator measured as a function
of $T_{s}$, the source temperature. The average total gain
$\left<G_{\Sigma,i}\right>=94.4\pm 0.2$ dB over the cavity resonance is
obtained from the linear fit (in red) to the data. This gain value
$\left<G_{\Sigma,i}\right>$ also includes frequency mixing
losses/amplification in the signal analyzer circuit part. The term
$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(\frac{hf_{i}}{2k_{b}T_{\text{preamp}}})$ sets the background for
the diagonal elements in the covariance matrix $4\mathbf{V}_{\text{off}}$.

Our system gain calibration procedure consists of a measurement of Johnson-
Nyquist noise spectral density emitted by a 50 $\Omega$ 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_{\text{HEMT}}$, and the noise of the room-temperature amplifiers
$T_{\text{RT}}$ multiplied by the system gain
$G_{\Sigma,i}=G_{\text{HEMT}_{i}}G_{\text{RT}_{i}}$:

$\frac{\left<I_{i}^{2}+Q_{i}^{2}\right>}{Z_{0}\Delta
f_{i}}=k_{b}G_{\Sigma,i}\left(T_{i}+T_{\text{HEMT}}+\frac{T_{\text{RT}}}{G_{\text{HEMT}}}\right).$
(33)

Here $i$ refers to the frequency of the spectral mode and $\Delta f_{i}$
refers to the bandwidth of the detection of quadratures $I_{i}$ and $Q_{i}$.
The total gain $G_{\Sigma,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 $\left<G_{\Sigma,i}\right>=94.4\pm
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(\hbar\omega/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\\{\tilde{\nu}_{k}^{*}\\}$=$\min\\{\tilde{\nu}_{k}^{*}\\}\pm
0.018$ for each partial bipartition. Random variations of the system
parameters were reduced by averaging the outcome by ten to twenty times.

###  VII.4 System parameter fitting

In order to determine coupling rates $\gamma$ and $\kappa$ 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 $\Phi_{DC}=0.383\,\Phi_{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 $\kappa$ and $\gamma$,
respectively. By fitting the measured resonance curve to the analytical
solution of the QLE ($\tilde{b}_{out}(\omega)/\tilde{b}_{in}(\omega)$),
derived for the linear case without any pump drive, we obtain the coupling
coefficients $\frac{\kappa}{2\pi}=4.44$ MHz and $\frac{\gamma}{2\pi}=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\tilde{a}^{\dagger}\tilde{a}\tilde{a}$ into account:

$\frac{\tilde{b}_{out}(\omega)}{\tilde{b}_{in}(\omega)}=1-\frac{\kappa}{(-i(\omega-\omega_{r})+\frac{\kappa+\gamma}{2})}.$
(34)

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(\omega_{\text{probe}}-\omega_{r},A)$ (Fig. 10) in the
cavity at large pump amplitudes, we obtain an estimate $K=0.54\omega_{r}$ for
the Kerr constant.

(a) 1

(b) 1

Figure 10: Plots (a) and (b) show gain coefficient measurement and QLE’s
simulation results (Eq. (26)) which used for fitting coupling and Kerr
constants. Vertical axis represents normalized pump amplitude
$A=\frac{\alpha}{\kappa+\gamma}$. Horizontal axis represents detuning between
probe signal frequency and resonance frequency
$\frac{\omega_{\text{probe}}-\omega_{r}}{2\pi}$. Pumping carried out in
degenerate mode, $\omega_{d}=2\omega_{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 $\kappa$.

###  VII.5 Cavity phase response

Optimal value of GME criterion, which governed by ”symmetric” covariance
matrix view in tripartite mode case, can be obtained with
$\\{\frac{\pi}{2};\frac{\pi}{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.

Figure 11: Cavity phase response given by fitted experimental parameters
$\frac{\kappa}{2\pi}=4.44$ MHz and $\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
$\Delta\varphi_{1}=\frac{\pi}{2}-\frac{\pi}{4}=\frac{\pi}{4}$ and
$\Delta\varphi_{2}=\frac{\pi}{2}+\frac{\pi}{4}=\frac{3\pi}{4}$ with
corresponding phase shift between pump tones $\frac{\pi}{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.

###  VII.6 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 $\tilde{a}$ as a vector of spectral modes:

$\tilde{a}=\\{{\tilde{a}_{1},\tilde{a}_{2},\tilde{a}_{3},\tilde{a}^{\dagger}_{1},\tilde{a}^{\dagger}_{2},\tilde{a}^{\dagger}_{3}\\}}^{T};$
(35)

where the creation $\tilde{a}_{i}^{\dagger}=\tilde{a}_{i}^{\dagger}(t)$ and
annihilation $\tilde{a}_{i}=\tilde{a}_{i}(t)$ operators are time-dependent.
After Fourier transform,

$\tilde{a}(\omega)=\\{{\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)\\}}^{T}.$
(36)

We define our spectral modes $a_{i}$ as
$\\{(-\frac{3\Delta}{2},-\frac{\Delta}{2});(-\frac{\Delta}{2},\frac{\Delta}{2});(\frac{\Delta}{2},\frac{3\Delta}{2})\\}$
according pump tone positions $\\{-\Delta,\Delta\\}$ (see Fig. 1). Similarly,
we define for the input and output modes $\hat{b}_{\text{in}/\text{out}}$:

$\displaystyle\tilde{b}_{\text{in/out}}=$
$\displaystyle\\{{\tilde{b}_{\text{in}1/\text{out}1},\tilde{b}_{\text{in}2/\text{out}2},\tilde{b}_{\text{in}3/\text{out}3}},$
(37)
$\displaystyle\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}\\}^{T}.$

The commutation relationships for the case of $N$ modes can be conveniently
expressed in matrix form. We use the common convention for
$\left[\tilde{a}_{i},\tilde{a}_{j}\right]$ from Ref. 58.

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 $\Delta\ll\kappa$; we also neglect internal dissipation
expressed by loss rate $\gamma$. 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 [23], we denote
$\tilde{a}_{i}(\omega)=\int\tilde{a}_{i}(t)e^{i\omega 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 $\tilde{a}^{\dagger}_{1,2}$:
$\int\tilde{a}^{\dagger}_{1}(t)e^{i\omega
t}e^{i\Delta_{1}t}dt=\tilde{a}^{\dagger}_{2}(-\omega)$;
$\int\tilde{a}^{\dagger}_{2}(t)e^{i\omega
t}e^{i\Delta_{2}t}dt=\tilde{a}^{\dagger}_{1}(-\omega)$, while for
$\tilde{a}^{\dagger}_{2,3}$: $\int\tilde{a}^{\dagger}_{2}(t)e^{i\omega
t}e^{-i\Delta_{2}t}dt=\tilde{a}^{\dagger}_{3}(-\omega)$;
$\int\tilde{a}^{\dagger}_{3}(t)e^{i\omega
t}e^{-i\Delta_{3}t}dt=\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,

$\displaystyle(-i(\omega-\Delta_{r})+$
$\displaystyle\frac{\kappa}{2})\tilde{a}(\omega)+i\alpha(\int\tilde{a}^{\dagger}(t)e^{i\omega
t}e^{-i\Delta_{d}t}dt+$ (38)
$\displaystyle\int\tilde{a}^{\dagger}(t)e^{i\omega
t}e^{i\Delta_{d}t}dt)=\sqrt{\kappa}\tilde{b}_{in}(\omega),$

the QLE yields the following system of linear equations:

$\displaystyle\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)$
$\displaystyle\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))$
$\displaystyle\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)$
$\displaystyle\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)$
$\displaystyle\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))$
$\displaystyle\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).$
(39)

We cast Eq. (39) into matrix form:

$\mathbf{M}\tilde{a}(\omega)=\sqrt{\kappa}\tilde{b}_{\text{in}}(\omega)$ (40)
$\mathbf{M}=\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}$ (41)

where $c_{1}=-i(\omega-\Delta_{r})+\frac{\kappa}{2}$ and
$c_{2}=-i(\omega+\Delta_{r})+\frac{\kappa}{2}$ Solving for the inverse of
$\hat{M}$ and using Eq. (27), we obtain

$\tilde{b}_{\text{out}}(\omega)=(\mathbf{I}-\kappa\mathbf{M}^{-1})\tilde{b}_{\text{in}}(\omega)$
(42)

for the outgoing radiation $\tilde{b}_{\text{out}}(\omega)$ in terms of
incoming waves $\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 $\tilde{a}$ with
equation
$\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 $\mathbf{I}$ and the multiplication factor $\kappa$ do not change the
final structure.

Assuming that the pump amplitude $\alpha$ is a real number and $c_{1}=c_{2}=c$
(zero detuning case), we have

$\displaystyle\mathbf{M}^{-1}=\frac{1}{c^{2}-2\alpha^{2}}\cdot$
$\displaystyle\begin{bmatrix}c-\frac{\alpha^{2}}{c}&0&\mathbf{\frac{\boldsymbol{\alpha}^{2}}{c}}&0&-i\alpha&0\\\
0&c&0&-i\alpha&0&-i\alpha\\\
\mathbf{\frac{\boldsymbol{\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{\boldsymbol{\alpha}^{2}}{c}}\\\
i\alpha&0&i\alpha&0&c&0\\\
0&i\alpha&0&\mathbf{\frac{\boldsymbol{\alpha}^{2}}{c}}&0&c-\frac{\alpha^{2}}{c}\end{bmatrix}$
(43)

This allows us to draw the generalized $\mathscr{\tilde{H}}$-graph for
describing the parametric interaction between the spectral modes, Fig. 2. The
off-diagonal beamsplitter elements proportional to $\alpha^{2}$ are set in
bold in Eq. (43).

Still, we want to construct the parametric interaction matrix
$\mathbf{S}^{-1}$ for quadrature vector operator $\tilde{r}$. Using a linear
operator matrix $\mathbf{K}$ to implement a change of basis

$\displaystyle\mathbf{K}=\frac{1}{2}\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},$ (44)

we obtain by a canonical transformation
$\mathbf{S}^{-1}=\sqrt{\kappa}\mathbf{K}\mathbf{M}^{-1}\mathbf{K}^{-1}$:

$\displaystyle\mathbf{S}^{-1}=\frac{\sqrt{\kappa}}{c^{2}-2\alpha^{2}}\cdot$
$\displaystyle\begin{bmatrix}c+\frac{\alpha^{2}}{c}&0&0&-\alpha&\mathbf{\frac{\boldsymbol{\alpha}^{2}}{c}}&0\\\
0&c+\frac{\alpha^{2}}{c}&-\alpha&0&0&\mathbf{\frac{\boldsymbol{\alpha}^{2}}{c}}\\\
0&-\alpha&c&0&0&-\alpha\\\ -\alpha&0&0&c&-\alpha&0\\\
\mathbf{\frac{\boldsymbol{\alpha}^{2}}{c}}&0&0&-\alpha&c+\frac{\alpha^{2}}{c}&0\\\
0&\mathbf{\frac{\boldsymbol{\alpha}^{2}}{c}}&-\alpha&0&0&c+\frac{\alpha^{2}}{c}\end{bmatrix}.$
(45)

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,
$\tilde{b}_{\text{in}}$ has a Gaussian covariance matrix of the form
$\mathbf{V}_{\text{in}}=\frac{1}{4}\mathbf{I}$. Consequently, the covariance
matrix of the cavity spectral modes $\tilde{a}_{i}$ can be represented as

$\mathbf{V}_{\text{a}}=\mathbf{S}^{-1}\mathbf{V}_{\text{in}}(\mathbf{S}^{-1})^{T}$
(46)

or, equivalently, for output modes $\tilde{b}_{\text{out}}$:
$\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 $\mathbf{V}_{\text{a}}$ and $\mathbf{V}_{\text{out}}$ can be
employed for studying the structure of parametric interactions between the
quadratures, because input-output 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 $\mathbf{M}$ due to
a phase shift is illustrated in Eq. (47), in which the phase of the first pump
has been rotated by $e^{\frac{i\pi}{2}}$.

$\hat{M}=\begin{bmatrix}c&0&0&0&-\boldsymbol{\alpha}&0\\\
0&c&0&-\boldsymbol{\alpha}&0&i\alpha\\\ 0&0&c&0&i\alpha&0\\\
0&-\boldsymbol{\alpha}&0&c&0&0\\\ -\boldsymbol{\alpha}&0&-i\alpha&0&c&0\\\
0&-i\alpha&0&0&0&c\end{bmatrix}.$ (47)

The elements affected by the rotation are indicated in bold in the matrix. The
elements in bold face indicate coupling between modes
$\tilde{a}_{1}(\omega)\leftrightarrow\tilde{a}_{2}(\omega)$ while the other
off-diagonal elements indicate squeezing across
$\tilde{a}_{2}(\omega)\leftrightarrow\tilde{a}_{3}(\omega)$. Note, that phase
rotation operates in opposite direction on rows related to
$\tilde{b}_{\text{in}}(\omega)$ and $\tilde{b}_{\text{in}}^{\dagger}(\omega)$.

The inversion of the rotated matrix $\mathbf{M}$ yields for the parametric
interaction matrix, where all the beamsplitter elements (in bold) have
acquired a $\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 $\mathbf{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.

$\displaystyle\mathbf{M}^{-1}=\frac{1}{c^{2}-2\alpha^{2}}\cdot$
$\displaystyle\begin{bmatrix}c-\frac{\alpha^{2}}{c}&0&\mathbf{\frac{\boldsymbol{i\alpha}^{2}}{c}}&0&\alpha&0\\\
0&c&0&\alpha&0&-i\alpha\\\
-\mathbf{\frac{\boldsymbol{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{\boldsymbol{i\alpha}^{2}}{c}}\\\
\alpha&0&i\alpha&0&c&0\\\
0&i\alpha&0&\mathbf{\frac{\boldsymbol{i\alpha}^{2}}{c}}&0&c-\frac{\alpha^{2}}{c}\end{bmatrix}$
(48)

The covariance matrix $\mathbf{V}_{\text{a}}$ obtained from matrix
$\mathbf{M}$ in Eq. (43) with zero pump phase shifts is given in Eq. (49). The
corresponding covariance matrix for $\pi/2$ phase rotation in the first pump
is displayed in Eq. (50). The matrix $\mathbf{V}_{\text{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^{\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 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.

$\mathbf{V}_{\text{a}}=\frac{\kappa}{4(c^{2}-2\alpha^{2})^{2}}\cdot\\\
\begin{bmatrix}-\alpha^{2}+\frac{2\alpha^{4}}{c^{2}}+c^{2}&0&0&-2\alpha
c&\boldsymbol{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&\boldsymbol{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\\\
\boldsymbol{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&\boldsymbol{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}$ (49)

$\mathbf{V}_{\text{a}}=\frac{\kappa}{4(c^{2}-2\alpha^{2})^{2}}\cdot\\\
\begin{bmatrix}-\alpha^{2}+\frac{2\alpha^{4}}{c^{2}}+c^{2}&0&\boldsymbol{2\alpha
c}&0&0&\boldsymbol{-3\alpha^{2}+\frac{2\alpha^{4}}{c^{2}}}\\\
0&-\alpha^{2}+\frac{2\alpha^{4}}{c^{2}}+c^{2}&0&\boldsymbol{-2\alpha
c}&\boldsymbol{3\alpha^{2}-\frac{2\alpha^{4}}{c^{2}}}&0\\\ \boldsymbol{2\alpha
c}&0&2\alpha^{2}+c^{2}&0&0&-2\alpha c\\\ 0&\boldsymbol{-2\alpha
c}&0&2\alpha^{2}+c^{2}&-2\alpha c&0\\\
0&\boldsymbol{3\alpha^{2}-\frac{2\alpha^{4}}{c^{2}}}&0&-2\alpha
c&-\alpha^{2}+\frac{2\alpha^{4}}{c^{2}}+c^{2}&0\\\
\boldsymbol{-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}$ (50)

$\mathbf{V}_{\text{a}}=\frac{\kappa}{4(c^{2}-2\alpha^{2})^{2}}\cdot\\\
\begin{bmatrix}-\alpha^{2}+\frac{2\alpha^{4}}{c^{2}}+c^{2}&0&2\alpha
c&0&\boldsymbol{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&\boldsymbol{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\\\
\boldsymbol{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&\boldsymbol{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}.$ (51)

4 mode case. We treat the case with three pump tones ($p=3$) in the same way
as we did with two pumps above. By taking
$\omega_{\Sigma}=\frac{\sum_{d=1}^{p}\omega_{d}}{p}=2\omega_{r}$, we have a
configuration of pump tones $\\{{-2\Delta,0,2\Delta\\}}$ with respect to
$2\omega_{r}$. We define our spectral modes around $\omega_{r}$ as
$\\{(-2\Delta,-\Delta);(-\Delta,0);(0,\Delta);(\Delta,2\Delta)\\}$
$\widehat{=}\\{\tilde{a}_{1}(\omega);\tilde{a}_{2}(\omega);\tilde{a}_{3}(\omega);\tilde{a}_{4}(\omega)\\}$.
Individual pump tones at $\omega_{1}$ and $\omega_{3}$ create TMS states
between neighboring spectral modes as in the 2 pump case above. However, the
middle pump creates TMS correlations between two pairs of spectral modes:
$\tilde{a}_{1}(\omega)\leftrightarrow\tilde{a}^{\dagger}_{4}(-\omega)$ and
$\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 $\mathscr{\tilde{H}}$-graph states (see
Fig. 2).

$\displaystyle\mathbf{M}=\begin{bmatrix}c&0&0&0&0&-\alpha&0&\boldsymbol{-\alpha}\\\
0&c&0&0&-\alpha&0&\boldsymbol{-\alpha}&0\\\
0&0&c&0&0&\boldsymbol{-\alpha}&0&-\alpha\\\
0&0&0&c&\boldsymbol{-\alpha}&0&-\alpha&0\\\
0&-\alpha&0&\boldsymbol{-\alpha}&c&0&0&0\\\
-\alpha&0&\boldsymbol{-\alpha}&0&0&c&0&0\\\
0&\boldsymbol{-\alpha}&0&-\alpha&0&0&c&0\\\
\boldsymbol{-\alpha}&0&-\alpha&0&0&0&0&c\\\ \end{bmatrix},$ (52)

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 $\mathbf{M}$ (Eq. (52)), where all $c$ coefficient are
equal since we have assumed
$\Delta_{r}=\omega_{r}-\frac{\omega_{\Sigma}}{2}=0$. The signs of $\alpha$’s
are governed by the choice of pump phases as $\\{{\alpha
e^{\frac{i\pi}{2}};\alpha e^{\frac{i\pi}{2}};\alpha e^{\frac{i\pi}{2}}\\}}$.
The correlations produced by the pump at $\omega_{2}=2\omega_{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 $\mathbf{M}^{-1}$ reveals the beamsplitter correlations
between $\tilde{a}_{1}(\omega)\leftrightarrow\tilde{a}_{3}(\omega)$ and
$\tilde{a}_{2}(\omega)\leftrightarrow\tilde{a}_{4}(\omega)$:

$\displaystyle\mathbf{M}^{-1}=\frac{1}{(c^{2}-4\alpha^{2})}\cdot$
$\displaystyle\begin{bmatrix}\begin{array}[]{cccccccc}c-\frac{2\alpha^{2}}{c}&0&\boldsymbol{\frac{2\alpha^{2}}{c}}&0&0&\alpha&0&\alpha\\\
0&c-\frac{2\alpha^{2}}{c}&0&\boldsymbol{\frac{2\alpha^{2}}{c}}&\alpha&0&\alpha&0\\\
\boldsymbol{\frac{2\alpha^{2}}{c}}&0&c-\frac{2\alpha^{2}}{c}&0&0&\alpha&0&\alpha\\\
0&\boldsymbol{\frac{2\alpha^{2}}{c}}&0&c-\frac{2\alpha^{2}}{c}&\alpha&0&\alpha&0\\\
0&\alpha&0&\alpha&c-\frac{2\alpha^{2}}{c}&0&\boldsymbol{\frac{2\alpha^{2}}{c}}&0\\\
\alpha&0&\alpha&0&0&c-\frac{2\alpha^{2}}{c}&0&\boldsymbol{\frac{2\alpha^{2}}{c}}\\\
0&\alpha&0&\alpha&\boldsymbol{\frac{2\alpha^{2}}{c}}&0&c-\frac{2\alpha^{2}}{c}&0\\\
\alpha&0&\alpha&0&0&\boldsymbol{\frac{2\alpha^{2}}{c}}&0&c-\frac{2\alpha^{2}}{c}\\\
\end{array}\end{bmatrix}$ (53)

The beamsplitter correlations are indicated in bold in this matrix
$\mathbf{M}^{-1}$. We see that there are two sequences of pump transformations
that yield BS correlations between modes
$\tilde{a}_{1}(\omega)\leftrightarrow\tilde{a}_{2}(\omega)$ and
$\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
$\mathscr{\tilde{H}}$-graphs is illustrated in Fig. 2.

The beamsplitter correlations are prominent also in the covariance matrix
$\mathbf{V}_{\text{a}}$ (see 46):

$\displaystyle\mathbf{V}_{\text{a}}=\frac{\kappa}{4(c^{2}-4\alpha^{2})^{2}}\cdot$
(54)

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
$\left(\tilde{x}_{1},\tilde{p}_{1},\dots\tilde{x}_{N},\tilde{p}_{N}\right)^{T}$.
So the BS correlations are between quadratures of
$\tilde{a}_{1}(\omega)\leftrightarrow\tilde{a}_{3}(\omega)$ and
$\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 $\\{{\alpha e^{\frac{-i\pi}{2}};\alpha e^{\frac{i\pi}{2}};\alpha
e^{\frac{i\pi}{2}}\\}}$ we are able to flip the sign of one minor diagonal
indicated by bold font in Eq. (55).

$\mathbf{M}=\begin{bmatrix}c&0&0&0&0&\boldsymbol{\alpha}&0&-\alpha\\\
0&c&0&0&\boldsymbol{\alpha}&0&-\alpha&0\\\ 0&0&c&0&0&-\alpha&0&-\alpha\\\
0&0&0&c&-\alpha&0&-\alpha&0\\\ 0&\boldsymbol{\alpha}&0&-\alpha&c&0&0&0\\\
\boldsymbol{\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}$ (55)

Interestingly, this choice of phases leads to full cancellation of the
beamsplitter correlation terms. This is seen in the structure of the inverse
matrix:

$\displaystyle\mathbf{M}^{-1}=\frac{1}{(c^{2}-2\alpha^{2})}\cdot$
$\displaystyle\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},$ (56)

which does not have any elements proportional to $\alpha^{2}$. Without any
beamsplitter correlations, Eq. (56) indicates a clear connection to a square-
lattice $\mathscr{\tilde{H}}$-graph.

Finally, the corresponding covariance matrix $\mathbf{V}_{\text{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^{\frac{i\pi}{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\times 8$ covariance matrix. In
particular, the influence of the beamsplitter correlations can be eliminated
from $\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. 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.

$\displaystyle\mathbf{V}_{\text{a}}=\frac{\kappa}{4(c^{2}-2\alpha^{2})^{2}}\cdot$
(57)