Squeezing and quantum state engineering with Josephson traveling wave amplifiers

We develop a quantum theory describing the input-output properties of Josephson traveling wave parametric amplifiers. This allows us to show how such a device can be used as a source of nonclassical radiation, and how dispersion engineering can be used to tailor gain profiles and squeezing spectra with attractive properties, ranging from genuinely broadband spectra to"squeezing combs"consisting of a number of discrete entangled quasimodes. The device's output field can be used to generate a multi-mode squeezed bath--a powerful resource for dissipative quantum state preparation. In particular, we show how it can be used to generate continuous variable cluster states that are universal for measurement based quantum computing. The favourable scaling properties of the preparation scheme makes it a promising path towards continuous variable quantum computing in the microwave regime.

We develop a quantum theory describing the input-output properties of Josephson traveling wave parametric amplifiers. This allows us to show how such a device can be used as a source of nonclassical radiation, and how dispersion engineering can be used to tailor gain profiles and squeezing spectra with attractive properties, ranging from genuinely broadband spectra to "squeezing combs" consisting of a number of discrete entangled quasimodes. The device's output field can be used to generate a multi-mode squeezed bath-a powerful resource for dissipative quantum state preparation. In particular, we show how it can be used to generate continuous variable cluster states that are universal for measurement based quantum computing. The favourable scaling properties of the preparation scheme makes it a promising path towards continuous variable quantum computing in the microwave regime.

I. INTRODUCTION
Superconducting microwave circuits can be used to behave as artificial atoms in engineered electromagnetic environments, where strong light-matter interaction is achieved by confining the electromagnetic field in microwave resonators [1] or one-dimensional waveguides [2]. This is in close analogy, respectively, with cavity and waveguide quantum electrodynamics with real atoms [3][4][5]. The flexibility offered by microwave engineering also allows experimentalists to go beyond the limits of conventional quantum optics in many ways. Examples include realizing light-matter coupling strengths that are unachievable with real atoms [6], and using nonlinear microwave resonators to simulate relativistic quantum effects [7].
A recent advancement in microwave quantum optics is the bottom-up design of nonlinear, one-dimensional metamaterials with strong photon-photon interactions and engineered dispersion relations [8][9][10]. The nonlinearity is provided by Josephson junctions embedded in a transmission line, with photon-photon interactions activated by a strong pump tone through a parametric fourwave mixing process. These devices have been dubbed Josephson traveling Wave Parametric Amplifiers (JTW-PAs) [9], and are analogous to one-dimensional χ (3) nonlinear crystals [11].
The development of JTWPAs is motivated by their potential use as amplifiers for readout of superconducting qubits. The extremely high measurement fidelity necessary for fault tolerant quantum computing requires phase preserving amplifiers with added noise near the fundamental quantum limit [12,13]. The JTWPA design is in this respect a very promising candidate. Early experimental realizations have shown impressive performance, and are already expected to have sufficient dynamic range * Electronic address: arne.loehre.grimsmo@usherbrooke.ca and bandwidth to read out several tens of superconducting qubits with a single device [9,10]. The key advantage to the JTWPA design is the operational bandwidth which is in the GHz range. This is in contrast to other nearquantum limited microwave amplifiers based on resonant cavity interactions, which typically have bandwidths of tens of MHz [14][15][16].
An amplifier operating near the quantum limit is, however, very different from a classical amplifier. Quantumlimited phase preserving amplification implies the presence of entanglement between the amplified signal and an "idler" signal in a two-mode squeezed state [12,17]. This motivates an alternative viewpoint on the JTWPA: Besides using the device to amplify a signal of interest, one can also view it as a broadband source of nonclassical radiation.
In this paper we show how the inherent flexibility in the bottom-up JTWPA construction allows for designing broadband squeezing spectra with attractive properties. In particular, we show how to tailor the squeezing spectrum such that some frequency ranges are unaffected by the nonlinear interaction. This is useful for example to avoid unwanted quantum heating [18] of systems placed at the JTWPA output.
We subsequently demonstrate how the JTWPA can be used as a resource for dissipative quantum state preparation, including resource states for universal measurement based quantum computing [19,20]. Dissipative quantum state preparation has over the last years emerged as an exciting alternative to preparation of entangled states using coherent Hamiltonian [21] or gate-based methods [22]. It has been shown that universal quantum computing can be achieved through dissipative processes alone [23], and in a similar vein that highly correlated states such as stabilizer states and projected entangled pair states can be created as stable steady states of dissipative processes [23,24]. The early theoretical proposals in Refs. [23,24], however, involve many-body dissipative interactions that are hard to realize in practice. As a consequence, searching for simpler schemes for dis-sipative state preparation that can be implemented in present day experiments has become an active area of research [25][26][27][28][29][30].
We show that broadband squeezed radiation, such as the radiation emitted by a JTWPA, is a particularly potent resource for dissipative quantum state preparation. The emitted radiation generates a broadband squeezed bath which can be used to cool quantum systems placed at the source's output into entangled states. This contrasts the amplifier mode of operation, where the systems of interest are located at the amplifier's input, and do not see the nonclassical radiation emitted by the device. We show that by engineering such a squeezed bath one can produce pairs of entangled qubits, as well as continuous variable cluster states that are universal for measurement based quantum computing. The preparation schemes are simple, requiring no Hamiltonian interactions or complicated reservoir engineering. For the case of JTWPAs as squeezing sources, the large bandwidth furthermore makes the process very hardware efficient, making this an attractive avenue for measurement based quantum computing with microwaves.
The purely dissipative nature of the preparation process distinguishes our proposal from similar approaches for generating cluster states in the optical regime [31][32][33][34][35][36]. A distinct advantage of a dissipative scheme is that it relaxes constraints on locality, which might allow for a more modular architecture that avoids spurious interactions and increases scalability [37].
Although we focus on JTPWAs as squeezing sources in this work, due to their design flexibility and large bandwidth, we emphasize that the dissipative quantum state preparation schemes we develop are relevant for any type of broadband squeezing source that can be integrated with coherent quantum systems, such as other types of traveling wave amplifiers [38,39], impedance engineered Josephson parametric amplifiers [40], squeezing sources based on reservoir engineering [41], or even the nonclassical radiation emitted by an ac-biased tunnel junction [42,43].

II. ASYMPTOTIC INPUT-OUTPUT THEORY
To describe the JTWPA's squeezing properties, we first need a quantized theory of its dynamics. Classical treatments of a JTWPA are presented in Refs. [8,10,44]. In the following we give a quantized Hamiltonian treatment of the nonlinear dynamics, taking into account dispersion and the continuum nature of the electromagnetic field in the waveguide.
The resulting theory is in general difficult to treat analytically due to time-ordering effects [45,46], and we therefore take a perturbative approach treating the nonlinearity to first order. An input-outpu relation linking the field entering the JTWPA to the emitted output field is derived in the usual scattering limit where the initial and final times of the problem are taken to minus and (a) plus infinity, respectively [45][46][47][48][49].
Equations of motion similar to those we derive here have been used previously by Caves and Crouch in a study of wideband traveling wave amplifier [50], where they were taken as operator versions of macroscopic Maxwell's equations for a nonlinear, homogeneous and dispersionless medium [51]. We here rigorously justify these equations by deriving them from a microscopic theory, taking into account the finite extent of the nonlinearity as well as dispersion, the latter stemming from both junctions' plasma oscillations, non-linear phasemodulation and engineered bandgaps in the medium.
The device we consider in this paper is depicted in Fig. 1. It consists of a series of identical coupled Josephson junctions with Josephson energies E J and junction capacitances C J . Each junction is coupled to ground by a passive, dissipationless element with impedance Z(ω), which is left arbitrary for now. By engineering Z(ω) one can modify the dispersion relation of waves propagating through the device as shown in Ref. [8]. We show in Sec. III how this can be used to tailor the squeezing properties of the output field leaving the device.
Realistic JTWPAs have several thousand junctions with a unit cell distance much smaller than the relevant wavelengths [9,10]. One can therefore approximate the device with a continuum description (formally taking the unit cell distance, a, to zero). We furthermore assume that the JTWPA is coupled to identical, semi-infinite and impedance matched transmission lines to the left and the right, as illustrated in placed by SQUIDs have recently been discussed [52,53]. We do not consider such modifications here, but the general approach we develop below can be used to formulate a quantum theory also in these cases. As shown in Appendix A, the position-dependent flux, φ(x) (in the Schrödinger picture), along a transmission line with a JTWPA section extending from x = 0 to x = z can in the continuum limit be expanded in terms of a set of left-and right-moving modes, where [â νω ,â † µω ] = δ νµ δ(ω − ω ) and the mode functions are given by Here, + (−) corresponds to ν = R (ν = L), k ω (x) = η ω (x)ω/v(x) is the wavevector, with η ω (x) the refractive index, and v(x) = 1/ c(x)l(x). The x-dependent parameters are defined such that they take one (constant) value inside the JTWPA section, and another value outside this section. c(x) is the capacitance to ground per unit cell, l(x) is the linear inductance of the transmission line. The only difference from the usual prescription for the quantized flux in a linear, homogeneous and dispersion free transmission line [54,55] is the x-dependent wavevector, which now takes a different form inside and outside the nonlinear section. Explicitly, the dispersion relation is found to be (see Appendix A) [8] k where z −1 (ω) = Z −1 (ω)/a is the admittance to ground per unit cell in the JTWPA section, and ω P is the junctions' plasma frequency. Implicit in the continuum description is that we are considering sufficiently low frequencies, such that the wavelengths are large compared to the unit cell distance,  3: Illustration of the disperion relation when Z(ω) (illustrated in the inset) describes a single resonant mode at a frequency ωr linearly coupled to the flux field in every unit cell. A bandgap opens up around the resonance frequency, close to 6 GHz in this example. The width of the bandgap is set by the coupling capacitance Cc shown in the inset.
a. Furthermore, plane wave solutions only exists when the right hand side of Eq. (3) is real. In practice, we are interested in frequencies ω 2 ω 2 P such that the dispersion due to the plasma oscillations of the junctions is relatively small. If, however, the admittance z −1 (ω) describes a linear element with a resonant mode, a bandgap opens around the resonance frequency for which no plane wave solutions exists. Physically such a resonant mode behaves as a "matter field" in the continuum limit, and the excitations of the systems resemble light-matter polaritons [49,56,57]. As long as we are away from any bandgap, however, these "matter fields" slave the photonic field,φ(x, t), and only modifies the dielectric properties of the medium, manifest in the dispersion relation Eq. (3). The behavior of the dispersion relation close to a bandgap is illustrated in Fig. 3.
As shown in detail in Appendix A a continuum limit Hamiltonian for the system can be written whereĤ 0 is a linear contribution containing all terms up to second order in the fields, andĤ 1 is a nonlinear contribution due to the Josephson junction potential. The linear Hamiltonian can be diagonalized in terms of the frequency modes introduced in Eq. (1), where we have omitted the zero-point energy.
For the nonlinear Hamiltonian, we systematically perform a series of approximations that are ultimately equivalent to those used in the classical treatment given in Refs. [8,10,44]. A quantized analog of the classical equation of motion found in previous work is shown to be a limiting case of a more general theory.
We assume the presence of a strong right-moving classical pump centered at a frequency Ω p with corresponding wavevector k p , and replaceâ Rω →â Rω + b(ω), with b(ω) a complex valued function centered at Ω p . The fields, a νω , are assumed to be sufficiently weak so that we can drop inĤ 1 terms that are smaller than second order in the pump. Also dropping fast rotating terms and the highly phase mismatched left moving field, this leads to the approximate Hamiltonian wherê describes cross-phase modulation due to the pump, and describes broadband squeezing. The dynamics of the classical pump is governed by a classical Hamiltonian which includes self-phase modulation, given in Eq. (A34). For notational convenience, we have defined the phase matching function [45,46] Φ(ω 1 , ω 2 , ω 3 , and a dimensionless pump amplitude, β(Ω), which as shown in Appendix A, can be written in terms of the ratio of the pump current to the Josephson junction critical current, where TreatingĤ 1 as a perturbation, it is natural to go to an interaction picture with respect toĤ 0 . The timeevolution operator for the problem in this picture iŝ whereĤ 1 (t) = exp iĤ 0 t Ĥ 1 exp −iĤ 0 t and T is the time-ordering operator. Solving the time-dynamics according to Eq. (11) is difficult in general [45,46]. A greatly simplified approximate theory can be derived, however, by 1) treatingĤ 1 as a perturbation to first order only, in which case the time-ordering in Eq. (11) can be dropped, and 2) taking the initial and final times to t 0 = −∞ and t 1 = ∞, respectively. The time-integral then gives rise to deltafunctions in frequency space, and we are left with approximate asymptotic evolution operator, or scattering matrix [49],Û Explicit expressions forK 1 for a general classical pump are given in Appendix A. We from now on focus on the monochromatic pump limit, taking b(ω) → b p δ(ω − Ω p ), with b p a c-number. In this limit we havê andK where we have defined β = β(Ω p ) and ∆k L (ω) here quantifies a phase-mismatch due to the linear dispersion in the JTWPA section. As we show below there is also an additional nonlinear contribution to the phase-mismatch that must be taken into account. Defining Heisenberg picture output fields,â out Rω = U †â RωÛ , we find the following input-output relation where the functions u(ω, z) and v(ω, z), defined in Appendix B, satisfy |u(ω, z)| 2 − |v(ω, z)| 2 = 1, and is the phase mismatch, including a nonlinear correction due to to the cross-and self-phase modulation of the pump. Eq. (18) is formally identical, up to a frequencydependent normalization of the wave amplitudes, to the classical solution derived in Refs. [8,10,44]. To summarize, this limiting equation is valid for weak nonlinearity and weak fields, only treating the nonlinear Hamiltonian H 1 to first order, a strong monochromatic classical pump at a frequency Ω p , and in an asymptotic large time limit where t 0 = −∞ and t 1 = ∞.
How can we interpret the asymptotic limit where the initial and final times are taken to minus and plus infinity, respectively? If we consider a situation where an initial wave packet is localized far away at x 0 at an early time t 0 0, this can be interpreted as a "scattering" limit, where we let the wave packet propagate through the nonlinearity and consider the asymptotic field at x z for a late time t 1 0 [48,49]. However, since the initial evolution before the wave packet enters the nonlinear section is governed byĤ 0 , it is trivial to propagate the wave packet forward towards the nonlinearity. The late evolution after the wave packet has left the nonlinear section is similarly trivial. We can therefore think ofâ Rω as a frequency domain input field entering the JTWPA andâ out Rω as the corresponding output field leaving the device. This is similar to the definition of input and output fields used in the description of damped quantum optical systems [48,58]. One should keep in mind, however, that the validity of this interpretation depends on the problem one is trying to solve: it is clearly not appropriate if, for example, the initial state of the field is delocalized over the nonlinear section.

III. ENGINEERING NONCLASSICAL RADIATION
The quantum input-output theory developed above allows us predict features of the JTPWA's output field, such as the device's gain profile and output field squeezing spectrum. In this section we show how output spectra can be tailored through dispersion engineering. We focus first on an ideal device and discuss the effect of loss in Sec. III B.

A. Ideal Squeezing Spectra
From Eq. (18), the JTWPA's amplitude gain is given by u(ω, z), and we define the power gain as G(ω, z) = |u(ω, z)| 2 [8,12,44]. This function grows exponentially with z for small phase mismatch, ∆k(ω) 0, but is only of order one if the phase mismatch is large (see Appendix B).
The squeezing of the JTWPA's output field is manifest in correlations between frequencies ω and 2Ω p − ω, symmetric around the pump frequency. It is convenient to define for the right moving field the thermal photon number the squeezing parameter and the squeezing spectrum [59] where we have defined quadratureŝ with fluctuations ∆Ŷ Rω =Ŷ Rω − Ŷ Rω , and θ is the squeezing angle, given through M R (ω, z) = |M R (ω, z)|e iθ . We emphasize that Eqs. (20)- (22) are defined exclusively in terms of the right-moving field. The left-moving field also contributes vacuum noise and might add to the total photon number, but will have zero squeezing parameter in the absence of left-moving pump fields. The squeezing spectrum is typically probed in experiments by heterodyne measurement of filtered field quadratures [42,[60][61][62]. We discuss this in some more detail in Appendix B. For a vacuum input field, where â R (ω, 0)â † R (ω , 0) = δ(ω − ω ) and all other second order moments vanish, it follows that , the maximum value allowed by the Heisenberg uncertainty relation and also imply quantum-limited amplification [12].
The gain and the squeezing at the output depends strongly on the phase mismatch ∆k(ω). The phase mismatch can however be compensated for by tuning Z(Ω p ), as this allows for tuning the pump wavevector k p = k Ωp according to Eq. (3). As was proposed theoretically in Ref. [8] and demonstrated experimentally in Refs. [9,10], it is possible to tune the phase mismatch to zero at the pump frequency, ∆k(Ω p ) 0, and greatly reduce it across the whole JTWPA bandwidth. This is done by placing LC (or transmission line) resonators with resonance frequency ω r0 Ω p regularly along the JTWPA transmission line. This technique is referred to as resonant phase matching (RPM) [8].
The effect of RPM on the gain and squeezing spectra is illustrated in Fig. 4 for a simulated device similar to what has been realized experimentally in Refs. [9,10]: The device length was chosen to be 2000 unit cells with characteristic impedance Z 0 = l/c = 50 Ω, critical current I c = (2π/Φ 0 )E J = 2.75 µA, dimensionless pump strength β = 0.125, junctions' plasma frequency Ω p /2π = 5.97 GHz, and pump frequency Ω 2 p /ω 2 P = 6.7 × 10 −3 . The green lines in Fig. 4 (a) show the gain profile and squeezing spectrum of the output field for the device without RPM, while the blue lines show results for an identical device where RPM has been used to tune ∆k(Ω p ) = 0. The circuit parameters for the LC resonator are C c = 10 fF, C r = 7.0 pF, L r = 100 pH, giving a resonance frequency of ω r0 /2π = 6.0 GHz. The corresponding impedances to ground in each unit cell, Z(ω), are illustrated schematically in Fig. 5 (a).  Two-mode squeezing has applications for entanglement generation [63,64], quantum teleportation [65], interferometry [66], creation of so-called quantum mechanics free subsystems [67], high-fidelity qubit readout [68,69] and logical operations [70], amongst others. A broadband squeezing source such as the JTWPA might have a great advantage for scalability, as tasks can be parallelized with many pairs of far-separated two-mode squeezed frequencies using a single device. It is, however, not necessarily desirable to have squeezing at all frequencies over the operational bandwidth as this might lead, e.g., to unwanted quantum heating [18]. This is the case, for example, for the qubit measurement scheme with Heisenberg limited scaling of the signal-to-noise ratio proposed in Ref. [69].
It is important in this scheme to not have high degrees of squeezing at the qubit frequency. The reason being that the qubit sees thermal noise with photon number N R (ω q , z), ω q being the qubit frequency-and even if the qubit is not directly coupled to the squeezing source, this leads to increased Purcell decay via the cavities (see Supplemental Material in Ref. [69]). This problem can be avoided with the JTWPA through dispersion engineering. In the following we show how to shape the squeezing spectrum to prohibit squeezing in certain frequency bands, and create spectra with a comb-like structure.
Building on the RPM technique, we consider placing additional resonances in each unit cell with resonance frequencies ω rk away from Ω p . This leads to a bandgap and a divergence in k(ω) close to each resonance ω rk , as illustrated in Fig. 3. The huge phase mismatch close to these resonances prohibits gain at ω ω rk and ω 2Ω p − ω rk , effectively punching two symmetric holes in the gain and squeezing spectra. This is illustrated by the orange lines in Fig. 4 (a), where a single additional resonace has been placed at ω r1 = 9.0 × 2π GHz. The parameters are otherwise as before, except that the second LC resonator is chosen to have twice the coupling capacitance, 2C c . This choice serves to illustrate how the width of the hole in the spectrum is determined by the coupling capacitance to the resonator, as is clearly seen when comparing the holes at ω r0 and ω r1 .
In Fig. 4 (b) we demonstrate how this technique can be used to engineer a "squeezing comb" where there is considerable gain and squeezing only for a discrete set of narrow quasimodes. With a larger number of closely spaced resonance frequencies-either using individual lumped LC circuits or the resonances of a multi-mode transmission line resonator-it is possible to have phase matching only over narrow frequency bandwiths. In Fig. 4 (b) we show the gain profile and squeezing spectrum where nineteen additional resonances at ω rk = ω r0 +k ×ω r0 /20, k = 1, 2, . . . , 19 has been used to create a squeezing comb with 38 quasimodes. Slightly different parameters were chosen for this device, to get similar gain and squeezing profiles as before: Z 0 = 14 Ω, I 0 = 2.75 µA, β = 0.069, while the additional coupling capacitances were chosen to be 3.0C c . The corresponding impedance to ground is illustrated in Fig. 5 For certain applications it might also be of interest to have a squeezing spectrum with a flatter profile than what is shown in Fig. 4. This can be achieved by suitably engineering the phase mismatch. In Fig. 6 we show a device where RPM has been used to tune ∆k(ω) = 0 for ω/2π 1.8 GHz, with the pump frequency close to to the resonance frequency at ω r0 /2π = 6 GHz. This leads to larger phase mismatch in the center region of the spectrum, close to the pump, giving the flatter profile shown in the figure. The simulated device otherwise has parameters Z 0 = 60 Ω, I 0 = 1.75 µA, β = 0.113.

B. Reduction in Squeezing Due to Loss
Internal loss in the JTWPA is likely to be a source of reduction in squeezing from the ideal results shown in Fig. 4. A simplified model for losses is a beam splitter with transmittance η(ω) placed after the JTPWA, with vacuum noise incident on the beam splitter's second input port [71]. This leads to a reduction in photon number, N R (ω, z) → |η(ω)|N R (ω, z), and squeezing parame- Taking η = η(ω) frequency independent for simplicity, this gives a reduction in squeezing, Note that distributed loss throughout the JTWPA can be taken into account through a simple phenomenological model [50], but this is beyond the scope of the present discussion. Fig. 7 shows the maximum squeezing level as a function of gain as the pump strength is ramped up. The parameters are otherwise identical to those used for the blue lines displayed in Fig. 4 (a). The solid lines show the maximally squeezed quadrature, while the dashed lines show the corresponding anti-squeezed quadrature, for three different values η = 0.75 (yellow), 0.99 (dark red) and 1.00 (blue). Note that the gain is also reduced by the loss, G(ω, z) = η|u(ω, z)| 2 , such that we have attenuation at zero pump power.
For a non-unity η, the squeezing level saturates with gain, while the anti-squeezed quadrature keeps growing proportionally. The maximal squeezing depends sensitively on η: while a quantum-limited device with η = 1 would produce more than 25 dB of squeezing at 20 dB of gain, a device with η = 0.75 only gives about 6.5 dB of squeezing for the same gain.

IV. PROBING THE OUTPUT
The examples discussed above demonstrate how the flexible JTWPA design allows for generating nonclassical light with interesting and useful squeezing spectra.
The squeezing spectrum can be found experimentally by measuring the variance of filtered two-mode quadratures (see Appendix B and, e.g., [42,[60][61][62]). However, this necessarily includes insertion loss and noise from subsequent parts of the amplification chain Ref. [9]. For a more direct probing of the JTWPA's performance we propose placing two superconducting qubits capacitively coupled to the transmission line at the output port.
For two off-resonant qubits with respective frequencies ω 1 = ω 2 , and ω 1 + ω 2 2Ω p , the qubits will be in uncorrelated thermally populated states. If, however, ω 1 +ω 2 = 2Ω p , the qubits become entangled and information about the JTWPA's squeezing spectrum is encoded in the joint two-qubit density matrix. This information can then be extracted by measuring qubit-qubit correlation functions.
Assuming for simplicity that the qubits are both located at the JTWPA output, x 0 > z, their reduced dynamics after tracing out the bath is governed by a Markovian master equation,ρ = Lρ. The form of L for the general case is given in Appendix C, while we here focus on the most interesting situation when the two qubits are tuned in with the squeezing interaction, such that ω 1 + ω 2 = 2Ω p , where ω m is the frequency of the mth qubit. We can then write the Lindbladian in the interaction picture is the decay rate of qubit m andσ − = |g e| (σ + = |e g|) is the qubit lowering (raising) operator. The Lindbladian has two contributions coming from the left-and the rightmoving field respectively. In general both fields can have non-zero thermal photon number N m,ν = N ν (ω m ) and squeezing parameter M ν = [M ν (ω 1 ) + M ν (ω 2 )] /2. If, on the other hand, the qubits are tuned out of resonance with the squeezing interaction, ω 1 + ω 2 2Ω p , the last line in Eq. (24) will be fast rotating and can be dropped in a rotating wave approximation (see Appendix C for more details).
Assuming for simplicity a single right-moving pump and a left-moving field in the vacuum state, we have that for ω 1 + ω 2 2Ω p , the steady state of the two qubits is the product state ρ = ρ 1 ⊗ ρ 2 , where ρ m is a thermal state with thermal population N R (ω m )/2 and inversion σ . On the other hand, for ω 1 + ω 2 = 2Ω p the qubits become entangled. Under the simplifying symmetric assumptions N R (ω m ) ≡ N and γ m ≡ γ we find that and in steady state, where M (ω i ) ≡ M . More general expressions are given in Appendix C. Hence, by measuring qubit-qubit correlation functions and single-qubit inversion using standard qubit readout protocols [72][73][74], one can map out the squeezing spectrum and the quantum efficiency of the JTWPA (the latter also requires knowledge of the thermal noise at the input, which could be probed in a similar way by a single qubit located at the input port, see [9] for a similar experiment). We can also turn this around and, rather than view the two qubits as a probe of the JTWPA's performance, view the JTWPA as a source of entanglement for the qubits. To achieve maximal degree of entanglement between the qubits, it is desirable to avoid the vacuum noise of the left-moving field. This can be achieved by squeezing the left-moving field with a separate JTWPA section, or more simply by operating the device in reflection mode, as illustrated in Fig. 8 (c).
Assuming ideal conditions where the qubits couple symmetrically to equally squeezed left-and right-moving , and ideal lossless squeezing, the steady state of the two qubits is the pure state (see Appendix C for more information) where θ is the squeezing angle. For large N , this pure state approaches a maximally entangled state with entan- Of practical importance is the steady state entanglement's dependence on the degree of loss, and the behaviour of the spectral gap of the Lindbladian in Eq. (24). The latter is important because it sets the time-scale for approaching the steady state. It is defined as ∆(L) = |Re λ 1 |, where λ 1 is the non-zero right-eigenvalue of L with real part closest to zero. In Fig. 9 we plot the steady state entanglement, quantified by the concurrence [75], and the spectral gap as a function of gain for different values of η (as defined in Sec. III B). These results show that the achievable entanglement is very sensitive to loss, but an upshot is that relatively modest gains are needed to achieve high degree of entanglement, which might facilitate creating devices with higher η. Furthermore, note that multiple pairs of qubits can be entangled using a single squeezing source. Due to the large bandwidth of the JTWPA several tens of entangled qubit pairs can likely be generated in this way using a single device.
Our scheme for probing the squeezing spectrum is similar to previous proposals for extracting information about single-mode squeezing through the resonance fluorescence emitted by a single atom [76,77]. The predictions of Refs. [76,77] were recently confirmed experimentally using a superconducting artificial atom coupled to the squeezed output field of a Josephson parametric amplifier [78]. Our scheme extends this to probing correlation between different frequency components of the squeezed radiation by going from a single to two qubits.

V. CONTINUOUS VARIABLE CLUSTER STATES
The two-qubit dynamics considered above demonstrates the JTWPA's potential for entanglement generation. By adding multiple pump tones, a single frequency can become entangled with multiple other "idler" frequencies in a multi-mode squeezed state, and complex patters of entanglement can emerge. Together with its broadband nature and the potential for dispersion engineering, this turns the JTWPA into a powerful resource for dissipative quantum state engineering.
As a demonstration of the JTWPA's potential as a source of nonclassical radiation, we show below how continuous variable (CV) cluster states can be generated through a dissipative and deterministic process, using the output field of multiple JTWPAs. Th cluster states are a powerful class of entangled many-body quantum states that are resource states for measurement based quantum computing. Given a universal cluster state, an algorithm is executed using only single-site measurements and classical feed forward on the state [19,20,[79][80][81].
A CV cluster state is defined with respect to a (weighted) simple graph G = (V, E), with V the set of vertices and E the set of edges. A CV quantum sys- is a bosonic annihilation (creation) operator, is associated to each vertex v. The ideal CV cluster state (with respect to G) is defined as the unique state |φ G satisfying [20,81,82] ŷ where N (v) is the neighborhood of v, i.e., all the vertices connected to v by an edge in E and a vw = a wv ∈ [−1, 1] is the weight of the edge {v, w} [96]. The operators L v ≡ y v − w∈N (v) a v,wxw are referred to as the nullifiers of |φ G . We can define an adjacency matrix A = [a vw ] for the graph, where a vw = 0 if there is no edge {v, w} ∈ E.
Since the adjacency matrix uniquely defines the graph, and vice versa, we use the symbol G to interchangeably refer to both the graph and its adjacency matrix in the following.
We focus here on a class of graphs, first studied in Refs. [31,32], satisfying two simplifying criteria: 1) The graph is bicolorable. This means that every vertex can be given one out of two colors, in such a way that every edge connects vertices of different colors (the square lattice is an example).
2) The graph's adjacency matrix is selfinverse, G = G −1 . The latter constraint has a simple geometric interpretation described in Ref. [32]. We show in Appendix D that for a graph G satisfying these critera, the Lindblad equationρ = L G ρ, with Lindbladian where M = N (N + 1) and S iM [A, B] is defined in Eq. (25), has a unique steady state |φ G (M ) that approaches |φ G as M → ∞. The existence of graphs satisfying all the listed criteria, with associated cluster states, |φ G , that are universal for quantum computing, was shown in Refs. [31,32].
In Ref. [34] Wang and coworkers showed how cluster states with graphs of the type considered here could be generated through Hamiltonian interactions between the modes of optical parametric oscillators (OPOs), followed by an interferometer combining modes from distinct OPOs. We adopt this scheme in the following, using JTWPAs (other types of broadband squeezing sources can also be used) in place of OPOs. The main difference between our proposal and that of Ref [34] and previous proposals [31,32] is that our scheme is purely dissipative: the CV modes of the cluster state never interact directly, but rather become entangled through absorption and stimulated emission of correlated photons from their environment. We focus primarily on a situation where the modes are embodied in multimode resonators, which is a particularly hardware efficient implementation. We emphasize, however, that due to the dissipative nature of the scheme, this is not a necessary constraint. The modes could in principle all be embodied in physically distinct and remote resonators, removing any constraints on locality. This is a distinct advantage of such a dissipative scheme.
Following Ref. [34], the modes of the cluster states are resonator modes with equally spaced frequencies ω m = ω 0 + m∆, where m is an integer, ω 0 is some frequency offset and ∆ the frequency separation. We require a number of degenerate modes for each frequency ω m : to create a D-dimensional cluster state requires a 2 × D-fold degeneracy per frequency. This can be achieved using 2 × D identical multi-mode resonators, as illustrated for D = 1 in Fig. 10. Each mode is a vertex in the cluster state graph, and as will become clear below, a set of degenerate modes can be thought of as a graph "macronode" [34]. It is convenient to relabel the frequencies with a "macronode index" M = (−1) m m.
We show in Appendix D that a master equation with Lindbladian of the form Eq. (30) is realized for a single resonator interacting with a bath generated by the output field of a JTWPA with a single pump frequency Ω p = ω 0 +p∆/2 where p = m+n for some choice of frequencies ω m = ω n . The graph is in this case a trivial graph consisting of a set of disjoint pairs of vertices connected by an edge, i.e., a set of two-mode cluster states which can be represented as G 0 = . . . The edges have weight +1, under the assumption of a quantum limited, flat squeezing spectrum M (ω) = iM = i N (N + 1) with N (ω) = N over the relevant bandwidth.
More complex and useful graphs can be constructed using these two-mode cluster states as basic building blocks [34]. Taking a number of JTW-PAs, each labelled by i and acting as a squeezing source independently generating a disjoint graph G i = as above, universal cluster states can be created by combining the output fields of the different sources on an interferometer. The action of the interferometer can be written as a graph represents an interferometer acting independently on each macronode M, i.e., each set of 2 × Dfold degenerate modes. R has to be orthogonal for the transformed graph to be self-inverse, G = G −1 , which we recall is one of the criteria for Eq. (30) to generate the corresponding cluster state. As shown in Ref. [34] this is the case if the 2D × 2D matrix H D is a Hadamard transformation H D = H ⊗D built from 2 × 2 Hadamard matrices Physically such a transformation can be realized by pairwise interfering the output fields of the JTWPAs on 50-50 beam splitters with beam splitter matrix as in Eq. (31). The network of beam splitters needed for the case D = 1 is illustrated in Fig. 10, for D = 2 in Fig. 11, and for higher dimensions in Ref. [34]. In Ref. [34] it was shown that graphs G constructed in this way can give rise to D-dimensional cluster states that are universal for measurement based quantum computing. Let us consider an example with D = 1 in some more detail to illustrate the basic principles, while referring the reader to [34] for more details. First, take two JTWPAs pumped individually with respective pump frequencies Ω i and Ω j , with i = −j = ∆M. On the macronode level, this gives exactly one edge between macronodes separated by |∆M|, as illustrated by the horizontal edges in Fig. 10  by Eq. (31), every node in each macronode becomes entangled with every node in the neighboring macronode, as illustrated by the diagonal edges in the figure. This gives a graph G with a linear structure, corresponding to a one-dimensional cluster state that is universal for single-mode quantum computation [32][33][34].
The scheme can straight-forwardly be scaled up to arbitrary D-dimensional cluster states using 2×D JTWPAs and the same number of beam splitter transformations as shown in Ref. [34]. D = 2 is sufficient for universal quantum computation; a possible setup of JTWPAs and resonators is illustrated in Fig. 11. As emphasized in Ref. [34], the relative ease of creating even higher dimensional cluster states is a very attractive property of the scheme. D = 3 might allow for error correction with high thresholds based on surface-code encodings [83], and D ≥ 4 might allow for simulating systems with topological self-correcting properties [84]. trum. We have shown how one can create multi-mode squeezed baths that can be used for dissipative quantum state preparation. By placing quantum systems at the output of a broadband squeezing source, the systems are cooled into non-trivial entangled states by interacting with a multi-mode squeezed vacuum. In particular we have shown how to prepare pairs of entangled qubits and continuous variable cluster states that are universal for quantum computing. In both cases the large bandwidth of the JTWPA makes the state preparation highly hardware efficient.
The ability to prepare cluster states demonstrates the universal power of broadband squeezing as a resource. We hope this motivates experimental efforts to demonstrate high degrees of squeezing over large bandwidths. It should also motivate a broader theoretical study of how squeezing sources such as the JTWPA can be used to generate quantum radiation with complex entanglement structures geared towards particular applications in quantum technology and information processing.  Fig. 1. Each junction has Josephson energy E J , junction capacitance C J , and is coupled to ground by an impedance Z(ω), describing a reactive and dissipationless element. We first treat the case of a single capacitance to ground, Z(ω) = 1/iωC, before considering the more general case where the impedance contains resonances.

Without Resonances
We write the Lagrangian of the JTWPA in terms of the node fluxes, φ n , following the standard lumped element approach [86], where ∆φ n = φ n+1 − φ n , Φ 0 = h/2e is the magnetic flux quantum, and we have expanded the Josephson junction cosine potential to fourth order.
In experimental realizations, the number of junctions is of the order of a few thousand, and the unit cell distance, a, is much smaller than the relevant microwave wavelengths. This justifies a continuum limit treatment of the device. Taking and defining the continuum parameters where ω P = (2π/Φ 0 ) E J /C J is the junctions' plasma frequency, we can formally take the continuum limit N → ∞, a → 0 such that the length of the device N a ≡ z is kept constant. A continuum Lagrangian can then be introduced The x-dependent parameters in the above expression are defined such that they take the values in Eq. (A5) for 0 < x < z, while outside this region we take c(x) = c 0 , l(x) = l 0 , ω P (x) = ∞ and γ = 0. Eq. (A6) thus represents a JTWPA section extending from x = 0 to x = z, sandwiched between two identical semi-infinite linear transmission line sections extending from x = −∞ to x = 0 and x = z to x = ∞, respectively. We furthermore assume that the sections are impedance matched, Z 0 = l 0 /c 0 = l/c, ensuring that there are no reflections at the boundaries.
The Euler-Lagrange equation found from Eq. (A6) is where we have left out the x and t dependence of the fields and the parameters for notational simplicity. It is useful to find stationary solutions of the form φ(x, t) = φ(x)e −iωt to the linear part of Eq. (A7), i.e., the left hand side only of this equation. These classical solutions serve to determine the spatial dependence of the field in the absence of nonlinearity and are useful for constructing quantized solutions to the full problem [87]. We find that solutions of the form φ(x, t) = Ae −iωt+ikωx satisfies the linear part of Eq. (A7) with a dispersion relation where v 1 = 1/ √ cl and v 0 = 1/ √ c 0 l 0 . Of course, the right hand side of Eq. (A8a) has to be positive for traveling wave solutions to exist. In practice, we are interested in frequencies ω 2 ω 2 P , such that the dispersion due to the junctions' plasma oscillations is relatively small.
We now wish to transform to a Hamiltonian. The canonical momentum to φ(x, t) is found from Eq. (A6) in the usual way, where L is the Lagrangian density. Note that this differs from the usual prescription, π = ∂L/∂(∂ t φ), due to the term proportional to 1/ω 2 P [49]. From this we can define which after an integration by parts and dropping a boundary term reads H = H 0 + H 1 with and Note that although we express the Hamiltonian density in terms of ∂ t φ at this stage, this should be considered a function of π. Quantization follows by the usual prescription of promoting the fields to operators φ(x, t) →φ(x, t), π(x, t) → π(x, t) and imposing the canonical commutation relation We attack the problem by first considering only the linear Hamiltonian H 0 , and diagonalizing this part by expandingφ(x, t) in a set of mode functions and inserting into the Hamiltonian, not taking the interaction H 1 into account for the moment. We will subsequently use these modes to construct the full-nonlinear Hamiltonian. We follow closely the treatment of a dielectric with an interface given by Santos and Loudon in Ref. [88] and write for the flux where [â νω ,â † µω ] = δ νµ δ(ω − ω ), and mode functions where the + sign is for ν = R and the − sign is for ν = L, with k ω (x) = η ω (x)ω/v(x) the wavevector. The refractive index, η ω (x), has a different value inside and outside the JTWPA section, according to Eq. (A8). Note that the refractive index should be real and symmetric, as we are not including absorption in the theory [88]. Using Eqs. (A9) and (A14) this gives for the canonical momentumπ(x, t), where we have defined the dielectric function Using Eq. (A8) this can also be expressed as (A18) Note that the canonical commutation relation, Eq. (A13), implies the following condition on the mode functions This equality was proven to hold true for canonical fields of the form Eqs. (A9) and (A14) in Ref. [88] for real and symmetric η ω (x) taking, as in Eq. (A8), different constant values in different sections of a dielectric with interfaces. As was pointed out by these authors, however, the theory is not entirely satisfactory as traveling wave solutions do not exist for all frequencies. As a result the integration over all frequencies is not strictly valid [88]. A more careful analysis has to take absorption into account leading to localized excitations of the junction degrees of freedom.
Using the mode decomposition ofφ(x, t) andπ(x, t), we can diagonalize the linear HamiltonianĤ 0 . This follows directly from the results from Ref. [88]. To that end, we first note that from Eq. (A11) we can write the differential equation (A20) After an integration by parts and dropping boundary terms, this reduces to which is exactly Equation (3.3) of Ref. [88] (with the identificationsÊ → −∂ tφ ,B → ∂ xφ ,D → −π, µ 0 → l). It immediately follows from the results there [Eqs. (3.3) to (3.6)] that, after performing the integration over space, the linear Hamiltonian becomeŝ where we have omitted the zero-point energy.

With Resonances
We now consider a general impedance Z(ω) describing a reactive, dissipationless element, coupling the node fluxes to ground in each unit cell. In general, we can represent this impedance with a circuit as illustrated in Fig. 12, and write for the inverse impedance (the admittance) [86,89], where we include a zero-frequency component represented by the capacitance C. Following identical steps . 12: Representation of a general lossless admittance, Z −1 (ω), in terms of a set of modes, ψnq, with resonance frequencies ωq = 1/ LqCq. Note that we assume the presence of a zero-frequency component represented by the capacitance C.
as before, the linear part of the continuum Lagrangian is now modified to where c q (x) = C q /a and l q (x) = L q /a for 0 < x < z are the capacitance and inductance per unit cell associated to the q th resonance, and ψ nq (t) → ψ q (x, t) is the continuum limit field for the resonance with frequency ω q . We take c q (x) = 0, l q (x) = ∞ outside the JTWPA section. Similarly the linear Hamiltonian now reads after an integration by parts and dropping a boundary term The field ψ q (x, t) plays the role of a "matter field," with a single resonance at ω q = 1/ L q C q coupling linearly to the "photonic field" φ(x, t). This type of light-matter coupling is a standard microscopic model, based on the Hopfield model [90], for dispersion in linear media [49,56,57].
The Euler-Lagrange equations for the fields are As before we seek stationary solutions of the form φ(x, t) = Ae −iωt+ikωx and ψ q (x, t) = B q e −iωt+ikωx .
The Euler-Lagrange equations can be used to find the dispersion relation of the medium, but also to relate the frequency components of the fields ψ q (x, t) to those of the field φ(x, t) [49]. Indeed, solutions are found for Eq. (A26) with and inside the JTWPA section and k 2 ω = ω 2 /v 2 0 elsewhere, as before. z −1 (ω) = Z −1 (ω)/a is here the admittance to ground per unit cell. Using Eq. (A14) for φ(x, t) and Eq. (A27) to relate the frequency components of the traveling wave solutions for the matter fields to those of the photonic field, we can write Performing steps analogous to those in the previous section again leads to a diagonal linear Hamiltonian, H 0 , as given in Eq. (A22). First we find a differential equation forĤ 0 , where we have used Eq. (A26) and performed a partial integration, dropping a boundary term. This is of the same form as Eq. (A21) if we define π(x, t) as in Eq. (A16) with a dielectric function which from Eq. (A28) reads (A32) With this, we can promote the fields to operators, and the canonical commutation relation as well as the diagonal form ofĤ 0 follows, again, from the results of Ref. [88]. Note that the numerator of the refractive index given by Eq. (A32) has the standard form of a Sellmeir expansion that one expects for a dispersive medium with material resonances at ω q [49]. The plasma oscillations of the junctions, however, play a somewhat different role from the "material" resonances due to the admittance Z −1 (ω). Indeed, a term of the type (∂ x ∂ t φ) 2 , as appearing in our Hamiltonian, is used in macroscopic models for dispersive dielectrics [49]. It is therefore interesting to note that we here have two distinct "types" of dispersion present at the same time-a term (∂ x ∂ t φ) 2 in the Hamiltonian, stemming from junction plasma oscillations, and resonances associated to harmonic oscillators coupled linearly to the fieldφ(x, t)-and both effects arise from a microscopic theory in our case.
It is of course important to keep in mind that traveling wave solutions only exists for frequencies such that the right hand side of Eq. (A28) is positive. Introducing resonances as illustrated in Fig. 5 introduces band gaps in the dispersion relation, as illustrated in Fig. 3. The theory developed here is valid away from any bandgap. A more general treatment would require keeping the "matter field" degrees of freedom and subsequently diagonalizing the linear Hamiltonian by introducing polariton fields [56,57].

Nonlinear Hamiltonian
Having diagonalized the linear Hamiltonian H 0 we now consider the nonlinear contribution We use the expansion ofφ(x, t) in frequency modesâ νω introduced above, and assume that the fieldsâ νω are small except for close to the single pump frequency. More precisely, we assume a strong right moving, classical pump centered at a frequency Ω p and corresponding wave number k p , and takeâ Rω →â Rω + b(ω), with b(ω) a c-number. After dropping fast-rotating terms, the highly phasemismatched left-moving field, and neglecting terms smaller than O[b(ω) 2 ] we arrive at Eqs. (6)- (8) in Sec. II. Note that there is also a contribution fromĤ 1 involving only the classical pump, which describes self-phase modulation. Here, the dimensionless pump amplitude is defined as where and we drop the x-argument on variables when there is no danger of confusion, as γ = 0 only inside the JTWPA section. The dimensionless pump amplitude can also be written in terms of the ratio of the pump current to the Josephson junction critical current, where I c = (2π/Φ 0 )E J is the critical current, and we have related the pump amplitude to the pump current through B(Ω) = I p (Ω)Z 0 /Ω where Z 0 = l/c is the characteristic impedance.
As explained in Sec. II, by defining an interaction picture evolution operator and taking the initial and final times to minus and plus infinity respectively, we can define an asymptotic evolution operator (to first order in describes cross-phase modulation due to the pump, describes broadband squeezing, and describes pump self-phase modulation. Finally, taking the monochromatic pump limit, b(ω) → b p δ(ω − Ω p ), with b p a c-number, we arrive at Eqs. (14) and (15), while the classical pump Hamiltonian is simply from which the squeezing spectrum can be computed easily. The squeezing spectrum is typically probed in experiments by heterodyne measurement of filtered field quadratures [42,[60][61][62]. We here give some more details on how the squeezing spectrum defined in Eq. (22) can be measured in this way. We first define filtered output fields [43,60] and refers to a narrowband low-pass filter capturing the experimental bandwidth [91]. We also define filtered quadraturesX Two-mode squeezing refers to the fluctuations in some of these two-mode quadratures being below the vacuum level. In the present case we find squeezing in the fluctuations where the overline denotes time-averaging. For a vacuum input field we find that the two-mode fluctuations are giving a filtered version of the squeezing spectrum defined in Eq. (22).

Appendix C: Squeezed Bath Engineering
In this appendix we derive a master equation for a collection of quantum systems interacting with multiple multimode squeezed baths, and give detailed results for the case of two qubits placed at the output of a JTWPA pumped by a single classical pump.

Master Equation for Multiple Systems and Multiple Squeezed Baths
We consider a number of quantum systems distributed along linear sections of a set of transmission lines. We label the transmission lines by the index i and the quantum systems coupled to the ith line by an index m i . The systems couple to the transmission line voltages, , at some position x i , via system operatorsĉ mi . These system operators are furthermore assumed to rotate with frequencies ω mi in the interaction picture. The interaction Hamiltonian in the interaction picture is then of the form are bath operators and κ mi describes the coupling strength of system m i to the ith bath.
We do not consider a situation where the systems are cascaded [92,93], i.e., none of the systems are driven by the output of any of the other. The different baths can, however, be correlated, with a defining set of correlation functions where N ji (ω) = N ij (ω) and M k ji (ω) = M k ij (ω). This generalizes Eq. (B10) and includes, for example, setups of the type illustrated in Figs. 10 and 11 where the output field of multiple JTWPAs are combined on beam splitters before being incident on the quantum systems.
Following the standard approach [94], we find the following master equation for the reduced system density matrix under the usual Born-Markov approximatioṅ where we have defined bath correlation functions with ω > 0 and the brackets refer to an expectation value with respect to the bath density matrix ρ B . At this stage we can invoke a rotating wave approximation (RWA). For example, we can approximate where we have used and dropped the small principal part. We find similar expressions for the other terms in Eq. (C4), leading to the interaction picture Master equatioṅ Note that this master equation includes dissipative interactions between systems even if they are not coupled to the same bath, i = j, if the baths are correlated according to Eq. (C3). As we will see in Appendix D, this is the type of interactions that give rise to the diagonal edges in the cluster state graph in Fig. 10. But first, we focus in the next subsection on the simple case of two qubits coupled to the left-and right-moving fields of a single transmission line.

Two Qubits
In this section we consider two qubits coupled to the left-and right moving field of a transmission line (in general both fields can be squeezed). We take both fields to have correlation functions as in Eq. (B10) without any cross correlations between the two fields (N ji , M ji ∝ δ ji in Eq. (C3)).
The qubits have free HamiltoniansĤ m = ω mσ where we have dropped off-resonant terms rotating at and γ m the decay rate of qubit m, which has contributions from both the left-and the right-moving field. N m,ν = N ν (ω m ) and M ν = [M ν (ω 1 ) + M ν (ω 2 )] /2 are, respectively, thermal photon number and squeezing parameters for each field. If the qubits are not tuned into resonance with the squeezing interaction induced by the pump fields, ω 1 + ω 2 2Ω ν , the last two lines of Eq. (C9) can be dropped in an RWA. In this case, the two qubits are effectively interacting with independent thermal baths, and will reach an uncorrelated thermal steady state, ρ ss = ρ 1 ⊗ρ 2 , where where N m,tot = N m,R + N m,L is the total thermal photon number at the frequency of qubit m. For a left moving field in the vacuum state, N L,m = 0, this gives the atomic inversion quoted in Sec. IV. As discussed in Sec. IV, when the qubits are tuned into resonance with the squeezing interaction, they become entangled and encode information about the squeezing spectrum. This can be exploited to use the qubits as a spectroscopic probe. Assuming a left moving field in the vacuum state, and taking ω 1 + ω 2 = 2Ω R , Eq. (C9) reduces to The steady state of Eq. (C12) can be found analytically if we assume that the qubits start in the ground state. We find (in the basis {|ee , |eg , |ge , |gg }) where This gives steady state expectation values from which both N m and M m , and thus the squeezing spectrum of the source, can be extracted. For the simplifying symmetric assumptions γ 1 = γ 2 ≡ γ and N 1 = N 2 ≡ N , this leads to Eqs. (26) and (27). Another interesting case is when both the left-and right moving fields are squeezed or, equivalently, the qubit couples to a single squeezed field. This could be realized, e.g., by ending the transmission line with a mirror as indicated in Fig. 8 (c). This gives higher degrees of entanglement between the qubits, as the vacuum noise of the left moving field otherwise reduces correlations between them. In particular, taking N m,ν = N/2, M ν = M/2 and γ 1 = γ 2 in Eq. (C9), and furthermore assuming |M | = N (N + 1), we find that the steady state is given by the pure state Eq. (28).

Appendix D: Generating cluster states
In this appendix we give more details regarding the dissipative generation of CV cluster states. Recall that the cluster state |φ G with respect to a graph G is the unique pure state satisfying where −1 ≤ a vw ≤ 1 is the weight of the (undirected) edge {v, w} and N (v) denotes the neighborhood of v. We can define an adjacency matrix, A = [a vw ] (with a vw = 0 if there is no edge {v, w} in the graph), and we use the graph G and its adjacency matrix A interchangeably when referring to the graph in the following. The general question of dissipatively preparing pure Gaussian states was addressed by Koga and Yamamoto in Ref. [95] and, in particular, they found a Lindblad master equationρ = L G (ε)ρ whose unique steady state approaces Eq. (D1), in the limit of infinite squeezing. A Lindbladian that achieves this is where κ is a decay rate setting the overall time-scale and the Lindblad operators are The steady state approaches Eq. (D1) in the limit ε → 0 + . In the special case of a bicolorable graph, meaning that every vertex can be given one out of two colors in such a way that every edge connects vertices of different colors, the Lindbladian Eq. (D2) is to first order in ε Recall that E is the set of all pairs of vertices {v, w} connected by an edge. Furthermore, we have defined E 2 as the set of all pairs of distinct vertices {w, w } such that there exists some v = w, w such that w, w ∈ N (v). In other words, E 2 is the set of all next-nearest neighbors. Note that the sum over edges in E is a sum over undirected edges, i.e., there is one term in the sum for every pair of vertices connected by an edge, and similarly for the sum over E 2 . We have furthermore defined The former can be recognized as the sum of the weights of all "2-paths" starting and ending at v, and the latter as the sum of the weights of all "2-paths" connecting the two vertices {w, w }. The weight of a 2-path is here defined to be the product of the weights of the two respective edges. Eq. (D4) takes a particularly simple form if the following two geometric conditions for the graph G are satisfied: 1. The sum of all 2-paths starting and ending at the same vertex add up to one: D v = 1.
2. The sum of all 2-paths connecting two different vertices cancel out: D w,w = 0.
Note that these two criteria are equivalent to saying that the adjacency matrix defining the graph is self-inverse, A = A −1 [32]. In this case, Or, by defining |M | = 1/2ε, and N through |M | = N (N + 1) we can write to first order in ε (D8) Next we show how Eq. (D8) can be generated through squeezed bath engineering. Remarkably, the three criteria we have assumed for the graph G, that it is bicolorable and satisfies the two conditions listed above, are exactly the criteria first studied by Menicucci, Flammia and Pfister in Refs. [31,32]. These authors showed that graphs satisfying these criteria corresponding to cluster states that are universal for quantum computing exists, and moreover how such cluster states can be prepared using a single optical parametric oscillator (OPO). Later, Wang and coworkers [34] presented an alternative implementation using a network of OPOs and beam splitters.
We adopt the scheme from Ref. [34] and use JTWPAs as squeezing sources in place of OPOs. The main difference between what we propose here and previous proposals is the dissipative nature of the interactions. While in Ref. [34] and previous work [31,32] the schemes were based on Hamiltonian interactions between internal cavity modes of OPOs, we are using dissipative interactions between resonator modes that are external to the squeezing sources. The CV modes of the cluster state never interact directly, but rather become entangled through absorption and stimulated emission of correlated photons from their environment. Even so, the graph constructions from Ref. [34] can be adopted almost directly for our purposes as well.
Following Ref. [34] the modes of the cluster states will be resonator modes with equally spaced frequencies ω m = ω 0 + m∆, where ω 0 is some frequency offset and ∆ the frequency separation. We require a number of degenerate modes for each frequency ω m : to create a D-dimensional cluster state requires a 2 × D-fold degeneracy per frequency. Each mode is a vertex in the cluster state graph, and as will become clear below, a set of degenerate modes can be thought of as a graph "macronode" [34]. It is convenient to relabel the frequencies ω m with a "macronode index" M = (−1) m m.
A simple implementation is to embody the modes of the cluster state in 2×D identical multi-mode resonators. Each resonator is assumed to couple to a single squeezed bath, which we label by an index i. We use the same number of squeezing sources, and each source is pumped by a frequency Ω i = ω 0 + i∆/2 where i = m + n for some choice of frequencies ω m = ω n . The key to creating useful cluster states is to correlate the different baths by combining the output fields of broadband squeezing sources using interferometers.
We first consider the simplest situation where each resonator is coupled to the output field of a single squeezing source, as illustrated in Fig. 13 where we assume that M j (ω) > 0 is real. With a JTPWA as a squeezing source, as long as the phase mismatch is small ∆k(ω) 0, the squeezing angle is frequency independent and M (ω) can thus be assumed real without loss of generality since any complex phase can be absorbed in a redefinition of the system operatorsĉ v . Using this in Eq. (C8) and doing an RWA leaves us with the Lindbladian where we have labeled the resonator modes by a vertex index v, N v = N j (ω v ), with ω v the frequency of mode v of the jth resonator, and the matrix M vw is non-zero with value M vw = M j (ω v ) only if ω v + ω w = 2Ω j and v and w both are associated to the jth resonator. This is of the form Eq. (D8) if we assume equal decay rates, κ v ≡ κ and N v ≡ N and we write M vw = |M |a vw with |M | = N (N + 1) and −1 ≤ a vw ≤ 1. The graph generated by Eq. (D10), defined by the adjacency matrix A ∼ M = [M vw ], is however not useful for measurement based quantum computing: it consists of a set of disjoint edges, as illustrated in Fig. 13 (c). Note that for the adjacency matrix of this graph to be self-inverse, A = A −1 , requires edges with weights a vw = 0, ±1, meaning a flat, quantum limited squeezing spectrum such that |M j (ω)| = |M i (ω )| = |M | = N (N + 1) within the relevant bandwidth. We assume in the following that M (ω) = M > 0, such that the graph consists of a collection of identical two-mode cluster states with egdge weights +1.
Note that Eq. (D10) can be brought into Lindblad form, withLindblad operators where we have used that A = A −1 and that we can write N = sinh(r) 2 , M = sinh(r) cosh(r).
There are, of course, no edges between vertices associated to modes embodied in different resonators at this stage, as they interact with independent, uncorrelated baths. The key to generate non-trivial cluster states is to combine the output fields of the squeezing sources on beam splitters before being incident on the resonators, as illustrated for the D = 1 case in Fig. 10 in Sec. V.
A general two-port beam splitter transformation for a pair of fields, givesb iω → t i (ω)b iω + r i (ω)b jω , where t i (ω) and r i (ω) are the transmittance and reflectance coefficients. The correlation functions in Eq. (D9) thus transform to b † jωb iω = δ ij N j δ(ω − ω ), where we have suppressed the ω argument on the functions on the right hand side for notational simplicity, and used that unitarity of the beam splitter requires t * j r i + r * j t i = 0. We see that the form of the Lindbladian Eq. (D8) is preserved, with a transformed squeezing matrix, M = [M vw ] → RM R T , where R = R ij represents the beam splitter matrix R ij acting on all pairs of degenerate modes belonging to resonators i and j.
Following Ref. [34] we take the beam splitter to be a Hadamard transformation such that the squeezing matrix where R = H represents a Hadamard matrix acting on pairs of degenerate modes. Since M ∼ A, the graph of course transforms in the same way. For the D = 1 case this gives the graph illustrated in Fig. 10. Note that the transformed adjacency matrix RAR T is still self-inverse since R is orthogonal [34].
This can be generalized in a straight forward manner to combine multiple squeezing sources. For 2 × D sources, we can construct a 2D × 2D Hadamard matrix by composing together 2 × 2 Hadamard matrices, H D = H ⊗D [34]. In practice this is done by combining the output fields pairwise on two-port beam splitters in a cascaded fashion, as illustrated for the case D = 2 in Fig. 11 and explained in more detail in Ref. [34]. The total transformation on the squeezing matrix is R = M H D , representing a 2D × 2D Hadamard transformation on each macronode M, i.e., each set of 2 × D degenerate modes.
As shown in Ref. [34], graphs defined by the adjacency matrix A ∼ M corresponding to D-dimensional cluster states can be generated in this way, where in particular the D = 2 case is universal for measurement based quantum computing.